contrib/gmp/mpn/generic/toom53_mul.c

   1 /* mpn_toom53_mul -- Multiply {ap,an} and {bp,bn} where an is nominally 5/3
   2    times as large as bn.  Or more accurately, (4/3)bn < an < (5/2)bn.
   3
   4    Contributed to the GNU project by Torbjorn Granlund and 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 2006, 2007, 2008 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 the GNU Lesser General Public License as published by
  19 the Free Software Foundation; either version 3 of the License, or (at your
  20 option) any later version.
  21
  22 The GNU MP Library is distributed in the hope that it will be useful, but
  23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  24 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  25 License for more details.
  26
  27 You should have received a copy of the GNU Lesser General Public License
  28 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  29
  30
  31 #include "gmp.h"
  32 #include "gmp-impl.h"
  33
  34 /* Evaluate in: 0, +1, -1, +2, -2, 1/2, +inf
  35
  36   <-s-><--n--><--n--><--n--><--n-->
  37    ___ ______ ______ ______ ______
  38   |a4_|___a3_|___a2_|___a1_|___a0_|
  39                |__b2|___b1_|___b0_|
  40                <-t--><--n--><--n-->
  41
  42   v0  =    a0                  *  b0          #    A(0)*B(0)
  43   v1  = (  a0+ a1+ a2+ a3+  a4)*( b0+ b1+ b2) #    A(1)*B(1)      ah  <= 4   bh <= 2
  44   vm1 = (  a0- a1+ a2- a3+  a4)*( b0- b1+ b2) #   A(-1)*B(-1)    |ah| <= 2   bh <= 1
  45   v2  = (  a0+2a1+4a2+8a3+16a4)*( b0+2b1+4b2) #    A(2)*B(2)      ah  <= 30  bh <= 6
  46   vm2 = (  a0-2a1+4a2-8a3+16a4)*( b0-2b1+4b2) #    A(2)*B(2)     -9<=ah<=20 -1<=bh<=4
  47   vh  = (16a0+8a1+4a2+2a3+  a4)*(4b0+2b1+ b2) #  A(1/2)*B(1/2)    ah  <= 30  bh <= 6
  48   vinf=                     a4 *          b2  #  A(inf)*B(inf)
  49 */
  50
  51 void
  52 mpn_toom53_mul (mp_ptr pp,
  53                 mp_srcptr ap, mp_size_t an,
  54                 mp_srcptr bp, mp_size_t bn,
  55                 mp_ptr scratch)
  56 {
  57   mp_size_t n, s, t;
  58   mp_limb_t cy;
  59   mp_ptr gp;
  60   mp_ptr as1, asm1, as2, asm2, ash;
  61   mp_ptr bs1, bsm1, bs2, bsm2, bsh;
  62   enum toom7_flags flags;
  63   TMP_DECL;
  64
  65 #define a0  ap
  66 #define a1  (ap + n)
  67 #define a2  (ap + 2*n)
  68 #define a3  (ap + 3*n)
  69 #define a4  (ap + 4*n)
  70 #define b0  bp
  71 #define b1  (bp + n)
  72 #define b2  (bp + 2*n)
  73
  74   n = 1 + (3 * an >= 5 * bn ? (an - 1) / (size_t) 5 : (bn - 1) / (size_t) 3);
  75
  76   s = an - 4 * n;
  77   t = bn - 2 * n;
  78
  79   ASSERT (0 < s && s <= n);
  80   ASSERT (0 < t && t <= n);
  81
  82   TMP_MARK;
  83
  84   as1  = TMP_SALLOC_LIMBS (n + 1);
  85   asm1 = TMP_SALLOC_LIMBS (n + 1);
  86   as2  = TMP_SALLOC_LIMBS (n + 1);
  87   asm2 = TMP_SALLOC_LIMBS (n + 1);
  88   ash  = TMP_SALLOC_LIMBS (n + 1);
  89
  90   bs1  = TMP_SALLOC_LIMBS (n + 1);
  91   bsm1 = TMP_SALLOC_LIMBS (n + 1);
  92   bs2  = TMP_SALLOC_LIMBS (n + 1);
  93   bsm2 = TMP_SALLOC_LIMBS (n + 1);
  94   bsh  = TMP_SALLOC_LIMBS (n + 1);
  95
  96   gp = pp;
  97
  98   /* Compute as1 and asm1.  */
  99   flags = toom7_w3_neg & mpn_toom_eval_pm1 (as1, asm1, 4, ap, n, s, gp);
 100
 101   /* Compute as2 and asm2. */
 102   flags |= toom7_w1_neg & mpn_toom_eval_pm2 (as2, asm2, 4, ap, n, s, gp);
 103
 104   /* Compute ash = 16 a0 + 8 a1 + 4 a2 + 2 a3 + a4
 105      = 2*(2*(2*(2*a0 + a1) + a2) + a3) + a4  */
 106 #if HAVE_NATIVE_mpn_addlsh1_n
 107   cy = mpn_addlsh1_n (ash, a1, a0, n);
 108   cy = 2*cy + mpn_addlsh1_n (ash, a2, ash, n);
 109   cy = 2*cy + mpn_addlsh1_n (ash, a3, ash, n);
 110   if (s < n)
 111     {
 112       mp_limb_t cy2;
 113       cy2 = mpn_addlsh1_n (ash, a4, ash, s);
 114       ash[n] = 2*cy + mpn_lshift (ash + s, ash + s, n - s, 1);
 115       MPN_INCR_U (ash + s, n+1-s, cy2);
 116     }
 117   else
 118     ash[n] = 2*cy + mpn_addlsh1_n (ash, a4, ash, n);
 119 #else
 120   cy = mpn_lshift (ash, a0, n, 1);
 121   cy += mpn_add_n (ash, ash, a1, n);
 122   cy = 2*cy + mpn_lshift (ash, ash, n, 1);
 123   cy += mpn_add_n (ash, ash, a2, n);
 124   cy = 2*cy + mpn_lshift (ash, ash, n, 1);
 125   cy += mpn_add_n (ash, ash, a3, n);
 126   cy = 2*cy + mpn_lshift (ash, ash, n, 1);
 127   ash[n] = cy + mpn_add (ash, ash, n, a4, s);
 128 #endif
 129
 130   /* Compute bs1 and bsm1.  */
 131   bs1[n] = mpn_add (bs1, b0, n, b2, t);         /* b0 + b2 */
 132 #if HAVE_NATIVE_mpn_add_n_sub_n
 133   if (bs1[n] == 0 && mpn_cmp (bs1, b1, n) < 0)
 134     {
 135       bs1[n] = mpn_add_n_sub_n (bs1, bsm1, b1, bs1, n) >> 1;
 136       bsm1[n] = 0;
 137       flags ^= toom7_w3_neg;
 138     }
 139   else
 140     {
 141       cy = mpn_add_n_sub_n (bs1, bsm1, bs1, b1, n);
 142       bsm1[n] = bs1[n] - (cy & 1);
 143       bs1[n] += (cy >> 1);
 144     }
 145 #else
 146   if (bs1[n] == 0 && mpn_cmp (bs1, b1, n) < 0)
 147     {
 148       mpn_sub_n (bsm1, b1, bs1, n);
 149       bsm1[n] = 0;
 150       flags ^= toom7_w3_neg;
 151     }
 152   else
 153     {
 154       bsm1[n] = bs1[n] - mpn_sub_n (bsm1, bs1, b1, n);
 155     }
 156   bs1[n] += mpn_add_n (bs1, bs1, b1, n);  /* b0+b1+b2 */
 157 #endif
 158
 159   /* Compute bs2 and bsm2. */
 160 #if HAVE_NATIVE_mpn_addlsh_n || HAVE_NATIVE_mpn_addlsh2_n
 161 #if HAVE_NATIVE_mpn_addlsh2_n
 162   cy = mpn_addlsh2_n (bs2, b0, b2, t);
 163 #else /* HAVE_NATIVE_mpn_addlsh_n */
 164   cy = mpn_addlsh_n (bs2, b0, b2, t, 2);
 165 #endif
 166   if (t < n)
 167     cy = mpn_add_1 (bs2 + t, b0 + t, n - t, cy);
 168   bs2[n] = cy;
 169 #else
 170   cy = mpn_lshift (gp, b2, t, 2);
 171   bs2[n] = mpn_add (bs2, b0, n, gp, t);
 172   MPN_INCR_U (bs2 + t, n+1-t, cy);
 173 #endif
 174
 175   gp[n] = mpn_lshift (gp, b1, n, 1);
 176
 177 #if HAVE_NATIVE_mpn_add_n_sub_n
 178   if (mpn_cmp (bs2, gp, n+1) < 0)
 179     {
 180       ASSERT_NOCARRY (mpn_add_n_sub_n (bs2, bsm2, gp, bs2, n+1));
 181       flags ^= toom7_w1_neg;
 182     }
 183   else
 184     {
 185       ASSERT_NOCARRY (mpn_add_n_sub_n (bs2, bsm2, bs2, gp, n+1));
 186     }
 187 #else
 188   if (mpn_cmp (bs2, gp, n+1) < 0)
 189     {
 190       ASSERT_NOCARRY (mpn_sub_n (bsm2, gp, bs2, n+1));
 191       flags ^= toom7_w1_neg;
 192     }
 193   else
 194     {
 195       ASSERT_NOCARRY (mpn_sub_n (bsm2, bs2, gp, n+1));
 196     }
 197   mpn_add_n (bs2, bs2, gp, n+1);
 198 #endif
 199
 200   /* Compute bsh = 4 b0 + 2 b1 + b0 = 2*(2*b0 + b1)+b0.  */
 201 #if HAVE_NATIVE_mpn_addlsh1_n
 202   cy = mpn_addlsh1_n (bsh, b1, b0, n);
 203   if (t < n)
 204     {
 205       mp_limb_t cy2;
 206       cy2 = mpn_addlsh1_n (bsh, b2, bsh, t);
 207       bsh[n] = 2*cy + mpn_lshift (bsh + t, bsh + t, n - t, 1);
 208       MPN_INCR_U (bsh + t, n+1-t, cy2);
 209     }
 210   else
 211     bsh[n] = 2*cy + mpn_addlsh1_n (bsh, b2, bsh, n);
 212 #else
 213   cy = mpn_lshift (bsh, b0, n, 1);
 214   cy += mpn_add_n (bsh, bsh, b1, n);
 215   cy = 2*cy + mpn_lshift (bsh, bsh, n, 1);
 216   bsh[n] = cy + mpn_add (bsh, bsh, n, b2, t);
 217 #endif
 218
 219   ASSERT (as1[n] <= 4);
 220   ASSERT (bs1[n] <= 2);
 221   ASSERT (asm1[n] <= 2);
 222   ASSERT (bsm1[n] <= 1);
 223   ASSERT (as2[n] <= 30);
 224   ASSERT (bs2[n] <= 6);
 225   ASSERT (asm2[n] <= 20);
 226   ASSERT (bsm2[n] <= 4);
 227   ASSERT (ash[n] <= 30);
 228   ASSERT (bsh[n] <= 6);
 229
 230 #define v0    pp                                /* 2n */
 231 #define v1    (pp + 2 * n)                      /* 2n+1 */
 232 #define vinf  (pp + 6 * n)                      /* s+t */
 233 #define v2    scratch                           /* 2n+1 */
 234 #define vm2   (scratch + 2 * n + 1)             /* 2n+1 */
 235 #define vh    (scratch + 4 * n + 2)             /* 2n+1 */
 236 #define vm1   (scratch + 6 * n + 3)             /* 2n+1 */
 237 #define scratch_out (scratch + 8 * n + 4)               /* 2n+1 */
 238   /* Total scratch need: 10*n+5 */
 239
 240   /* Must be in allocation order, as they overwrite one limb beyond
 241    * 2n+1. */
 242   mpn_mul_n (v2, as2, bs2, n + 1);              /* v2, 2n+1 limbs */
 243   mpn_mul_n (vm2, asm2, bsm2, n + 1);           /* vm2, 2n+1 limbs */
 244   mpn_mul_n (vh, ash, bsh, n + 1);              /* vh, 2n+1 limbs */
 245
 246   /* vm1, 2n+1 limbs */
 247 #ifdef SMALLER_RECURSION
 248   mpn_mul_n (vm1, asm1, bsm1, n);
 249   if (asm1[n] == 1)
 250     {
 251       cy = bsm1[n] + mpn_add_n (vm1 + n, vm1 + n, bsm1, n);
 252     }
 253   else if (asm1[n] == 2)
 254     {
 255 #if HAVE_NATIVE_mpn_addlsh1_n
 256       cy = 2 * bsm1[n] + mpn_addlsh1_n (vm1 + n, vm1 + n, bsm1, n);
 257 #else
 258       cy = 2 * bsm1[n] + mpn_addmul_1 (vm1 + n, bsm1, n, CNST_LIMB(2));
 259 #endif
 260     }
 261   else
 262     cy = 0;
 263   if (bsm1[n] != 0)
 264     cy += mpn_add_n (vm1 + n, vm1 + n, asm1, n);
 265   vm1[2 * n] = cy;
 266 #else /* SMALLER_RECURSION */
 267   vm1[2 * n] = 0;
 268   mpn_mul_n (vm1, asm1, bsm1, n + ((asm1[n] | bsm1[n]) != 0));
 269 #endif /* SMALLER_RECURSION */
 270
 271   /* v1, 2n+1 limbs */
 272 #ifdef SMALLER_RECURSION
 273   mpn_mul_n (v1, as1, bs1, n);
 274   if (as1[n] == 1)
 275     {
 276       cy = bs1[n] + mpn_add_n (v1 + n, v1 + n, bs1, n);
 277     }
 278   else if (as1[n] == 2)
 279     {
 280 #if HAVE_NATIVE_mpn_addlsh1_n
 281       cy = 2 * bs1[n] + mpn_addlsh1_n (v1 + n, v1 + n, bs1, n);
 282 #else
 283       cy = 2 * bs1[n] + mpn_addmul_1 (v1 + n, bs1, n, CNST_LIMB(2));
 284 #endif
 285     }
 286   else if (as1[n] != 0)
 287     {
 288       cy = as1[n] * bs1[n] + mpn_addmul_1 (v1 + n, bs1, n, as1[n]);
 289     }
 290   else
 291     cy = 0;
 292   if (bs1[n] == 1)
 293     {
 294       cy += mpn_add_n (v1 + n, v1 + n, as1, n);
 295     }
 296   else if (bs1[n] == 2)
 297     {
 298 #if HAVE_NATIVE_mpn_addlsh1_n
 299       cy += mpn_addlsh1_n (v1 + n, v1 + n, as1, n);
 300 #else
 301       cy += mpn_addmul_1 (v1 + n, as1, n, CNST_LIMB(2));
 302 #endif
 303     }
 304   v1[2 * n] = cy;
 305 #else /* SMALLER_RECURSION */
 306   v1[2 * n] = 0;
 307   mpn_mul_n (v1, as1, bs1, n + ((as1[n] | bs1[n]) != 0));
 308 #endif /* SMALLER_RECURSION */
 309
 310   mpn_mul_n (v0, a0, b0, n);                    /* v0, 2n limbs */
 311
 312   /* vinf, s+t limbs */
 313   if (s > t)  mpn_mul (vinf, a4, s, b2, t);
 314   else        mpn_mul (vinf, b2, t, a4, s);
 315
 316   mpn_toom_interpolate_7pts (pp, n, flags, vm2, vm1, v2, vh, s + t,
 317                              scratch_out);
 318
 319   TMP_FREE;
 320 }