Re: [Help-smalltalk] Re: [PATCH] Fraction and Integer asFloat: answer nearest floating point

[Paolo Bonzini]

> Yes, our posts did cross... your trailing zero must be "my" lowBit - 1, 
> isn't it?
> I called it lowBit in reference to highBit.
> And i should not say "my" because I am sure there is already a lowBit 
> implemented in Squeak, probably in other dialects too.

Ok, I changed my trailingZeros to lowBit, found a bug in detecting the 
number of digits necessary to print a number, fixed the lex.c bug with 
denormals, and made Float>>#raisedToInteger: more accurate by computing 
the reciprocal (if necessary) at the end of the computation rather than 
at the beginning.  :-)

Paolo

2007-01-29  Paolo Bonzini  <address@hidden>

        * kernel/LargeInt.st: Rename #trailingZeros to #lowBit.
        * kernel/SmallInt.st: Rename #trailingZeros to #lowBit.
        * kernel/Integer.st: Rename #trailingZeros to #lowBit.

        * kernel/Float.st: Fix bug in printing exact floating-point numbers.
        Support gradual underflow in #raisedToInteger:.

        * tests/floatmath.st: Add accuracy tests.
        * tests/floatmath.ok: Adjust.

        * libgst/lex.c: Rename ipowl to mul_powl, support gradual underflow.

--- orig/kernel/Float.st
+++ mod/kernel/Float.st
@@ -132,6 +132,41 @@ negated
 integerPart
     "Return the receiver's integer part"
     ^self - self fractionPart
+!
+
+raisedToInteger: anInteger
+    "Return self raised to the anInteger-th power"
+
+    | exp adjustExp val mant |
+
+    "Some special cases first"
+    anInteger isInteger ifFalse: [ SystemExceptions.WrongClass signalOn: 
anInteger mustBe: Integer ].
+    anInteger = 0 ifTrue: [ ^self unity ].
+    anInteger = 1 ifTrue: [ ^self ].
+
+    "Avoid overflow when the result is denormal and we would have an
+     unrepresentable intermediate result for its reciprocal."
+    adjustExp := self exponent.
+    exp := anInteger abs.
+    (anInteger > 0 or: [ (adjustExp + 1) * exp < self class emax ])
+       ifTrue: [
+           mant := self.
+           adjustExp := 0 ]
+       ifFalse: [
+           mant := self timesTwoPower: 0 - adjustExp.
+           adjustExp := adjustExp * anInteger ].
+
+    "Fire the big loop."
+    val := mant
+        raisedToInteger: exp
+        withCache: ((Array new: (255 min: exp))
+                        at: 1 put: mant;
+                        yourself).
+
+(mant->exp->val->adjustExp) printNl.
+    anInteger < 0 ifTrue: [ val := val reciprocal ].
+    adjustExp = 0 ifFalse: [ val := val timesTwoPower: adjustExp ].
+    ^val
 ! !
 
 
@@ -371,7 +406,7 @@ printOn: aStream special: whatToPrintArr
     "Private - Print a decimal representation of the receiver on aStream,
      printing one of the three elements of whatToPrintArray if it is
      infinity, negative infinity, or a NaN"
-    | me exponential small num den gcd
+    | me exponential small num numLog den denLog gcd
       intFactor precision int rounding digits digitStream exponent
       dotPrinted |
 
@@ -415,11 +450,17 @@ printOn: aStream special: whatToPrintArr
 
     "Get the first `me class decimalDigits' base-10 digits of num // den,
      appropriately rounded"
-    intFactor := 10 raisedToInteger: (den ceilingLog: 10).
-    rounding := (10 raisedToInteger: (num ceilingLog: 10)
-                                    - me class decimalDigits) + 1.
+    numLog := num ceilingLog: 10.
+    denLog := den ceilingLog: 10.
+    numLog > me class decimalDigits
+       ifTrue: [
+           intFactor := 10 raisedToInteger: denLog.
+           rounding := 10 raisedToInteger: numLog - me class decimalDigits ]
+       ifFalse: [
+           intFactor := 10 raisedToInteger: me class decimalDigits - denLog.
+           rounding := 0 ].
 
