source: project/release/5/mathh/tags/4.2.0/mathh.scm @ 36246

Last change on this file since 36246 was 36246, checked in by kon, 5 weeks ago

add modf*, rel 4.2.0

File size: 12.5 KB
Line 
1;;;; mathh.scm  -*- Scheme -*-
2;;;; Kon Lovett, Jul '18
3;;;; Peter Bex, Aug '15 (argvector work!)
4;;;; Kon Lovett, Mar '09 - '10, '15, Sep '17
5;;;; Kon Lovett & John Cowen, '07 - '08
6
7;;; Provides access to ISO C math functions in <math.h>
8;;; that are not defined by the Chicken core.
9
10;; Issues
11;;
12;; - Some public names do not match the C name: fpmod (fmod), bessel-j0 (j0), ...
13;;
14;; - Windows Visual C++ 2005 deprecates hypot, j0, j1, jn, y0, y1, yn
15;; & renames _hypot, _j0, etc.
16;;
17;; - Windows does not provide log2, log1p, lgamma, gamma, scalbn,
18;; acosh, asinh, atanh, erf, erfc, signbit, cbrt.
19;; Some are implemented (poorly) here: log2, log1p, erf, erfc, scalbn, signbit.
20;; These are now defunct log2, log1p, erf, erfc, scalbn
21;;
22;; - 'gamma' is deprecated in favor of 'tgamma' but not available
23;; yet on some platforms, so we use 'gamma' for now.
24;;
25;; - http://www.davidegrayson.com/windev/msys2/ : doesn't define the _WIN32
26;; preprocessor macro by default, but it will if you provide the -mwin32 options
27
28#>
29#include <math.h>
30
31#if defined(__sun__) || defined(__sun)
32# include <sunmath.h> /* log2 */
33#endif
34
35#if 0 /*defined(_WIN32)*/
36/* log2, log1p, erf, erfc, scalbn were not originally provided */
37static double log2( double x )
38{
39# ifndef M_LN2
40# define M_LN2 0.693147180559945309417232121458176568
41# endif
42  return log( x ) / M_LN2;
43}
44
45static double log1p( double x )
46{
47  /* very imprecise */
48  return log( 1.0 + x );
49}
50
51/* from a Python implementation by John D. Cook */
52static double erf( double x )
53{
54# define A1 0.254829592
55# define A2 -0.284496736
56# define A3 1.421413741
57# define A4 -1.453152027
58# define A5 1.061405429
59# define P  0.3275911
60
61  double abs_x = fabs( x ); /* handle singage */
62
63  /* per Abramowitz & Stegun 7.1.26 */
64  double t = 1.0 / (1.0 + P * abs_x);
65  double y = 1.0 - (((((A5 * t + A4) * t) + A3) * t + A2) * t + A1) * t * exp( -abs_x * abs_x );
66
67  return _copysign( y, x ); /* handle signage */
68
69# undef P
70# undef A5
71# undef A4
72# undef A3
73# undef A2
74# undef A1
75# undef P
76}
77
78static double erfc( double x )
79{
80  /* very imprecise */
81  return 1.0 - erf( x );
82}
83
84static double scalbn( double x, int n )
85{
86  /* not efficient */
87  return ldexp( x, n );
88}
89#endif
90<#
91
92(module mathh
93
94(;export
95  fpclassify fpclass
96  signbit
97  copysign nextafter
98  fpmod
99  cbrt
100  hypot
101  log10 log2 log1p
102  log-with-base log/base
103  modf modf* frexp
104  ldexp scalbn
105  cosh sinh tanh
106  acosh asinh atanh
107  gamma tgamma lgamma
108  bessel-j0 bessel-j1 bessel-jn
109  bessel-y0 bessel-y1 bessel-yn
110  erf erfc)
111
112(import scheme
113  (chicken base)
114  (chicken flonum)
115  (chicken syntax)
116  (chicken type)
117  (chicken foreign))
118
119;;; Unimplemented Support
120
121(define-inline (%unimplemented-error name)
122  (error name (##core#immutable '"this function is unavailable on this platform")) )
123
124#; ;UNUSED
125(define-syntax define-unimplemented
126  (syntax-rules ()
127    ((define-unimplemented ?name)
128     (define (?name . _) (%unimplemented-error ?name) ) ) ) )
129
130(define-syntax lambda-unimplemented
131  (syntax-rules ()
132    ((lambda-unimplemented ?name)
133     (lambda _ (%unimplemented-error ?name) ) ) ) )
134
135;;;
136
137;; Bessel functions of the 1st kind
138
139(: bessel-j0 (float --> float))
140;
141(define bessel-j0
142  (cond-expand
143    (windows  (foreign-lambda double "_j0" double) )
144    (else     (foreign-lambda double "j0" double) ) ) )
145
146(: bessel-j1 (float --> float))
147;
148(define bessel-j1
149  (cond-expand
150    (windows  (foreign-lambda double "_j1" double) )
151    (else     (foreign-lambda double "j1" double) ) ) )
152
153(: bessel-jn (fixnum float --> float))
154;
155(define bessel-jn
156  (cond-expand
157    (windows  (foreign-lambda double "_jn" int double) )
158    (else     (foreign-lambda double "jn" int double) ) ) )
159
160;; Bessel functions of the 2nd kind
161
162(: bessel-y0 (float --> float))
163;
164(define bessel-y0
165  (cond-expand
166    (windows  (foreign-lambda double "_y0" double) )
167    (else     (foreign-lambda double "y0" double) ) ) )
168
169(: bessel-y1 (float --> float))
170;
171(define bessel-y1
172  (cond-expand
173    (windows  (foreign-lambda double "_y1" double) )
174    (else     (foreign-lambda double "y1" double) ) ) )
175
176(: bessel-yn (fixnum float --> float))
177;
178(define bessel-yn
179  (cond-expand
180    (windows  (foreign-lambda double "_yn" int double) )
181    (else     (foreign-lambda double "yn" int double) ) ) )
182
183;; Error functions
184
185(: erf (float --> float))
186;
187(define erf (foreign-lambda double "erf" double))
188
189(: erfc (float --> float))
190;
191(define erfc (foreign-lambda double "erfc" double))
192
193;; Hyperbolic functions
194
195(: cosh (float --> float))
196;
197(define cosh (foreign-lambda double "cosh" double))
198
199(: sinh (float --> float))
200;
201(define sinh (foreign-lambda double "sinh" double))
202
203(: tanh (float --> float))
204;
205(define tanh (foreign-lambda double "tanh" double))
206
207;; Inverse Hyperbolic functions
208
209(: acosh (float --> float))
210;
211(define acosh
212  (cond-expand
213    (windows  (lambda-unimplemented 'acosh) )
214    (else     (foreign-lambda double "acosh" double)) ) )
215
216(: asinh (float --> float))
217;
218(define asinh
219  (cond-expand
220    (windows  (lambda-unimplemented 'asinh) )
221    (else     (foreign-lambda double "asinh" double)) ) )
222
223(: atanh (float --> float))
224;
225(define atanh
226  (cond-expand
227    (windows  (lambda-unimplemented 'atanh) )
228    (else     (foreign-lambda double "atanh" double)) ) )
229
230;; Euclidean distance function
231
232(: hypot (float float --> float))
233;
234(define hypot (foreign-lambda double "hypot" double double) )
235
236;; Gamma function
237
238(: gamma (float --> float))
239;
240(define gamma
241  (cond-expand
242    (windows  (lambda-unimplemented 'gamma) )
243    (linux    (foreign-lambda double "tgamma" double) )
244    (macosx   (foreign-lambda double "tgamma" double) )
245    (else     (foreign-lambda double "gamma" double) ) ) )
246
247(define tgamma gamma)
248
249;; Ln Abs Gamma function
250
251(: lgamma (float --> float))
252;
253(define lgamma
254  (cond-expand
255    (windows  (lambda-unimplemented 'lgamma) )
256    (else     (foreign-lambda double "lgamma" double) ) ) )
257
258;; Base-10 logarithm
259
260(: log10 (float --> float))
261;
262(define log10 (foreign-lambda double "log10" double))
263
264;; Base-2 logarithm
265
266(: log2 (float --> float))
267;
268(define log2 (foreign-lambda double "log2" double) )
269
270;; Natural logarithm of 1+x accurate for very small x
271
272(: log1p (float --> float))
273;
274(define log1p (foreign-lambda double "log1p" double) )
275
276;; Compute x * 2**n
277
278(: ldexp (float integer --> float))
279;
280(define ldexp (foreign-lambda double "ldexp" double integer))
281
282;; Efficiently compute x * 2**n
283
284(: scalbn (float integer --> float))
285;
286(define scalbn (foreign-lambda double "scalbn" double integer) )
287
288;; Log function for base n
289
290(: *log-with-base (fixnum --> (procedure (float) float)))
291;
292(define (*log-with-base b)
293  (let ((lnb (log b)))
294    (lambda (n)
295      ((foreign-lambda* double ((double x) (double lnb))
296        "C_return( log( x ) / lnb );") n lnb)) ) )
297
298(: log-with-base (fixnum --> (procedure (float) float)))
299;
300(define (log-with-base b)
301        (case b
302          ((2)  log2 )
303    ((10) log10 )
304    (else (*log-with-base b) ) ) )
305
306(define log/base log-with-base)
307
308;; Flonum remainder
309
310(: fpmod (float float --> float))
311;
312(define fpmod (foreign-lambda double "fmod" double double))
313
314;; Return integer & fraction (as multiple values) of a flonum
315
316(: modf (float --> float float))
317;
318(define modf (foreign-primitive ((double x)) "
319  double ipart;
320  double result  = modf( x, &ipart );
321  C_word* values = C_alloc( 2 * C_SIZEOF_FLONUM );
322  C_word value1  = C_flonum( &values, ipart );
323  C_word value2  = C_flonum( &values, result );
324  C_word av[4] = { C_SCHEME_UNDEFINED, C_k, value1, value2 };
325  C_values( 4, av );
326  ") )
327
328(: modf* (float --> integer float))
329;
330(define (modf* x)
331  (let-values (((i f) (modf x)))
332    (values (inexact->exact i) f) ) )
333
334;; Return mantissa & exponent (as multiple values) of a flonum
335
336(: frexp (float --> float fixnum))
337;
338(define frexp (foreign-primitive ((double x)) "
339  int exp;
340  double result = frexp( x, &exp );
341  C_word* values = C_alloc( C_SIZEOF_FLONUM );
342  C_word value1 = C_flonum( &values, result );
343  C_word value2 = C_fix( exp );
344  C_word av[4] = { C_SCHEME_UNDEFINED, C_k, value1, value2 };
345  C_values( 4, av );
346  ") )
347
348;; Returns arg1 with same sign as arg2
349
350(: copysign (float float --> float))
351;
352(define copysign
353  (cond-expand
354    (windows  (foreign-lambda double "_copysign" double double))
355    (else     (foreign-lambda double "copysign" double double)) ) )
356
357;; Increments/decrements arg1 in the direction of arg2
358
359(: nextafter (float float --> float))
360;
361(define nextafter
362  (cond-expand
363    (windows  (foreign-lambda double "_nextafter" double double))
364    (else     (foreign-lambda double "nextafter" double double)) ) )
365
366;; #t when negative, #f otherwise
367
368(: signbit (float --> boolean))
369;
370(define signbit
371  (cond-expand
372    (windows
373      (foreign-lambda* bool ((double n)) "C_return( _copysign( 1.0, n ) < 0 );"))
374    (else
375      (lambda (num)
376        (or
377          (fp= -0.0 num)
378          ((foreign-lambda bool "signbit" double) num) ) ) ) ) )
379
380;; Cube Root
381
382(: cbrt (float --> float))
383;
384(define cbrt
385  (cond-expand
386    (windows  (lambda-unimplemented 'cbrt))
387    (else     (foreign-lambda double "cbrt" double)) ) )
388
389;; Returns a symbol denoting the kind of floating-point number.
390
391(: fpclass (float --> symbol))
392;
393(define fpclass
394  (cond-expand
395    ((and windows (not cygwin))
396      (foreign-lambda* symbol ((double x)) "
397        char *name;
398        switch (_fpclass( x )) {
399        case _FPCLASS_SNAN:
400          name = \"signaling-nan\";
401          break;
402        case _FPCLASS_QNAN:
403          name = \"quiet-nan\";
404          break;
405        case _FPCLASS_NINF:
406          name = \"negative-infinite\";
407          break;
408        case _FPCLASS_NN:
409          name = \"negative-normal\";
410          break;
411        case _FPCLASS_ND:
412          name = \"negative-subnormal\";
413          break;
414        case _FPCLASS_NZ:
415          name = \"negative-zero\";
416          break;
417        case _FPCLASS_PZ:
418          name = \"positive-zero\";
419          break;
420        case _FPCLASS_PD:
421          name = \"positive-subnormal\";
422          break;
423        case _FPCLASS_PN:
424          name = \"positive-normal\";
425          break;
426        case _FPCLASS_PINF:
427          name = \"positive-infinite\";
428          break;
429        default:
430          name = \"unclassified\";
431          break;
432        }
433        C_return( name );") )
434    (else
435      (lambda (num)
436        ;
437        (define C_fpclassify
438          (foreign-lambda* symbol ((double x)) "
439            char *name;
440            switch (fpclassify( x )) {
441            case FP_INFINITE:
442              name = x < 0 ? \"negative-infinite\" : \"positive-infinite\";
443              break;
444            case FP_NAN:
445              /*FIXME A quiet nan can be distinguished by bit inspection*/
446              name = \"signaling-nan\";
447              break;
448            case FP_NORMAL:
449              name = x < 0 ? \"negative-normal\" : \"positive-normal\";
450              break;
451            case FP_SUBNORMAL:
452              name = x < 0 ? \"negative-subnormal\" : \"positive-subnormal\";
453              break;
454            case FP_ZERO:
455              name = signbit( x ) ? \"negative-zero\" : \"positive-zero\";
456              break;
457            default:
458              name = \"unclassified\";
459              break;
460            }
461            C_return( name );"))
462        ;
463        (if (fp= -0.0 num)
464          'negative-zero
465          (C_fpclassify num) ) ) ) ) )
466
467;; Returns a symbol denoting the kind of floating-point number.
468
469(: fpclassify (float --> symbol))
470;
471(define fpclassify
472  (cond-expand
473    ((and windows (not cygwin))
474      (foreign-lambda* symbol ((double x)) "
475        char *name;
476        switch (_fpclass( x )) {
477        case _FPCLASS_SNAN:
478        case _FPCLASS_QNAN:
479          name = \"nan\";
480          break;
481        case _FPCLASS_NINF:
482        case _FPCLASS_PINF:
483          name = \"infinite\";
484          break;
485        case _FPCLASS_NN:
486        case _FPCLASS_PN:
487          name = \"normal\";
488          break;
489        case _FPCLASS_ND:
490        case _FPCLASS_PD:
491          name = \"subnormal\";
492          break;
493        case _FPCLASS_NZ:
494        case _FPCLASS_PZ:
495          name = \"zero\";
496          break;
497        default:
498          name = \"unclassified\";
499          break;
500        }
501        C_return( name );") )
502    (else
503      (foreign-lambda* symbol ((double x)) "
504        char *name;
505        switch (fpclassify( x )) {
506        case FP_INFINITE:
507          name = \"infinite\";
508          break;
509        case FP_NAN:
510          name = \"nan\";
511          break;
512        case FP_NORMAL:
513          name = \"normal\";
514          break;
515        case FP_SUBNORMAL:
516          name = \"subnormal\";
517          break;
518        case FP_ZERO:
519          name = \"zero\";
520          break;
521        default:
522          name = \"unclassified\";
523          break;
524        }
525        C_return( name );") ) ) )
526
527) ;module mathh
Note: See TracBrowser for help on using the repository browser.