mpn/generic/mullo_n.c

1 /* mpn_mullo_n -- multiply two n-limb numbers and return the low n limbs
25 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
32 #define MULLO_BASECASE_THRESHOLD 0	/* never use mpn_mul_basecase */
58   Compute the least significant half of the product {xy,n}*{yp,n}, or
59   formally {rp,n} = {xy,n}*{yp,n} Mod (B^n).
63   are used to obtain the final result. The more natural strategy is to
68   split n = n1 + n2, where an = n1 <= n2 = (1-a)n; i.e. 0 < a <= 1/2.
71   given size ML(n) is a fraction of the cost of a full product with
72   same size M(n), and the cost M(n)=n^e for some exponent 1 < e <= 2;
75   ML(n) = 2*ML(an) + M((1-a)n) => k*M(n) = 2*k*M(n)*a^e + M(n)*(1-a)^e
112   ^ n2 ^n1^
114   For an actual implementation, the assumption that M(n)=n^e is
115   incorrect, as a consequence also the assumption that ML(n)=k*M(n)
129 mpn_mullo_n_itch (mp_size_t n)
131   return 2*n;
135     mpn_dc_mullo_n requires a scratch space of 2*n limbs at tp.
139 mpn_dc_mullo_n (mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t n, mp_ptr tp)
141   mp_size_t n2, n1;
142   ASSERT (n >= 2);
143   ASSERT (! MPN_OVERLAP_P (rp, n, xp, n));
144   ASSERT (! MPN_OVERLAP_P (rp, n, yp, n));
145   ASSERT (MPN_SAME_OR_SEPARATE2_P(rp, n, tp, 2*n));
149   /* We need fractional approximation of the value 0 < a <= 1/2
152   if (MAYBE_range_basecase && BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD*36/(36-11)))
153     n1 = n >> 1;
154   else if (MAYBE_range_toom22 && BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD*36/(36-11)))
155     n1 = n * 11 / (size_t) 36;	/* n1 ~= n*(1-.694...) */
156   else if (BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD*40/(40-9)))
157     n1 = n * 9 / (size_t) 40;	/* n1 ~= n*(1-.775...) */
158   else if (BELOW_THRESHOLD (n, MUL_TOOM8H_THRESHOLD*10/9))
159     n1 = n * 7 / (size_t) 39;	/* n1 ~= n*(1-.821...) */
160   /* n1 = n * 4 / (size_t) 31;	// n1 ~= n*(1-.871...) [TOOM66] */
162     n1 = n / (size_t) 10;		/* n1 ~= n*(1-.899...) [TOOM88] */
164   n2 = n - n1;
166   /* Split as x = x1 2^(n2 GMP_NUMB_BITS) + x0,
167 	      y = y1 2^(n2 GMP_NUMB_BITS) + y0 */
170   mpn_mul_n (tp, xp, yp, n2);
171   MPN_COPY (rp, tp, n2);
173   /* x1 * y0 * 2^(n2 GMP_NUMB_BITS) */
174   if (BELOW_THRESHOLD (n1, MULLO_BASECASE_THRESHOLD))
175     mpn_mul_basecase (tp + n, xp + n2, n1, yp, n1);
176   else if (BELOW_THRESHOLD (n1, MULLO_DC_THRESHOLD))
177     mpn_mullo_basecase (tp + n, xp + n2, yp, n1);
179     mpn_dc_mullo_n (tp + n, xp + n2, yp, n1, tp + n);
180   mpn_add_n (rp + n2, tp + n2, tp + n, n1);
182   /* x0 * y1 * 2^(n2 GMP_NUMB_BITS) */
183   if (BELOW_THRESHOLD (n1, MULLO_BASECASE_THRESHOLD))
184     mpn_mul_basecase (tp + n, xp, n1, yp + n2, n1);
185   else if (BELOW_THRESHOLD (n1, MULLO_DC_THRESHOLD))
186     mpn_mullo_basecase (tp + n, xp, yp + n2, n1);
188     mpn_dc_mullo_n (tp + n, xp, yp + n2, n1, tp + n);
189   mpn_add_n (rp + n2, rp + n2, tp + n, n1);
205 mpn_mullo_n (mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t n)
207   ASSERT (n >= 1);
208   ASSERT (! MPN_OVERLAP_P (rp, n, xp, n));
209   ASSERT (! MPN_OVERLAP_P (rp, n, yp, n));
211   if (BELOW_THRESHOLD (n, MULLO_BASECASE_THRESHOLD))
215       mpn_mul_basecase (tp, xp, n, yp, n);
216       MPN_COPY (rp, tp, n);
218   else if (BELOW_THRESHOLD (n, MULLO_DC_THRESHOLD))
220       mpn_mullo_basecase (rp, xp, yp, n);
227       tp = TMP_ALLOC_LIMBS (mpn_mullo_n_itch (n));
228       if (BELOW_THRESHOLD (n, MULLO_MUL_N_THRESHOLD))
230 	  mpn_dc_mullo_n (rp, xp, yp, n, tp);
234 	  /* For really large operands, use plain mpn_mul_n but throw away upper n
237 	  mpn_fft_mul (tp, xp, n, yp, n);
239 	  mpn_mul_n (tp, xp, yp, n);
241 	  MPN_COPY (rp, tp, n);