[Axiom-developer] 20070811.03.tpd.patch applied to silver

[daly]

diff --git a/changelog b/changelog
index 6b8cd6f..bede2a3 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,6 @@
+20070811 sxw src/interp/i-output.boot use digits-by-radix
+20070811 sxw src/interp/sys-pkg.lisp add digits-by-radix export
+20070811 sxw src/interp/vmlisp.lisp add digits-by-radix function
 20070811 gxv src/interp/metalex.lisp remove trace commands
 20070811 tpd src/input/Makefile add mathml.input
 20070811 tpd src/input/mathml.input write test cases
diff --git a/src/interp/i-output.boot.pamphlet 
b/src/interp/i-output.boot.pamphlet
index a1a84ab..b1066aa 100644
--- a/src/interp/i-output.boot.pamphlet
+++ b/src/interp/i-output.boot.pamphlet
@@ -9,25 +9,6 @@
 \eject
 \tableofcontents
 \eject
-\section{GCL\_log10\_bug}
-In some versions of GCL the LOG10 function returns improperly rounded values.
-The symptom is:
-\begin{verbatim}
-(24) -> [1000]
-   (24)  [100]
-\end{verbatim}
-The common lisp failure can be shown with:
-\begin{verbatim}
-(25) -> )lisp (log10 1000)
-Value = 2.9999999999999996
-\end{verbatim}
-This previous boot code was:
-\begin{verbatim}
-    u < MOST_-POSITIVE_-LONG_-FLOAT => 1+negative+FLOOR LOG10 u
-\end{verbatim}
-and should be restored when the GCL bug is fixed.
-<<GCLlog10bug>>=
-    u < MOST_-POSITIVE_-LONG_-FLOAT => 1+negative+FLOOR ((LOG10 u) + 0.0000001)
 @ 
 \section{License}
 <<license>>=
@@ -887,18 +868,12 @@ WIDTH u ==
     u.0="%" and ((u.1 = char 'b) or (u.1 = char 'd)) => 1
     #u
   INTEGERP u => 
+    u = 0 => 1
     if (u < 1) then 
       negative := 1
-      u := -u
     else
       negative := 0
-    -- Try and be fairly exact for smallish integers:
-    u = 0 => 1
-<<GCLlog10bug>>
-    -- Rough guess: integer-length returns log2 rounded up, so divide it by
-    -- roughly log2(10). This should return an over-estimate, but for objects
-    -- this big does it matter?
-    FLOOR(INTEGER_-LENGTH(u)/3.3)
+    DIGITS_-BY_-RADIX(u, 10) + negative
   atom u => # atom2String u
   putWidth u is [[.,:n],:.] => n
   THROW('outputFailure,'outputFailure)
diff --git a/src/interp/sys-pkg.lisp.pamphlet b/src/interp/sys-pkg.lisp.pamphlet
index eaa76b7..9a56b7c 100644
--- a/src/interp/sys-pkg.lisp.pamphlet
+++ b/src/interp/sys-pkg.lisp.pamphlet
@@ -337,7 +337,7 @@ provides support for compiler code.
 <<GCL.DEFINE-MACRO>>
 <<GCL.MEMQ>>
 <<GCL.PNAME>>
-         VMLISP::PUT
+         VMLISP::PUT VMLISP::DIGITS-BY-RADIX
         VMLISP::QVELT-1 VMLISP::QSETVELT-1 vmlisp::throw-protect
         VMLISP::|directoryp| VMLISP::EQCAR
         VMLISP::DEFIOSTREAM VMLISP::RDEFIOSTREAM VMLISP::MLAMBDA
diff --git a/src/interp/vmlisp.lisp.pamphlet b/src/interp/vmlisp.lisp.pamphlet
index 7247d7a..2478e04 100644
--- a/src/interp/vmlisp.lisp.pamphlet
+++ b/src/interp/vmlisp.lisp.pamphlet
@@ -95,6 +95,62 @@ Contributed by Juergen Weiss.
 (defun get-current-directory ()
   (namestring (truename "")))
 
+@ 
+\section{The digits-by-radix function} 
+The purpose of the following function is to calculate the number of
+digits in the radix $B$ expansion of an arbitrary Lisp integer $n$.
+The width of an integer can be determined rapidly when the radix is a
+power of two, otherwise an approach based on successive divisions is
+used.
+
+<<digits-by-radix>>=
+(defun digits-by-radix (n &optional (radix 10))
+  (flet (<<power-of-two-width>>
+         <<iterative-width>>)
+    (assert (>= radix 2) (radix) 
+            "Bad radix ~D < 2 given to DIGITS-BY-RADIX." radix)
+    (setq n (abs n))
+    (cond
+      ((zerop n) (values 1))
+      ((zerop (logand radix (1- radix))) (power-of-two-width n radix))
+      (t (iterative-width n radix)))))
+
+@ When the radix $B$ is of the form $2^b$, $b$ bits are needed to
+represent one radix $B$ digit. The radix $B$ width of $n$ is obtained
+by dividing the width of the binary representation of $n$ by $b$, and
+incrementing the result when the remainder is non-zero.
+
+<<power-of-two-width>>=
+ (power-of-two-width (n radix)
+   (let ((bits (integer-length n))
+         (radix-bits (1- (integer-length radix))))
+     (multiple-value-bind (quo rem) (floor bits radix-bits)
+       (if (zerop rem) quo (1+ quo)))))
+
+@ When the radix is not a power of two, we choose a power $p$ of the
+radix $B$ and use $B^p$ as a divisor.  Each division counts as $p$
+digits in the radix $B$ expansion.  The power, bound to the variable
+[[digits]] below, is chosen so that $B^p <$
+\texttt{most-positive-long-float}. This allows use of [[log]] to
+compute $p$ without concern for floating point overflow.  Once a
+quotient is produced which is smaller than the divisor, we complete
+the calculation by repeated divisions using the radix itself.
+
+<<iterative-width>>=
+ (iterative-width (n radix)
+   (multiple-value-bind (q width)
+       (let* ((target (if (< n most-positive-long-float)
+                          (values n)
+                          (values most-positive-long-float)))
+              (digits (let ((d (floor (log target radix))))
+                        (if (zerop d) 1 d)))
+              (div (expt radix digits)))
+         (loop for q = n then (floor q div)
+               until (< q div) sum digits into width
+               finally (return (values q width))))
+     (+ width (loop for r = q then (floor r radix)
+                    until (zerop r) count t))))
+
 @
 \section{License}
 <<license>>=
@@ -133,6 +189,7 @@ Contributed by Juergen Weiss.
 <<*>>=
 <<license>>
 
+
 ;      VM LISP EMULATION PACKAGE
 ;      Lars Ericson, Barry Trager, Martial Schor, tim daly, LVMCL, et al
 ;      IBM Thomas J. Watson Research Center
@@ -945,6 +1002,8 @@ Contributed by Juergen Weiss.
 
 ; 12.0 Operations on Numbers
 
+<<digits-by-radix>>
+
 ; 12.1 Conversion
 
 (define-function 'FIX #'truncate)