/* _mpn_mul -- Multiply two natural numbers. Copyright (C) 1991 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 General Public License as published by the Free Software Foundation; either version 2, 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 General Public License for more details. You should have received a copy of the GNU General Public License along with the GNU MP Library; see the file COPYING. If not, write to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */ #include #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" /* disabled for now, _mpn_mul_classic has vanished somehow... */ #if defined(DEBUG) && 0 #define MPN_MUL_VERIFY(res_ptr,res_size,op1_ptr,op1_size,op2_ptr,op2_size) \ _mpn_mul_verify (res_ptr, res_size, op1_ptr, op1_size, op2_ptr, op2_size) #include static void _mpn_mul_verify (res_ptr, res_size, op1_ptr, op1_size, op2_ptr, op2_size) mp_ptr res_ptr, op1_ptr, op2_ptr; mp_size res_size, op1_size, op2_size; { mp_ptr tmp_ptr; mp_size tmp_size; tmp_ptr = alloca ((op1_size + op2_size) * BYTES_PER_MP_LIMB); /* debug only */ if (op1_size >= op2_size) tmp_size = _mpn_mul_classic (tmp_ptr, op1_ptr, op1_size, op2_ptr, op2_size); else tmp_size = _mpn_mul_classic (tmp_ptr, op2_ptr, op2_size, op1_ptr, op1_size); if (tmp_size != res_size || _mpn_cmp (res_ptr, res_size, tmp_ptr, tmp_size) != 0) { fprintf (stderr, "GNU MP internal error: Wrong result in _mpn_mul.\n"); fprintf (stderr, "op1{%d} = ", op1_size); _mpn_print (op1_ptr, op1_size); fprintf (stderr, "op2{%d} = ", op2_size); _mpn_print (op2_ptr, op2_size); abort(); } } _mpn_print (mp_ptr p, mp_size size) { mp_size i; for (i = size - 1; i >= 0; i--) printf ("%08X", p[i]); puts (""); } #else #define MPN_MUL_VERIFY(a,b,c,d,e,f) #endif /* Multiply the natural numbers u (pointed to by UP, with USIZE limbs) and v (pointed to by VP, with VSIZE limbs), and store the result at PRODP. USIZE + VSIZE limbs are always stored, but if the input operands are normalized, the return value will reflect the true result size (which is either USIZE + VSIZE, or USIZE + VSIZE -1). NOTE: The space pointed to by PRODP is overwritten before finished with U and V, so overlap is an error. Argument constraints: 1. USIZE >= VSIZE. 2. PRODP != UP and PRODP != VP, i.e. the destination must be distinct from the multiplier and the multiplicand. */ #define KARATSUBA_THRESHOLD 100000000 /* Ridiculous large, but we will never do alloca ! */ /* If KARATSUBA_THRESHOLD is not already defined, define it to a value which is good on most machines. */ #ifndef KARATSUBA_THRESHOLD #define KARATSUBA_THRESHOLD 8 #endif /* The code can't handle KARATSUBA_THRESHOLD smaller than 4. */ #if KARATSUBA_THRESHOLD < 4 #undef KARATSUBA_THRESHOLD #define KARATSUBA_THRESHOLD 4 #endif mp_size _mpn_mul (mp_ptr prodp, mp_srcptr up, mp_size usize, mp_srcptr vp, mp_size vsize) { mp_size n; mp_size prod_size; mp_limb cy; if (vsize < KARATSUBA_THRESHOLD) { /* Handle simple cases with traditional multiplication. This is the most critical code of the entire function. All multiplies rely on this, both small and huge. Small ones arrive here immediately. Huge ones arrive here as this is the base case for the recursive algorithm below. */ Int i, j; mp_limb prod_low, prod_high; mp_limb cy_limb; mp_limb v_limb, c; if (vsize == 0) return 0; /* Offset UP and PRODP so that the inner loop can be faster. */ up += usize; prodp += usize; /* Multiply by the first limb in V separately, as the result can be stored (not added) to PROD. We also avoid a loop for zeroing. */ v_limb = vp[0]; cy_limb = 0; j = -usize; do { umul_ppmm (prod_high, prod_low, up[j], v_limb); add_ssaaaa (cy_limb, prodp[j], prod_high, prod_low, 0, cy_limb); j++; } while (j < 0); prodp[j] = cy_limb; prodp++; /* For each iteration in the outer loop, multiply one limb from U with one limb from V, and add it to PROD. */ for (i = 1; i < vsize; i++) { v_limb = vp[i]; cy_limb = 0; j = -usize; /* Inner loops. Simulate the carry flag by jumping between these loops. The first is used when there was no carry in the previois iteration; the second when there was carry. */ do { umul_ppmm (prod_high, prod_low, up[j], v_limb); add_ssaaaa (cy_limb, prod_low, prod_high, prod_low, 0, cy_limb); c = prodp[j]; prod_low += c; prodp[j] = prod_low; if (prod_low < c) goto cy_loop; ncy_loop: j++; } while (j < 0); prodp[j] = cy_limb; prodp++; continue; do { umul_ppmm (prod_high, prod_low, up[j], v_limb); add_ssaaaa (cy_limb, prod_low, prod_high, prod_low, 0, cy_limb); c = prodp[j]; prod_low += c + 1; prodp[j] = prod_low; if (prod_low > c) goto ncy_loop; cy_loop: j++; } while (j < 0); cy_limb += 1; prodp[j] = cy_limb; prodp++; } return usize + vsize - (cy_limb == 0); } n = (usize + 1) / 2; /* Is USIZE larger than 1.5 times VSIZE? Avoid Karatsuba's algorithm. */ if (2 * usize > 3 * vsize) { /* If U has at least twice as many limbs as V. Split U in two pieces, U1 and U0, such that U = U0 + U1*(2**BITS_PER_MP_LIMB)**N, and recursively multiply the two pieces separately with V. */ mp_size u0_size; mp_ptr tmp; mp_size tmp_size; /* V1 (the high part of V) is zero. */ /* Calculate the length of U0. It is normally equal to n, but of course not for sure. */ for (u0_size = n; u0_size > 0 && up[u0_size - 1] == 0; u0_size--) ; /* Perform (U0 * V). */ if (u0_size >= vsize) prod_size = _mpn_mul (prodp, up, u0_size, vp, vsize); else prod_size = _mpn_mul (prodp, vp, vsize, up, u0_size); MPN_MUL_VERIFY (prodp, prod_size, up, u0_size, vp, vsize); /* We have to zero-extend the lower partial product to n limbs, since the mpn_add some lines below expect the first n limbs to be well defined. (This is normally a no-op. It may do something when U1 has many leading 0 limbs.) */ while (prod_size < n) prodp[prod_size++] = 0; tmp = (mp_ptr) alloca ((size_t)((usize + vsize - n) * BYTES_PER_MP_LIMB)); /* Perform (U1 * V). Make sure the first source argument to _mpn_mul is not less than the second source argument. */ if (vsize <= usize - n) tmp_size = _mpn_mul (tmp, up + n, usize - n, vp, vsize); else tmp_size = _mpn_mul (tmp, vp, vsize, up + n, usize - n); MPN_MUL_VERIFY (tmp, tmp_size, up + n, usize - n, vp, vsize); /* In this addition hides a potentially large copying of TMP. */ if (prod_size - n >= tmp_size) cy = _mpn_add (prodp + n, prodp + n, prod_size - n, tmp, tmp_size); else cy = _mpn_add (prodp + n, tmp, tmp_size, prodp + n, prod_size - n); if (cy) abort (); /* prodp[prod_size] = cy; */ return tmp_size + n; } else { /* Karatsuba's divide-and-conquer algorithm. Split U in two pieces, U1 and U0, such that U = U0 + U1*(B**n), and V in V1 and V0, such that V = V0 + V1*(B**n). UV is then computed recursively using the identity 2n n n n UV = (B + B )U V + B (U -U )(V -V ) + (B + 1)U V 1 1 1 0 0 1 0 0 Where B = 2**BITS_PER_MP_LIMB. */ /* It's possible to decrease the temporary allocation by using the prodp area for temporary storage of the middle term, and doing that recursive multiplication first. (Do this later.) */ mp_size u0_size; mp_size v0_size; mp_size u0v0_size; mp_size u1v1_size; mp_ptr temp; mp_size temp_size; mp_size utem_size; mp_size vtem_size; mp_ptr ptem; mp_size ptem_size; int negflg; mp_ptr pp; pp = (mp_ptr) alloca ((size_t)(4 * n * BYTES_PER_MP_LIMB)); /* Calculate the lengths of U0 and V0. They are normally equal to n, but of course not for sure. */ for (u0_size = n; u0_size > 0 && up[u0_size - 1] == 0; u0_size--) ; for (v0_size = n; v0_size > 0 && vp[v0_size - 1] == 0; v0_size--) ; /*** 1. PROD]2n..0] := U0 x V0 (Recursive call to _mpn_mul may NOT overwrite input operands.) ________________ ________________ |________________||____U0 x V0_____| */ if (u0_size >= v0_size) u0v0_size = _mpn_mul (pp, up, u0_size, vp, v0_size); else u0v0_size = _mpn_mul (pp, vp, v0_size, up, u0_size); MPN_MUL_VERIFY (pp, u0v0_size, up, u0_size, vp, v0_size); /* Zero-extend to 2n limbs. */ while (u0v0_size < 2 * n) pp[u0v0_size++] = 0; /*** 2. PROD]4n..2n] := U1 x V1 (Recursive call to _mpn_mul may NOT overwrite input operands.) ________________ ________________ |_____U1 x V1____||____U0 x V0_____| */ u1v1_size = _mpn_mul (pp + 2*n, up + n, usize - n, vp + n, vsize - n); MPN_MUL_VERIFY (pp + 2*n, u1v1_size, up + n, usize - n, vp + n, vsize - n); prod_size = 2 * n + u1v1_size; /*** 3. PTEM]2n..0] := (U1-U0) x (V0-V1) (Recursive call to _mpn_mul may overwrite input operands.) ________________ |_(U1-U0)(V0-V1)_| */ temp = (mp_ptr) alloca ((size_t)((2 * n + 1) * BYTES_PER_MP_LIMB)); if (_mpn_cmp (up + n, usize - n, up, u0_size) >= 0) { utem_size = usize - n + _mpn_sub (temp, up + n, usize - n, up, u0_size); negflg = 0; } else { utem_size = u0_size + _mpn_sub (temp, up, u0_size, up + n, usize - n); negflg = 1; } if (_mpn_cmp (vp + n, vsize - n, vp, v0_size) >= 0) { vtem_size = vsize - n + _mpn_sub (temp + n, vp + n, vsize - n, vp, v0_size); negflg ^= 1; } else { vtem_size = v0_size + _mpn_sub (temp + n, vp, v0_size, vp + n, vsize - n); /* No change of NEGFLG. */ } ptem = (mp_ptr) alloca ((size_t)(2 * n * BYTES_PER_MP_LIMB)); if (utem_size >= vtem_size) ptem_size = _mpn_mul (ptem, temp, utem_size, temp + n, vtem_size); else ptem_size = _mpn_mul (ptem, temp + n, vtem_size, temp, utem_size); MPN_MUL_VERIFY (ptem, ptem_size, temp, utem_size, temp + n, vtem_size); /*** 4. TEMP]2n..0] := PROD]2n..0] + PROD]4n..2n] ________________ |_____U1 x V1____| ________________ |_____U0_x_V0____| */ cy = _mpn_add (temp, pp, 2*n, pp + 2*n, u1v1_size); if (cy != 0) { temp[2*n] = cy; temp_size = 2*n + 1; } else { /* Normalize temp. pp[2*n-1] might have been zero in the _mpn_add call above, and thus temp might be unnormalized. */ for (temp_size = 2*n; temp_size > 0 && temp[temp_size - 1] == 0; temp_size--) ; } if (prod_size - n >= temp_size) cy = _mpn_add (pp + n, pp + n, prod_size - n, temp, temp_size); else { /* This is a weird special case that should not happen (often)! */ cy = _mpn_add (pp + n, temp, temp_size, pp + n, prod_size - n); prod_size = temp_size + n; } if (cy != 0) { pp[prod_size] = cy; prod_size++; } #ifdef DEBUG if (prod_size > 4 * n) abort(); #endif if (negflg) prod_size = prod_size + _mpn_sub (pp + n, pp + n, prod_size - n, ptem, ptem_size); else { if (prod_size - n < ptem_size) abort(); cy = _mpn_add (pp + n, pp + n, prod_size - n, ptem, ptem_size); if (cy != 0) { pp[prod_size] = cy; prod_size++; #ifdef DEBUG if (prod_size > 4 * n) abort(); #endif } } MPN_COPY (prodp, pp, prod_size); return prod_size; } }