source: project/release/4/mathh/trunk/mathh.scm @ 14004

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

Added routines

File size: 9.1 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  signbit copysign nextafter
55  cbrt
56  fpclassify
57  fpclass)
58
59(import scheme chicken)
60
61;;; Unimplemented Support
62
63#; ;UNUSED
64(define-syntax define-unimplemented
65  (syntax-rules ()
66    ((_ ?name)
67     (define (?name . _) (%unimplemented-error ?name) ) ) ) )
68
69(define-syntax lambda-unimplemented
70  (syntax-rules ()
71    ((_ ?name)
72     (lambda _ (%unimplemented-error ?name) ) ) ) )
73
74;;;
75
76;; Bessel functions of the 1st kind
77
78(define bessel-j0 (foreign-lambda double "j0" double))
79(define bessel-j1 (foreign-lambda double "j1" double))
80(define bessel-jn (foreign-lambda double "jn" int double))
81
82;; Bessel functions of the 2nd kind
83
84(define bessel-y0 (foreign-lambda double "y0" double))
85(define bessel-y1 (foreign-lambda double "y1" double))
86(define bessel-yn (foreign-lambda double "yn" int double))
87
88;; Hyperbolic functions
89
90(define cosh (foreign-lambda double "cosh" double))
91(define sinh (foreign-lambda double "sinh" double))
92(define tanh (foreign-lambda double "tanh" double))
93
94;; Euclidean distance function
95
96(define hypot (foreign-lambda double "hypot" double double))
97
98;; Gamma function
99
100(define gamma
101  (cond-expand
102    (windows  (lambda-unimplemented 'gamma) )
103    (linux    (foreign-lambda double "tgamma" double) )
104    (macosx   (foreign-lambda double "tgamma" double) )
105    (else     (foreign-lambda double "gamma" double) ) ) )
106
107(define tgamma gamma)
108
109;; Ln Abs Gamma function
110
111(define lgamma
112  (cond-expand
113    (windows  (lambda-unimplemented 'lgamma) )
114    (else     (foreign-lambda double "lgamma" double) ) ) )
115
116;; Base-10 logarithm
117
118(define log10 (foreign-lambda double "log10" double))
119
120;; Base-2 logarithm
121
122(define log2
123  (cond-expand
124    ((or linux macosx bsd)
125      (foreign-lambda double "log2" double) )
126    (else
127      (foreign-lambda* double ((double x)) "
128        #ifndef M_LN2
129        #define #define M_LN2 0.693147180559945309417232121458176568
130        #endif
131        return( log( x ) / M_LN2 );")
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        return( (log( 2.0 * M_PI * x) / 2.0) + (x * log( x / M_E )) );") ) ) )
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)) "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)) "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)) "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 arg1 with same sign as arg2
203
204(define copysign
205  (cond-expand
206    (windows (foreign-lambda double "_copysign" double double))
207    (else (foreign-lambda double "copysign" double double)) ) )
208
209;; Increments/decrements arg1 in the direction of arg2
210
211(define nextafter
212  (cond-expand
213    (windows (foreign-lambda double "_nextafter" double double))
214    (else (foreign-lambda double "nextafter" double double)) ) )
215
216;; #t when negative, #f otherwise
217
218(define signbit
219  (cond-expand
220    (windows (foreign-lambda* bool ((double n)) "return( _copysign( 1.0, n ) < 0 );"))
221    (else (foreign-lambda bool "signbit" double)) ) )
222
223;; Cube Root
224
225(define cbrt
226  (cond-expand
227    (windows (lambda-unimplemented 'cbrt))
228    (else (foreign-lambda double "cbrt" double)) ) )
229
230;; Returns a symbol denoting the kind of floating-point number.
231
232(define fpclass
233  (cond-expand
234    ((and windows (not cygwin))
235      (foreign-lambda* symbol ((double x)) "
236        char *name;
237        switch (_fpclass( x )) {
238        case _FPCLASS_SNAN:
239          name = \"signaling-nan\";
240          break;
241        case _FPCLASS_QNAN:
242          name = \"quiet-nan\";
243          break;
244        case _FPCLASS_NINF:
245          name = \"negative-infinite\";
246          break;
247        case _FPCLASS_NN:
248          name = \"negative-normal\";
249          break;
250        case _FPCLASS_ND:
251          name = \"negative-subnormal\";
252          break;
253        case _FPCLASS_NZ:
254          name = \"negative-zero\";
255          break;
256        case _FPCLASS_PZ:
257          name = \"positive-zero\";
258          break;
259        case _FPCLASS_PD:
260          name = \"positive-subnormal\";
261          break;
262        case _FPCLASS_PN:
263          name = \"positive-normal\";
264          break;
265        case _FPCLASS_PINF:
266          name = \"positive-infinite\";
267          break;
268        default:
269          name = \"unclassified\";
270          break;
271        }
272        return( name );") )
273    (else
274      (foreign-lambda* symbol ((double x)) "
275        char *name;
276        switch (fpclassify( x )) {
277        case FP_INFINITE:
278          name = x < 0 ? \"negative-infinite\" : \"positive-infinite\";
279          break;
280        case FP_NAN:
281          /*FIXME A quiet nan can be distinguished by bit inspection*/
282          name = \"signaling-nan\";
283          break;
284        case FP_NORMAL:
285          name = x < 0 ? \"negative-normal\" : \"positive-normal\";
286          break;
287        case FP_SUBNORMAL:
288          name = x < 0 ? \"negative-subnormal\" : \"positive-subnormal\";
289          break;
290        case FP_ZERO:
291          name = signbit( x ) ? \"negative-zero\" : \"positive-zero\";
292          break;
293        default:
294          name = \"unclassified\";
295          break;
296        }
297        return( name );") ) ) )
298
299;; Returns a symbol denoting the kind of floating-point number.
300
301(define fpclassify
302  (cond-expand
303    ((and windows (not cygwin))
304      (foreign-lambda* symbol ((double x)) "
305        char *name;
306        switch (_fpclass( x )) {
307        case _FPCLASS_SNAN:
308        case _FPCLASS_QNAN:
309          name = \"nan\";
310          break;
311        case _FPCLASS_NINF:
312        case _FPCLASS_PINF:
313          name = \"infinite\";
314          break;
315        case _FPCLASS_NN:
316        case _FPCLASS_PN:
317          name = \"normal\";
318          break;
319        case _FPCLASS_ND:
320        case _FPCLASS_PD:
321          name = \"subnormal\";
322          break;
323        case _FPCLASS_NZ:
324        case _FPCLASS_PZ:
325          name = \"zero\";
326          break;
327        default:
328          name = \"unclassified\";
329          break;
330        }
331        return( name );") )
332    (else
333      (foreign-lambda* symbol ((double x)) "
334        char *name;
335        switch (fpclassify( x )) {
336        case FP_INFINITE:
337          name = \"infinite\";
338          break;
339        case FP_NAN:
340          name = \"nan\";
341          break;
342        case FP_NORMAL:
343          name = \"normal\";
344          break;
345        case FP_SUBNORMAL:
346          name = \"subnormal\";
347          break;
348        case FP_ZERO:
349          name = \"zero\";
350          break;
351        default:
352          name = \"unclassified\";
353          break;
354        }
355        return( name );") ) ) )
356
357) ;module mathh
Note: See TracBrowser for help on using the repository browser.