/src/gmp/mpn/toom43_mul.c

Source
/* mpn_toom43_mul -- Multiply {ap,an} and {bp,bn} where an is nominally 4/3
   times as large as bn.  Or more accurately, bn < an < 2 bn.

   Contributed to the GNU project by Marco Bodrato.

   The idea of applying toom to unbalanced multiplication is due to Marco
   Bodrato and Alberto Zanoni.

   THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
   SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.

Copyright 2009, 2020 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of either:

  * the GNU Lesser General Public License as published by the Free
    Software Foundation; either version 3 of the License, or (at your
    option) any later version.

or

  * the GNU General Public License as published by the Free Software
    Foundation; either version 2 of the License, or (at your option) any
    later version.

or both in parallel, as here.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received copies of the GNU General Public License and the
GNU Lesser General Public License along with the GNU MP Library.  If not,
see https://www.gnu.org/licenses/.  */


#include "gmp-impl.h"

/* Evaluate in: -2, -1, 0, +1, +2, +inf

  <-s-><--n--><--n--><--n-->
   ___ ______ ______ ______
  |a3_|___a2_|___a1_|___a0_|
  |_b2_|___b1_|___b0_|
  <-t--><--n--><--n-->

  v0  =  a0             * b0          #   A(0)*B(0)
  v1  = (a0+ a1+ a2+ a3)*(b0+ b1+ b2) #   A(1)*B(1)      ah  <= 3  bh <= 2
  vm1 = (a0- a1+ a2- a3)*(b0- b1+ b2) #  A(-1)*B(-1)    |ah| <= 1 |bh|<= 1
  v2  = (a0+2a1+4a2+8a3)*(b0+2b1+4b2) #   A(2)*B(2)      ah  <= 14 bh <= 6
  vm2 = (a0-2a1+4a2-8a3)*(b0-2b1+4b2) #  A(-2)*B(-2)    |ah| <= 9 |bh|<= 4
  vinf=              a3 *         b2  # A(inf)*B(inf)
*/

