source: project/release/3/mathh/trunk/mathh.scm @ 8510

Last change on this file since 8510 was 8510, checked in by Kon Lovett, 12 years ago

Save w/ canonical dir form.

File size: 4.5 KB
Line 
1;;;; mathh.scm
2;;;; Provides access to ISO C math functions in <math.h>
3;;;; that are not defined by Chicken
4;;;; Public domain
5
6;; Issues
7;;
8;; - Windows Visual C++ 2005 deprecates
9;; hypot, j0, j1, jn, y0, y1, yn
10;; & renames
11;; _hypot, _j0, etc.
12;;
13;; - Windows does not provide
14;; log2, log1p, lgamma, tgamma. scalbn
15;;
16;; - 'gamma' is deprecated in favor of 'tgamma' but not available
17;; yet on some platforms, so we use 'gamma' for now.
18;;
19;; - Solaris log2 in sunmath.h
20
21;;; Scheme-level declarations
22(declare
23        (usual-integrations)
24        (no-procedure-checks-for-usual-bindings)
25        (unused
26          unimplemented-error)
27        (bound-to-procedure
28                bessel-j0 bessel-j1 bessel-jn
29                bessel-y0 bessel-y1 bessel-yn
30                cosh sinh tanh
31                hypot
32                gamma tgamma lgamma
33                log10 log2 log1p
34                fpmod modf ldexp scalbn frexp)
35        (export
36                bessel-j0 bessel-j1 bessel-jn
37                bessel-y0 bessel-y1 bessel-yn
38                cosh sinh tanh
39                hypot
40                gamma tgamma lgamma
41                log10 log2 make-log/base log1p
42                fpmod modf ldexp scalbn frexp))
43
44;;; Include proper C header
45#>
46#include <math.h>
47#include <float.h>
48<#
49
50;;; Unimplemented Support
51
52;from posixwin.scm
53
54(define (unimplemented-error name)
55  (error name (##core#immutable '"this function is not available on this platform")) )
56
57(define-macro (define-unimplemented name)
58  `(define (,name . _) (unimplemented-error ',name)) )
59
60(define-macro (lambda-unimplemented name)
61  `(lambda _ (unimplemented-error ',name)) )
62
63;;; Bessel functions of the 1st kind
64(define bessel-j0 (foreign-lambda double j0 double))
65(define bessel-j1 (foreign-lambda double j1 double))
66(define bessel-jn (foreign-lambda double jn int double))
67
68;;; Bessel functions of the 2nd kind
69(define bessel-y0 (foreign-lambda double y0 double))
70(define bessel-y1 (foreign-lambda double y1 double))
71(define bessel-yn (foreign-lambda double yn int double))
72
73;;; Hyperbolic functions
74(define cosh (foreign-lambda double cosh double))
75(define sinh (foreign-lambda double sinh double))
76(define tanh (foreign-lambda double tanh double))
77
78;;; Euclidean distance function
79(define hypot (foreign-lambda double hypot double double))
80
81;;; Gamma function
82(define gamma
83  (cond-expand
84    [windows
85      (lambda-unimplemented 'gamma)]
86    [linux
87      (foreign-lambda double "tgamma" double)]
88    [macosx
89      (foreign-lambda double "tgamma" double)]
90    [else
91      (foreign-lambda double "gamma" double)]))
92
93(define tgamma gamma)
94
95;;; Ln Abs Gamma function
96(define lgamma
97  (cond-expand
98    [windows
99      (lambda-unimplemented 'lgamma)]
100    [else
101      (foreign-lambda double lgamma double)]))
102
103;;; Base-10 logarithm
104(define log10 (foreign-lambda double log10 double))
105
106;;; Base-2 logarithm
107(define log2
108  (cond-expand
109    [(or linux macosx bsd)
110      (foreign-lambda double "log2" double)]
111    [else
112      (foreign-lambda* double ((double x)) "return (log(x) / M_LN2);")]))
113
114;; Log function for base n
115(define (make-log/base b)
116        (when (fixnum? b)
117                (set! b (exact->inexact b)))
118        (unless (flonum? b)
119                (error 'make-log/base "invalid flonum" b))
120        (cond
121                [(= 2.0 b)
122                        log2]
123                [(= 10.0 b)
124                        log10]
125                [else
126                        (let ([lnb (log b)]
127                                          [logf
128                                            (foreign-lambda* double ((double x) (double lnb))
129                                              "return (log(x) / lnb);")])
130                                (lambda (n)
131                                        (logf n lnb)))]) )
132
133;;; Natural logarithm of 1+x accurate for very small x
134(define log1p
135  (cond-expand
136    [windows ; potentially inaccurate but ...
137      (foreign-lambda* double ((double x))
138        "return (log(1.0 + x));")]
139    [else
140                  (foreign-lambda double log1p double)]))
141
142;;; Flonum remainder
143(define fpmod (foreign-lambda double fmod double double))
144
145;;; Return integer, fraction (as multiple values) of a flonum
146(define modf (foreign-primitive ((double x))
147        "double ipart;
148         double result  = modf(x, &ipart);
149         C_word* values = C_alloc(2 * C_SIZEOF_FLONUM);
150         C_word value1  = C_flonum(&values, ipart);
151         C_word value2  = C_flonum(&values, result);
152         C_values(4, C_SCHEME_UNDEFINED, C_k, value1, value2);
153"))
154
155;;; Compute x * 2**n
156(define ldexp (foreign-lambda double ldexp double integer))
157
158;;; Efficiently compute x * 2**n
159(define scalbn
160  (cond-expand
161    [windows ; not efficient but ...
162      (foreign-lambda* double ((double x) (integer n))
163        "return (ldexp(x, n));")]
164    [else
165      (foreign-lambda double scalbn double integer)]) )
166
167;;; Return mantissa, exponent (as multiple values) of a flonum
168(define frexp (foreign-primitive ((double x))
169        "int exp;
170         double result = frexp(x, &exp);
171         C_word* values = C_alloc(C_SIZEOF_FLONUM);
172         C_word value1 = C_flonum(&values, result);
173         C_word value2 = C_fix(exp);
174         C_values(4, C_SCHEME_UNDEFINED, C_k, value1, value2);
175"))
Note: See TracBrowser for help on using the repository browser.