Skip to content

Commit c8c2e62

Browse files
rreusserkgryte
authored andcommitted
Add Fresnel integral
* Add Fresnel integral * Add equation to readme * Add additional Fresnel integral examples * Add equation to readme * Fix incorrect value in example output * Remove unneeded blank line * Remove unused fixture generation code * Fix labels in tests * Fix test labels and comment format * Remove equation file * Fix wording * Add Cephes method and copyright * Fix incorrect values in repl docs * Add a note on scaling * Fix PR review issues * Align equations * Align equations * Fix formatting issues * Move scipy benchmark
1 parent c3684cd commit c8c2e62

19 files changed

Lines changed: 1406 additions & 0 deletions

File tree

docs/links/database.json

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1798,6 +1798,17 @@
17981798
"frechet"
17991799
]
18001800
},
1801+
"https://en.wikipedia.org/wiki/Fresnel_integral": {
1802+
"id": "fresnel-integral",
1803+
"description": "Wikipedia entry for Fresnel integral.",
1804+
"short_url": "",
1805+
"keywords": [
1806+
"math",
1807+
"special",
1808+
"fresnel",
1809+
"integral"
1810+
]
1811+
},
18011812
"https://en.wikipedia.org/wiki/Function_composition_%28computer_science%29": {
18021813
"id": "function-composition",
18031814
"description": "Wikipedia entry for function composition.",
Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
# fresnel
2+
3+
> Compute the [Fresnel integrals][fresnel-integral] S(x) and C(x).
4+
5+
<section class="intro">
6+
7+
The [Fresnel integrals][fresnel-integral] are defined as
8+
9+
<!-- <equation class="equation" label="eq:fresnel_integrals" align="center" raw="\begin{align} S(x) &= \int_0^x \sin\left(\frac{\pi}{2} t^2\right)\,\mathrm{d}t, \\ C(x) &= \int_0^x \cos\left(\frac{\pi}{2} t^2\right)\,\mathrm{d}t. \end{align}" alt="Fresnel integral"> -->
10+
11+
<!-- </equation> -->
12+
13+
Some sources define the Fresnel integrals using t<sup>2</sup> for the argument of the sine and cosine. To get these functions, multiply the computed integrals by `√(π/2)` and multiply the argument `x` by `√(2/π)`.
14+
15+
</section>
16+
17+
<!-- /.intro -->
18+
19+
<section class="usage">
20+
21+
## Usage
22+
23+
```javascript
24+
var fresnel = require( '@stdlib/math/base/special/fresnel' );
25+
```
26+
27+
#### fresnel( \[out,] x )
28+
29+
Simultaneously computes the [Fresnel integrals][fresnel-integral] S(x) and C(x).
30+
31+
```javascript
32+
var v = fresnel( 0.0 );
33+
// returns [ ~0.0, ~0.0 ]
34+
35+
v = fresnel( 1.0 );
36+
// returns [ ~0.438, ~0.780 ]
37+
38+
v = fresnel( Infinity );
39+
// returns [ ~0.5, ~0.5 ]
40+
41+
v = fresnel( -Infinity );
42+
// returns [ ~-0.5, ~-0.5 ]
43+
44+
v = fresnel( NaN );
45+
// returns [ NaN, NaN ]
46+
```
47+
48+
By default, the function returns the S(x) and C(x) as a two-element `array`. To avoid extra memory allocation, the function supports providing an output (destination) object.
49+
50+
```javascript
51+
var out = new Float64Array( 2 );
52+
53+
var v = fresnel( out, 0.0 );
54+
// return <Float64Array>[ ~0.0, ~0.0 ]
55+
56+
var bool = ( v === out );
57+
// returns true
58+
```
59+
60+
</section>
61+
62+
<!-- /.usage -->
63+
64+
<section class="examples">
65+
66+
## Examples
67+
68+
```javascript
69+
var linspace = require( '@stdlib/math/utils/linspace' );
70+
var fresnel = require( '@stdlib/math/base/special/fresnel' );
71+
72+
var x = linspace( 0.0, 10.0, 100 );
73+
var i;
74+
75+
for ( i = 0; i < x.length; i++ ) {
76+
console.log( fresnel( x[ i ] ) );
77+
}
78+
```
79+
80+
</section>
81+
82+
<!-- /.examples -->
83+
84+
<section class="links">
85+
86+
[fresnel-integral]: https://en.wikipedia.org/wiki/Fresnel_integral
87+
88+
</section>
89+
90+
<!-- /.links -->
Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
'use strict';
2+
3+
// MODULES //
4+
5+
var bench = require( '@stdlib/bench' );
6+
var randu = require( '@stdlib/math/base/random/randu' );
7+
var isnan = require( '@stdlib/math/base/assert/is-nan' );
8+
var pkg = require( './../package.json' ).name;
9+
var fresnel = require( './../lib' );
10+
11+
12+
// MAIN //
13+
14+
bench( pkg, function benchmark( b ) {
15+
var x;
16+
var y;
17+
var i;
18+
19+
b.tic();
20+
for ( i = 0; i < b.iterations; i++ ) {
21+
x = ( randu()*20.0 ) - 10.0;
22+
y = fresnel( x );
23+
if ( isnan( y[0] ) || isnan( y[1] ) ) {
24+
b.fail( 'should not return NaN' );
25+
}
26+
}
27+
b.toc();
28+
if ( isnan( y[0] ) || isnan( y[1] ) ) {
29+
b.fail( 'should not return NaN' );
30+
}
31+
b.pass( 'benchmark finished' );
32+
b.end();
33+
});
34+
35+
bench( pkg+'::in-place', function benchmark( b ) {
36+
var x;
37+
var y;
38+
var i;
39+
40+
y = [ 0.0, 0.0 ];
41+
b.tic();
42+
for ( i = 0; i < b.iterations; i++ ) {
43+
x = ( randu()*20.0 ) - 10.0;
44+
fresnel( y, x );
45+
if ( isnan( y[0] ) || isnan( y[1] ) ) {
46+
b.fail( 'should not return NaN' );
47+
}
48+
}
49+
b.toc();
50+
if ( isnan( y[0] ) || isnan( y[1] ) ) {
51+
b.fail( 'should not return NaN' );
52+
}
53+
b.pass( 'benchmark finished' );
54+
b.end();
55+
});
Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
2+
# VARIABLES #
3+
4+
ifndef VERBOSE
5+
QUIET := @
6+
endif
7+
8+
# Specify the path to Cephes:
9+
CEPHES ?=
10+
11+
# Specify a list of Cephes source files:
12+
CEPHES_SRC ?=
13+
14+
# Determine the OS:
15+
#
16+
# [1]: https://en.wikipedia.org/wiki/Uname#Examples
17+
# [2]: http://stackoverflow.com/a/27776822/2225624
18+
OS ?= $(shell uname)
19+
ifneq (, $(findstring MINGW,$(OS)))
20+
OS := WINNT
21+
else
22+
ifneq (, $(findstring MSYS,$(OS)))
23+
OS := WINNT
24+
else
25+
ifneq (, $(findstring CYGWIN,$(OS)))
26+
OS := WINNT
27+
endif
28+
endif
29+
endif
30+
31+
# Define the program used for compiling C source files:
32+
ifdef C_COMPILER
33+
CC := $(C_COMPILER)
34+
else
35+
CC := gcc
36+
endif
37+
38+
# Define the command-line options when compiling C files:
39+
CFLAGS ?= \
40+
-std=c99 \
41+
-O3 \
42+
-Wall \
43+
-pedantic
44+
45+
# Determine whether to generate [position independent code][1]:
46+
#
47+
# [1]: https://gcc.gnu.org/onlinedocs/gcc/Code-Gen-Options.html#Code-Gen-Options
48+
# [2]: http://stackoverflow.com/questions/5311515/gcc-fpic-option
49+
ifeq ($(OS), WINNT)
50+
fPIC ?=
51+
else
52+
fPIC ?= -fPIC
53+
endif
54+
55+
# List of C targets:
56+
c_targets := benchmark.out
57+
58+
59+
# TARGETS #
60+
61+
# Default target.
62+
#
63+
# This target is the default target.
64+
65+
all: $(c_targets)
66+
67+
.PHONY: all
68+
69+
70+
# Compile C source.
71+
#
72+
# This target compiles C source files.
73+
74+
$(c_targets): %.out: %.c
75+
$(QUIET) $(CC) $(CFLAGS) $(fPIC) -o $@ $(CEPHES_SRC) $< -lm
76+
77+
78+
# Run a benchmark.
79+
#
80+
# This target runs a benchmark.
81+
82+
run: $(c_targets)
83+
$(QUIET) ./$<
84+
85+
.PHONY: run
86+
87+
88+
# Perform clean-up.
89+
#
90+
# This target removes generated files.
91+
92+
clean:
93+
$(QUIET) -rm -f *.o *.out
94+
95+
.PHONY: clean
Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
/**
2+
* Benchmark Cephes `fresnl`.
3+
*/
4+
#include <stdlib.h>
5+
#include <stdio.h>
6+
#include <math.h>
7+
#include <sys/time.h>
8+
9+
#define NAME "fresnel"
10+
#define ITERATIONS 1000000
11+
#define REPEATS 3
12+
13+
/**
14+
* Define prototypes for external functions.
15+
*/
16+
extern int fresnl( double x, double *S, double *C );
17+
18+
/**
19+
* Prints the TAP version.
20+
*/
21+
void print_version() {
22+
printf( "TAP version 13\n" );
23+
}
24+
25+
/**
26+
* Prints the TAP summary.
27+
*
28+
* @param total total number of tests
29+
* @param passing total number of passing tests
30+
*/
31+
void print_summary( int total, int passing ) {
32+
printf( "#\n" );
33+
printf( "1..%d\n", total ); // TAP plan
34+
printf( "# total %d\n", total );
35+
printf( "# pass %d\n", passing );
36+
printf( "#\n" );
37+
printf( "# ok\n" );
38+
}
39+
40+
/**
41+
* Prints benchmarks results.
42+
*
43+
* @param elapsed elapsed time in seconds
44+
*/
45+
void print_results( double elapsed ) {
46+
double rate = (double)ITERATIONS / elapsed;
47+
printf( " ---\n" );
48+
printf( " iterations: %d\n", ITERATIONS );
49+
printf( " elapsed: %0.9f\n", elapsed );
50+
printf( " rate: %0.9f\n", rate );
51+
printf( " ...\n" );
52+
}
53+
54+
/**
55+
* Returns a clock time.
56+
*
57+
* @return clock time
58+
*/
59+
double tic() {
60+
struct timeval now;
61+
gettimeofday( &now, NULL );
62+
return (double)now.tv_sec + (double)now.tv_usec/1.0e6;
63+
}
64+
65+
/**
66+
* Generates a random double on the interval [0,1].
67+
*
68+
* @return random double
69+
*/
70+
double rand_double() {
71+
int r = rand();
72+
return (double)r / ( (double)RAND_MAX + 1.0 );
73+
}
74+
75+
/**
76+
* Runs a benchmark.
77+
*
78+
* @return elapsed time in seconds
79+
*/
80+
double benchmark() {
81+
double elapsed;
82+
double x;
83+
double S;
84+
double C;
85+
double t;
86+
int i;
87+
88+
t = tic();
89+
for ( i = 0; i < ITERATIONS; i++ ) {
90+
x = ( 20.0*rand_double() ) - 10.0;
91+
fresnl( x, &S, &C );
92+
if ( S != S || C != C ) {
93+
printf( "should not return NaN\n" );
94+
break;
95+
}
96+
}
97+
elapsed = tic() - t;
98+
if ( S != S || C != C ) {
99+
printf( "should not return NaN\n" );
100+
}
101+
return elapsed;
102+
}
103+
104+
/**
105+
* Main execution sequence.
106+
*/
107+
int main( void ) {
108+
double elapsed;
109+
int i;
110+
111+
// Use the current time to seed the random number generator:
112+
srand( time( NULL ) );
113+
114+
print_version();
115+
for ( i = 0; i < REPEATS; i++ ) {
116+
printf( "# c::cephes::%s\n", NAME );
117+
elapsed = benchmark();
118+
print_results( elapsed );
119+
printf( "ok %d benchmark finished\n", i+1 );
120+
}
121+
print_summary( REPEATS, REPEATS );
122+
}

0 commit comments

Comments
 (0)