void
mpn_toom43_mul (mp_ptr pp,
    mp_srcptr ap, mp_size_t an,
    mp_srcptr bp, mp_size_t bn, mp_ptr scratch)
{
  mp_size_t n, s, t;
  enum toom6_flags flags;
  mp_limb_t cy;

#define a0  ap
#define a1  (ap + n)
#define a2  (ap + 2 * n)
#define a3  (ap + 3 * n)
#define b0  bp
#define b1  (bp + n)
#define b2  (bp + 2 * n)

  n = 1 + (3 * an >= 4 * bn ? (an - 1) >> 2 : (bn - 1) / (size_t) 3);

  s = an - 3 * n;
  t = bn - 2 * n;

  ASSERT (0 < s && s <= n);
  ASSERT (0 < t && t <= n);

  /* This is true whenever an >= 25 or bn >= 19, I think. It
     guarantees that we can fit 5 values of size n+1 in the product
     area. */
  ASSERT (s+t >= 5);

#define v0    pp        /* 2n */
#define vm1   (scratch)        /* 2n+1 */
#define v1    (pp + 2*n)      /* 2n+1 */
#define vm2   (scratch + 2 * n + 1)    /* 2n+1 */
#define v2    (scratch + 4 * n + 2)    /* 2n+1 */
#define vinf  (pp + 5 * n)      /* s+t */
#define bs1    pp        /* n+1 */
#define bsm1  (scratch + 2 * n + 2)    /* n+1 */
#define asm1  (scratch + 3 * n + 3)    /* n+1 */
#define asm2  (scratch + 4 * n + 4)    /* n+1 */
#define bsm2  (pp + n + 1)      /* n+1 */
#define bs2   (pp + 2 * n + 2)      /* n+1 */
#define as2   (pp + 3 * n + 3)      /* n+1 */
#define as1   (pp + 4 * n + 4)      /* n+1 */

  /* Total sccratch need is 6 * n + 3 + 1; we allocate one extra
     limb, because products will overwrite 2n+2 limbs. */

#define a0a2  scratch
#define b0b2  scratch
#define a1a3  asm1
#define b1d   bsm1

  /* Compute as2 and asm2.  */
  flags = (enum toom6_flags) (toom6_vm2_neg & mpn_toom_eval_dgr3_pm2 (as2, asm2, ap, n, s, a1a3));

  /* Compute bs2 and bsm2.  */
  b1d[n] = mpn_lshift (b1d, b1, n, 1);      /*       2b1      */
#if HAVE_NATIVE_mpn_addlsh2_n
  cy = mpn_addlsh2_n (b0b2, b0, b2, t);     /*  4b2      + b0 */
#else
  cy  = mpn_lshift (b0b2, b2, t, 2);      /*  4b2           */
  cy += mpn_add_n (b0b2, b0b2, b0, t);     /*  4b2      + b0 */
#endif
  if (t != n)
    cy = mpn_add_1 (b0b2 + t, b0 + t, n - t, cy);
  b0b2[n] = cy;

#if HAVE_NATIVE_mpn_add_n_sub_n
  if (mpn_cmp (b0b2, b1d, n+1) < 0)
    {
      mpn_add_n_sub_n (bs2, bsm2, b1d, b0b2, n+1);
      flags = (enum toom6_flags) (flags ^ toom6_vm2_neg);
    }
  else
    {
      mpn_add_n_sub_n (bs2, bsm2, b0b2, b1d, n+1);
    }
#else
  mpn_add_n (bs2, b0b2, b1d, n+1);
  if (mpn_cmp (b0b2, b1d, n+1) < 0)
    {
      mpn_sub_n (bsm2, b1d, b0b2, n+1);
      flags = (enum toom6_flags) (flags ^ toom6_vm2_neg);
    }
  else
    {
      mpn_sub_n (bsm2, b0b2, b1d, n+1);
    }
#endif

  /* Compute as1 and asm1.  */
  flags = (enum toom6_flags) (flags ^ (toom6_vm1_neg & mpn_toom_eval_dgr3_pm1 (as1, asm1, ap, n, s, a0a2)));

  /* Compute bs1 and bsm1.  */
  bsm1[n] = mpn_add (bsm1, b0, n, b2, t);
#if HAVE_NATIVE_mpn_add_n_sub_n
  if (bsm1[n] == 0 && mpn_cmp (bsm1, b1, n) < 0)
    {
      cy = mpn_add_n_sub_n (bs1, bsm1, b1, bsm1, n);
      bs1[n] = cy >> 1;
      flags = (enum toom6_flags) (flags ^ toom6_vm1_neg);
    }
  else
    {
      cy = mpn_add_n_sub_n (bs1, bsm1, bsm1, b1, n);
      bs1[n] = bsm1[n] + (cy >> 1);
      bsm1[n]-= cy & 1;
    }
#else
  bs1[n] = bsm1[n] + mpn_add_n (bs1, bsm1, b1, n);
  if (bsm1[n] == 0 && mpn_cmp (bsm1, b1, n) < 0)
    {
      mpn_sub_n (bsm1, b1, bsm1, n);
      flags = (enum toom6_flags) (flags ^ toom6_vm1_neg);
    }
  else
    {
      bsm1[n] -= mpn_sub_n (bsm1, bsm1, b1, n);
    }
#endif

  ASSERT (as1[n] <= 3);
  ASSERT (bs1[n] <= 2);
  ASSERT (asm1[n] <= 1);
  ASSERT (bsm1[n] <= 1);
  ASSERT (as2[n] <=14);
  ASSERT (bs2[n] <= 6);
  ASSERT (asm2[n] <= 9);
  ASSERT (bsm2[n] <= 4);

  /* vm1, 2n+1 limbs */
  vm1[2*n] = 0;
  mpn_mul_n (vm1, asm1, bsm1, n + (asm1[n] | bsm1[n]));  /* W4 */

  /* vm2, 2n+1 limbs */
  mpn_mul_n (vm2, asm2, bsm2, n+1);  /* W2 */

  /* v2, 2n+1 limbs */
  mpn_mul_n (v2, as2, bs2, n+1);  /* W1 */

  /* v1, 2n+1 limbs */
  mpn_mul_n (v1, as1, bs1, n+1);  /* W3 */

  /* vinf, s+t limbs */   /* W0 */
  if (s > t)  mpn_mul (vinf, a3, s, b2, t);
  else        mpn_mul (vinf, b2, t, a3, s);

  /* v0, 2n limbs */
  mpn_mul_n (v0, ap, bp, n);  /* W5 */

  mpn_toom_interpolate_6pts (pp, n, flags, vm1, vm2, v2, t + s);

#undef v0
#undef vm1
#undef v1
#undef vm2
#undef v2
#undef vinf
#undef bs1
#undef bs2
#undef bsm1
#undef bsm2
#undef asm1
#undef asm2
/* #undef as1 */
/* #undef as2 */
#undef a0a2
#undef b0b2
#undef a1a3
#undef b1d
#undef a0
#undef a1
#undef a2
#undef a3
#undef b0
#undef b1
#undef b2
}

