--- /dev/null
+/* mpn_mul_n and helper function -- Multiply/square natural numbers.
+
+ THE HELPER FUNCTIONS IN THIS FILE (meaning everything except mpn_mul_n) ARE
+ INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY SAFE TO REACH THEM THROUGH
+ DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE
+ OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
+
+Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003,
+2005 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 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.
+
+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 Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+
+/* Multiplies using 3 half-sized mults and so on recursively.
+ * p[0..2*n-1] := product of a[0..n-1] and b[0..n-1].
+ * No overlap of p[...] with a[...] or b[...].
+ * ws is workspace.
+ */
+
+void
+mpn_kara_mul_n (mp_ptr p, mp_srcptr a, mp_srcptr b, mp_size_t n, mp_ptr ws)
+{
+ mp_limb_t w, w0, w1;
+ mp_size_t n2;
+ mp_srcptr x, y;
+ mp_size_t i;
+ int sign;
+
+ n2 = n >> 1;
+ ASSERT (n2 > 0);
+
+ if ((n & 1) != 0)
+ {
+ /* Odd length. */
+ mp_size_t n1, n3, nm1;
+
+ n3 = n - n2;
+
+ sign = 0;
+ w = a[n2];
+ if (w != 0)
+ w -= mpn_sub_n (p, a, a + n3, n2);
+ else
+ {
+ i = n2;
+ do
+ {
+ --i;
+ w0 = a[i];
+ w1 = a[n3 + i];
+ }
+ while (w0 == w1 && i != 0);
+ if (w0 < w1)
+ {
+ x = a + n3;
+ y = a;
+ sign = ~0;
+ }
+ else
+ {
+ x = a;
+ y = a + n3;
+ }
+ mpn_sub_n (p, x, y, n2);
+ }
+ p[n2] = w;
+
+ w = b[n2];
+ if (w != 0)
+ w -= mpn_sub_n (p + n3, b, b + n3, n2);
+ else
+ {
+ i = n2;
+ do
+ {
+ --i;
+ w0 = b[i];
+ w1 = b[n3 + i];
+ }
+ while (w0 == w1 && i != 0);
+ if (w0 < w1)
+ {
+ x = b + n3;
+ y = b;
+ sign = ~sign;
+ }
+ else
+ {
+ x = b;
+ y = b + n3;
+ }
+ mpn_sub_n (p + n3, x, y, n2);
+ }
+ p[n] = w;
+
+ n1 = n + 1;
+ if (n2 < MUL_KARATSUBA_THRESHOLD)
+ {
+ if (n3 < MUL_KARATSUBA_THRESHOLD)
+ {
+ mpn_mul_basecase (ws, p, n3, p + n3, n3);
+ mpn_mul_basecase (p, a, n3, b, n3);
+ }
+ else
+ {
+ mpn_kara_mul_n (ws, p, p + n3, n3, ws + n1);
+ mpn_kara_mul_n (p, a, b, n3, ws + n1);
+ }
+ mpn_mul_basecase (p + n1, a + n3, n2, b + n3, n2);
+ }
+ else
+ {
+ mpn_kara_mul_n (ws, p, p + n3, n3, ws + n1);
+ mpn_kara_mul_n (p, a, b, n3, ws + n1);
+ mpn_kara_mul_n (p + n1, a + n3, b + n3, n2, ws + n1);
+ }
+
+ if (sign)
+ mpn_add_n (ws, p, ws, n1);
+ else
+ mpn_sub_n (ws, p, ws, n1);
+
+ nm1 = n - 1;
+ if (mpn_add_n (ws, p + n1, ws, nm1))
+ {
+ mp_limb_t x = (ws[nm1] + 1) & GMP_NUMB_MASK;
+ ws[nm1] = x;
+ if (x == 0)
+ ws[n] = (ws[n] + 1) & GMP_NUMB_MASK;
+ }
+ if (mpn_add_n (p + n3, p + n3, ws, n1))
+ {
+ mpn_incr_u (p + n1 + n3, 1);
+ }
+ }
+ else
+ {
+ /* Even length. */
+ i = n2;
+ do
+ {
+ --i;
+ w0 = a[i];
+ w1 = a[n2 + i];
+ }
+ while (w0 == w1 && i != 0);
+ sign = 0;
+ if (w0 < w1)
+ {
+ x = a + n2;
+ y = a;
+ sign = ~0;
+ }
+ else
+ {
+ x = a;
+ y = a + n2;
+ }
+ mpn_sub_n (p, x, y, n2);
+
+ i = n2;
+ do
+ {
+ --i;
+ w0 = b[i];
+ w1 = b[n2 + i];
+ }
+ while (w0 == w1 && i != 0);
+ if (w0 < w1)
+ {
+ x = b + n2;
+ y = b;
+ sign = ~sign;
+ }
+ else
+ {
+ x = b;
+ y = b + n2;
+ }
+ mpn_sub_n (p + n2, x, y, n2);
+
+ /* Pointwise products. */
+ if (n2 < MUL_KARATSUBA_THRESHOLD)
+ {
+ mpn_mul_basecase (ws, p, n2, p + n2, n2);
+ mpn_mul_basecase (p, a, n2, b, n2);
+ mpn_mul_basecase (p + n, a + n2, n2, b + n2, n2);
+ }
+ else
+ {
+ mpn_kara_mul_n (ws, p, p + n2, n2, ws + n);
+ mpn_kara_mul_n (p, a, b, n2, ws + n);
+ mpn_kara_mul_n (p + n, a + n2, b + n2, n2, ws + n);
+ }
+
+ /* Interpolate. */
+ if (sign)
+ w = mpn_add_n (ws, p, ws, n);
+ else
+ w = -mpn_sub_n (ws, p, ws, n);
+ w += mpn_add_n (ws, p + n, ws, n);
+ w += mpn_add_n (p + n2, p + n2, ws, n);
+ MPN_INCR_U (p + n2 + n, 2 * n - (n2 + n), w);
+ }
+}
+
+void
+mpn_kara_sqr_n (mp_ptr p, mp_srcptr a, mp_size_t n, mp_ptr ws)
+{
+ mp_limb_t w, w0, w1;
+ mp_size_t n2;
+ mp_srcptr x, y;
+ mp_size_t i;
+
+ n2 = n >> 1;
+ ASSERT (n2 > 0);
+
+ if ((n & 1) != 0)
+ {
+ /* Odd length. */
+ mp_size_t n1, n3, nm1;
+
+ n3 = n - n2;
+
+ w = a[n2];
+ if (w != 0)
+ w -= mpn_sub_n (p, a, a + n3, n2);
+ else
+ {
+ i = n2;
+ do
+ {
+ --i;
+ w0 = a[i];
+ w1 = a[n3 + i];
+ }
+ while (w0 == w1 && i != 0);
+ if (w0 < w1)
+ {
+ x = a + n3;
+ y = a;
+ }
+ else
+ {
+ x = a;
+ y = a + n3;
+ }
+ mpn_sub_n (p, x, y, n2);
+ }
+ p[n2] = w;
+
+ n1 = n + 1;
+
+ /* n2 is always either n3 or n3-1 so maybe the two sets of tests here
+ could be combined. But that's not important, since the tests will
+ take a miniscule amount of time compared to the function calls. */
+ if (BELOW_THRESHOLD (n3, SQR_BASECASE_THRESHOLD))
+ {
+ mpn_mul_basecase (ws, p, n3, p, n3);
+ mpn_mul_basecase (p, a, n3, a, n3);
+ }
+ else if (BELOW_THRESHOLD (n3, SQR_KARATSUBA_THRESHOLD))
+ {
+ mpn_sqr_basecase (ws, p, n3);
+ mpn_sqr_basecase (p, a, n3);
+ }
+ else
+ {
+ mpn_kara_sqr_n (ws, p, n3, ws + n1); /* (x-y)^2 */
+ mpn_kara_sqr_n (p, a, n3, ws + n1); /* x^2 */
+ }
+ if (BELOW_THRESHOLD (n2, SQR_BASECASE_THRESHOLD))
+ mpn_mul_basecase (p + n1, a + n3, n2, a + n3, n2);
+ else if (BELOW_THRESHOLD (n2, SQR_KARATSUBA_THRESHOLD))
+ mpn_sqr_basecase (p + n1, a + n3, n2);
+ else
+ mpn_kara_sqr_n (p + n1, a + n3, n2, ws + n1); /* y^2 */
+
+ /* Since x^2+y^2-(x-y)^2 = 2xy >= 0 there's no need to track the
+ borrow from mpn_sub_n. If it occurs then it'll be cancelled by a
+ carry from ws[n]. Further, since 2xy fits in n1 limbs there won't
+ be any carry out of ws[n] other than cancelling that borrow. */
+
+ mpn_sub_n (ws, p, ws, n1); /* x^2-(x-y)^2 */
+
+ nm1 = n - 1;
+ if (mpn_add_n (ws, p + n1, ws, nm1)) /* x^2+y^2-(x-y)^2 = 2xy */
+ {
+ mp_limb_t x = (ws[nm1] + 1) & GMP_NUMB_MASK;
+ ws[nm1] = x;
+ if (x == 0)
+ ws[n] = (ws[n] + 1) & GMP_NUMB_MASK;
+ }
+ if (mpn_add_n (p + n3, p + n3, ws, n1))
+ {
+ mpn_incr_u (p + n1 + n3, 1);
+ }
+ }
+ else
+ {
+ /* Even length. */
+ i = n2;
+ do
+ {
+ --i;
+ w0 = a[i];
+ w1 = a[n2 + i];
+ }
+ while (w0 == w1 && i != 0);
+ if (w0 < w1)
+ {
+ x = a + n2;
+ y = a;
+ }
+ else
+ {
+ x = a;
+ y = a + n2;
+ }
+ mpn_sub_n (p, x, y, n2);
+
+ /* Pointwise products. */
+ if (BELOW_THRESHOLD (n2, SQR_BASECASE_THRESHOLD))
+ {
+ mpn_mul_basecase (ws, p, n2, p, n2);
+ mpn_mul_basecase (p, a, n2, a, n2);
+ mpn_mul_basecase (p + n, a + n2, n2, a + n2, n2);
+ }
+ else if (BELOW_THRESHOLD (n2, SQR_KARATSUBA_THRESHOLD))
+ {
+ mpn_sqr_basecase (ws, p, n2);
+ mpn_sqr_basecase (p, a, n2);
+ mpn_sqr_basecase (p + n, a + n2, n2);
+ }
+ else
+ {
+ mpn_kara_sqr_n (ws, p, n2, ws + n);
+ mpn_kara_sqr_n (p, a, n2, ws + n);
+ mpn_kara_sqr_n (p + n, a + n2, n2, ws + n);
+ }
+
+ /* Interpolate. */
+ w = -mpn_sub_n (ws, p, ws, n);
+ w += mpn_add_n (ws, p + n, ws, n);
+ w += mpn_add_n (p + n2, p + n2, ws, n);
+ MPN_INCR_U (p + n2 + n, 2 * n - (n2 + n), w);
+ }
+}
+
+/******************************************************************************
+ * *
+ * Toom 3-way multiplication and squaring *
+ * *
+ *****************************************************************************/
+
+/* Starts from:
+ {v0,2k} (stored in {c,2k})
+ {vm1,2k+1} (which sign is sa, and absolute value is stored in {vm1,2k+1})
+ {v1,2k+1} (stored in {c+2k,2k+1})
+ {v2,2k+1}
+ {vinf,twor} (stored in {c+4k,twor}, except the first limb, saved in vinf0)
+
+ ws is temporary space, and should have at least twor limbs.
+
+ put in {c, 2n} where n = 2k+twor the value of {v0,2k} (already in place)
+ + B^k * {tm1, 2k+1}
+ + B^(2k) * {t1, 2k+1}
+ + B^(3k) * {t2, 2k+1}
+ + B^(4k) * {vinf,twor} (high twor-1 limbs already in place)
+ where {t1, 2k+1} = ({v1, 2k+1} + sa * {vm1, 2k+1}- 2*{v0,2k})/2-*{vinf,twor}
+ {t2, 2k+1} = (3*({v1,2k+1}-{v0,2k})-sa*{vm1,2k+1}+{v2,2k+1})/6-2*{vinf,twor}
+ {tm1,2k+1} = ({v1,2k+1}-sa*{vm1,2k+1}/2-{t2,2k+1}
+
+ Exact sequence described in a comment in mpn_toom3_mul_n.
+ mpn_toom3_mul_n() and mpn_toom3_sqr_n() implement steps 1-2.
+ mpn_toom_interpolate_5pts() implements steps 3-4.
+
+ Reference: What About Toom-Cook Matrices Optimality? Marco Bodrato
+ and Alberto Zanoni, October 19, 2006, http://bodrato.it/papers/#CIVV2006
+
+ ************* saved note ****************
+ Think about:
+
+ The evaluated point a-b+c stands a good chance of having a zero carry
+ limb, a+b+c would have a 1/4 chance, and 4*a+2*b+c a 1/8 chance, roughly.
+ Perhaps this could be tested and stripped. Doing so before recursing
+ would be better than stripping at the start of mpn_toom3_mul_n/sqr_n,
+ since then the recursion could be based on the new size. Although in
+ truth the kara vs toom3 crossover is never so exact that one limb either
+ way makes a difference.
+
+ A small value like 1 or 2 for the carry could perhaps also be handled
+ with an add_n or addlsh1_n. Would that be faster than an extra limb on a
+ (recursed) multiply/square?
+*/
+
+#define TOOM3_MUL_REC(p, a, b, n, ws) \
+ do { \
+ if (MUL_TOOM3_THRESHOLD / 3 < MUL_KARATSUBA_THRESHOLD \
+ && BELOW_THRESHOLD (n, MUL_KARATSUBA_THRESHOLD)) \
+ mpn_mul_basecase (p, a, n, b, n); \
+ else if (BELOW_THRESHOLD (n, MUL_TOOM3_THRESHOLD)) \
+ mpn_kara_mul_n (p, a, b, n, ws); \
+ else \
+ mpn_toom3_mul_n (p, a, b, n, ws); \
+ } while (0)
+
+#define TOOM3_SQR_REC(p, a, n, ws) \
+ do { \
+ if (SQR_TOOM3_THRESHOLD / 3 < SQR_BASECASE_THRESHOLD \
+ && BELOW_THRESHOLD (n, SQR_BASECASE_THRESHOLD)) \
+ mpn_mul_basecase (p, a, n, a, n); \
+ else if (SQR_TOOM3_THRESHOLD / 3 < SQR_KARATSUBA_THRESHOLD \
+ && BELOW_THRESHOLD (n, SQR_KARATSUBA_THRESHOLD)) \
+ mpn_sqr_basecase (p, a, n); \
+ else if (BELOW_THRESHOLD (n, SQR_TOOM3_THRESHOLD)) \
+ mpn_kara_sqr_n (p, a, n, ws); \
+ else \
+ mpn_toom3_sqr_n (p, a, n, ws); \
+ } while (0)
+
+/* The necessary temporary space T(n) satisfies T(n)=0 for n < THRESHOLD,
+ and T(n) <= max(2n+2, 6k+3, 4k+3+T(k+1)) otherwise, where k = ceil(n/3).
+
+ Assuming T(n) >= 2n, 6k+3 <= 4k+3+T(k+1).
+ Similarly, 2n+2 <= 6k+2 <= 4k+3+T(k+1).
+
+ With T(n) = 2n+S(n), this simplifies to S(n) <= 9 + S(k+1).
+ Since THRESHOLD >= 17, we have n/(k+1) >= 19/8
+ thus S(n) <= S(n/(19/8)) + 9 thus S(n) <= 9*log(n)/log(19/8) <= 8*log2(n).
+*/
+
+void
+mpn_toom3_mul_n (mp_ptr c, mp_srcptr a, mp_srcptr b, mp_size_t n, mp_ptr t)
+{
+ mp_size_t k, k1, kk1, r, twok, twor;
+ mp_limb_t cy, cc, saved, vinf0;
+ mp_ptr trec;
+ int sa, sb;
+ mp_ptr c1, c2, c3, c4, c5;
+
+ ASSERT(GMP_NUMB_BITS >= 6);
+ ASSERT(n >= 17); /* so that r <> 0 and 5k+3 <= 2n */
+
+ /*
+ The algorithm is the following:
+
+ 0. k = ceil(n/3), r = n - 2k, B = 2^(GMP_NUMB_BITS), t = B^k
+ 1. split a and b in three parts each a0, a1, a2 and b0, b1, b2
+ with a0, a1, b0, b1 of k limbs, and a2, b2 of r limbs
+ 2. Evaluation: vm1 may be negative, the other can not.
+ v0 <- a0*b0
+ v1 <- (a0+a1+a2)*(b0+b1+b2)
+ v2 <- (a0+2*a1+4*a2)*(b0+2*b1+4*b2)
+ vm1 <- (a0-a1+a2)*(b0-b1+b2)
+ vinf <- a2*b2
+ 3. Interpolation: every result is positive, all divisions are exact
+ t2 <- (v2 - vm1)/3
+ tm1 <- (v1 - vm1)/2
+ t1 <- (v1 - v0)
+ t2 <- (t2 - t1)/2
+ t1 <- (t1 - tm1 - vinf)
+ t2 <- (t2 - 2*vinf)
+ tm1 <- (tm1 - t2)
+ 4. result is c0+c1*t+c2*t^2+c3*t^3+c4*t^4 where
+ c0 <- v0
+ c1 <- tm1
+ c2 <- t1
+ c3 <- t2
+ c4 <- vinf
+ */
+
+ k = (n + 2) / 3; /* ceil(n/3) */
+ twok = 2 * k;
+ k1 = k + 1;
+ kk1 = k + k1;
+ r = n - twok; /* last chunk */
+ twor = 2 * r;
+
+ c1 = c + k;
+ c2 = c1 + k;
+ c3 = c2 + k;
+ c4 = c3 + k;
+ c5 = c4 + k;
+
+ trec = t + 4 * k + 3; /* trec = v2 + (2k+2) */
+
+ /* put a0+a2 in {c, k+1}, and b0+b2 in {c+4k+2, k+1};
+ put a0+a1+a2 in {t, k+1} and b0+b1+b2 in {t+k+1,k+1}
+ [????requires 5k+3 <= 2n, ie. n >= 9] */
+ cy = mpn_add_n (c, a, a + twok, r);
+ cc = mpn_add_n (c4 + 2, b, b + twok, r);
+ if (r < k)
+ {
+ __GMPN_ADD_1 (cy, c + r, a + r, k - r, cy);
+ __GMPN_ADD_1 (cc, c4 + 2 + r, b + r, k - r, cc);
+ }
+
+ /* Put in {t, k+1} the sum
+ * (a_0+a_2) - stored in {c, k+1} -
+ * +
+ * a_1 - stored in {a+k, k} */
+ t[k] = (c1[0] = cy) + mpn_add_n (t, c, a + k, k);
+ /* ^ ^
+ * carry of a_0 + a_2 carry of (a_0+a_2) + a_1
+
+ */
+
+ /* Put in {t+k+1, k+1} the sum of the two values
+ * (b_0+b_2) - stored in {c1+1, k+1} -
+ * +
+ * b_1 - stored in {b+k, k} */
+ t[kk1] = (c5[3] = cc) + mpn_add_n (t + k1, c4 + 2, b + k, k);
+ /* ^ ^
+ * carry of b_0 + b_2 carry of (b_0+b_2) + b_1 */
+
+#define v2 (t+2*k+1)
+
+ /* compute v1 := (a0+a1+a2)*(b0+b1+b2) in {t, 2k+1};
+ since v1 < 9*B^(2k), v1 uses only 2k+1 words if GMP_NUMB_BITS >= 4 */
+ TOOM3_MUL_REC (c2, t, t + k1, k1, trec);
+
+ /* c c2 c4 t
+ {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
+ v1 */
+
+ /* put |a0-a1+a2| in {c, k+1} and |b0-b1+b2| in {c+4k+2,k+1} */
+ /* (They're already there, actually) */
+
+ /* sa = sign(a0-a1+a2) */
+ sa = (cy != 0) ? 1 : mpn_cmp (c, a + k, k);
+ c[k] = (sa >= 0) ? cy - mpn_sub_n (c, c, a + k, k)
+ : mpn_sub_n (c, a + k, c, k);
+
+ sb = (cc != 0) ? 1 : mpn_cmp (c4 + 2, b + k, k);
+ c5[2] = (sb >= 0) ? cc - mpn_sub_n (c4 + 2, c4 + 2, b + k, k)
+ : mpn_sub_n (c4 + 2, b + k, c4 + 2, k);
+ sa *= sb; /* sign of vm1 */
+
+ /* compute vm1 := (a0-a1+a2)*(b0-b1+b2) in {t, 2k+1};
+ since |vm1| < 4*B^(2k), vm1 uses only 2k+1 limbs */
+ TOOM3_MUL_REC (t, c, c4 + 2, k1, trec);
+
+ /* {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
+ v1 vm1
+ */
+
+ /* compute a0+2a1+4a2 in {c, k+1} and b0+2b1+4b2 in {c+4k+2, k+1}
+ [requires 5k+3 <= 2n, i.e. n >= 17] */
+#ifdef HAVE_NATIVE_mpn_addlsh1_n
+ c1[0] = mpn_addlsh1_n (c, a + k, a + twok, r);
+ c5[2] = mpn_addlsh1_n (c4 + 2, b + k, b + twok, r);
+ if (r < k)
+ {
+ __GMPN_ADD_1 (c1[0], c + r, a + k + r, k - r, c1[0]);
+ __GMPN_ADD_1 (c5[2], c4 + 2 + r, b + k + r, k - r, c5[2]);
+ }
+ c1[0] = 2 * c1[0] + mpn_addlsh1_n (c, a, c, k);
+ c5[2] = 2 * c5[2] + mpn_addlsh1_n (c4 + 2, b, c4 + 2, k);
+#else
+ c[r] = mpn_lshift (c, a + twok, r, 1);
+ c4[r + 2] = mpn_lshift (c4 + 2, b + twok, r, 1);
+ if (r < k)
+ {
+ MPN_ZERO(c + r + 1, k - r);
+ MPN_ZERO(c4 + r + 3, k - r);
+ }
+ c1[0] += mpn_add_n (c, c, a + k, k);
+ c5[2] += mpn_add_n (c4 + 2, c4 + 2, b + k, k);
+ mpn_lshift (c, c, k1, 1);
+ mpn_lshift (c4 + 2, c4 + 2, k1, 1);
+ c1[0] += mpn_add_n (c, c, a, k);
+ c5[2] += mpn_add_n (c4 + 2, c4 + 2, b, k);
+#endif
+
+ /* compute v2 := (a0+2a1+4a2)*(b0+2b1+4b2) in {t+2k+1, 2k+1}
+ v2 < 49*B^k so v2 uses at most 2k+1 limbs if GMP_NUMB_BITS >= 6 */
+ TOOM3_MUL_REC (v2, c, c4 + 2, k1, trec);
+
+ /* {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
+ v1 vm1 v2
+ */
+
+ /* compute v0 := a0*b0 in {c, 2k} */
+ TOOM3_MUL_REC (c, a, b, k, trec);
+
+ /* {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
+ v0 v1 vm1 v2 */
+
+ /* compute vinf := a2*b2 in {t+4k+2, 2r}: in {c4, 2r} */
+
+ saved = c4[0]; /* Remember v1's highest byte (will be overwritten). */
+ TOOM3_MUL_REC (c4, a + twok, b + twok, r, trec); /* Overwrites c4[0]. */
+ vinf0 = c4[0]; /* Remember vinf's lowest byte (will be overwritten).*/
+ c4[0] = saved; /* Overwriting. Now v1 value is correct. */
+
+ /* {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
+ v0 v1 vinf[1..] vm1 v2 */
+
+ mpn_toom_interpolate_5pts (c, v2, t, k, 2*r, sa, vinf0, trec);
+
+#undef v2
+}
+
+void
+mpn_toom3_sqr_n (mp_ptr c, mp_srcptr a, mp_size_t n, mp_ptr t)
+{
+ mp_size_t k, k1, kk1, r, twok, twor;
+ mp_limb_t cy, saved, vinf0;
+ mp_ptr trec;
+ int sa;
+ mp_ptr c1, c2, c3, c4;
+
+ ASSERT(GMP_NUMB_BITS >= 6);
+ ASSERT(n >= 17); /* so that r <> 0 and 5k+3 <= 2n */
+
+ /* the algorithm is the same as mpn_toom3_mul_n, with b=a */
+
+ k = (n + 2) / 3; /* ceil(n/3) */
+ twok = 2 * k;
+ k1 = k + 1;
+ kk1 = k + k1;
+ r = n - twok; /* last chunk */
+ twor = 2 * r;
+
+ c1 = c + k;
+ c2 = c1 + k;
+ c3 = c2 + k;
+ c4 = c3 + k;
+
+ trec = t + 4 * k + 3; /* trec = v2 + (2k+2) */
+
+ cy = mpn_add_n (c, a, a + twok, r);
+ if (r < k)
+ __GMPN_ADD_1 (cy, c + r, a + r, k - r, cy);
+ t[k] = (c1[0] = cy) + mpn_add_n (t, c, a + k, k);
+
+#define v2 (t+2*k+1)
+
+ TOOM3_SQR_REC (c2, t, k1, trec);
+
+ sa = (cy != 0) ? 1 : mpn_cmp (c, a + k, k);
+ c[k] = (sa >= 0) ? cy - mpn_sub_n (c, c, a + k, k)
+ : mpn_sub_n (c, a + k, c, k);
+
+ TOOM3_SQR_REC (t, c, k1, trec);
+
+#ifdef HAVE_NATIVE_mpn_addlsh1_n
+ c1[0] = mpn_addlsh1_n (c, a + k, a + twok, r);
+ if (r < k)
+ __GMPN_ADD_1 (c1[0], c + r, a + k + r, k - r, c1[0]);
+ c1[0] = 2 * c1[0] + mpn_addlsh1_n (c, a, c, k);
+#else
+ c[r] = mpn_lshift (c, a + twok, r, 1);
+ if (r < k)
+ MPN_ZERO(c + r + 1, k - r);
+ c1[0] += mpn_add_n (c, c, a + k, k);
+ mpn_lshift (c, c, k1, 1);
+ c1[0] += mpn_add_n (c, c, a, k);
+#endif
+
+ TOOM3_SQR_REC (v2, c, k1, trec);
+
+ TOOM3_SQR_REC (c, a, k, trec);
+
+ saved = c4[0];
+ TOOM3_SQR_REC (c4, a + twok, r, trec);
+ vinf0 = c4[0];
+ c4[0] = saved;
+
+ mpn_toom_interpolate_5pts (c, v2, t, k, 2*r, 1, vinf0, trec);
+
+#undef v2
+}
+
+void
+mpn_mul_n (mp_ptr p, mp_srcptr a, mp_srcptr b, mp_size_t n)
+{
+ ASSERT (n >= 1);
+ ASSERT (! MPN_OVERLAP_P (p, 2 * n, a, n));
+ ASSERT (! MPN_OVERLAP_P (p, 2 * n, b, n));
+
+ if (BELOW_THRESHOLD (n, MUL_KARATSUBA_THRESHOLD))
+ {
+ mpn_mul_basecase (p, a, n, b, n);
+ }
+ else if (BELOW_THRESHOLD (n, MUL_TOOM3_THRESHOLD))
+ {
+ /* Allocate workspace of fixed size on stack: fast! */
+ mp_limb_t ws[MPN_KARA_MUL_N_TSIZE (MUL_TOOM3_THRESHOLD_LIMIT-1)];
+ ASSERT (MUL_TOOM3_THRESHOLD <= MUL_TOOM3_THRESHOLD_LIMIT);
+ mpn_kara_mul_n (p, a, b, n, ws);
+ }
+ else if (BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD))
+ {
+ mp_ptr ws;
+ TMP_SDECL;
+ TMP_SMARK;
+ ws = TMP_SALLOC_LIMBS (MPN_TOOM3_MUL_N_TSIZE (n));
+ mpn_toom3_mul_n (p, a, b, n, ws);
+ TMP_SFREE;
+ }
+#if WANT_FFT || TUNE_PROGRAM_BUILD
+ else if (BELOW_THRESHOLD (n, MUL_FFT_THRESHOLD))
+#else
+ else if (BELOW_THRESHOLD (n, MPN_TOOM44_MAX_N))
+#endif
+ {
+ mp_ptr ws;
+ TMP_SDECL;
+ TMP_SMARK;
+ ws = TMP_SALLOC_LIMBS (mpn_toom44_mul_itch (n, n));
+ mpn_toom44_mul (p, a, n, b, n, ws);
+ TMP_SFREE;
+ }
+ else
+#if WANT_FFT || TUNE_PROGRAM_BUILD
+ {
+ /* The current FFT code allocates its own space. That should probably
+ change. */
+ mpn_mul_fft_full (p, a, n, b, n);
+ }
+#else
+ {
+ /* Toom4 for large operands. */
+ mp_ptr ws;
+ TMP_DECL;
+ TMP_MARK;
+ ws = TMP_BALLOC_LIMBS (mpn_toom44_mul_itch (n, n));
+ mpn_toom44_mul (p, a, n, b, n, ws);
+ TMP_FREE;
+ }
+#endif
+}
+
+void
+mpn_sqr_n (mp_ptr p, mp_srcptr a, mp_size_t n)
+{
+ ASSERT (n >= 1);
+ ASSERT (! MPN_OVERLAP_P (p, 2 * n, a, n));
+
+#if 0
+ /* FIXME: Can this be removed? */
+ if (n == 0)
+ return;
+#endif
+
+ if (BELOW_THRESHOLD (n, SQR_BASECASE_THRESHOLD))
+ { /* mul_basecase is faster than sqr_basecase on small sizes sometimes */
+ mpn_mul_basecase (p, a, n, a, n);
+ }
+ else if (BELOW_THRESHOLD (n, SQR_KARATSUBA_THRESHOLD))
+ {
+ mpn_sqr_basecase (p, a, n);
+ }
+ else if (BELOW_THRESHOLD (n, SQR_TOOM3_THRESHOLD))
+ {
+ /* Allocate workspace of fixed size on stack: fast! */
+ mp_limb_t ws[MPN_KARA_SQR_N_TSIZE (SQR_TOOM3_THRESHOLD_LIMIT-1)];
+ ASSERT (SQR_TOOM3_THRESHOLD <= SQR_TOOM3_THRESHOLD_LIMIT);
+ mpn_kara_sqr_n (p, a, n, ws);
+ }
+ else if (BELOW_THRESHOLD (n, SQR_TOOM4_THRESHOLD))
+ {
+ mp_ptr ws;
+ TMP_SDECL;
+ TMP_SMARK;
+ ws = TMP_SALLOC_LIMBS (MPN_TOOM3_SQR_N_TSIZE (n));
+ mpn_toom3_sqr_n (p, a, n, ws);
+ TMP_SFREE;
+ }
+#if WANT_FFT || TUNE_PROGRAM_BUILD
+ else if (BELOW_THRESHOLD (n, SQR_FFT_THRESHOLD))
+#else
+ else if (BELOW_THRESHOLD (n, MPN_TOOM44_MAX_N))
+#endif
+ {
+ mp_ptr ws;
+ TMP_SDECL;
+ TMP_SMARK;
+ ws = TMP_SALLOC_LIMBS (mpn_toom4_sqr_itch (n));
+ mpn_toom4_sqr (p, a, n, ws);
+ TMP_SFREE;
+ }
+ else
+#if WANT_FFT || TUNE_PROGRAM_BUILD
+ {
+ /* The current FFT code allocates its own space. That should probably
+ change. */
+ mpn_mul_fft_full (p, a, n, a, n);
+ }
+#else
+ {
+ /* Toom4 for large operands. */
+ mp_ptr ws;
+ TMP_DECL;
+ TMP_MARK;
+ ws = TMP_BALLOC_LIMBS (mpn_toom4_sqr_itch (n));
+ mpn_toom4_sqr (p, a, n, ws);
+ TMP_FREE;
+ }
+#endif
+}
projects / msp430-gcc.git / blobdiff