Coverage Report

Created: 2023-02-22 06:39

/src/gmp-6.2.1/mpn/mullo_n.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_mullo_n -- multiply two n-limb numbers and return the low n limbs
2
   of their products.
3
4
   Contributed to the GNU project by Torbjorn Granlund and Marco Bodrato.
5
6
   THIS IS (FOR NOW) AN INTERNAL FUNCTION.  IT IS ONLY SAFE TO REACH THIS
7
   FUNCTION THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST GUARANTEED
8
   THAT IT'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
9
10
Copyright 2004, 2005, 2009, 2010, 2012 Free Software Foundation, Inc.
11
12
This file is part of the GNU MP Library.
13
14
The GNU MP Library is free software; you can redistribute it and/or modify
15
it under the terms of either:
16
17
  * the GNU Lesser General Public License as published by the Free
18
    Software Foundation; either version 3 of the License, or (at your
19
    option) any later version.
20
21
or
22
23
  * the GNU General Public License as published by the Free Software
24
    Foundation; either version 2 of the License, or (at your option) any
25
    later version.
26
27
or both in parallel, as here.
28
29
The GNU MP Library is distributed in the hope that it will be useful, but
30
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
32
for more details.
33
34
You should have received copies of the GNU General Public License and the
35
GNU Lesser General Public License along with the GNU MP Library.  If not,
36
see https://www.gnu.org/licenses/.  */
37
38
#include "gmp-impl.h"
39
40
41
#if TUNE_PROGRAM_BUILD || WANT_FAT_BINARY
42
#define MAYBE_range_basecase 1
43
#define MAYBE_range_toom22   1
44
#else
45
#define MAYBE_range_basecase                                           \
46
546
  ((MULLO_DC_THRESHOLD == 0 ? MULLO_BASECASE_THRESHOLD : MULLO_DC_THRESHOLD) < MUL_TOOM22_THRESHOLD*36/(36-11))
47
#define MAYBE_range_toom22                                             \
48
546
  ((MULLO_DC_THRESHOLD == 0 ? MULLO_BASECASE_THRESHOLD : MULLO_DC_THRESHOLD) < MUL_TOOM33_THRESHOLD*36/(36-11) )