Coverage Report

Created: 2026-03-31 06:37

Line	Count	Source
1		/* mpn_toom43_mul -- Multiply {ap,an} and {bp,bn} where an is nominally 4/3
2		times as large as bn. Or more accurately, bn < an < 2 bn.
3
4		Contributed to the GNU project by Marco Bodrato.
5
6		The idea of applying toom to unbalanced multiplication is due to Marco
7		Bodrato and Alberto Zanoni.
8
9		THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
10		SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
11		GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
12
13		Copyright 2009, 2020 Free Software Foundation, Inc.
14
15		This file is part of the GNU MP Library.
16
17		The GNU MP Library is free software; you can redistribute it and/or modify
18		it under the terms of either:
19
20		* the GNU Lesser General Public License as published by the Free
21		Software Foundation; either version 3 of the License, or (at your
22		option) any later version.
23
24		or
25
26		* the GNU General Public License as published by the Free Software
27		Foundation; either version 2 of the License, or (at your option) any
28		later version.
29
30		or both in parallel, as here.
31
32		The GNU MP Library is distributed in the hope that it will be useful, but
33		WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
34		or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
35		for more details.
36
37		You should have received copies of the GNU General Public License and the
38		GNU Lesser General Public License along with the GNU MP Library. If not,
39		see https://www.gnu.org/licenses/. */
40
41
42		#include "gmp-impl.h"
43
44		/* Evaluate in: -2, -1, 0, +1, +2, +inf
45
46		<-s-><--n--><--n--><--n-->
47		___ ______ ______ ______
48		\|a3_\|___a2_\|___a1_\|___a0_\|
49		\|_b2_\|___b1_\|___b0_\|
50		<-t--><--n--><--n-->
51
52		v0 = a0 * b0 # A(0)*B(0)
53		v1 = (a0+ a1+ a2+ a3)(b0+ b1+ b2) # A(1)B(1) ah <= 3 bh <= 2
54		vm1 = (a0- a1+ a2- a3)(b0- b1+ b2) # A(-1)B(-1) \|ah\| <= 1 \|bh\|<= 1
55		v2 = (a0+2a1+4a2+8a3)(b0+2b1+4b2) # A(2)B(2) ah <= 14 bh <= 6
56		vm2 = (a0-2a1+4a2-8a3)(b0-2b1+4b2) # A(-2)B(-2) \|ah\| <= 9 \|bh\|<= 4
57		vinf= a3 * b2 # A(inf)*B(inf)
58		*/
59
60		void
61		mpn_toom43_mul (mp_ptr pp,
62		mp_srcptr ap, mp_size_t an,
63		mp_srcptr bp, mp_size_t bn, mp_ptr scratch)
64	445	{
65	445	mp_size_t n, s, t;
66	445	enum toom6_flags flags;
67	445	mp_limb_t cy;
68
69	445	#define a0 ap
70	445	#define a1 (ap + n)
71	445	#define a2 (ap + 2 * n)
72	445	#define a3 (ap + 3 * n)
73	1.22k	#define b0 bp
74	1.73k	#define b1 (bp + n)
75	1.33k	#define b2 (bp + 2 * n)
76
77	445	n = 1 + (3 * an >= 4 * bn ? (an - 1) >> 2 : (bn - 1) / (size_t) 3);
78
79	445	s = an - 3 * n;
80	445	t = bn - 2 * n;
81
82	445	ASSERT (0 < s && s <= n);
83	445	ASSERT (0 < t && t <= n);
84
85		/* This is true whenever an >= 25 or bn >= 19, I think. It
86		guarantees that we can fit 5 values of size n+1 in the product
87		area. */
88	445	ASSERT (s+t >= 5);
89
90	445	#define v0 pp /* 2n */
91	1.33k	#define vm1 (scratch) /* 2n+1 */
92	445	#define v1 (pp + 2n) / 2n+1 */
93	890	#define vm2 (scratch + 2 * n + 1) /* 2n+1 */
94	890	#define v2 (scratch + 4 * n + 2) /* 2n+1 */
95	445	#define vinf (pp + 5 * n) /* s+t */
96	1.33k	#define bs1 pp /* n+1 */
97	6.87k	#define bsm1 (scratch + 2 * n + 2) /* n+1 */
98	1.78k	#define asm1 (scratch + 3 * n + 3) /* n+1 */
99	890	#define asm2 (scratch + 4 * n + 4) /* n+1 */
100	890	#define bsm2 (pp + n + 1) /* n+1 */
101	890	#define bs2 (pp + 2 * n + 2) /* n+1 */
102	890	#define as2 (pp + 3 * n + 3) /* n+1 */
103	890	#define as1 (pp + 4 * n + 4) /* n+1 */
104
105		/* Total sccratch need is 6 * n + 3 + 1; we allocate one extra
106		limb, because products will overwrite 2n+2 limbs. */
107
108	445	#define a0a2 scratch
109	3.44k	#define b0b2 scratch
110	445	#define a1a3 asm1
111	2.22k	#define b1d bsm1
112
113		/* Compute as2 and asm2. */
114	445	flags = (enum toom6_flags) (toom6_vm2_neg & mpn_toom_eval_dgr3_pm2 (as2, asm2, ap, n, s, a1a3));
115
116		/* Compute bs2 and bsm2. */
117	445	b1d[n] = mpn_lshift (b1d, b1, n, 1); /* 2b1 */
118		#if HAVE_NATIVE_mpn_addlsh2_n
119		cy = mpn_addlsh2_n (b0b2, b0, b2, t); /* 4b2 + b0 */
120		#else
121	445	cy = mpn_lshift (b0b2, b2, t, 2); /* 4b2 */
122	445	cy += mpn_add_n (b0b2, b0b2, b0, t); /* 4b2 + b0 */
123	445	#endif
124	445	if (t != n)
125	334	cy = mpn_add_1 (b0b2 + t, b0 + t, n - t, cy);
126	445	b0b2[n] = cy;
127
128		#if HAVE_NATIVE_mpn_add_n_sub_n
129		if (mpn_cmp (b0b2, b1d, n+1) < 0)
130		{
131		mpn_add_n_sub_n (bs2, bsm2, b1d, b0b2, n+1);
132		flags = (enum toom6_flags) (flags ^ toom6_vm2_neg);
133		}
134		else
135		{
136		mpn_add_n_sub_n (bs2, bsm2, b0b2, b1d, n+1);
137		}
138		#else
139	445	mpn_add_n (bs2, b0b2, b1d, n+1);
140	445	if (mpn_cmp (b0b2, b1d, n+1) < 0)
141	256	{
142	256	mpn_sub_n (bsm2, b1d, b0b2, n+1);
143	256	flags = (enum toom6_flags) (flags ^ toom6_vm2_neg);
144	256	}
145	189	else
146	189	{
147	189	mpn_sub_n (bsm2, b0b2, b1d, n+1);
148	189	}
149	445	#endif
150
151		/* Compute as1 and asm1. */
152	445	flags = (enum toom6_flags) (flags ^ (toom6_vm1_neg & mpn_toom_eval_dgr3_pm1 (as1, asm1, ap, n, s, a0a2)));
153
154		/* Compute bs1 and bsm1. */
155	445	bsm1[n] = mpn_add (bsm1, b0, n, b2, t);
156		#if HAVE_NATIVE_mpn_add_n_sub_n
157		if (bsm1[n] == 0 && mpn_cmp (bsm1, b1, n) < 0)
158		{
159		cy = mpn_add_n_sub_n (bs1, bsm1, b1, bsm1, n);
160		bs1[n] = cy >> 1;
161		flags = (enum toom6_flags) (flags ^ toom6_vm1_neg);
162		}
163		else
164		{
165		cy = mpn_add_n_sub_n (bs1, bsm1, bsm1, b1, n);
166		bs1[n] = bsm1[n] + (cy >> 1);
167		bsm1[n]-= cy & 1;
168		}
169		#else
170	445	bs1[n] = bsm1[n] + mpn_add_n (bs1, bsm1, b1, n);
171	445	if (bsm1[n] == 0 && mpn_cmp (bsm1, b1, n) < 0)
172	196	{
173	196	mpn_sub_n (bsm1, b1, bsm1, n);
174	196	flags = (enum toom6_flags) (flags ^ toom6_vm1_neg);
175	196	}
176	249	else
177	249	{
178	249	bsm1[n] -= mpn_sub_n (bsm1, bsm1, b1, n);
179	249	}
180	445	#endif
181
182	445	ASSERT (as1[n] <= 3);
183	445	ASSERT (bs1[n] <= 2);
184	445	ASSERT (asm1[n] <= 1);
185	445	ASSERT (bsm1[n] <= 1);
186	445	ASSERT (as2[n] <=14);
187	445	ASSERT (bs2[n] <= 6);
188	445	ASSERT (asm2[n] <= 9);
189	445	ASSERT (bsm2[n] <= 4);
190
191		/* vm1, 2n+1 limbs */
192	445	vm1[2*n] = 0;
193	445	mpn_mul_n (vm1, asm1, bsm1, n + (asm1[n] \| bsm1[n])); /* W4 */
194
195		/* vm2, 2n+1 limbs */
196	445	mpn_mul_n (vm2, asm2, bsm2, n+1); /* W2 */
197
198		/* v2, 2n+1 limbs */
199	445	mpn_mul_n (v2, as2, bs2, n+1); /* W1 */
200
201		/* v1, 2n+1 limbs */
202	445	mpn_mul_n (v1, as1, bs1, n+1); /* W3 */
203
204		/* vinf, s+t limbs / / W0 */
205	445	if (s > t) mpn_mul (vinf, a3, s, b2, t);
206	221	else mpn_mul (vinf, b2, t, a3, s);
207
208		/* v0, 2n limbs */
209	445	mpn_mul_n (v0, ap, bp, n); /* W5 */
210
211	445	mpn_toom_interpolate_6pts (pp, n, flags, vm1, vm2, v2, t + s);
212
213	445	#undef v0
214	445	#undef vm1
215	445	#undef v1
216	445	#undef vm2
217	445	#undef v2
218	445	#undef vinf
219	445	#undef bs1
220	445	#undef bs2
221	445	#undef bsm1
222	445	#undef bsm2
223	445	#undef asm1
224	445	#undef asm2
225		/* #undef as1 */
226		/* #undef as2 */
227	445	#undef a0a2
228	445	#undef b0b2
229	445	#undef a1a3
230	445	#undef b1d
231	445	#undef a0
232	445	#undef a1
233	445	#undef a2
234	445	#undef a3
235	445	#undef b0
236	445	#undef b1
237	445	#undef b2
238	445	}