618251 - Math.pow accumulates error

Status: RESOLVED WONTFIX

(Tamarin Graveyard :: Library, defect)

Description

[Lars T Hansen]

From Shannon Hickey:

  var oneWay:Number = Math.pow(10, 305);
  var theOther:Number = 1e+305;
  trace(oneWay);
  trace(theOther);
  trace(oneWay == theOther);
  trace(oneWay - theOther);

  Results:

  1.00000000000000e+305
  1e+305
  false
  5.84718840683999e+289

  Result line 1 is unexpected.
  Result line 3 is odd.
  Result line 4 is completely crazy.

My analysis:

  Primarily this looks like an inaccuracy in our code for Math.pow.  You
  will notice that the difference between the two numbers is in the 16th
  digit (305-289=16), which cannot be represented reliably by a Number:
  Number only has 15 decimal digits of precision, generally.  The pow()
  code does specialize for integer exponents but there could be inaccuracies
  creeping in around the extremes of the floating-point range (which
  extends to e+308) - I'd have to dig deeper to understand precisely 
  what's going on, though.

  The reason the result of Math.pow() is being written as 
  1.00000000000000e+305 and not (eg) 1.00000000000001e+305, which is 
  presumably the more accurate representation, is a separate bug in 
  the number formatter.

For integer exponents our Math.pow implementation uses a standard doubling algorithm.  It appears that that doubling algorithm accumulates error.  We ought to be able to do better.  That said, 1e305 isn't exactly representable in a Number (it has 708 leading non-zero digits in the binary representation) so it's just a question of how good our approximation is.

Comment 1

[Lars T Hansen]

The problem is the precision loss in the doubling of 'base' in our algorithm.

The following two procedures are Scheme realizations of the algorithm used in Tamarin for an integer exponent.  The first computes e4 using a floating-point accumulator for 'value' but an exact accumulator for 'base'; the second computes e5 using a floating-point accumulator also for 'base' (which is equivalent to what Tamarin does).

(define e4
  (let loop ((value 1.0) (base 10) (e 305))
    (cond ((= e 0) 
	   value)
	  ((not (zero? (remainder e 2)))
	   (loop (* value base) base (- e 1)))
	  (else
	   (loop value (* base base) (quotient e 2))))))

(define e5
  (let loop ((value 1.0) (base 10.0) (e 305))
    (cond ((= e 0) 
	   value)
	  ((not (zero? (remainder e 2)))
	   (loop (* value base) base (- e 1)))
	  (else
	   (loop value (* base base) (quotient e 2))))))

The error accumulation becomes clear when we print the values:

> e4
1.0e305

> e5
1.0000000000000005e305

Comment 2

[Lars T Hansen]

It seems plausible that a table of precisely represented powers of ten within a reasonable range could be used to both speed up the exponentiation and to get better precision.

Comment 3

[Lars T Hansen]

(In reply to comment #2)
> It seems plausible that a table of precisely represented powers of ten within a
> reasonable range could be used to both speed up the exponentiation and to get
> better precision.

Well, that was not very insightful, given that there are many bases for which we would like to have good precision.

Comment 4

[Lars T Hansen]

http://hal.inria.fr/docs/00/16/46/07/PDF/prod.pdf might be relevant.

Comment 5

[Sylvestre Ledru [:Sylvestre]]

Tamarin isn't maintained anymore. WONTFIX remaining bugs.