| ... |
... |
@@ -24,6 +24,18 @@ |
|
24
|
24
|
(defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
|
|
25
|
25
|
;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
|
|
26
|
26
|
|
|
|
27
|
+(defconstant +log-2d0+
|
|
|
28
|
+ (log 2d0)
|
|
|
29
|
+ "Natural log of 2d0")
|
|
|
30
|
+
|
|
|
31
|
+(defconstant +log-2w0+
|
|
|
32
|
+ (log 2w0)
|
|
|
33
|
+ "Natural log of 2w0")
|
|
|
34
|
+
|
|
|
35
|
+(defconstant +log2-e+
|
|
|
36
|
+ (log (exp 1d0) 2d0)
|
|
|
37
|
+ "Base 2 logarithm of e")
|
|
|
38
|
+
|
|
27
|
39
|
;;; Make these INLINE, since the call to C is at least as compact as a Lisp
|
|
28
|
40
|
;;; call, and saves number consing to boot.
|
|
29
|
41
|
;;;
|
| ... |
... |
@@ -488,14 +500,14 @@ |
|
488
|
500
|
complex)
|
|
489
|
501
|
(if (and (zerop base) (plusp (realpart power)))
|
|
490
|
502
|
(* base power)
|
|
491
|
|
- (exp (* power (* (log2 base 1w0) (log 2w0))))))
|
|
|
503
|
+ (exp (* power (* (log2 base 1w0) +log-2w0+)))))
|
|
492
|
504
|
(((foreach fixnum (or bignum ratio) single-float double-float)
|
|
493
|
505
|
(foreach (complex double-float)))
|
|
494
|
506
|
;; Result should have double-float accuracy. Use log2 in
|
|
495
|
507
|
;; case the base won't fit in a double-float.
|
|
496
|
508
|
(if (and (zerop base) (plusp (realpart power)))
|
|
497
|
509
|
(* base power)
|
|
498
|
|
- (exp (* power (* (log2 base) (log 2d0))))))
|
|
|
510
|
+ (exp (* power (* (log2 base) +log-2d0+)))))
|
|
499
|
511
|
((double-float
|
|
500
|
512
|
(foreach (complex rational) (complex single-float)))
|
|
501
|
513
|
(if (and (zerop base) (plusp (realpart power)))
|
| ... |
... |
@@ -508,7 +520,7 @@ |
|
508
|
520
|
;; in case the base won't fit in a double-float.
|
|
509
|
521
|
(if (and (zerop base) (plusp (realpart power)))
|
|
510
|
522
|
(* base power)
|
|
511
|
|
- (exp (* power (* (log2 base 1w0) (log 2w0))))))
|
|
|
523
|
+ (exp (* power (* (log2 base 1w0) +log-2w0+)))))
|
|
512
|
524
|
(((foreach fixnum (or bignum ratio) single-float)
|
|
513
|
525
|
(foreach (complex rational) (complex single-float)))
|
|
514
|
526
|
(if (and (zerop base) (plusp (realpart power)))
|
| ... |
... |
@@ -566,18 +578,18 @@ |
|
566
|
578
|
;; log(2), with the precision specified by the type of F
|
|
567
|
579
|
(number-dispatch ((f real))
|
|
568
|
580
|
((double-float)
|
|
569
|
|
- #.(log 2d0))
|
|
|
581
|
+ +log-2d0+)
|
|
570
|
582
|
#+double-double
|
|
571
|
583
|
((double-double-float)
|
|
572
|
|
- #.(log 2w0))))
|
|
|
584
|
+ +log-2w0+)))
|
|
573
|
585
|
(log-2-pi (f)
|
|
574
|
586
|
;; log(pi), with the precision specified by the type of F
|
|
575
|
587
|
(number-dispatch ((f real))
|
|
576
|
588
|
((double-float)
|
|
577
|
|
- #.(/ pi (log 2d0)))
|
|
|
589
|
+ #.(/ pi +log-2d0+))
|
|
578
|
590
|
#+double-double
|
|
579
|
591
|
((double-double-float)
|
|
580
|
|
- #.(/ dd-pi (log 2w0)))))
|
|
|
592
|
+ #.(/ dd-pi +log-2w0+))))
|
|
581
|
593
|
(log1p (x)
|
|
582
|
594
|
;; log(1+x), with the precision specified by the type of
|
|
583
|
595
|
;; X
|
| ... |
... |
@@ -827,7 +839,7 @@ |
|
827
|
839
|
(if (minusp number)
|
|
828
|
840
|
(complex (coerce (log (- number)) 'single-float)
|
|
829
|
841
|
(coerce pi 'single-float))
|
|
830
|
|
- (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float)))
|
|
|
842
|
+ (coerce (/ (log2 number) +log2-e+) 'single-float)))
|
|
831
|
843
|
((ratio)
|
|
832
|
844
|
(if (minusp number)
|
|
833
|
845
|
(complex (coerce (log (- number)) 'single-float)
|
| ... |
... |
@@ -852,7 +864,7 @@ |
|
852
|
864
|
(coerce (%log1p (coerce (- number 1) 'double-float))
|
|
853
|
865
|
'single-float)
|
|
854
|
866
|
(coerce (/ (- (log2 top) (log2 bot))
|
|
855
|
|
- #.(log (exp 1d0) 2d0))
|
|
|
867
|
+ +log2-e+)
|
|
856
|
868
|
'single-float)))))
|
|
857
|
869
|
(((foreach single-float double-float))
|
|
858
|
870
|
;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
|
| ... |
... |
@@ -1393,7 +1405,7 @@ This is for use with J /= 0 only when |z| is huge." |
|
1393
|
1405
|
(let ((t0 #.(/ 1 (sqrt 2.0d0)))
|
|
1394
|
1406
|
(t1 1.2d0)
|
|
1395
|
1407
|
(t2 3d0)
|
|
1396
|
|
- (ln2 #.(log 2d0))
|
|
|
1408
|
+ (ln2 +log-2d0+)
|
|
1397
|
1409
|
(x (float (realpart z) 1.0d0))
|
|
1398
|
1410
|
(y (float (imagpart z) 1.0d0)))
|
|
1399
|
1411
|
(multiple-value-bind (rho k)
|