source: project/release/4/mathh/tags/2.0.0/mathh.scm @ 13570

Last change on this file since 13570 was 13570, checked in by Kon Lovett, 11 years ago

Release.

File size: 8.4 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.