Changeset 15627


Ignore:
Timestamp:
Jan 31, 2013, 5:03:29 PM (6 years ago)
Author:
gb
Message:

Use (large) parts of dfindlay's "reference implementation" to fix
problems noted in ticket:869 in the trunk.

todo: check other issues (e.g., sin/cos/tab of ratios, etc.)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/source/level-0/l0-float.lisp

    r15553 r15627  
    2323  (require "NUMBER-MACROS")
    2424  (require :number-case-macro)
     25  (defconstant single-float-pi (coerce pi 'single-float))
     26  (defconstant double-float-half-pi (/ pi 2))
     27  (defconstant single-float-half-pi (coerce (/ pi 2) 'single-float))
     28  (defconstant single-float-log2 0.6931472)                ; (log 2)
     29  (defconstant double-float-log2 0.6931471805599453d0)     ; (log 2.0d0)
     30  (defconstant double-float-log2^23 15.942385152878742d0)  ; (log (expt 2 23))
    2531)
    2632
     
    729735
    730736
    731 
     737;;; Multiply by i (with additional optional scale factor)
     738;;; Does the "right thing" with minus zeroes (see CLTL2)
     739(defun i* (number &optional (scale 1))
     740  (complex (* (- scale) (imagpart number))
     741           (* scale (realpart number))))
     742
     743;;; complex atanh
     744(defun %complex-atanh (z)
     745  (let* ((rx (realpart z))
     746         (ix (imagpart z))
     747         (sign (typecase rx
     748                 (double-float (%double-float-sign rx))
     749                 (short-float (%short-float-sign rx))
     750                 (t (minusp rx))))
     751         (x rx)
     752         (y ix)
     753         (y1 (abs y))
     754         ra ia)
     755    ;; following code requires non-negative x
     756    (when sign
     757      (setf x (- x))
     758      (setf y (- y)))
     759    (cond ((> (max x y1) 1.8014399e+16)
     760           ;; large value escape
     761           (setq ra (if (> x y1)
     762                      (let ((r (/ y x)))
     763                        (/ (/ x) (1+ (* r r))))
     764                      (let ((r (/ x y)))
     765                        (/ (/ r y) (1+ (* r r))))))
     766           (setq ia (typecase y
     767                      (double-float (float-sign y double-float-half-pi))
     768                      (single-float (float-sign y single-float-half-pi))
     769                      (t (if (minusp y) #.(- single-float-half-pi) single-float-half-pi)))))
     770          ((= x 1)
     771           (setq ra (if (< y1 1e-9)
     772                      (/ (log-e (/ 2 y1)) 2)
     773                      (/ (log1+ (/ 4 (* y y))) 4)))
     774           (setq ia (/ (atan (/ 2 y) -1) 2)))
     775          (t
     776           (let ((r2 (+ (* x x) (* y y))))
     777             (setq ra (/ (log1+ (/ (* -4 x) (1+ (+ (* 2 x) r2)))) -4))
     778             (setq ia (/ (atan (* 2 y) (- 1 r2)) 2)))))
     779    ;; fixup signs, with special case for real arguments
     780    (cond (sign
     781           (setq ra (- ra))
     782           (when (typep z 'complex)
     783             (setq ia (- ia))))
     784          (t
     785           (unless (typep z 'complex)
     786             (setq ia (- ia)))))
     787    (complex ra ia)))
    732788
    733789(defun atan (y &optional (x nil x-p))
    734790  "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
    735   (if x-p
    736     (if (or (typep x 'double-float)
    737             (typep y 'double-float))
    738       (with-stack-double-floats ((dy y)
    739                                  (dx x))
    740         (%df-atan2 dy dx))
    741       #+32-bit-target
    742       (target::with-stack-short-floats ((sy y)
    743                                 (sx x))
    744         (%sf-atan2! sy sx))
    745       #+64-bit-target
    746       (%sf-atan2 (%short-float y) (%short-float x)))
    747     (if (typep y 'complex)
    748       (let ((iy (complex (- (imagpart y)) (realpart y))))
    749         (/ (- (log (+ 1 iy)) (log (- 1 iy)))
    750            #c(0 2)))
    751       (if (typep y 'double-float)
    752         (%double-float-atan! y (%make-dfloat))
    753         #+32-bit-target
    754         (target::with-stack-short-floats ((sy y))
    755           (%single-float-atan! sy (%make-sfloat)))
    756         #+64-bit-target
    757         (%single-float-atan (%short-float y))
    758         ))))
     791  (cond (x-p
     792         (cond ((or (typep x 'double-float)
     793                    (typep y 'double-float))
     794                (with-stack-double-floats ((dy y)
     795                                           (dx x))
     796                  (%df-atan2 dy dx)))
     797               (t
     798                (when (and (rationalp x) (rationalp y))
     799                  ;; rescale arguments so that the maximum absolute value is 1
     800                  (let ((x1 (abs x)) (y1 (abs y)))
     801                    (cond ((> y1 x1)
     802                           (setf x (/ x y1))
     803                           (setf y (signum y)))
     804                          ((not (zerop x))
     805                           (setf y (/ y x1))
     806                           (setf x (signum x))))))
     807                #+32-bit-target
     808                (target::with-stack-short-floats ((sy y)
     809                                                  (sx x))
     810                  (%sf-atan2! sy sx))
     811                #+64-bit-target
     812                (%sf-atan2 (%short-float y) (%short-float x)))))
     813        ((typep y 'double-float)
     814         (%double-float-atan! y (%make-dfloat)))
     815        ((typep y 'single-float)
     816         #+32-bit-target
     817         (%single-float-atan! y (%make-sfloat))
     818         #+64-bit-target
     819         (%single-float-atan y))
     820        ((typep y 'rational)
     821         (cond ((<= (abs y) most-positive-short-float)
     822                #+32-bit-target
     823                (target::with-stack-short-floats ((sy y))
     824                  (%single-float-atan! sy (%make-sfloat)))
     825                #+64-bit-target
     826                (%single-float-atan (%short-float y)))
     827               ((minusp y)
     828                #.(- single-float-half-pi))
     829               (t
     830                single-float-half-pi)))
     831        (t
     832         (let ((r (realpart y))
     833               (i (imagpart y)))
     834           (if (zerop i)
     835             (complex (atan r) i)
     836             (i* (%complex-atanh (complex (- i) r)) -1))))))
    759837
    760838
     
    766844      (if (zerop x)
    767845        (report-bad-arg x '(not (satisfies zerop) ))
    768         (if (floatp x) (float 0.0d0 x) 0))
    769       (/ (log-e x) (log-e b)))
     846        ;; ** CORRECT THE CONTAGION for complex args **
     847        (+ 0 (* x b)))
     848      ;; do the float/rational contagion before the division
     849      ;; but take care with negative zeroes
     850      (let ((x1 (realpart x))
     851            (b1 (realpart b)))
     852        (if (and (typep x1 'float)
     853                 (typep b1 'float))
     854          (/ (log-e (* x (float 1.0 b1)))
     855             (log-e (* b (float 1.0 x1))))
     856          (let ((r (/ (cond ((typep x 'rational)
     857                             (%rational-log x 1.0d0))
     858                            ((typep x1 'rational)
     859                             (%rational-complex-log x 1.0d0))
     860                            (t
     861                             (log-e (* x 1.0d0))))
     862                      (cond ((typep b 'rational)
     863                             (%rational-log b 1.0d0))
     864                            ((typep b1 'rational)
     865                             (%rational-complex-log b 1.0d0))
     866                            (t
     867                             (log-e (* b 1.0d0)))))))
     868            (cond ((or (typep x1 'double-float)
     869                       (typep b1 'double-float))
     870                   r)
     871                  ((complexp r)
     872                   (complex (%short-float (realpart r))
     873                            (%short-float (imagpart r))))
     874                  (t
     875                   (%short-float r)))))))
    770876    (log-e x)))
    771877
     878
     879
    772880(defun log-e (x)
    773   (cond
    774     ((bignump x)
    775      (if (minusp x)
    776        (complex (log-e (- x)) pi)
    777        (let* ((base1 3)
    778               (guess (floor (1- (integer-length x))
    779                             (log base1 2)))
    780               (guess1 (* guess (log-e base1))))
    781          (+ guess1 (log-e (/ x (expt base1 guess)))))))
    782     ((and (ratiop x) 
    783           (or (> x most-positive-short-float)
    784               (< x most-negative-short-float)))
    785      (- (log-e (%numerator x)) (log-e (%denominator x))))
     881   (cond
     882     ((typep x 'double-float)
     883      (if (%double-float-sign x)
     884        (with-stack-double-floats ((dx x))
     885          (complex (%double-float-log! (%%double-float-abs! dx dx) (%make-dfloat)) pi))
     886        (%double-float-log! x (%make-dfloat))))
     887    ((typep x 'short-float)
     888     #+32-bit-target
     889     (if (%short-float-sign x)
     890       (target::with-stack-short-floats ((sx x))
     891         (complex (%single-float-log! (%%short-float-abs! sx sx) (%make-sfloat))
     892                  single-float-pi))
     893       (%single-float-log! x (%make-sfloat)))
     894     #+64-bit-target
     895     (if (%short-float-sign x)
     896       (complex (%single-float-log (%short-float-abs (%short-float x)))
     897                single-float-pi)
     898       (%single-float-log (%short-float x))))
    786899    ((typep x 'complex)
    787      (complex (log-e (abs x)) (phase x)))
    788     ((typep x 'double-float)
    789      (with-stack-double-floats ((dx x))
    790        (if (minusp x)
    791          (complex (%double-float-log! (%%double-float-abs! dx dx) (%make-dfloat)) pi)
    792          (%double-float-log! dx (%make-dfloat)))))
     900     (if (typep (realpart x) 'rational)
     901       (%rational-complex-log x 1.0s0)
     902       ;; take care that intermediate results do not overflow/underflow:
     903       ;; pre-scale argument by an appropriate power of two;
     904       ;; we only need to scale for very large and very small values -
     905       ;;  hence the various 'magic' numbers (values may not be too
     906       ;;  critical but do depend on the sqrt update to fix abs's operation)
     907       (let ((m (max (abs (realpart x))
     908                     (abs (imagpart x))))
     909             (log-s 0)
     910             (s 1))
     911         (if (typep m 'short-float)
     912           (let ((expon (- (%short-float-exp m) IEEE-single-float-bias)))
     913             (cond ((> expon 126)
     914                    (setq log-s double-float-log2^23)
     915                    (setq s #.(ash 1 23)))
     916                   ((< expon -124)
     917                    (setq log-s #.(- double-float-log2^23))
     918                    (setq s #.(/ 1.0s0 (ash 1 23))))))
     919           (let ((expon (- (%double-float-exp m) IEEE-double-float-bias)))
     920             (cond ((> expon 1022)
     921                    (setq log-s double-float-log2^23)
     922                    (setq s #.(ash 1 23)))
     923                   ((< expon -1020)
     924                    (setq log-s #.(- double-float-log2^23))
     925                    (setq s #.(/ 1.0d0 (ash 1 23)))))))
     926         (if (eql s 1)
     927           (complex (log-abs x) (phase x))
     928           (let ((temp (log-abs (/ x s))))
     929             (complex (float (+ log-s temp) temp) (phase x)))))))
    793930    (t
    794      #+32-bit-target
    795      (target::with-stack-short-floats ((sx x))
    796        (if (minusp x)
    797          (complex (%single-float-log! (%%short-float-abs! sx sx) (%make-sfloat))
    798                   #.(coerce pi 'short-float))
    799          (%single-float-log! sx (%make-sfloat))))
    800      #+64-bit-target
    801      (if (minusp x)
    802        (complex (%single-float-log (%short-float-abs (%short-float x))) #.(coerce pi 'single-float))
    803        (%single-float-log (%short-float x)))
    804      )))
    805 
    806 
     931     (%rational-log x 1.0s0))))
     932
     933;;; helpers for rational log
     934(defun %rational-log (x prototype)
     935  (cond ((minusp x)
     936         (complex (%rational-log (- x) prototype)
     937                  (if (typep prototype 'short-float)
     938                    single-float-pi
     939                    pi)))
     940        ((bignump x)
     941         ;(let* ((base1 3)
     942         ;       (guess (floor (1- (integer-length x))
     943         ;                     (log base1 2)))
     944         ;       (guess1 (* guess (log-e base1))))
     945         ;  (+ guess1 (log-e (/ x (expt base1 guess)))))
     946         ; Using base1 = 2 is *much* faster. Was there a reason for 3?
     947         (let* ((guess (1- (integer-length x)))
     948                (guess1 (* guess double-float-log2)))
     949           (float (+ guess1 (log-e (float (/ x (ash 1 guess)) 1.0d0))) prototype)))
     950        ((and (ratiop x)
     951              ;; Rational arguments near +1 can be specified with great precision: don't lose it
     952              (cond ((< 0.5 x 3)
     953                     (log1+ (float (- x 1) prototype)))
     954                    (
     955                     ;; Escape out small values as well as large
     956                     (or (> x most-positive-short-float)
     957                         (< x least-positive-normalized-short-float))
     958                     ;; avoid loss of precision due to subtracting logs of numerator and denominator
     959                     (let* ((n (%numerator x))
     960                            (d (%denominator x))
     961                            (sn (1- (integer-length n)))
     962                            (sd (1- (integer-length d))))
     963                       (float (+ (* (- sn sd) double-float-log2)
     964                                 (- (log1+ (float (1- (/ n (ash 1 sn))) 1.0d0))
     965                                    (log1+ (float (1- (/ d (ash 1 sd))) 1.0d0))))
     966                              prototype))))))
     967        ((typep prototype 'short-float)
     968         #+32-bit-target
     969         (target::with-stack-short-floats ((sx x))
     970           (%single-float-log! sx (%make-sfloat)))
     971         #+64-bit-target
     972         (%single-float-log (%short-float x)))
     973        (t
     974         (with-stack-double-floats ((dx x))
     975           (%double-float-log! dx (%make-dfloat))))))
     976
     977;;; (log1+ x) = (log (1+ x))
     978;;; but is much more precise near x = 0
     979(defun log1+ (x)
     980  ;;(cond ((typep x 'complex)
     981  ;;      (let ((r (realpart x))
     982  ;;            (i (imagpart x)))
     983  ;;        (if (and (< (abs r) 0.5)
     984  ;;                 (< (abs i) 3))
     985  ;;          (let* ((n (+ (* r (+ 2 r)) (* i i)))
     986  ;;                 (d (1+ (sqrt (1+ n)))))
     987  ;;            (complex (log1+ (/ n d)) (atan i (1+ r))))
     988  ;;         (log (1+ x)))))
     989  ;;     (t
     990  (if (and (typep x 'ratio)
     991           (< -0.5 x 2))
     992    (setq x (%short-float x)))
     993  (let ((y (1+ x)))
     994    (if (eql y x)
     995      (log-e y)
     996      (let ((e (1- y)))
     997        (if (zerop e)
     998          (* x 1.0)
     999          (- (log-e y) (/ (- e x) y)))))))
     1000
     1001;;; helper for complex log
     1002;;; uses more careful approach when (abs x) is near 1
     1003(defun log-abs (x)
     1004  (let ((a (abs x)))
     1005    (if (< 0.5 a 3)
     1006      (let* ((r (realpart x))
     1007             (i (imagpart x))
     1008             (n (if (> (abs r) (abs i))
     1009                  (+ (* (1+ r) (1- r)) (* i i))
     1010                  (+ (* r r) (* (1+ i) (1- i))))))
     1011        (log1+ (/ n (1+ a))))
     1012      (log-e a))))
     1013
     1014(defun %rational-complex-log (x prototype &aux ra ia)
     1015  (let* ((rx (realpart x))
     1016         (ix (imagpart x))
     1017         (x (abs rx))
     1018         (y (abs ix)))
     1019    (if (> y x)
     1020      (let ((r (float (/ rx y) 1.0d0)))
     1021        (setq ra (+ (%rational-log y 1.0d0)
     1022                    (/ (log1+ (* r r)) 2)))
     1023        (setq ia (atan (if (minusp ix) -1.0d0 1.0d0) r)))
     1024      (let ((r (float (/ ix x) 1.0d0)))
     1025        (setq ra (+ (%rational-log x 1.0d0)
     1026                    (/ (log1+ (* r r)) 2)))
     1027        (setq ia (atan r (if (minusp rx) -1.0d0 1.0d0)))))
     1028    (if (typep prototype 'short-float)
     1029      (complex (%short-float ra) (%short-float ia))
     1030      (complex ra ia))))
    8071031
    8081032(defun exp (x)
     
    8541078  (cond ((zerop x) x)
    8551079        ((complexp x)
    856          (let* ((i (imagpart x)))
    857            (if (zerop i)                ; has to be a float
    858              (let* ((zero (if (typep i 'double-float) 0.0d0 0.0f0))
    859                     (r (realpart x)))
    860                (if (< r zero)
    861                  (%make-complex zero (float-sign i (sqrt (- r))))
    862                  (%make-complex (abs (sqrt r)) (float-sign i zero))))
    863              (* (sqrt (abs x)) (cis (/ (phase x) 2))))))
     1080         (let ((rx (realpart x))
     1081               (ix (imagpart x)))
     1082           (cond ((rationalp rx)
     1083                  (if (zerop rx)
     1084                    (let ((s (sqrt (/ (abs ix) 2))))
     1085                      (complex s (if (minusp ix) (- s) s)))
     1086                    (let* ((s (+ (* rx rx) (* ix ix)))
     1087                           (d (if (ratiop s)
     1088                                (/ (isqrt (%numerator s))
     1089                                   (isqrt (%denominator s)))
     1090                                (isqrt s))))
     1091                      (unless (eql s (* d d))
     1092                        (setf d (%double-float-hypot (%double-float rx)
     1093                                                     (%double-float ix))))
     1094                      (cond ((minusp rx)
     1095                             (setq b (sqrt (/ (- d rx) 2)))
     1096                             (when (minusp ix)
     1097                               (setq b (- b)))
     1098                             (setq a (/ ix (+ b b))))
     1099                            (t
     1100                             (setq a (sqrt (/ (+ rx d) 2)))
     1101                             (setq b (/ ix (+ a a)))))
     1102                      (if (rationalp a)
     1103                        (complex a b)
     1104                        (complex (%short-float a) (%short-float b))))))
     1105                 ((minusp rx)
     1106                  (if (zerop ix)
     1107                    (complex 0 (float-sign ix (sqrt (- rx))))
     1108                    (let ((shift (cond ((< rx -1) -3)
     1109                                       ((and (> rx -5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25)
     1110                                       (t -1))))
     1111                      (setq rx (scale-float rx shift))
     1112                      (let ((s (fsqrt (- (abs (complex rx (scale-float ix shift))) rx))))
     1113                        (setq b (scale-float s (ash (- -1 shift) -1)))
     1114                        (when (minusp ix)
     1115                          (setq b (- b)))
     1116                        (setq a (/ ix (scale-float b 1)))
     1117                        (complex a b)))))
     1118                 (t
     1119                  (if (zerop ix)
     1120                    (complex (sqrt rx) ix)
     1121                    (let ((shift (cond ((> rx 1) -3)
     1122                                       ((and (< rx 5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25)
     1123                                       (t -1))))
     1124                      (setq rx (scale-float rx shift))
     1125                      (let ((s (fsqrt (+ rx (abs (complex rx (scale-float ix shift)))))))
     1126                        (setq a (scale-float s (ash (- -1 shift) -1)))
     1127                        (setq b (/ ix (scale-float a 1)))
     1128                        (complex a b))))))))
    8641129        ((minusp x) (complex 0 (sqrt (- x))))
    8651130        ((floatp x)
     
    8741139         (/ a b))         
    8751140        (t
    876          #+32-bit-target
    877          (target::with-stack-short-floats ((f1))
    878            (fsqrt (%short-float x f1)))
    879          #+64-bit-target
    880          (fsqrt (%short-float x)))))
     1141         (float (fsqrt (float x 0.0d0)) 1.0s0))))
    8811142
    8821143
Note: See TracChangeset for help on using the changeset viewer.