1 | ;;;; mathh.scm |
---|
2 | ;;;; Kon Lovett & John Cowen, '07 - '08 |
---|
3 | ;;;; Kon Lovett, Mar '09 |
---|
4 | |
---|
5 | ;;; Provides access to ISO C math functions in <math.h> |
---|
6 | ;;; that are not defined by the Chicken core. |
---|
7 | |
---|
8 | ;; Issues |
---|
9 | ;; |
---|
10 | ;; - Windows Visual C++ 2005 deprecates hypot, j0, j1, jn, y0, y1, yn |
---|
11 | ;; & renames _hypot, _j0, etc. |
---|
12 | ;; |
---|
13 | ;; - Windows does not provide log2, log1p, lgamma, tgamma. scalbn. |
---|
14 | ;; |
---|
15 | ;; - 'gamma' is deprecated in favor of 'tgamma' but not available |
---|
16 | ;; yet on some platforms, so we use 'gamma' for now. |
---|
17 | ;; |
---|
18 | ;; - Solaris log2 in <sunmath.h>. |
---|
19 | |
---|
20 | |
---|
21 | ;;; Prelude |
---|
22 | |
---|
23 | (declare |
---|
24 | (usual-integrations) |
---|
25 | (inline) |
---|
26 | (local) |
---|
27 | (number-type generic) |
---|
28 | (no-procedure-checks) |
---|
29 | (no-bound-checks) |
---|
30 | (bound-to-procedure |
---|
31 | ##sys#check-inexact ) ) |
---|
32 | |
---|
33 | #> |
---|
34 | #include <math.h> |
---|
35 | #include <float.h> |
---|
36 | <# |
---|
37 | |
---|
38 | |
---|
39 | ;;; Unimplemented Support |
---|
40 | |
---|
41 | (define-inline (%unimplemented-error name) |
---|
42 | (error name (##core#immutable '"this function is unavailable on this platform")) ) |
---|
43 | |
---|
44 | |
---|
45 | ;;; Mathh |
---|
46 | |
---|
47 | (module mathh (;export |
---|
48 | bessel-j0 bessel-j1 bessel-jn |
---|
49 | bessel-y0 bessel-y1 bessel-yn |
---|
50 | cosh sinh tanh |
---|
51 | hypot |
---|
52 | gamma tgamma lgamma |
---|
53 | log10 log2 make-log/base log1p |
---|
54 | fpmod modf ldexp scalbn frexp |
---|
55 | fpclassify |
---|
56 | fpclass) |
---|
57 | |
---|
58 | (import scheme chicken) |
---|
59 | |
---|
60 | ;;; Unimplemented Support |
---|
61 | |
---|
62 | #; ;UNUSED |
---|
63 | (define-syntax define-unimplemented |
---|
64 | (syntax-rules () |
---|
65 | [(_ ?name) |
---|
66 | (define (?name . _) (%unimplemented-error ?name) ) ] ) ) |
---|
67 | |
---|
68 | (define-syntax lambda-unimplemented |
---|
69 | (syntax-rules () |
---|
70 | [(_ ?name) |
---|
71 | (lambda _ (%unimplemented-error ?name) ) ] ) ) |
---|
72 | |
---|
73 | ;;; |
---|
74 | |
---|
75 | ;; Bessel functions of the 1st kind |
---|
76 | |
---|
77 | (define bessel-j0 (foreign-lambda double "j0" double)) |
---|
78 | (define bessel-j1 (foreign-lambda double "j1" double)) |
---|
79 | (define bessel-jn (foreign-lambda double "jn" int double)) |
---|
80 | |
---|
81 | ;; Bessel functions of the 2nd kind |
---|
82 | |
---|
83 | (define bessel-y0 (foreign-lambda double "y0" double)) |
---|
84 | (define bessel-y1 (foreign-lambda double "y1" double)) |
---|
85 | (define bessel-yn (foreign-lambda double "yn" int double)) |
---|
86 | |
---|
87 | ;; Hyperbolic functions |
---|
88 | |
---|
89 | (define cosh (foreign-lambda double "cosh" double)) |
---|
90 | (define sinh (foreign-lambda double "sinh" double)) |
---|
91 | (define tanh (foreign-lambda double "tanh" double)) |
---|
92 | |
---|
93 | ;; Euclidean distance function |
---|
94 | |
---|
95 | (define hypot (foreign-lambda double "hypot" double double)) |
---|
96 | |
---|
97 | ;; Gamma function |
---|
98 | |
---|
99 | (define gamma |
---|
100 | (cond-expand |
---|
101 | [windows (lambda-unimplemented 'gamma) ] |
---|
102 | [linux (foreign-lambda double "tgamma" double) ] |
---|
103 | [macosx (foreign-lambda double "tgamma" double) ] |
---|
104 | [else (foreign-lambda double "gamma" double) ] ) ) |
---|
105 | |
---|
106 | (define tgamma gamma) |
---|
107 | |
---|
108 | ;; Ln Abs Gamma function |
---|
109 | |
---|
110 | (define lgamma |
---|
111 | (cond-expand |
---|
112 | [windows (lambda-unimplemented 'lgamma) ] |
---|
113 | [else (foreign-lambda double "lgamma" double) ] ) ) |
---|
114 | |
---|
115 | ;; Base-10 logarithm |
---|
116 | |
---|
117 | (define log10 (foreign-lambda double "log10" double)) |
---|
118 | |
---|
119 | ;; Base-2 logarithm |
---|
120 | |
---|
121 | (define log2 |
---|
122 | (cond-expand |
---|
123 | [(or linux macosx bsd) |
---|
124 | (foreign-lambda double "log2" double) ] |
---|
125 | [else |
---|
126 | (foreign-lambda* double ((double x)) " |
---|
127 | #ifndef M_LN2 |
---|
128 | #define #define M_LN2 0.693147180559945309417232121458176568 |
---|
129 | #endif |
---|
130 | C_return( log( x ) / M_LN2 ); |
---|
131 | ") |
---|
132 | #; |
---|
133 | (foreign-lambda* double ((double x)) " |
---|
134 | #ifndef M_PI |
---|
135 | # define M_PI 3.14159265358979323846264338327950288 |
---|
136 | #endif |
---|
137 | #ifndef M_E |
---|
138 | # define M_E 2.71828182845904523536028747135266250 |
---|
139 | #endif |
---|
140 | C_return( (log( 2.0 * M_PI * x) / 2.0) + (x * log( x / M_E )) ); |
---|
141 | ") ] ) ) |
---|
142 | |
---|
143 | ;; Natural logarithm of 1+x accurate for very small x |
---|
144 | |
---|
145 | (define log1p |
---|
146 | (cond-expand |
---|
147 | [windows ; potentially inaccurate but ... |
---|
148 | (foreign-lambda* double ((double x)) "C_return( log( 1.0 + x ) );") ] |
---|
149 | [else |
---|
150 | (foreign-lambda double "log1p" double) ] ) ) |
---|
151 | |
---|
152 | ;; Compute x * 2**n |
---|
153 | |
---|
154 | (define ldexp (foreign-lambda double "ldexp" double integer)) |
---|
155 | |
---|
156 | ;; Efficiently compute x * 2**n |
---|
157 | |
---|
158 | (define scalbn |
---|
159 | (cond-expand |
---|
160 | [windows ; not efficient but ... |
---|
161 | (foreign-lambda* double ((double x) (integer n)) "C_return( ldexp( x, n ) );")] |
---|
162 | [else |
---|
163 | (foreign-lambda double "scalbn" double integer)]) ) |
---|
164 | |
---|
165 | ;; Log function for base n |
---|
166 | |
---|
167 | (define (make-log/base b) |
---|
168 | (when (fixnum? b) (set! b (exact->inexact b))) |
---|
169 | (##sys#check-inexact b 'make-log/base) |
---|
170 | (cond [(fp= 2.0 b) log2 ] |
---|
171 | [(fp= 10.0 b) log10 ] |
---|
172 | [else |
---|
173 | (let ([lnb (log b)]) |
---|
174 | (lambda (n) |
---|
175 | ((foreign-lambda* double ((double x) (double lnb)) "C_return( log( x ) / lnb );") n lnb)) ) ] ) ) |
---|
176 | |
---|
177 | ;; Flonum remainder |
---|
178 | |
---|
179 | (define fpmod (foreign-lambda double "fmod" double double)) |
---|
180 | |
---|
181 | ;; Return integer & fraction (as multiple values) of a flonum |
---|
182 | |
---|
183 | (define modf (foreign-primitive ((double x)) " |
---|
184 | double ipart; |
---|
185 | double result = modf( x, &ipart ); |
---|
186 | C_word* values = C_alloc( 2 * C_SIZEOF_FLONUM ); |
---|
187 | C_word value1 = C_flonum( &values, ipart ); |
---|
188 | C_word value2 = C_flonum( &values, result ); |
---|
189 | C_values( 4, C_SCHEME_UNDEFINED, C_k, value1, value2 ); |
---|
190 | ") ) |
---|
191 | |
---|
192 | ;; Return mantissa & exponent (as multiple values) of a flonum |
---|
193 | |
---|
194 | (define frexp (foreign-primitive ((double x)) " |
---|
195 | int exp; |
---|
196 | double result = frexp( x, &exp ); |
---|
197 | C_word* values = C_alloc( C_SIZEOF_FLONUM ); |
---|
198 | C_word value1 = C_flonum( &values, result ); |
---|
199 | C_word value2 = C_fix( exp ); |
---|
200 | C_values( 4, C_SCHEME_UNDEFINED, C_k, value1, value2 ); |
---|
201 | ") ) |
---|
202 | |
---|
203 | ;; Returns a symbol denoting the kind of floating-point number. |
---|
204 | |
---|
205 | (define fpclass |
---|
206 | (cond-expand |
---|
207 | [(and windows (not cygwin)) |
---|
208 | (foreign-lambda* symbol ((double x)) " |
---|
209 | char *name; |
---|
210 | switch (_fpclass( x )) { |
---|
211 | case _FPCLASS_SNAN: |
---|
212 | name = \"signaling-nan\"; |
---|
213 | break; |
---|
214 | case _FPCLASS_QNAN: |
---|
215 | name = \"quiet-nan\"; |
---|
216 | break; |
---|
217 | case _FPCLASS_NINF: |
---|
218 | name = \"negative-infinite\"; |
---|
219 | break; |
---|
220 | case _FPCLASS_NN: |
---|
221 | name = \"negative-normal\"; |
---|
222 | break; |
---|
223 | case _FPCLASS_ND: |
---|
224 | name = \"negative-subnormal\"; |
---|
225 | break; |
---|
226 | case _FPCLASS_NZ: |
---|
227 | name = \"negative-zero\"; |
---|
228 | break; |
---|
229 | case _FPCLASS_PZ: |
---|
230 | name = \"positive-zero\"; |
---|
231 | break; |
---|
232 | case _FPCLASS_PD: |
---|
233 | name = \"positive-subnormal\"; |
---|
234 | break; |
---|
235 | case _FPCLASS_PN: |
---|
236 | name = \"positive-normal\"; |
---|
237 | break; |
---|
238 | case _FPCLASS_PINF: |
---|
239 | name = \"positive-infinite\"; |
---|
240 | break; |
---|
241 | default: |
---|
242 | name = \"unclassified\"; |
---|
243 | break; |
---|
244 | } |
---|
245 | C_return( name ); |
---|
246 | ") ] |
---|
247 | [else |
---|
248 | (foreign-lambda* symbol ((double x)) " |
---|
249 | char *name; |
---|
250 | switch (fpclassify( x )) { |
---|
251 | case FP_INFINITE: |
---|
252 | name = x < 0 ? \"negative-infinite\" : \"positive-infinite\"; |
---|
253 | break; |
---|
254 | case FP_NAN: |
---|
255 | /*FIXME A quiet nan can be distinguished by bit inspection*/ |
---|
256 | name = \"signaling-nan\"; |
---|
257 | break; |
---|
258 | case FP_NORMAL: |
---|
259 | name = x < 0 ? \"negative-normal\" : \"positive-normal\"; |
---|
260 | break; |
---|
261 | case FP_SUBNORMAL: |
---|
262 | name = x < 0 ? \"negative-subnormal\" : \"positive-subnormal\"; |
---|
263 | break; |
---|
264 | case FP_ZERO: |
---|
265 | name = x == -0.0 ? \"negative-zero\" : \"positive-zero\"; |
---|
266 | break; |
---|
267 | default: |
---|
268 | name = \"unclassified\"; |
---|
269 | break; |
---|
270 | } |
---|
271 | C_return( name ); |
---|
272 | ") ] ) ) |
---|
273 | |
---|
274 | ;; Returns a symbol denoting the kind of floating-point number. |
---|
275 | |
---|
276 | (define fpclassify |
---|
277 | (cond-expand |
---|
278 | [(and windows (not cygwin)) |
---|
279 | (foreign-lambda* symbol ((double x)) " |
---|
280 | char *name; |
---|
281 | switch (_fpclass( x )) { |
---|
282 | case _FPCLASS_SNAN: |
---|
283 | case _FPCLASS_QNAN: |
---|
284 | name = \"nan\"; |
---|
285 | break; |
---|
286 | case _FPCLASS_NINF: |
---|
287 | case _FPCLASS_PINF: |
---|
288 | name = \"infinite\"; |
---|
289 | break; |
---|
290 | case _FPCLASS_NN: |
---|
291 | case _FPCLASS_PN: |
---|
292 | name = \"normal\"; |
---|
293 | break; |
---|
294 | case _FPCLASS_ND: |
---|
295 | case _FPCLASS_PD: |
---|
296 | name = \"subnormal\"; |
---|
297 | break; |
---|
298 | case _FPCLASS_NZ: |
---|
299 | case _FPCLASS_PZ: |
---|
300 | name = \"zero\"; |
---|
301 | break; |
---|
302 | default: |
---|
303 | name = \"unclassified\"; |
---|
304 | break; |
---|
305 | } |
---|
306 | C_return( name ); |
---|
307 | ") ] |
---|
308 | [else |
---|
309 | (foreign-lambda* symbol ((double x)) " |
---|
310 | char *name; |
---|
311 | switch (fpclassify( x )) { |
---|
312 | case FP_INFINITE: |
---|
313 | name = \"infinite\"; |
---|
314 | break; |
---|
315 | case FP_NAN: |
---|
316 | name = \"nan\"; |
---|
317 | break; |
---|
318 | case FP_NORMAL: |
---|
319 | name = \"normal\"; |
---|
320 | break; |
---|
321 | case FP_SUBNORMAL: |
---|
322 | name = \"subnormal\"; |
---|
323 | break; |
---|
324 | case FP_ZERO: |
---|
325 | name = \"zero\"; |
---|
326 | break; |
---|
327 | default: |
---|
328 | name = \"unclassified\"; |
---|
329 | break; |
---|
330 | } |
---|
331 | C_return( name ); |
---|
332 | ") ] ) ) |
---|
333 | |
---|
334 | ) ;module mathh |
---|