1;;; 5_3.ms
2;;; Copyright 1984-2017 Cisco Systems, Inc.
3;;;
4;;; Licensed under the Apache License, Version 2.0 (the "License");
5;;; you may not use this file except in compliance with the License.
6;;; You may obtain a copy of the License at
7;;;
8;;; http://www.apache.org/licenses/LICENSE-2.0
9;;;
10;;; Unless required by applicable law or agreed to in writing, software
11;;; distributed under the License is distributed on an "AS IS" BASIS,
12;;; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13;;; See the License for the specific language governing permissions and
14;;; limitations under the License.
15
16(define <<
17   (case-lambda
18      [(x y)
19       (and (flonum? x)
20            (flonum? y)
21            (if (and (fl= x 0.0) (fl= y 0.0))
22                (fl< (fl/ 1.0 x) (fl/ 1.0 y))
23                (fl< x y)))]
24      [(x y z)
25       (and (<< x y) (<< y z))]))
26
27(mat inexact
28   (== (inexact 0) +0.0)
29   (== (inexact #e+1e-400) +0.0)
30   (== (inexact #e-1e-400) -0.0)
31   (== (inexact #e+1e+400) +inf.0)
32   (== (inexact #e-1e+400) -inf.0)
33   (== (inexact #e+1e-5000) +0.0)
34   (== (inexact #e-1e-5000) -0.0)
35   (== (inexact #e+1e+5000) +inf.0)
36   (== (inexact #e-1e+5000) -inf.0)
37
38  ; make sure inexact rounds to even whenever exactly half way to next
39  ; (assuming 52-bit mantissa + hidden bit)
40  ; ratios
41   (fl= (inexact (+ (ash 1 52) 0/2)) #x10000000000000.0)
42   (fl= (inexact (+ (ash 1 52) 1/2)) #x10000000000000.0)
43   (fl= (inexact (+ (ash 1 52) 2/2)) #x10000000000001.0)
44   (fl= (inexact (+ (ash 1 52) 3/2)) #x10000000000002.0)
45   (fl= (inexact (+ (ash 1 52) 4/2)) #x10000000000002.0)
46   (fl= (inexact (+ (ash 1 52) 5/2)) #x10000000000002.0)
47  ; integers
48   (fl= (inexact (* (+ (ash 1 52) 0/2) 2)) #x20000000000000.0)
49   (fl= (inexact (* (+ (ash 1 52) 1/2) 2)) #x20000000000000.0)
50   (fl= (inexact (* (+ (ash 1 52) 2/2) 2)) #x20000000000002.0)
51   (fl= (inexact (* (+ (ash 1 52) 3/2) 2)) #x20000000000004.0)
52   (fl= (inexact (* (+ (ash 1 52) 4/2) 2)) #x20000000000004.0)
53   (fl= (inexact (* (+ (ash 1 52) 5/2) 2)) #x20000000000004.0)
54   (fl= (inexact (ash (* (+ (ash 1 52) 0/2) 2) 40)) #x200000000000000000000000.0)
55   (fl= (inexact (ash (* (+ (ash 1 52) 1/2) 2) 40)) #x200000000000000000000000.0)
56   (fl= (inexact (ash (* (+ (ash 1 52) 2/2) 2) 40)) #x200000000000020000000000.0)
57   (fl= (inexact (ash (* (+ (ash 1 52) 3/2) 2) 40)) #x200000000000040000000000.0)
58   (fl= (inexact (ash (* (+ (ash 1 52) 4/2) 2) 40)) #x200000000000040000000000.0)
59   (fl= (inexact (ash (* (+ (ash 1 52) 5/2) 2) 40)) #x200000000000040000000000.0)
60  ; make sure inexact rounds up when more than half way to next
61  ; (assuming 52-bit mantissa + hidden bit)
62  ; ratios
63   (fl= (inexact (+ (ash 1 52) 0/2 1/4)) #x10000000000000.0)
64   (fl= (inexact (+ (ash 1 52) 1/2 1/4)) #x10000000000001.0)
65   (fl= (inexact (+ (ash 1 52) 2/2 1/4)) #x10000000000001.0)
66   (fl= (inexact (+ (ash 1 52) 3/2 1/4)) #x10000000000002.0)
67   (fl= (inexact (+ (ash 1 52) 4/2 1/4)) #x10000000000002.0)
68   (fl= (inexact (+ (ash 1 52) 5/2 1/4)) #x10000000000003.0)
69   (fl= (inexact (+ (ash 1 52) 1/2 1/8)) #x10000000000001.0)
70   (fl= (inexact (+ (ash 1 52) 3/2 1/8)) #x10000000000002.0)
71   (fl= (inexact (+ (ash 1 52) 1/2 (expt 2 -80))) #x10000000000001.0)
72   (fl= (inexact (+ (ash 1 52) 3/2 (expt 2 -80))) #x10000000000002.0)
73  ; integers
74   (fl= (inexact (* (+ (ash 1 52) 0/2 1/4) 4)) #x40000000000000.0)
75   (fl= (inexact (* (+ (ash 1 52) 1/2 1/4) 4)) #x40000000000004.0)
76   (fl= (inexact (* (+ (ash 1 52) 2/2 1/4) 4)) #x40000000000004.0)
77   (fl= (inexact (* (+ (ash 1 52) 3/2 1/4) 4)) #x40000000000008.0)
78   (fl= (inexact (* (+ (ash 1 52) 4/2 1/4) 4)) #x40000000000008.0)
79   (fl= (inexact (* (+ (ash 1 52) 5/2 1/4) 4)) #x4000000000000C.0)
80   (fl= (inexact (* (+ (ash 1 52) 1/2 1/8) 8)) #x80000000000008.0)
81   (fl= (inexact (* (+ (ash 1 52) 3/2 1/8) 8)) #x80000000000010.0)
82   (fl= (inexact (* (+ (ash 1 52) 1/2 (expt 2 -80)) (expt 2 80)))
83        #x1000000000000100000000000000000000.0)
84   (fl= (inexact (* (+ (ash 1 52) 3/2 (expt 2 -80)) (expt 2 80)))
85        #x1000000000000200000000000000000000.0)
86   ; verify fix for incorrect input of 2.2250738585072011e-308 reported by leppie
87   ; 2.2250738585072011e-308 falls right on the edge between normalized and denormalized numbers,
88   ; and should not be rounded up to a normalized number
89   (equal?
90     (number->string (string->number "2.2250738585072011e-308"))
91     "2.225073858507201e-308|52")
92   (equal?
93     (decode-float (string->number "2.2250738585072011e-308"))
94     '#(#b1111111111111111111111111111111111111111111111111111 -1074 1))
95   ; similar case in binary...
96   (equal?
97     (decode-float (string->number "#b1.111111111111111111111111111111111111111111111111111011e-1111111111"))
98     '#(#b1111111111111111111111111111111111111111111111111111 -1074 1))
99   (equal?
100     (number->string (string->number "#b1.111111111111111111111111111111111111111111111111111011e-1111111111"))
101     "2.225073858507201e-308|52")
102   ; slightly higher number should be rounded up
103   (equal?
104     (number->string (string->number "2.2250738585072012e-308"))
105     "2.2250738585072014e-308")
106   (equal?
107     (number->string (string->number "#b1.111111111111111111111111111111111111111111111111111100e-1111111111"))
108     "2.2250738585072014e-308")
109)
110
111(mat exact
112   (error? (exact (nan)))
113   (error? (exact +inf.0))
114   (error? (exact -inf.0))
115   (eq? (exact +0.0) 0)
116   (eq? (exact -0.0) 0)
117)
118
119(mat ==
120   (== 1.0 1.0)
121   (== -1.0 -1.0)
122   (not (== -1.0 +1.0))
123   (not (== +1.0 -1.0))
124   (== 0.0 0.0)
125   (== -0.0 -0.0)
126   (not (== -0.0 +0.0))
127   (not (== +0.0 -0.0))
128   (== +inf.0 +inf.0)
129   (== -inf.0 -inf.0)
130   (not (== -inf.0 +inf.0))
131   (not (== +inf.0 -inf.0))
132   (== (nan) (nan))
133   (not (== +inf.0 (nan)))
134   (not (== (nan) -inf.0))
135   (not (== 0.0 0.0-0.0i))
136   (== +e +e)
137   (== -e -e)
138   (not (== +e +0.0))
139   (not (== -e -0.0))
140 )
141
142(mat <<
143   (<< -1.0 1.0)
144   (not (<< +1.0 -1.0))
145   (not (<< 0.0 0.0))
146   (<< -0.0 +0.0)
147   (not (<< +0.0 -0.0))
148   (<< -inf.0 +inf.0)
149   (not (<< +inf.0 -inf.0))
150   (not (<< (nan) (nan)))
151   (not (<< (nan) +0.0))
152   (not (<< +0.0 (nan)))
153   (<< -e +0.0 +e)
154   (<< -e -0.0 +e)
155   (not (<< +e +e))
156   (not (<< -e -e))
157 )
158
159(mat fl=
160   (let ((n (read (open-input-string "+nan.0"))))
161     (not (fl= n n)))
162   (not (fl= (nan)))
163   (not (fl= (nan) +inf.0))
164   (not (fl= (nan) -inf.0))
165   (not (fl= (nan) (nan)))
166   (not (fl= (nan) 0.0))
167   (fl= +inf.0 +inf.0)
168   (fl= -inf.0 -inf.0)
169   (not (fl= -inf.0 +inf.0))
170   (fl= +0.0 -0.0)
171 )
172
173(mat fl<
174   (not (fl< (nan)))
175   (not (fl< (nan) (nan)))
176   (not (fl< (nan) 0.0))
177   (not (fl< 0.0 (nan)))
178   (fl< -inf.0 0.0)
179 )
180
181(mat fl>
182   (not (fl> (nan)))
183   (not (fl> (nan) (nan)))
184   (not (fl> (nan) 0.0))
185   (not (fl> 0.0 (nan)))
186   (fl> +inf.0 -inf.0)
187   (fl> +inf.0 0.0)
188   (not (fl> +0.0 -0.0))
189 )
190
191(mat fl<=
192   (not (fl<= (nan)))
193   (not (fl<= (nan) (nan)))
194   (not (fl<= (nan) 0.0))
195   (not (fl<= 0.0 (nan)))
196 )
197
198(mat fl>=
199   (not (fl>= (nan)))
200   (not (fl>= (nan) (nan)))
201   (not (fl>= (nan) 0.0))
202   (not (fl>= 0.0 (nan)))
203 )
204
205(mat fl-
206   (== (fl- +0.0) -0.0)
207   (== (fl- -0.0) +0.0)
208   (== (fl- +inf.0) -inf.0)
209   (== (fl- -inf.0) +inf.0)
210   (== (fl- (nan)) (nan))
211   (== (fl- -0.0 -0.0) +0.0)
212   (== (fl- +0.0 -0.0) +0.0)
213   (== (fl- -0.0 +0.0) -0.0)
214   (== (fl- +0.0 +0.0) +0.0)
215   (andmap
216     (lambda (a)
217       (andmap
218         (lambda (b)
219           (andmap
220             (lambda (c) (== (fl- a b c) (fl- (fl- a b) c)))
221             '(0.0 -0.0)))
222         '(0.0 -0.0)))
223     '(0.0 -0.0))
224   (let ()
225     (define-syntax ff
226       (syntax-rules ()
227         [(_ k1 k2) (lambda (x) (eqv? (fl- k1 x k2) (fl- (fl- k1 x) k2)))]))
228     (andmap
229       (lambda (p) (and (p +0.0) (p -0.0)))
230       (list (ff +0.0 +0.0) (ff +0.0 -0.0) (ff -0.0 +0.0) (ff -0.0 -0.0))))
231   (error? (fl- 3.0 5.4 'a))
232   (error? (fl- 'a 3.0 5.4))
233   (error? (fl- 3.0 'a 5.4))
234   (== (fl- 5.0 4.0 3.0 2.0) -4.0)
235   (== (fl- 5.0 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0) -2.0)
236   (begin
237     (define ($fl-f x y) (fl- -0.0 x y))
238     (procedure? $fl-f))
239   (== ($fl-f 3.0 4.0) -7.0)
240   (== (fl- 1e30 1e30 7.0) -7.0)
241 )
242
243(mat +
244  ; just in case we're ever tempted to combine nested generic arithmetic operators...
245  (begin
246    (define f1a (lambda (x) (= (+ x 2) (+ (+ x 1) 1))))
247    (define f1b (lambda (x) (= (+ (+ x 1) 1) x)))
248    (define f2 (lambda (x) (= (- (+ x 1e308) 1e308) +inf.0)))
249    #t)
250  (f1a 0)
251  (not (f1a (inexact (expt 2 53))))
252  (not (f1b 0))
253  (f1b (inexact (expt 2 53)))
254  (not (f2 (inexact 0)))
255  (f2 +inf.0)
256  (not (f2 +nan.0))
257  (f2 1e308)
258 )
259
260(mat -
261   (== (- +0.0) -0.0)
262   (== (- -0.0) +0.0)
263   (== (- +inf.0) -inf.0)
264   (== (- -inf.0) +inf.0)
265   (== (- (nan)) (nan))
266   (== (- -0.0 -0.0) +0.0)
267   (== (- +0.0 -0.0) +0.0)
268   (== (- -0.0 +0.0) -0.0)
269   (== (- +0.0 +0.0) +0.0)
270   (andmap
271     (lambda (a)
272       (andmap
273         (lambda (b)
274           (andmap
275             (lambda (c) (== (- a b c) (- (- a b) c)))
276             '(0.0 -0.0)))
277         '(0.0 -0.0)))
278     '(0.0 -0.0))
279   (error? (- 3.0 5.4 'a))
280   (error? (- 'a 3.0 5.4))
281   (error? (- 3.0 'a 5.4))
282   (== (- 1e30 1e30 7.0) -7.0)
283   (begin
284     (define $ieee-foo
285       (lambda (x)
286         (- x 1e30 7.0)))
287     #t)
288   (== ($ieee-foo 1e30) -7.0)
289 )
290
291(mat fl+
292   (== (fl+ -0.0 -0.0) -0.0)
293   (== (fl+ +0.0 -0.0) +0.0)
294   (== (fl+ -0.0 +0.0) +0.0)
295   (== (fl+ +0.0 +0.0) +0.0)
296 )
297
298(mat fl*
299   (== (fl* -1.0 +0.0) -0.0)
300   (== (fl* -1.0 -0.0) +0.0)
301   (== (fl* +1.0 +0.0) +0.0)
302   (== (fl* +1.0 -0.0) -0.0)
303 )
304
305(mat fl/
306   (== (fl/ +0.0) +inf.0)
307   (== (fl/ -0.0) -inf.0)
308   (== (fl/ +inf.0) +0.0)
309   (== (fl/ -inf.0) -0.0)
310   (== (fl/ (nan)) (nan))
311   (== (fl/ +1.0 +0.0) +inf.0)
312   (== (fl/ +1.0 -0.0) -inf.0)
313   (== (fl/ -1.0 +0.0) -inf.0)
314   (== (fl/ -1.0 -0.0) +inf.0)
315   (== (fl/ +0.0 +0.0) (nan))
316   (== (fl/ +0.0 -0.0) (nan))
317   (== (fl/ -0.0 +0.0) (nan))
318   (== (fl/ -0.0 -0.0) (nan))
319   (andmap
320     (lambda (a)
321       (andmap
322         (lambda (b)
323           (andmap
324             (lambda (c) (== (fl/ a b c) (fl/ (fl/ a b) c)))
325             '(1e300 1e250)))
326         '(1e300 1e250)))
327     '(1e300 1e250))
328   (error? (fl/ 3.0 5.4 'a))
329   (error? (fl/ 'a 3.0 5.4))
330   (error? (fl/ 3.0 'a 5.4))
331   (== (fl/ 16.0 2.0 -2.0 2.0) -2.0)
332   (== (fl/ 16.0 2.0 -2.0 2.0 4.0 1.0 -1.0) 0.5)
333   (== (fl/ 1e300 1e300 1e300) 1e-300)
334 )
335
336(mat /
337   (== (/ +0.0) +inf.0)
338   (== (/ -0.0) -inf.0)
339   (== (/ +inf.0) +0.0)
340   (== (/ -inf.0) -0.0)
341   (== (/ (nan)) (nan))
342   (== (/ +1.0 +0.0) +inf.0)
343   (== (/ +1.0 -0.0) -inf.0)
344   (== (/ -1.0 +0.0) -inf.0)
345   (== (/ -1.0 -0.0) +inf.0)
346   (== (/ +0.0 +0.0) (nan))
347   (== (/ +0.0 -0.0) (nan))
348   (== (/ -0.0 +0.0) (nan))
349   (== (/ -0.0 -0.0) (nan))
350   (andmap
351     (lambda (a)
352       (andmap
353         (lambda (b)
354           (andmap
355             (lambda (c) (== (/ a b c) (/ (/ a b) c)))
356             '(1e300 1e250)))
357         '(1e300 1e250)))
358     '(1e300 1e250))
359   (error? (/ 3.0 5.4 'a))
360   (error? (/ 'a 3.0 5.4))
361   (error? (/ 3.0 'a 5.4))
362   (== (fl/ 1e300 1e300 1e300) 1e-300)
363 )
364
365(mat expt
366   (== (expt +0.0 +0.0) +1.0)
367   (== (expt -0.0 +0.0) +1.0)
368   (== (expt +0.0 -0.0) +1.0)
369   (== (expt -0.0 -0.0) +1.0)
370   (== (expt +1.0 +0.0) +1.0)
371   (== (expt -1.0 +0.0) +1.0)
372   (== (expt +0.0 +1.0) +0.0)
373   (== (expt -0.0 +1.0) -0.0)
374   (== (expt -0.0 +2.0) +0.0)
375   (== (expt -0.0 +3.0) -0.0)
376   (== (expt +inf.0 +0.0) +1.0)
377   (== (expt +inf.0 +1.0) +inf.0)
378   (== (expt -inf.0 +0.0) +1.0)
379   (== (expt -inf.0 +1.0) -inf.0)
380   (== (expt +inf.0 +inf.0) +inf.0)
381   (== (expt +inf.0 -inf.0) +0.0)
382   (== (expt -inf.0 +inf.0) +inf.0)
383   (== (expt -inf.0 -inf.0) +0.0)
384   (== (expt +inf.0 +.5) +inf.0)
385   (== (expt (nan) +.5) (nan))
386   (== (expt +.5 (nan)) (nan))
387   (== (expt (nan) (nan)) (nan))
388   (== (expt (nan) +0.0) +1.0)
389   (== (expt +0.0 (nan)) (nan))
390   (== (expt +0.0 (nan)) (nan))
391   (== (expt +inf.0+2i 2) +inf.0+0.0i)
392   (== (let ([n (expt 2 32)]) (expt 2 (make-rectangular n n))) -inf.0+inf.0i)
393 )
394
395(mat magnitude
396   (== (magnitude -0.0) 0.0)
397   (== (magnitude 0.0) 0.0)
398   (== (magnitude 0.0-0.0i) 0.0)
399   (== (magnitude -1.0) 1.0)
400   (== (magnitude 1.0) 1.0)
401   (== (magnitude 0.0+1.0i) 1.0)
402   (== (magnitude +inf.0) +inf.0)
403   (== (magnitude -inf.0) +inf.0)
404   (== (magnitude +inf.0+inf.0i) +inf.0)
405   (== (magnitude +inf.0+2.0i) +inf.0)
406   (== (magnitude +2.0+inf.0i) +inf.0)
407   (== (magnitude (nan)) (nan))
408   (== (magnitude (make-rectangular (nan) (nan))) (nan))
409   (== (magnitude (make-rectangular +0.0 (nan))) (nan))
410   (<< +0.0 (magnitude (make-rectangular +e +e)))
411   (<< +0.0 (magnitude (make-rectangular -e -e)))
412 )
413
414(mat sqrt
415   ; from Kahan
416   (== (sqrt -0.0) +0.0+0.0i)
417   (== (sqrt -4.0) +0.0+2.0i)
418   (== (sqrt -inf.0) +0.0+inf.0i)
419   (== (sqrt 0.0+inf.0i) +inf.0+inf.0i)
420   (== (sqrt 4.0+inf.0i) +inf.0+inf.0i)
421   (== (sqrt +inf.0+inf.0i) +inf.0+inf.0i)
422   (== (sqrt -0.0+inf.0i) +inf.0+inf.0i)
423   (== (sqrt -4.0+inf.0i) +inf.0+inf.0i)
424   (== (sqrt -inf.0+inf.0i) +inf.0+inf.0i)
425   (== (sqrt 0.0-inf.0i) +inf.0-inf.0i)
426   (== (sqrt 4.0-inf.0i) +inf.0-inf.0i)
427   (== (sqrt +inf.0-inf.0i) +inf.0-inf.0i)
428   (== (sqrt -0.0-inf.0i) +inf.0-inf.0i)
429   (== (sqrt -4.0-inf.0i) +inf.0-inf.0i)
430   (== (sqrt -inf.0-inf.0i) +inf.0-inf.0i)
431   (== (sqrt (make-rectangular (nan) +0.0)) (make-rectangular (nan)(nan)))
432   (== (sqrt (make-rectangular 0.0 (nan))) (make-rectangular (nan) (nan)))
433   (== (sqrt (make-rectangular (nan) (nan))) (make-rectangular (nan) (nan)))
434   (== (sqrt +inf.0+0.0i) +inf.0+0.0i)
435   (== (sqrt +inf.0+4.0i) +inf.0+0.0i)
436   (== (sqrt +inf.0-0.0i) +inf.0-0.0i)
437   (== (sqrt +inf.0-4.0i) +inf.0-0.0i)
438   (== (sqrt (make-rectangular +inf.0 (nan))) (make-rectangular +inf.0 (nan)))
439   (== (sqrt -inf.0+0.0i) +0.0+inf.0i)
440   (== (sqrt -inf.0+4.0i) +0.0+inf.0i)
441   (== (sqrt -inf.0-0.0i) +0.0-inf.0i)
442   (== (sqrt -inf.0-4.0i) +0.0-inf.0i)
443   (let ([z (sqrt (make-rectangular -inf.0 (nan)))])
444      (and (== (real-part z) (nan)) (== (abs (imag-part z)) +inf.0)))
445   ; others
446   (== (sqrt +0.0) +0.0)
447   (== (sqrt +1.0) +1.0)
448   (== (sqrt +4.0) +2.0)
449   (== (sqrt +inf.0) +inf.0)
450   (== (sqrt +0.0+0.0i) +0.0+0.0i)
451   (== (sqrt +1.0+0.0i) +1.0+0.0i)
452   (== (sqrt +4.0+0.0i) +2.0+0.0i)
453   (== (sqrt +inf.0+0.0i) +inf.0+0.0i)
454   (== (sqrt -0.0+0.0i) +0.0+0.0i)
455   (== (sqrt -1.0+0.0i) +0.0+1.0i)
456   (== (sqrt -4.0+0.0i) +0.0+2.0i)
457   (== (sqrt -inf.0+0.0i) +0.0+inf.0i)
458   (== (sqrt -0.0-0.0i) +0.0-0.0i)
459   (== (sqrt -1.0-0.0i) +0.0-1.0i)
460   (== (sqrt -inf.0-0.0i) +0.0-inf.0i)
461   (== (sqrt +0.0-0.0i) +0.0-0.0i)
462   (== (sqrt +1.0-0.0i) +1.0-0.0i)
463   (== (sqrt +inf.0-0.0i) +inf.0-0.0i)
464   (== (sqrt (nan)) (nan))
465 )
466
467(mat exp
468   (== (exp +0.0) +1.0)
469   (== (exp -0.0) +1.0)
470   (== (exp +inf.0) +inf.0)
471   (== (exp -inf.0) +0.0)
472   (== (exp (nan)) (nan))
473   (== (exp +0.0+0.0i) +1.0+0.0i)
474   (== (exp -0.0-0.0i) +1.0-0.0i)
475  ; if exp treats x+0.0i the same as x:
476   (== (exp +inf.0+0.0i) +inf.0+0.0i)
477  ; otherwise:
478   #;(== (exp +inf.0+0.0i) +inf.0+nan.0i)
479   (== (exp +inf.0-0.0i) +inf.0-0.0i)
480   (== (exp -inf.0+0.0i) 0.0+0.0i)
481   (== (exp -inf.0-0.0i) 0.0-0.0i)
482  ; if exp treats x+0.0i the same as x:
483   (== (exp (make-rectangular (nan) +0.0)) (make-rectangular (nan) +0.0))
484  ; otherwise:
485   #;(== (exp (make-rectangular (nan) +0.0)) (make-rectangular (nan) (nan)))
486  ; if exp treats x+0.0i the same as x:
487   (== (exp (make-rectangular (nan) -0.0)) (make-rectangular (nan) -0.0))
488  ; otherwise:
489   #;(== (exp (make-rectangular (nan) -0.0)) (make-rectangular (nan) (nan)))
490   (~= (exp 700.0+.75i) 7.421023049046266e303+6.913398801654868e303i)
491   (~= (exp 700.0-.75i) 7.421023049046266e303-6.913398801654868e303i)
492   (== (exp 800.0+.75i) +inf.0+inf.0i)
493   (== (exp 800.0-.75i) +inf.0-inf.0i)
494   (== (exp 800.0+1e-200i) +inf.0+2.7263745721125063e147i)
495   (== (exp 800.0-1e-200i) +inf.0-2.7263745721125063e147i)
496   (== (exp +inf.0+1.0i) +inf.0+inf.0i)
497   (== (exp +inf.0+2.0i) -inf.0+inf.0i)
498   (== (exp +inf.0+3.0i) -inf.0+inf.0i)
499   (== (exp +inf.0+4.0i) -inf.0-inf.0i)
500   (== (exp +inf.0+123.0i) -inf.0-inf.0i)
501 )
502
503(mat log
504   (== (log 0.0) -inf.0)
505   (== (log 1.0) 0.0)
506   (== (log +inf.0) +inf.0)
507
508   (== (log -0.0) (make-rectangular -inf.0 +pi))
509   (== (log -1.0) (make-rectangular 0.0 +pi))
510   (== (log -inf.0) (make-rectangular +inf.0 +pi))
511
512   (== (log +1.0i) (make-rectangular 0.0 +pi/2))
513   (== (log -1.0i) (make-rectangular 0.0 -pi/2))
514
515   (== (log -0.0+0.0i) (make-rectangular -inf.0 +pi))
516   (== (log -0.0-0.0i) (make-rectangular -inf.0 -pi))
517   (== (log +0.0+0.0i) -inf.0+0.0i)
518   (== (log +0.0-0.0i) -inf.0-0.0i)
519
520   (== (log +1.0+0.0i) 0.0+0.0i)
521   (== (log -1.0+0.0i) (make-rectangular 0.0 +pi))
522   (== (log +1.0-0.0i) 0.0-0.0i)
523   (== (log -1.0-0.0i) (make-rectangular 0.0 -pi))
524 )
525
526(mat fllog
527   (== (log 0.0) -inf.0)
528   (== (log 1.0) 0.0)
529   (== (log +inf.0) +inf.0)
530
531   (== (log -0.0) (make-rectangular -inf.0 +pi))
532   (== (log -1.0) (make-rectangular 0.0 +pi))
533   (== (log -inf.0) (make-rectangular +inf.0 +pi))
534
535   (== (log +1.0i) (make-rectangular 0.0 +pi/2))
536   (== (log -1.0i) (make-rectangular 0.0 -pi/2))
537
538   (== (log -0.0+0.0i) (make-rectangular -inf.0 +pi))
539   (== (log -0.0-0.0i) (make-rectangular -inf.0 -pi))
540   (== (log +0.0+0.0i) -inf.0+0.0i)
541   (== (log +0.0-0.0i) -inf.0-0.0i)
542
543   (== (log +1.0+0.0i) 0.0+0.0i)
544   (== (log -1.0+0.0i) (make-rectangular 0.0 +pi))
545   (== (log +1.0-0.0i) 0.0-0.0i)
546   (== (log -1.0-0.0i) (make-rectangular 0.0 -pi))
547)
548
549(mat sin
550   (== (sin +0.0) +0.0)
551   (== (sin -0.0) -0.0)
552   (== (sin +inf.0) (nan))
553   (== (sin -inf.0) (nan))
554   (== (sin (nan)) (nan))
555 )
556
557(mat cos
558   (== (cos +0.0) +1.0)
559   (== (cos -0.0) +1.0)
560   (== (cos +inf.0) (nan))
561   (== (cos -inf.0) (nan))
562   (== (cos (nan)) (nan))
563 )
564
565(mat tan
566   (== (tan +0.0) +0.0)
567   (== (tan -0.0) -0.0)
568   (== (tan +inf.0) (nan))
569   (== (tan -inf.0) (nan))
570   (== (tan (nan)) (nan))
571   (== (tan -0.0+0.0i) -0.0+0.0i)
572 )
573
574(mat asin
575   (== (asin +0.0) +0.0)
576   (== (asin -0.0) -0.0)
577   (== (asin +1.0) +pi/2)
578   (== (asin -1.0) -pi/2)
579   (== (asin (nan)) (nan))
580   (== (asin -0.0+0.0i) -0.0+0.0i)
581 )
582
583(mat acos
584   (== (acos +1.0) +0.0)
585   (== (acos -1.0) +pi)
586   (== (acos +0.0) +pi/2)
587   (== (acos -0.0) +pi/2)
588   (== (acos (nan)) (nan))
589 )
590
591(mat atan
592   ; cases from Steele (CLtL)
593   (== (atan +0.0 +e) +0.0)
594   (== (atan +0.0 +inf.0) +0.0)
595   (<< +0.0 (atan +e +e) +pi/2)
596   (<< +0.0 (atan +inf.0 +inf.0) +pi/2)
597   (== (atan +e +0.0) +pi/2)
598   (== (atan +inf.0 +0.0) +pi/2)
599   (== (atan +e -0.0) +pi/2)
600   (== (atan +inf.0 -0.0) +pi/2)
601   (<< +pi/2 (atan +e -e) +pi)
602   (<< +pi/2 (atan +inf.0 -inf.0) +pi)
603   (== (atan +0.0 -e) +pi)
604   (== (atan +0.0 -inf.0) +pi)
605   (== (atan -0.0 -e) -pi)             ; Steele erroneously says +pi
606   (== (atan -0.0 -inf.0) -pi)         ; Steele erroneously says +pi
607   (<< -pi (atan -e -e) -pi/2)
608   (<< -pi (atan -inf.0 -inf.0) -pi/2)
609   (== (atan -e +0.0) -pi/2)
610   (== (atan -e -0.0) -pi/2)
611   (== (atan -inf.0 +0.0) -pi/2)
612   (== (atan -inf.0 -0.0) -pi/2)
613   (<< -pi/2 (atan -e +e) -0.0)
614   (<< -pi/2 (atan -inf.0 +inf.0) -0.0)
615   (== (atan -0.0 +e) -0.0)
616   (== (atan -0.0 +inf.0) -0.0)
617   (== (atan +0.0 +0.0) +0.0)
618   (== (atan -0.0 +0.0) -0.0)
619   (== (atan +0.0 -0.0) +pi)
620   (== (atan -0.0 -0.0) -pi)
621
622   (== (atan -inf.0) -pi/2)
623   (== (atan +inf.0) +pi/2)
624   (== (atan +0.0) +0.0)
625   (== (atan -0.0) -0.0)
626   (if (memq (machine-type) '(i3qnx ti3qnx))
627       (~= (atan +1.0) +pi/4)
628       (== (atan +1.0) +pi/4))
629   (== (atan -1.0) -pi/4)
630   (== (atan (nan)) (nan))
631   (== (atan -0.0+0.0i) -0.0+0.0i)
632)
633
634(mat sinh
635   (== (sinh 0.0) 0.0)
636   (== (sinh -0.0) -0.0)
637   (== (sinh +inf.0) +inf.0)
638   (== (sinh -inf.0) -inf.0)
639   (== (sinh (nan)) (nan))
640   (== (sinh -0.0+0.0i) -0.0+0.0i)
641 )
642
643(mat cosh
644   (== (cosh 0.0) 1.0)
645   (== (cosh -0.0) 1.0)
646   (== (cosh +inf.0) +inf.0)
647   (== (cosh -inf.0) +inf.0)
648   (== (cosh (nan)) (nan))
649 )
650
651(mat tanh
652   (== (tanh 0.0) 0.0)
653   (== (tanh -0.0) -0.0)
654   (== (tanh +inf.0) +1.0)
655   (== (tanh -inf.0) -1.0)
656   (== (tanh (nan)) (nan))
657   (== (tanh -0.0+0.0i) -0.0+0.0i)
658 )
659
660(mat asinh
661   (== (asinh 0.0) 0.0)
662   (== (asinh -0.0) -0.0)
663   (== (asinh +inf.0) +inf.0)
664   (== (asinh -inf.0) -inf.0)
665   (== (asinh (nan)) (nan))
666   (== (asinh -0.0+0.0i) -0.0+0.0i)
667 )
668
669(mat acosh
670   (== (acosh 1.0) 0.0)
671   (== (acosh +inf.0) +inf.0)
672   (== (acosh (nan)) (nan))
673 )
674
675(mat atanh
676   (== (atanh 0.0) 0.0)
677   (== (atanh -0.0) -0.0)
678   (== (atanh +1.0) +inf.0)
679   (== (atanh -1.0) -inf.0)
680   (== (atanh (nan)) (nan))
681   (== (atanh -0.0+0.0i) -0.0+0.0i)
682   (== (atanh -0.0+0.0i) -0.0+0.0i)
683 )
684
685(mat flonum->fixnum
686   (error? (flonum->fixnum +inf.0))
687   (error? (flonum->fixnum -inf.0))
688   (error? (flonum->fixnum (nan)))
689   (eq? (flonum->fixnum -0.0) 0)
690 )
691
692(mat fllp
693  (error? (fllp 3))
694  (eqv? (fllp 0.0) 0)
695  (eqv? (fllp 1.0) 2046)
696  (eqv? (fllp -1.0) 2046)
697  (eqv? (fllp 1.5) 2047)
698  (eqv? (fllp -1.5) 2047)
699  (and (memv (fllp +nan.0) '(4094 4095)) #t)
700  (eqv? (fllp +inf.0) 4094)
701  (eqv? (fllp -inf.0) 4094)
702  (eqv?
703    (fllp #b1.1111111111111111111111111111111111111111111111111111e1111111111)
704    4093)
705  (eqv? (fllp #b1.0e-1111111110) 2)
706  (or (eqv? #b.1e-1111111110 0.0)
707      (eqv? (fllp #b.1e-1111111110) 1))
708  (eqv? (fllp #b.01e-1111111110) 0)
709 )
710
711(mat fp-output
712  (equal? (number->string 1e23) "1e23")
713  (equal? (number->string 4.450147717014403e-308) "4.450147717014403e-308")
714  (equal? (number->string 1.1665795231290236e-302) "1.1665795231290236e-302")
715 ; fp printing algorithm always rounds up on ties
716  (equal? (number->string 3.6954879760742188e-6) "3.6954879760742188e-6")
717  (equal? (number->string 5.629499534213123e14) "5.629499534213123e14")
718 )
719
720(mat string->number
721  (equal? (string->number "+1e-400") +0.0)
722  (equal? (string->number "-1e-400") -0.0)
723  (equal? (string->number "+1e+400") +inf.0)
724  (equal? (string->number "-1e+400") -inf.0)
725  (equal? (string->number "+1e-5000") +0.0)
726  (equal? (string->number "-1e-5000") -0.0)
727  (equal? (string->number "+1e+5000") +inf.0)
728  (equal? (string->number "-1e+5000") -inf.0)
729  (equal? (string->number "+1e-50000") +0.0)
730  (equal? (string->number "-1e-50000") -0.0)
731  (equal? (string->number "+1e+50000") +inf.0)
732  (equal? (string->number "-1e+50000") -inf.0)
733  (equal? (string->number "+1e-500000") +0.0)
734  (equal? (string->number "-1e-500000") -0.0)
735  (equal? (string->number "+1e+500000") +inf.0)
736  (equal? (string->number "-1e+500000") -inf.0)
737
738  (equal? (string->number "#b1e-10000110010") 5e-324)
739  (equal? (string->number "5e-324") 5e-324)
740  (equal? (string->number "#b-1e-10000110010") -5e-324)
741  (equal? (string->number "-5e-324") -5e-324)
742  (equal? (string->number "#b1e-10000110010") (inexact (* 5 (expt 10 -324))))
743  (equal? (string->number "5e-324") (inexact (* 5 (expt 10 -324))))
744  (equal? (string->number "#b-1e-10000110010") (inexact (* -5 (expt 10 -324))))
745  (equal? (string->number "-5e-324") (inexact (* -5 (expt 10 -324))))
746
747  (if (memq (machine-type) '(a6nt ta6nt)) ; tolerably inaccurate
748      (equal? (string->number "#b1e-10000110100") 0.0)
749      (equal? (string->number "#b1e-10000110011") 0.0))
750  (if (memq (machine-type) '(a6nt ta6nt)) ; tolerably inaccurate
751      (equal? (string->number "#b-1e-10000110100") -0.0)
752      (equal? (string->number "#b-1e-10000110011") -0.0))
753  (equal? (string->number "5e-325") 0.0)
754  (equal? (string->number "-5e-325") -0.0)
755
756  (equal? (string->number "1.7976931348623157e308") 1.7976931348623157e308)
757  (equal? (string->number "-1.7976931348623157e308") -1.7976931348623157e308)
758  (equal? (string->number "#b1.1111111111111111111111111111111111111111111111111111e1111111111") 1.7976931348623157e308)
759  (equal? (string->number "#b-1.1111111111111111111111111111111111111111111111111111e1111111111") -1.7976931348623157e308)
760  (equal? (string->number "1.7976931348623157e308") (inexact (* 9007199254740991 (expt 2 971))))
761  (equal? (string->number "-1.7976931348623157e308") (inexact (* -9007199254740991 (expt 2 971))))
762  (equal? (string->number "#b1.1111111111111111111111111111111111111111111111111111e1111111111") (inexact (* 9007199254740991 (expt 2 971))))
763  (equal? (string->number "#b-1.1111111111111111111111111111111111111111111111111111e1111111111") (inexact (* -9007199254740991 (expt 2 971))))
764
765  (equal? (string->number "#b1.11111111111111111111111111111111111111111111111111111e1111111111") +inf.0)
766  (equal? (string->number "#b-1.11111111111111111111111111111111111111111111111111111e1111111111") -inf.0)
767  (equal? (string->number "1.7976931348623159e308") +inf.0)
768  (equal? (string->number "-1.7976931348623159e308") -inf.0)
769
770  (equal? (string->number "1e100000000000000000000") +inf.0)
771  (equal? (string->number "-1e100000000000000000000") -inf.0)
772  (equal? (string->number "1e-100000000000000000000") 0.0)
773  (equal? (string->number "-1e-100000000000000000000") -0.0)
774)
775