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

Last change on this file since 13657 was 13657, 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  (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
Note: See TracBrowser for help on using the repository browser.