[racket-dev] Sublinear functions of superfloat numbers

From: Robby Findler (robby at eecs.northwestern.edu)
Date: Sun Jul 1 20:17:12 EDT 2012

Oh, that helps a lot, thanks.

Is it not the problem that this:

  (sqrt (expt 10 402))

has no Racket number that could represent it? It is too big for a
float (I think?), it isn't a rational, there isn't anything else...?
No?

Robby

On Sun, Jul 1, 2012 at 7:00 PM, Neil Toronto <neil.toronto at gmail.com> wrote:
> 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 ⊥


Posted on the dev mailing list.