# Diff of /src/code/irrat.lisp

revision 1.34 by toy, Tue Dec 31 14:28:08 2002 UTC revision 1.35 by toy, Wed Jan 8 23:28:47 2003 UTC
# Line 388  Line 388
388           (if (minusp number)           (if (minusp number)
389               (complex (coerce (log (- number)) 'single-float)               (complex (coerce (log (- number)) 'single-float)
390                        (coerce pi 'single-float))                        (coerce pi 'single-float))
391               (coerce (- (log (numerator number))               ;; What happens when the ratio is close to 1?  We need to
392                          (log (denominator number)))               ;; be careful to preserve accuracy.
393                       'single-float)))               (let ((top (numerator number))
394                       (bot (denominator number)))
395                   ;; If the number of bits in the numerator and
396                   ;; denominator are different, just use the fact
397                   ;; log(x/y) = log(x) - log(y).  However, if they have
398                   ;; the same number of bits, implying the quotient is
399                   ;; near one, we use log1p(x) = log(1+x).  Since the
400                   ;; number is rational, we don't lose precision
401                   ;; subtracting 1 from it, and converting it to
402                   ;; double-float is accurate.
403                   (if (= (integer-length top)
404                                   (integer-length bot))
405                       (coerce (%log1p (coerce (- number 1) 'double-float))
406                               'single-float)
407                       (coerce (- (log2 top) (log2 bot)) 'single-float)))))
408          (((foreach single-float double-float))          (((foreach single-float double-float))
409           ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?           ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
410           ;; Since this doesn't seem to be an implementation issue           ;; Since this doesn't seem to be an implementation issue
# Line 990  Z may be any number, but the result is a Line 1004  Z may be any number, but the result is a
1004  ;;  ;;
1005  ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)  ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1006  ;;  ;;
1007  ;;.and these two expressions are equal if and only if arg conj z =  ;; and these two expressions are equal if and only if arg conj z =
1008  ;;-arg z, which is clearly true for all z.  ;; -arg z, which is clearly true for all z.
1009
1010  (defun complex-acos (z)  (defun complex-acos (z)
1011    "Compute acos z = pi/2 - asin z    "Compute acos z = pi/2 - asin z

Legend:
 Removed from v.1.34 changed lines Added in v.1.35