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