49
#endif
50
51
/*  THINK: The DC strategy uses different constants in different Toom's
52
   ranges. Something smoother?
53
*/
54
55
/*
56
  Compute the least significant half of the product {xy,n}*{yp,n}, or
57
  formally {rp,n} = {xy,n}*{yp,n} Mod (B^n).
58
59
  Above the given threshold, the Divide and Conquer strategy is used.
60
  The operands are split in two, and a full product plus two mullo
61
  are used to obtain the final result. The more natural strategy is to
62
  split in two halves, but this is far from optimal when a
63
  sub-quadratic multiplication is used.
64
65
  Mulders suggests an unbalanced split in favour of the full product,
66
  split n = n1 + n2, where an = n1 <= n2 = (1-a)n; i.e. 0 < a <= 1/2.
67
68
  To compute the value of a, we assume that the cost of mullo for a
69
  given size ML(n) is a fraction of the cost of a full product with
70
  same size M(n), and the cost M(n)=n^e for some exponent 1 < e <= 2;
71
  then we can write:
72
73
  ML(n) = 2*ML(an) + M((1-a)n) => k*M(n) = 2*k*M(n)*a^e + M(n)*(1-a)^e
74
75
  Given a value for e, want to minimise the value of k, i.e. the
76
  function k=(1-a)^e/(1-2*a^e).
77
78
  With e=2, the exponent for schoolbook multiplication, the minimum is
79
  given by the values a=1-a=1/2.
80
81
  With e=log(3)/log(2), the exponent for Karatsuba (aka toom22),
82
  Mulders compute (1-a) = 0.694... and we approximate a with 11/36.
83
84
  Other possible approximations follow:
85
  e=log(5)/log(3) [Toom-3] -> a ~= 9/40
86
  e=log(7)/log(4) [Toom-4] -> a ~= 7/39
87
  e=log(11)/log(6) [Toom-6] -> a ~= 1/8
88
  e=log(15)/log(8) [Toom-8] -> a ~= 1/10
89
90
  The values above where obtained with the following trivial commands
91
  in the gp-pari shell:
92
93
fun(e,a)=(1-a)^e/(1-2*a^e)
94
mul(a,b,c)={local(m,x,p);if(b-c<1/10000,(b+c)/2,m=1;x=b;forstep(p=c,b,(b-c)/8,if(fun(a,p)<m,m=fun(a,p);x=p));mul(a,(b+x)/2,(c+x)/2))}
95
contfracpnqn(contfrac(mul(log(2*2-1)/log(2),1/2,0),5))
96
contfracpnqn(contfrac(mul(log(3*2-1)/log(3),1/2,0),5))
97
contfracpnqn(contfrac(mul(log(4*2-1)/log(4),1/2,0),5))
98
contfracpnqn(contfrac(mul(log(6*2-1)/log(6),1/2,0),3))
99
contfracpnqn(contfrac(mul(log(8*2-1)/log(8),1/2,0),3))
100
101
  ,
102
  |\
103
  | \
104
  +----,
105
  |    |
106
  |    |
107
  |    |\
108
  |    | \
109
  +----+--`
110
  ^ n2 ^n1^
111
112
  For an actual implementation, the assumption that M(n)=n^e is
113
  incorrect, as a consequence also the assumption that ML(n)=k*M(n)
114
  with a constant k is wrong.
115
116
  But theory suggest us two things:
117
  - the best the multiplication product is (lower e), the more k
118
    approaches 1, and a approaches 0.
119
120
  - A value for a smaller than optimal is probably less bad than a
121
    bigger one: e.g. let e=log(3)/log(2), a=0.3058_ the optimal
122
    value, and k(a)=0.808_ the mul/mullo speed ratio. We get
123
    k(a+1/6)=0.929_ but k(a-1/6)=0.865_.
124
*/
125
126
static mp_size_t
127
mpn_mullo_n_itch (mp_size_t n)
128
546
{
129
546
  return 2*n;
130
546
}
131
132
/*
133
    mpn_dc_mullo_n requires a scratch space of 2*n limbs at tp.
134
    It accepts tp == rp.
135
*/
136
static void
137
mpn_dc_mullo_n (mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t n, mp_ptr tp)
138
273
{
139
273
  mp_size_t n2, n1;
140
273
  ASSERT (n >= 2);
141
273
  ASSERT (! MPN_OVERLAP_P (rp, n, xp, n));
142
273
  ASSERT (! MPN_OVERLAP_P (rp, n, yp, n));
143
273
  ASSERT (MPN_SAME_OR_SEPARATE2_P(rp, n, tp, 2*n));
144
145
  /* Divide-and-conquer */
146
147
  /* We need fractional approximation of the value 0 < a <= 1/2
148
     giving the minimum in the function k=(1-a)^e/(1-2*a^e).
149
  */
150
273
  if (MAYBE_range_basecase && BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD*36/(36-11)))
151
0
    n1 = n >> 1;
152
273
  else if (MAYBE_range_toom22 && BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD*36/(36-11)))
153
273
    n1 = n * 11 / (size_t) 36;  /* n1 ~= n*(1-.694...) */
154
0
  else if (BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD*40/(40-9)))
155
0
    n1 = n * 9 / (size_t) 40; /* n1 ~= n*(1-.775...) */
156
0
  else if (BELOW_THRESHOLD (n, MUL_TOOM8H_THRESHOLD*10/9))
157
0
    n1 = n * 7 / (size_t) 39; /* n1 ~= n*(1-.821...) */
158
  /* n1 = n * 4 / (size_t) 31;  // n1 ~= n*(1-.871...) [TOOM66] */
159
0
  else
160
0
    n1 = n / (size_t) 10;   /* n1 ~= n*(1-.899...) [TOOM88] */
161
162
273
  n2 = n - n1;
163
164
  /* Split as x = x1 2^(n2 GMP_NUMB_BITS) + x0,
165
        y = y1 2^(n2 GMP_NUMB_BITS) + y0 */
166
167
  /* x0 * y0 */
168
273
  mpn_mul_n (tp, xp, yp, n2);
