mul_fft, cleaning up some details of the code

[Paul Zimmermann]

       Dear Marco,

> Date: Sun, 06 Mar 2022 11:14:49 +0100
> From: Marco Bodrato <bodrato at 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