From: Kaz K. <kky...@gm...> - 2009-01-17 11:09:28
|
Here is a crazy function I hacked up (defun bignum-log (x &optional (b 10)) "logarithm function for CLISP that tries to handle large numbers when the answer is real, and the LOG function blows up. This uses the transformation log x = log x / y + log y. A suitable y is found which is a power of the base that is close to magnitude to x. CLISP can deal with logs of small rations between large bignums better than logs of bignum integers." ;; cases we don't handle (when (or (minusp x) (minusp b) (not (integerp x)) (not (integerp b))) (log x b)) ;; Try LOG function first, but don't bother ;; if X is ``obviously'' too large by inspection. (let* ((obviously-too-large (load-time-value (expt 10 50))) (try-log (if (not (> x obviously-too-large)) (ignore-errors (log x b))))) (if try-log try-log (loop with y = (expt b 100) with lowest-acceptable = (/ x (expt b 4)) with low-probe = 0 with high-probe = nil do (flet ((bisect() (let* ((delta (/ high-probe low-probe)) (digits (log delta b))) (* low-probe (expt b (truncate digits 2)))))) (cond ((< y lowest-acceptable) (when (> y low-probe) (setf low-probe y)) (if high-probe (setf y (bisect)) (setf y (* y y)))) ((> y x) (when (or (null high-probe) (< y high-probe)) (setf high-probe y)) (setf y (bisect))) (t (return (+ (log (/ x y) b) (log y b)))))))))) |