169
273
  MPN_COPY (rp, tp, n2);
170
171
  /* x1 * y0 * 2^(n2 GMP_NUMB_BITS) */
172
273
  if (BELOW_THRESHOLD (n1, MULLO_BASECASE_THRESHOLD))
173
0
    mpn_mul_basecase (tp + n, xp + n2, n1, yp, n1);
174
273
  else if (BELOW_THRESHOLD (n1, MULLO_DC_THRESHOLD))
175
273
    mpn_mullo_basecase (tp + n, xp + n2, yp, n1);
176
0
  else
177
0
    mpn_dc_mullo_n (tp + n, xp + n2, yp, n1, tp + n);
178
273
  mpn_add_n (rp + n2, tp + n2, tp + n, n1);
179
180
  /* x0 * y1 * 2^(n2 GMP_NUMB_BITS) */
181
273
  if (BELOW_THRESHOLD (n1, MULLO_BASECASE_THRESHOLD))
182
0
    mpn_mul_basecase (tp + n, xp, n1, yp + n2, n1);
183
273
  else if (BELOW_THRESHOLD (n1, MULLO_DC_THRESHOLD))
184
273
    mpn_mullo_basecase (tp + n, xp, yp + n2, n1);
185
0
  else
186
0
    mpn_dc_mullo_n (tp + n, xp, yp + n2, n1, tp + n);
187
273
  mpn_add_n (rp + n2, rp + n2, tp + n, n1);
188
273
}
189
190
/* Avoid zero allocations when MULLO_BASECASE_THRESHOLD is 0.  */
191
#define MUL_BASECASE_ALLOC \
192
 (MULLO_BASECASE_THRESHOLD_LIMIT == 0 ? 1 : 2*MULLO_BASECASE_THRESHOLD_LIMIT)
193
194
/* FIXME: This function should accept a temporary area; dc_mullow_n
195
   accepts a pointer tp, and handle the case tp == rp, do the same here.
196
   Maybe recombine the two functions.
197
   THINK: If mpn_mul_basecase is always faster than mpn_mullo_basecase
198
    (typically thanks to mpn_addmul_2) should we unconditionally use
199
    mpn_mul_n?
200
*/
201
202
void
203
mpn_mullo_n (mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t n)
204
422
{
205
422
  ASSERT (n >= 1);
206
422
  ASSERT (! MPN_OVERLAP_P (rp, n, xp, n));
207
422
  ASSERT (! MPN_OVERLAP_P (rp, n, yp, n));
208
209
422
  if (BELOW_THRESHOLD (n, MULLO_BASECASE_THRESHOLD))
210
0
    {
211
      /* Allocate workspace of fixed size on stack: fast! */
212
0
      mp_limb_t tp[MUL_BASECASE_ALLOC];
213
0
      mpn_mul_basecase (tp, xp, n, yp, n);
214
0
      MPN_COPY (rp, tp, n);
215
0
    }
216
422
  else if (BELOW_THRESHOLD (n, MULLO_DC_THRESHOLD))
217
149
    {
218
149
      mpn_mullo_basecase (rp, xp, yp, n);
219
149
    }
220
273
  else
221
273
    {
222
273
      mp_ptr tp;
223
273
      TMP_DECL;
224
273
      TMP_MARK;
225
273
      tp = TMP_ALLOC_LIMBS (mpn_mullo_n_itch (n));
226
273
      if (BELOW_THRESHOLD (n, MULLO_MUL_N_THRESHOLD))
227
273
  {
228
273
    mpn_dc_mullo_n (rp, xp, yp, n, tp);
229
273
  }
230
0
      else
231
0
  {
232
    /* For really large operands, use plain mpn_mul_n but throw away upper n
233
       limbs of result.  */
234
0
#if !TUNE_PROGRAM_BUILD && (MULLO_MUL_N_THRESHOLD > MUL_FFT_THRESHOLD)
235
0
    mpn_fft_mul (tp, xp, n, yp, n);
236
#else
237
    mpn_mul_n (tp, xp, yp, n);
238
#endif
239
0
    MPN_COPY (rp, tp, n);
240
0
  }
241
273
      TMP_FREE;
242
273
    }
243
422
}