/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: 2023-09-25 06:33

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