-    int := ((num * intFactor) + ((den // 2) * rounding)) // den.
+    int := ((num * intFactor) + ((den // 2) * (rounding + 1))) // den.
     digits := int printString.
     digits size > me class decimalDigits
         ifTrue: [ digits := digits copyFrom: 1 to: me class decimalDigits ].


--- orig/kernel/Integer.st
+++ mod/kernel/Integer.st
@@ -143,9 +143,8 @@ noMask: anInteger
     ^(self bitAnd: anInteger) = 0
 !
 
-trailingZeros
-    "Return the index of the lowest order 1 bit of the receiver,
-     minus 1."
+lowBit
+    "Return the index of the lowest order 1 bit of the receiver."
     self subclassResponsibility
 !
 


--- orig/kernel/LargeInt.st
+++ mod/kernel/LargeInt.st
@@ -493,15 +493,15 @@ negated
 
 !LargeInteger methodsFor: 'bit operations'!
 
-trailingZeros
-
-    "Answer the number of trailing zero bits in the receiver"
+lowBit
+    "Return the index of the lowest order 1 bit of the receiver."
 
     | each |
     1 to: self size do: [ :index |
         (each := self digitAt: index) = 0
-            ifFalse: [ ^index * 8 - 8 + (TrailingZeros at: each) ].
+            ifFalse: [ ^index * 8 - 7 + (TrailingZeros at: each) ].
     ].
-    ^self highBit - 1
+    ^self highBit
 !
 
 bitAnd: aNumber
@@ -1000,7 +1000,7 @@ asFloat: characterization
 
     "Apply IEEE 754 round to nearest even default rounding mode"
     carry := self bitAt: nTruncatedBits.
-    (carry = 1 and: [mantissa odd or: [self trailingZeros >= nTruncatedBits]]) 
+    (carry = 1 and: [mantissa odd or: [self lowBit > nTruncatedBits]])
         ifTrue: [mantissa := mantissa + 1].
     ^(characterization coerce: mantissa) timesTwoPower: exponent!
 


--- orig/kernel/SmallInt.st
+++ mod/kernel/SmallInt.st
@@ -94,14 +94,16 @@ generality
 
 !SmallInteger methodsFor: 'bit arithmetic'!
 
-trailingZeros
-    "Return the index of the trailing zero bits in the receiver"
+lowBit
+    "Return the index of the lowest order 1 bit of the receiver."
 
     | n bit |
     self = 0 ifTrue: [ ^0 ].
 
     n := self.
-    bit := 0.
+    "The result is 1-based, but we start from 2 to compensate with the
+     subtraction in the final line."
+    bit := 2.
 
     (n bitAnd: 16r3FFFFFFF) = 0 ifTrue: [
        bit := bit + 30. n := n bitShift: -30 ].
@@ -110,7 +112,7 @@ trailingZeros
     (n bitAnd: 16rFF) = 0 ifTrue: [ bit := bit + 8. n := n bitShift: -8 ].
     (n bitAnd: 16rF) = 0 ifTrue: [ bit := bit + 4. n := n bitShift: -4 ].
     (n bitAnd: 16r3) = 0 ifTrue: [ bit := bit + 2. n := n bitShift: -2 ].
-    ^bit + (1 - (n bitAnd: 16r1))!
+    ^bit - (n bitAnd: 16r1)!
 
 highBit
     "Return the index of the highest order 1 bit of the receiver"


--- orig/libgst/lex.c
+++ mod/libgst/lex.c
@@ -167,9 +167,10 @@ static int scan_symbol (int c,
 static int digit_to_int (int c,
                         int base);
 
-/* Raise BASE to the N-th power.  */
-static inline long double ipowl (long double base,
-                                int n);
+/* Raise BASE to the N-th power and multiply MANT by the result.  */
+static inline long double mul_powl (long double mant,
+                                   long double base,
+                                   int n);
 
 #ifdef LEXDEBUG
 static void print_token (int token,
@@ -649,22 +650,33 @@ scan_ident (int c,
 
 
 long double
-ipowl (long double base, 
-       int n)
+mul_powl (long double mant, 
+         long double base,
+         int n)
 {
-  int k = 1;
   long double result = 1.0;
-  while (n)
+  int n_abs = n < 0 ? -n : n;
+  int exp;
+  int k = 1;
+
+  /* Restrict base to the range from 1.0 to 2.0.  */
+  base = frexpl (base, &exp);
+  base *= 2;
+  exp--;
+
+  while (n_abs)
     {
-      if (n & k)
+      if (n_abs & k)
         {
           result *= base;
-          n ^= k;
+          n_abs ^= k;
         }
       base *= base;
       k <<= 1;
     }
-  return result;
+
+  mant = (n < 0) ? mant / result : mant * result;
+  return ldexpl (mant, exp * n);
 }
 
 int
@@ -836,10 +848,8 @@ scan_number (int c,
   else
     _gst_unread_char (ic);
 
-  if (exponent > 0)
-    num *= ipowl ((long double) base, exponent);
-  else
-    num /= ipowl ((long double) base, -exponent);
+  if (exponent != 0)
+    num = mul_powl (num, (long double) base, exponent);
 
   if (isNegative)
     num = -num;


--- orig/tests/floatmath.ok
+++ mod/tests/floatmath.ok
@@ -220,6 +220,10 @@ true->1.2345e16
 returned value is 1.23450e+16
 
 Execution begins...
+true->1.25
+returned value is 1.25000
+
+Execution begins...
 true->10.0
 returned value is 10.0000
 


--- orig/tests/floatmath.st
+++ mod/tests/floatmath.st
@@ -172,6 +172,7 @@ test
 1e16 test!
 1.2345d16 test!
 1.2345e16 test!
+1.25 test!
 10.0 test!
 (20 - 2.2325251) test!
 0.12345 test!