Anomaly in mpn_sqrtrem and mpn_rottrem

[Niels Möller]

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.