ni...@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. _______________________________________________ gmp-devel mailing list gmp-devel@gmplib.org https://gmplib.org/mailman/listinfo/gmp-devel