Here is the "simple algorithm" in Scheme for you to experiment with where I
added the discussed test cases. The code should be standard R6RS except
for Chez's printf in the test procedure. The structure and names follow
the original code proposed for GNU MP.
PS: Around the time when I started this thread here, I asked about
mpq_get_nearest_d on the mailing list [1]. If you know people from GNU MP,
maybe you can bring it to their attention (if you think it is a good idea).
;;; THE CODE
(define approximate
(lambda (x p)
(assert (and (exact? x) (rational? x)))
(assert (and (exact? p) (integer? p) (positive? p)))
(let ([a (numerator x)]
[b (denominator x)])
(let ([sa (bitwise-length a)]
[sb (bitwise-length b)])
(let-values
([(aa bb)
(if (>= sa sb)
(values a (bitwise-arithmetic-shift-left b (- sa sb)))
(values (bitwise-arithmetic-shift-left a (- sb sa)) b))])
(let-values
([(bb sa)
(if (>= aa bb)
(values (bitwise-arithmetic-shift-left bb 1) (+ sa 1))
(values bb sa))])
(let ([aa (bitwise-arithmetic-shift-left aa (+ p 1))]
[sb (+ sb p 1)])
(let-values ([(aa bb) (div-and-mod aa bb)])
(let
([aa
(cond
[(not (bitwise-bit-set? aa 0))
aa]
[(not (zero? bb))
(if (positive? aa)
(+ aa 1)
(- aa 1))]
[(not (bitwise-bit-set? aa 1))
(if (positive? aa)
(- aa 1)
(+ aa 1))]
[(positive? aa)
(+ aa 1)]
[else
(- aa 1)])])
(* aa (expt 2 (- sa sb))))))))))))
(define test
(lambda (x p)
(printf "(approximate ~s ~s) => ~s~%"
x p (approximate x p))))
(test #e0.1 1)
(test 99999999999999983222784 54)
(test 99999999999999983222784 53)
(test 99999999999999983222783 54)
(test 99999999999999983222783 53)
Am Mi., 4. Sept. 2024 um 08:11 Uhr schrieb Marc Nieper-Wißkirchen <
[email protected]>:
> Am Mi., 4. Sept. 2024 um 03:30 Uhr schrieb Bradley Lucier <
> [email protected]>:
>
>> I *really* don't want to get into this discussion right now. HOWEVER ...
>>
>> Perhaps it would be good for Marc to show the current Scheme
>> implementation of his idea.
>>
>
> The algorithm is basically the simple one from [1] (without a hardcoded
> mantissa width), but you can use any other algorithm that rounds a rational
> number to p significant digits. Proper Scheme code will follow in some
> subsequent post.
>
> I like to make two more comments that I should have made in my previous
> post but forgot:
>
> (1) Everything I proved about Claim 2 did not take the existence of
> subnormal numbers into account as this wasn't relevant to the actual
> discussion.
>
> (2) The more important comment is this one, where I try to explain again
> in different words why "double rounding" is not the right term: When an
> implementation supporting explicitly given mantissa widths writes a
> binary-floating point number and it will choose the shortest decimal syntax
> x in printing x|p, the actual x has no meaning in itself. It only has a
> meaning together with the p. For example, an implementation that outputs a
> number object in the subnormal range for that implementation, will print
> something like, say, x|50, encoding the 50 significant binary digits in a
> decimal number x that is as short as possible so that x|p gives the actual
> subnormal number. It is not important to try to approximate x as well as
> possible.
>
> Marc
>
> --
>
> [1] https://gmplib.org/list-archives/gmp-devel/2013-April/003223.html
>
>
>> Brad
>>
>
