(in-package :cl-user)
;;; This code is written in cl-user, not ccl
;;; The code for new-sqrt, new-log and new-atan follows that for CCL's sqrt, log
;;; and atan (in l0-float.lisp) in terms of the case structure; but where the branch
;;; is not changed the CCL function is called (rather than copying the code in CCL's
;;; branch, which probably wouldn't work here). I hope it shouldn't be too difficult
;;; to cut and paste the appropriate code bits into CCL/s code.
;;; Use of new %hypot is assumed:
;;; (but this version isn't as efficient as the *real* one)
(defun ccl::%hypot (x y)
(let ((a (abs x))
(b (abs y)))
(when (> a b)
(psetq a b b a))
(if (= b 0.0d0)
b
(let ((c (/ a b)))
(* b (sqrt (+ 1.0d0 (* c c))))))))
(defmacro single-float-sqrt (x)
;; usual expression for single float sqrt
;; (referred to several times in new-sqrt)
`(sqrt ,x))
;; fix sqrt to return values for rationals whose square root can be represented
;; as single float, even though they can't be themselves.
;; also take care to avoid intermediate overflow for some complex arguments
(defun new-sqrt (x)
"Return the square root of NUMBER."
(cond ((zerop x) x)
((complexp x)
;; take care that intermediate results do not overflow/underflow:
;; pre-scale argument by an appropriate power of two
;; (there may be more efficient ways of doing this; and strictly
;; if doesn't need to be done for 'reasonably sized' arguments)
(let ((m (max (abs (realpart x))
(abs (imagpart x))))
(s 1)
(s2 1))
(typecase m
(integer (let ((expon (floor (integer-length m) 2)))
(setq s (ash 1 expon))
(setq s2 (ash s expon))))
(ratio (let ((expon (floor (- (integer-length (numerator m))
(integer-length (denominator m)))
2)))
(setq s (expt 2 expon))
(setq s2 (* s s))))
(float (let ((expon (truncate (nth-value 1 (decode-float m)) 2)))
(if (minusp expon)
(incf expon)
(decf expon))
(setq s (scale-float (float 1 m) expon))
(setq s2 (scale-float s expon)))))
(let ((z (/ x s2)))
(* s (* (sqrt (abs z)) (cis (/ (phase z) 2)))))))
((minusp x) (complex 0 (new-sqrt (- x))))
((floatp x)
;; usual floating point sqrt
(sqrt x))
((integerp x)
(let* ((a (isqrt x))
(test (* a a)))
(cond ((eql x test)
a)
((<= x most-positive-single-float)
(single-float-sqrt x))
(t
;; handle case where (sqrt x) can be represented as single float
;; but x can't
(let ((error (/ x test)))
(* (float a 1.0)
(single-float-sqrt error)))))))
((ratiop x)
;; similar to integer case - be careful when (sqrt x) can be represented
;; as single float but x can't
(let* ((n (numerator x))
(a (isqrt n))
(test1 (* a a)))
(unless (eql n test1)
(when (and (<= x most-positive-single-float)
(>= x least-positive-normalized-single-float))
(return-from new-sqrt (single-float-sqrt x))))
(let* ((d (denominator x))
(b (isqrt d))
(test2 (* b b)))
(cond ((and (eql d test2)
(eql n test1))
(/ a b))
((and (<= x most-positive-single-float)
(>= x least-positive-normalized-single-float))
(single-float-sqrt x))
(t
(let ((error (/ x (/ test1 test2))))
(* (float (/ a b) 1.0)
(single-float-sqrt error))))))))
(t
(single-float-sqrt x))))
;; makes log work for small rationals (it already worked for big ones)
;; correct wrong return type for negative bignums
;; fix float/rational contagion in two argument case
;; propagate rational behaviour to complex rationals as well
(defun new-log (x &optional (b nil b-p))
"Return the logarithm of NUMBER in the base BASE, which defaults to e."
(if b-p
(if (zerop b)
(if (zerop x)
;(report-bad-arg x '(not (satisfies zerop) ))
(error "NEW-LOG: zero arguments")
;(if (floatp x) (float 0.0d0 x) 0)) ** CORRECT THE CONTAGION for complex args **
; (is this the best way to get a zero of the right type?)
(+ 0 (* x b)))
;; ** Don't we have to do the float/rational contagion before the division?
;; (is this the best way to do it?)
(/ (new-log-e (+ x (* b 0)))
(new-log-e (+ b (* x 0)))))
(new-log-e x)))
;;; constant required: will need to be 'magic' number somewhere
(defconstant +log2+ #.(log 2.0d0)) ; = 0.6931471805599453D0
(defun new-log-e (x)
(cond
((bignump x)
(if (minusp x)
(complex (new-log-e (- x)) #.(coerce pi 'single-float)) ; ** FIX RETURN TYPE ERROR **
(let* ((base1 3)
(guess (floor (1- (integer-length x))
(log base1 2)))
(guess1 (* guess (new-log-e base1))))
(+ guess1 (new-log-e (/ x (expt base1 guess)))))))
((and (ratiop x) ; Escape out small values as well as large
(if (minusp x)
(or (> x least-negative-normalized-short-float)
(< x most-negative-short-float))
(or (> x most-positive-short-float)
(< x least-positive-normalized-short-float))))
(- (new-log-e (numerator x)) (new-log-e (denominator x))))
((typep x 'complex)
;; take care that intermediate results do not overflow/underflow:
;; pre-scale argument by an appropriate power of two
;; we only need to scale for very large and very small values -
;; hence the various 'magic' numbers (values may not be too
;; critical but do depend on the sqrt update to fix abs's operation)
(let ((m (max (abs (realpart x))
(abs (imagpart x))))
(log-s 0)
(s 1))
(typecase m
(integer (let ((expon (integer-length m)))
(when (> expon 126)
(setq log-s (float (* expon +log2+) 1.0))
(setq s (ash 1 expon)))))
(ratio (let ((expon (- (integer-length (numerator m))
(integer-length (denominator m)))))
(when (or (> expon 125)
(< expon -122))
(setq log-s (float (* expon +log2+) 1.0))
(setq s (expt 2 expon)))))
(single-float (let ((expon (nth-value 1 (decode-float m))))
(when (or (> expon 126)
(< expon -124))
(if (minusp expon)
(incf expon)
(decf expon))
(setq log-s (float (* expon +log2+) 1.0)) ; assumes (float-radix m) = 2
(setq s (scale-float 1.0 expon)))))
(double-float (let ((expon (nth-value 1 (decode-float m))))
(when (or (> expon 1022)
(< expon -1020))
(if (minusp expon)
(incf expon)
(decf expon))
(setq log-s (* expon +log2+)) ; assumes (float-radix m) = 2
(setq s (scale-float 1.0d0 expon))))))
(if (eql s 1)
(complex (new-log-e (abs x)) (phase x))
(let ((z (/ x s)))
(complex (+ log-s (new-log-e (abs z))) (phase z))))))
((typep x 'double-float)
;; usual double float log-e
(log x))
(t
;; usual single float log-e
(log x)
)))
;; atan needs to be more careful about small and large rational arguments
;; (this will make phase as careful too for complex rationals)
(defun new-atan (y &optional (x nil x-p))
"Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
(cond (x-p
(when (or (rationalp x) (rationalp y))
(unless (zerop y) ; don't tread on special zero handling
;; rescale arguments so that the maximum absolute value is 1
(cond ((> (abs y) (abs x))
(setf x (/ x (abs y)))
(setf y (signum y)))
(t
(setf y (/ y (abs x)))
(setf x (signum x))))))
;; usual expression
(atan y x))
(t
;; usual expressions
;; [but (sqrt -1) in the complex case could be just #C(0 1)]
(atan y))))
;;; FIXES
#|
(sqrt (expt 10 47)) => 3.1622778E+23
(sqrt (/ (expt 10 47) 3)) => 1.8257418E+23
(sqrt (complex (expt 10 46) (expt 10 47))) => #C(2.3505187E+23 2.12719E+23)
(sqrt (complex most-positive-short-float most-positive-short-float)) => #C(2.0267142E+19 8.394926E+18)
(log (expt 10 -66)) => -151.97063
(log (- (expt 10 66))) => #C(151.97063 3.1415927) [corrects type]
(log (complex (expt 10 65) (expt 10 66))) => #C(151.9756 1.4711276)
(log (complex (expt 10 -65) (expt 10 -66))) => #C(-149.66307 0.09966865)
(log 8.0d0 2) => 3.0D0 [corrects accuracy]
(log #C(0.0 1.0) 0) => #C(0.0 0.0) [corrects type]
(atan (expt 10 46) (expt 10 47)) => 0.09966865
(atan (expt 10 -46) (expt 10 -47)) => 1.4711276 [corrects loss of accuracy]
|#