Dear Marco,
> Date: Sun, 06 Mar 2022 11:14:49 +0100
> From: Marco Bodrato <bodr...@mail.dm.unipi.it>
>
> Ciao,
>
> I looked into the code published by Samuel Vivien.
> I realised that in mul_fft there are a lot of _add_1 and _sub_1. At
> least some of them can be easily replaced by _INCR_U or _DECR_U...
>
> Specifically I'd focus into a suspicious piece of code, shared by both
> our current code and Vivien's.
>
> The function mpn_mul_fft_decompose has a branch "if (dif > Kl)", that
> ends with the following lines:
>
> if (cy >= 0)
> cy = mpn_add_1 (tmp, tmp, Kl, cy);
> else
> cy = mpn_sub_1 (tmp, tmp, Kl, -cy);
>
> Can those lines be correct? Can the value in cy be used after this code
> both if it contains the carry or if it contains the borrow of the
> operation on tmp?
thank you, I agree with your analysis. Computations in that part are done
modulo B^Kl+1, thus any carry (resp. borrow) at tmp[Kl] should be subtracted
(resp. added) at tmp[0]. Here are some comments added to the code (from
GMP 6.2.1) and a fix:
--- mul_fft.c.orig 2022-03-07 09:56:28.133700601 +0100
+++ mul_fft.c 2022-03-07 10:12:20.127270744 +0100
@@ -712,15 +712,28 @@
cy += mpn_sub (tmp, tmp, Kl, n, dif);
else
cy -= mpn_add (tmp, tmp, Kl, n, dif);
+ /* cy is the borrow at tmp[Kl], thus we should subtract
+ cy at tmp+Kl, or equivalently add cy at tmp, since
+ B^Kl = -1 */
if (cy >= 0)
cy = mpn_add_1 (tmp, tmp, Kl, cy);
+ /* cy is now the carry at tmp[Kl] */
else
+ {
cy = mpn_sub_1 (tmp, tmp, Kl, -cy);
+ /* cy is now the borrow at tmp[Kl] */
+ if (cy)
+ cy = mpn_add_1 (tmp, tmp, Kl, cy);
+ /* cy is now the carry at tmp[Kl] */
+ }
}
else /* dif <= Kl, i.e. nl <= 2 * Kl */
{
cy = mpn_sub (tmp, n, Kl, n + Kl, dif);
+ /* cy is the borrow at tmp[Kl], thus we should add
+ cy at tmp[0] */
cy = mpn_add_1 (tmp, tmp, Kl, cy);
+ /* cy is now the carry at tmp[Kl] */
}
tmp[Kl] = cy;
nl = Kl + 1;
> I believe this is an error, I'd change the last line into
> cy = 1 - mpn_sub_1 (tmp, tmp, Kl, -cy - 1);
> and I propose the attached patch changing this and some other _1
> functions.
Yes this seems ok. Assume we are in the case cy < 0.
Since cy is the borrow at tmp[Kl], we should add cy at tmp[0],
thus (since cy < 0) subtract -cy (say b=-cy) at tmp[0].
Since b > 0 and B^Kl+1 = 0 it is equivalent to subtract b-1 at
tmp[0] and add 1 at tmp[Kl].
> This is not really needed to solve a bug.
> The comment before the mpn_mul_fft_decompose function says "We must have
> nl <= 2*K*l", this means that we should never have ((dif = nl - Kl) >
> Kl), and the code in that branch should never be used.
> According to the coverage analisys, the code is not explored by the
> tests:
> https://gmplib.org/devel/lcov/shell/var/tmp/lcov/gmp/mpn/mul_fft.c.gcov#691
> The "bug" is never triggered.
>
> Maybe the branch could be used if someone re-enables mpn_mul_fft_full
> and uses it for very unbalanced (more than 4x1) operands?
yes it could be used in that case.
Best regards,
Paul
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
