[racket-dev] Sublinear functions of superfloat numbers

How about more words and examples?

"Argument reduction" is using function properties to reduce the 
magnitude of arguments to make computation more tractable or more accurate.

I'll bet `log' uses this property:

     (log (sqrt x)) = (log (expt x 1/2)) = (* 1/2 (log x))

This form is nice for doing algebraic manipulations; not so much for 
computation. So multiply the outer sides by 2:

     (* 2 (log (sqrt x))) = (log x)

You could define your own log function like this:

     (define (my-log x)
       (* 2 (log (sqrt x))))

If `sqrt' could compute square roots of rational numbers larger than 
+max.0 (about 2^1024), then `my-log' could compute logs for those as 
well. But it can't.

Here's a version of log that does reduces its argument with 
`integer-sqrt', which computes truncated square roots of bigints:

(require racket/flonum
          (only-in unstable/flonum +max.0))

(define (real-log x)
   (cond [(x . = . +inf.0)   +inf.0]
         [(x . <= . +max.0)  (fllog (real->double-flonum x))]
         [(x . > . +max.0)
          (let loop ([x  (exact-round x)])
            (cond [(x . > . +max.0)  (* 2.0 (loop (integer-sqrt x)))]
                  [else  (fllog (real->double-flonum x))]))]
         [else  +nan.0]))

Computing (real-log #e1e800242) takes exactly the same amount of time as 
(log #e1e800242), and gives the same answer as well. (And I just 
discovered that that's how it's implemented in number.c.)

Square root (for large real numbers) is much simpler. Floating-point 
numbers above 2^52 are all integers, so bigints above 2^1024 can be 
thought of as floating-point numbers with at least 1024 - 52 = 972 bits 
of precision. That means `integer-sqrt' will do the job perfectly, 
despite the fact that it always returns integers.

(define (real-sqrt x)
   (cond [(x . = . +inf.0)   +inf.0]
         [(x . <= . +max.0)  (flsqrt (real->double-flonum x))]
         [(x . > . +max.0)   (real->double-flonum
                              (integer-sqrt (exact-round x)))]
         [else  +nan.0]))

But sine is much harder because it's periodic with an irrational period. 
It would look something like this, but with a rational approximation of 
pi whose precision depends on the magnitude of the argument:

(define (real-sin x)
   (cond ...
         [(x . > . +max.0)
          (let ([x  (- x (truncate (/ x (* 2 pi))))])
            (flsin (real->double-flonum x)))]
         ...))

On 07/01/2012 04:04 PM, Robby Findler wrote:
> 3. Can you explain the issue again, using smaller words? (I think I
> understand the first example, but then I'm lost.)
>
> Robby
>
> On Sun, Jul 1, 2012 at 5:02 PM, Matthias Felleisen <matthias at ccs.neu.edu> wrote:
>>
>> 1. What's the computational cost of such changes?

The additional cost when applying `sqrt' and `sin' to numbers in typical 
ranges would be small: the cost of wrapping a kernel function and 
checking the size of arguments.

>> 2. What is the impact on TR?

None that I can tell. But TR would remove the additional cost when it 
could prove the arguments to `sqrt' and `sin' are Float.

----

I think I've been trying to come up with a general rule for when 
argument reduction is necessary. But I can't, because there isn't one. 
For example, this is infinite:

     (log (make-rectangular #e1e400 1))

even though the actual answer is representable as a Float-Complex. 
Apparently, argument reduction only happens to reals, and only in a few 
functions. (But I'm glad it does, because the plot library uses `log' to 
format huge numbers.)

Neil ⊥