/src/gmp-6.2.1/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 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      */
  cy  = mpn_lshift (b0b2, b2, t, 2);      /*  4b2           */
  cy += mpn_add_n (b0b2, b0b2, b0, t);     /*  4b2      + b0 */
  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 */
  mpn_mul_n (vm1, asm1, bsm1, n+1);  /* 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: 2025-04-11 06:45

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 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	345	{
65	345	mp_size_t n, s, t;
66	345	enum toom6_flags flags;
67	345	mp_limb_t cy;
68
69	345	#define a0 ap
70	345	#define a1 (ap + n)
71	345	#define a2 (ap + 2 * n)
72	345	#define a3 (ap + 3 * n)
73	948	#define b0 bp
74	1.37k	#define b1 (bp + n)
75	1.03k	#define b2 (bp + 2 * n)
76
77	345	n = 1 + (3 * an >= 4 * bn ? (an - 1) >> 2 : (bn - 1) / (size_t) 3);
78
79	345	s = an - 3 * n;
80	345	t = bn - 2 * n;
81
82	345	ASSERT (0 < s && s <= n);
83	345	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	345	ASSERT (s+t >= 5);
89
90	345	#define v0 pp /* 2n */
91	690	#define vm1 (scratch) /* 2n+1 */
92	345	#define v1 (pp + 2n) / 2n+1 */
93	690	#define vm2 (scratch + 2 * n + 1) /* 2n+1 */
94	690	#define v2 (scratch + 4 * n + 2) /* 2n+1 */
95	345	#define vinf (pp + 5 * n) /* s+t */
96	1.03k	#define bs1 pp /* n+1 */
97	4.98k	#define bsm1 (scratch + 2 * n + 2) /* n+1 */
98	1.03k	#define asm1 (scratch + 3 * n + 3) /* n+1 */
99	690	#define asm2 (scratch + 4 * n + 4) /* n+1 */
100	690	#define bsm2 (pp + n + 1) /* n+1 */
101	690	#define bs2 (pp + 2 * n + 2) /* n+1 */
102	690	#define as2 (pp + 3 * n + 3) /* n+1 */
103	690	#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	345	#define a0a2 scratch
109	2.67k	#define b0b2 scratch
110	345	#define a1a3 asm1
111	1.72k	#define b1d bsm1
112
113		/* Compute as2 and asm2. */
114	345	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	345	b1d[n] = mpn_lshift (b1d, b1, n, 1); /* 2b1 */
118	345	cy = mpn_lshift (b0b2, b2, t, 2); /* 4b2 */
119	345	cy += mpn_add_n (b0b2, b0b2, b0, t); /* 4b2 + b0 */
120	345	if (t != n)
121	258	cy = mpn_add_1 (b0b2 + t, b0 + t, n - t, cy);
122	345	b0b2[n] = cy;
123
124		#if HAVE_NATIVE_mpn_add_n_sub_n
125		if (mpn_cmp (b0b2, b1d, n+1) < 0)
126		{
127		mpn_add_n_sub_n (bs2, bsm2, b1d, b0b2, n+1);
128		flags = (enum toom6_flags) (flags ^ toom6_vm2_neg);
129		}
130		else
131		{
132		mpn_add_n_sub_n (bs2, bsm2, b0b2, b1d, n+1);
133		}
134		#else
135	345	mpn_add_n (bs2, b0b2, b1d, n+1);
136	345	if (mpn_cmp (b0b2, b1d, n+1) < 0)
137	249	{
138	249	mpn_sub_n (bsm2, b1d, b0b2, n+1);
139	249	flags = (enum toom6_flags) (flags ^ toom6_vm2_neg);
140	249	}
141	96	else
142	96	{
143	96	mpn_sub_n (bsm2, b0b2, b1d, n+1);
144	96	}
145	345	#endif
146
147		/* Compute as1 and asm1. */
148	345	flags = (enum toom6_flags) (flags ^ (toom6_vm1_neg & mpn_toom_eval_dgr3_pm1 (as1, asm1, ap, n, s, a0a2)));
149
150		/* Compute bs1 and bsm1. */
151	345	bsm1[n] = mpn_add (bsm1, b0, n, b2, t);
152		#if HAVE_NATIVE_mpn_add_n_sub_n
153		if (bsm1[n] == 0 && mpn_cmp (bsm1, b1, n) < 0)
154		{
155		cy = mpn_add_n_sub_n (bs1, bsm1, b1, bsm1, n);
156		bs1[n] = cy >> 1;
157		flags = (enum toom6_flags) (flags ^ toom6_vm1_neg);
158		}
159		else
160		{
161		cy = mpn_add_n_sub_n (bs1, bsm1, bsm1, b1, n);
162		bs1[n] = bsm1[n] + (cy >> 1);
163		bsm1[n]-= cy & 1;
164		}
165		#else
166	345	bs1[n] = bsm1[n] + mpn_add_n (bs1, bsm1, b1, n);
167	345	if (bsm1[n] == 0 && mpn_cmp (bsm1, b1, n) < 0)
168	178	{
169	178	mpn_sub_n (bsm1, b1, bsm1, n);
170	178	flags = (enum toom6_flags) (flags ^ toom6_vm1_neg);
171	178	}
172	167	else
173	167	{
174	167	bsm1[n] -= mpn_sub_n (bsm1, bsm1, b1, n);
175	167	}
176	345	#endif
177
178	345	ASSERT (as1[n] <= 3);
179	345	ASSERT (bs1[n] <= 2);
180	345	ASSERT (asm1[n] <= 1);
181	345	ASSERT (bsm1[n] <= 1);
182	345	ASSERT (as2[n] <=14);
183	345	ASSERT (bs2[n] <= 6);
184	345	ASSERT (asm2[n] <= 9);
185	345	ASSERT (bsm2[n] <= 4);
186
187		/* vm1, 2n+1 limbs */
188	345	mpn_mul_n (vm1, asm1, bsm1, n+1); /* W4 */
189
190		/* vm2, 2n+1 limbs */
191	345	mpn_mul_n (vm2, asm2, bsm2, n+1); /* W2 */
192
193		/* v2, 2n+1 limbs */
194	345	mpn_mul_n (v2, as2, bs2, n+1); /* W1 */
195
196		/* v1, 2n+1 limbs */
197	345	mpn_mul_n (v1, as1, bs1, n+1); /* W3 */
198
199		/* vinf, s+t limbs / / W0 */
200	345	if (s > t) mpn_mul (vinf, a3, s, b2, t);
201	307	else mpn_mul (vinf, b2, t, a3, s);
202
203		/* v0, 2n limbs */
204	345	mpn_mul_n (v0, ap, bp, n); /* W5 */
205
206	345	mpn_toom_interpolate_6pts (pp, n, flags, vm1, vm2, v2, t + s);
207
208	345	#undef v0
209	345	#undef vm1
210	345	#undef v1
211	345	#undef vm2
212	345	#undef v2
213	345	#undef vinf
214	345	#undef bs1
215	345	#undef bs2
216	345	#undef bsm1
217	345	#undef bsm2
218	345	#undef asm1
219	345	#undef asm2
220		/* #undef as1 */
221		/* #undef as2 */
222	345	#undef a0a2
223	345	#undef b0b2
224	345	#undef a1a3
225	345	#undef b1d
226	345	#undef a0
227	345	#undef a1
228	345	#undef a2
229	345	#undef a3
230	345	#undef b0
231	345	#undef b1
232	345	#undef b2
233	345	}