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

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

Save.

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