source/libs/gmp/gmp-src/mpn/generic/mul_basecase.c

   1 /* mpn_mul_basecase -- Internal routine to multiply two natural numbers
   2    of length m and n.
   3
   4    THIS IS AN INTERNAL FUNCTION WITH A MUTABLE INTERFACE.  IT IS ONLY
   5    SAFE TO REACH THIS FUNCTION THROUGH DOCUMENTED INTERFACES.
   6
   7 Copyright 1991-1994, 1996, 1997, 2000-2002 Free Software Foundation, Inc.
   8
   9 This file is part of the GNU MP Library.
  10
  11 The GNU MP Library is free software; you can redistribute it and/or modify
  12 it under the terms of either:
  13
  14   * the GNU Lesser General Public License as published by the Free
  15     Software Foundation; either version 3 of the License, or (at your
  16     option) any later version.
  17
  18 or
  19
  20   * the GNU General Public License as published by the Free Software
  21     Foundation; either version 2 of the License, or (at your option) any
  22     later version.
  23
  24 or both in parallel, as here.
  25
  26 The GNU MP Library is distributed in the hope that it will be useful, but
  27 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  28 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  29 for more details.
  30
  31 You should have received copies of the GNU General Public License and the
  32 GNU Lesser General Public License along with the GNU MP Library.  If not,
  33 see https://www.gnu.org/licenses/.  */
  34
  35 #include "gmp.h"
  36 #include "gmp-impl.h"
  37
  38
  39 /* Multiply {up,usize} by {vp,vsize} and write the result to
  40    {prodp,usize+vsize}.  Must have usize>=vsize.
  41
  42    Note that prodp gets usize+vsize limbs stored, even if the actual result
  43    only needs usize+vsize-1.
  44
  45    There's no good reason to call here with vsize>=MUL_TOOM22_THRESHOLD.
  46    Currently this is allowed, but it might not be in the future.
  47
  48    This is the most critical code for multiplication.  All multiplies rely
  49    on this, both small and huge.  Small ones arrive here immediately, huge
  50    ones arrive here as this is the base case for Karatsuba's recursive
  51    algorithm.  */
  52
  53 void
  54 mpn_mul_basecase (mp_ptr rp,
  55                   mp_srcptr up, mp_size_t un,
  56                   mp_srcptr vp, mp_size_t vn)
  57 {
  58   ASSERT (un >= vn);
  59   ASSERT (vn >= 1);
  60   ASSERT (! MPN_OVERLAP_P (rp, un+vn, up, un));
  61   ASSERT (! MPN_OVERLAP_P (rp, un+vn, vp, vn));
  62
  63   /* We first multiply by the low order limb (or depending on optional function
  64      availability, limbs).  This result can be stored, not added, to rp.  We
  65      also avoid a loop for zeroing this way.  */
  66
  67 #if HAVE_NATIVE_mpn_mul_2
  68   if (vn >= 2)
  69     {
  70       rp[un + 1] = mpn_mul_2 (rp, up, un, vp);
  71       rp += 2, vp += 2, vn -= 2;
  72     }
  73   else
  74     {
  75       rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
  76       return;
  77     }
  78 #else
  79   rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
  80   rp += 1, vp += 1, vn -= 1;
  81 #endif
  82
  83   /* Now accumulate the product of up[] and the next higher limb (or depending
  84      on optional function availability, limbs) from vp[].  */
  85
  86 #define MAX_LEFT MP_SIZE_T_MAX  /* Used to simplify loops into if statements */
  87
  88
  89 #if HAVE_NATIVE_mpn_addmul_6
  90   while (vn >= 6)
  91     {
  92       rp[un + 6 - 1] = mpn_addmul_6 (rp, up, un, vp);
  93       if (MAX_LEFT == 6)
  94         return;
  95       rp += 6, vp += 6, vn -= 6;
  96       if (MAX_LEFT < 2 * 6)
  97         break;
  98     }
  99 #undef MAX_LEFT
 100 #define MAX_LEFT (6 - 1)
 101 #endif
 102
 103 #if HAVE_NATIVE_mpn_addmul_5
 104   while (vn >= 5)
 105     {
 106       rp[un + 5 - 1] = mpn_addmul_5 (rp, up, un, vp);
 107       if (MAX_LEFT == 5)
 108         return;
 109       rp += 5, vp += 5, vn -= 5;
 110       if (MAX_LEFT < 2 * 5)
 111         break;
 112     }
 113 #undef MAX_LEFT
 114 #define MAX_LEFT (5 - 1)
 115 #endif
 116
 117 #if HAVE_NATIVE_mpn_addmul_4
 118   while (vn >= 4)
 119     {
 120       rp[un + 4 - 1] = mpn_addmul_4 (rp, up, un, vp);
 121       if (MAX_LEFT == 4)
 122         return;
 123       rp += 4, vp += 4, vn -= 4;
 124       if (MAX_LEFT < 2 * 4)
 125         break;
 126     }
 127 #undef MAX_LEFT
 128 #define MAX_LEFT (4 - 1)
 129 #endif
 130
 131 #if HAVE_NATIVE_mpn_addmul_3
 132   while (vn >= 3)
 133     {
 134       rp[un + 3 - 1] = mpn_addmul_3 (rp, up, un, vp);
 135       if (MAX_LEFT == 3)
 136         return;
 137       rp += 3, vp += 3, vn -= 3;
 138       if (MAX_LEFT < 2 * 3)
 139         break;
 140     }
 141 #undef MAX_LEFT
 142 #define MAX_LEFT (3 - 1)
 143 #endif
 144
 145 #if HAVE_NATIVE_mpn_addmul_2
 146   while (vn >= 2)
 147     {
 148       rp[un + 2 - 1] = mpn_addmul_2 (rp, up, un, vp);
 149       if (MAX_LEFT == 2)
 150         return;
 151       rp += 2, vp += 2, vn -= 2;
 152       if (MAX_LEFT < 2 * 2)
 153         break;
 154     }
 155 #undef MAX_LEFT
 156 #define MAX_LEFT (2 - 1)
 157 #endif
 158
 159   while (vn >= 1)
 160     {
 161       rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
 162       if (MAX_LEFT == 1)
 163         return;
 164       rp += 1, vp += 1, vn -= 1;
 165     }
 166 }