Anomaly in mpn_sqrtrem and mpn_rottrem
Niels Möller
nisse at lysator.liu.se
Thu Jul 9 20:56:02 UTC 2015
nisse at lysator.liu.se (Niels Möller) writes:
> sqrt_nm1 (n, A):
>
> Input: A = {a_{n-1}, ..., a_0}, n >= 2
> Output: X = {x_{n-2}, ..., x_0} \approx sqrt(B^{n-2} A)
>
> if (n == 2)
> x_0 <-- floor (sqrt ({a_1, a_0}))
> return X = {x_0}
>
> k <-- floor (n/2)
> H <-- sqrt_nm1 (k+1, {a_{n-1},...,a_{n-1-k}}) // k limbs
>
> if (n odd) // n = 2k+1
> E <-- B H^2 - A
> X <-- B^k H - floor (B^{k-1} E / 2H)
This seems to give a much worse error than I expected. Let's focus on
the case n == 3, and B = 2^64. Also switch sign of the residual. We get
k = 1, so this boils down to
x_1 = floor(sqrt({a_2, a_1}))
E = A - B x_1^2
x_0 = floor (E / 2 x_1)
Here's a failing case from my tests
A = 1fe000000000 f8000000000003ff ffe1fffc00ffffff
We're expecting to compute something close to floor (sqrt (A B)), which
is
X = 5a552d07350bac a6d3c0a5a9b624a1
We get x_1 = 5a552d07350bac, good so far. We next compute the residual,
E = 75bbe6c23fc86f ffe1fffc00ffffff
Dividn by 2 x_1, we get finally
x_0 = a6d3c0a5a9b6253b
This differs from the correct value by 0x9a, where I hoped for an error
of at most one or two.
I suspect that the problem is the "high half" ought to include half of
the output bits (plus a little margin). And in the failing example, the
correct square root is 119 bits, but the "high half" x_1 is only 55
bits. And even with a normalized 192-bit input, we'd get only 64 out of
128 bits, which seems a bit too tight.
So back to the drawing board.
Regards,
/Niels
--
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.
More information about the gmp-devel
mailing list