Annotation of OpenXM_contrib/gmp/mpz/tests/t-mul.c, Revision 1.1
1.1 ! maekawa 1: /* Test mpz_add, mpz_cmp, mpz_cmp_ui, mpz_divmod, mpz_mul.
! 2:
! 3: Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
! 4:
! 5: This file is part of the GNU MP Library.
! 6:
! 7: The GNU MP Library is free software; you can redistribute it and/or modify
! 8: it under the terms of the GNU Library General Public License as published by
! 9: the Free Software Foundation; either version 2 of the License, or (at your
! 10: option) any later version.
! 11:
! 12: The GNU MP Library is distributed in the hope that it will be useful, but
! 13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Library General Public License
! 18: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 20: MA 02111-1307, USA. */
! 21:
! 22: #include <stdio.h>
! 23: #include "gmp.h"
! 24: #include "gmp-impl.h"
! 25: #include "longlong.h"
! 26: #include "urandom.h"
! 27:
! 28: void debug_mp ();
! 29: mp_size_t _mpn_mul_classic ();
! 30: void mpz_refmul ();
! 31:
! 32: #ifndef SIZE
! 33: #define SIZE 128
! 34: #endif
! 35:
! 36: main (argc, argv)
! 37: int argc;
! 38: char **argv;
! 39: {
! 40: mpz_t multiplier, multiplicand;
! 41: mpz_t product, ref_product;
! 42: mpz_t quotient, remainder;
! 43: mp_size_t multiplier_size, multiplicand_size;
! 44: int i;
! 45: int reps = 10000;
! 46:
! 47: if (argc == 2)
! 48: reps = atoi (argv[1]);
! 49:
! 50: mpz_init (multiplier);
! 51: mpz_init (multiplicand);
! 52: mpz_init (product);
! 53: mpz_init (ref_product);
! 54: mpz_init (quotient);
! 55: mpz_init (remainder);
! 56:
! 57: for (i = 0; i < reps; i++)
! 58: {
! 59: multiplier_size = urandom () % SIZE - SIZE/2;
! 60: multiplicand_size = urandom () % SIZE - SIZE/2;
! 61:
! 62: mpz_random2 (multiplier, multiplier_size);
! 63: mpz_random2 (multiplicand, multiplicand_size);
! 64:
! 65: mpz_mul (product, multiplier, multiplicand);
! 66: mpz_refmul (ref_product, multiplier, multiplicand);
! 67: if (mpz_cmp_ui (multiplicand, 0) != 0)
! 68: mpz_divmod (quotient, remainder, product, multiplicand);
! 69:
! 70: if (mpz_cmp (product, ref_product))
! 71: dump_abort (multiplier, multiplicand);
! 72:
! 73: if (mpz_cmp_ui (multiplicand, 0) != 0)
! 74: if (mpz_cmp_ui (remainder, 0) || mpz_cmp (quotient, multiplier))
! 75: dump_abort (multiplier, multiplicand);
! 76: }
! 77:
! 78: exit (0);
! 79: }
! 80:
! 81: void
! 82: mpz_refmul (w, u, v)
! 83: mpz_t w;
! 84: const mpz_t u;
! 85: const mpz_t v;
! 86: {
! 87: mp_size_t usize = u->_mp_size;
! 88: mp_size_t vsize = v->_mp_size;
! 89: mp_size_t wsize;
! 90: mp_size_t sign_product;
! 91: mp_ptr up, vp;
! 92: mp_ptr wp;
! 93: mp_ptr free_me = NULL;
! 94: size_t free_me_size;
! 95: TMP_DECL (marker);
! 96:
! 97: TMP_MARK (marker);
! 98: sign_product = usize ^ vsize;
! 99: usize = ABS (usize);
! 100: vsize = ABS (vsize);
! 101:
! 102: if (usize < vsize)
! 103: {
! 104: /* Swap U and V. */
! 105: {const __mpz_struct *t = u; u = v; v = t;}
! 106: {mp_size_t t = usize; usize = vsize; vsize = t;}
! 107: }
! 108:
! 109: up = u->_mp_d;
! 110: vp = v->_mp_d;
! 111: wp = w->_mp_d;
! 112:
! 113: /* Ensure W has space enough to store the result. */
! 114: wsize = usize + vsize;
! 115: if (w->_mp_alloc < wsize)
! 116: {
! 117: if (wp == up || wp == vp)
! 118: {
! 119: free_me = wp;
! 120: free_me_size = w->_mp_alloc;
! 121: }
! 122: else
! 123: (*_mp_free_func) (wp, w->_mp_alloc * BYTES_PER_MP_LIMB);
! 124:
! 125: w->_mp_alloc = wsize;
! 126: wp = (mp_ptr) (*_mp_allocate_func) (wsize * BYTES_PER_MP_LIMB);
! 127: w->_mp_d = wp;
! 128: }
! 129: else
! 130: {
! 131: /* Make U and V not overlap with W. */
! 132: if (wp == up)
! 133: {
! 134: /* W and U are identical. Allocate temporary space for U. */
! 135: up = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
! 136: /* Is V identical too? Keep it identical with U. */
! 137: if (wp == vp)
! 138: vp = up;
! 139: /* Copy to the temporary space. */
! 140: MPN_COPY (up, wp, usize);
! 141: }
! 142: else if (wp == vp)
! 143: {
! 144: /* W and V are identical. Allocate temporary space for V. */
! 145: vp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
! 146: /* Copy to the temporary space. */
! 147: MPN_COPY (vp, wp, vsize);
! 148: }
! 149: }
! 150:
! 151: wsize = _mpn_mul_classic (wp, up, usize, vp, vsize);
! 152: w->_mp_size = sign_product < 0 ? -wsize : wsize;
! 153: if (free_me != NULL)
! 154: (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
! 155:
! 156: TMP_FREE (marker);
! 157: }
! 158:
! 159: mp_size_t
! 160: _mpn_mul_classic (prodp, up, usize, vp, vsize)
! 161: mp_ptr prodp;
! 162: mp_srcptr up;
! 163: mp_size_t usize;
! 164: mp_srcptr vp;
! 165: mp_size_t vsize;
! 166: {
! 167: mp_size_t i, j;
! 168: mp_limb_t prod_low, prod_high;
! 169: mp_limb_t cy_dig;
! 170: mp_limb_t v_limb, c;
! 171:
! 172: if (vsize == 0)
! 173: return 0;
! 174:
! 175: /* Offset UP and PRODP so that the inner loop can be faster. */
! 176: up += usize;
! 177: prodp += usize;
! 178:
! 179: /* Multiply by the first limb in V separately, as the result can
! 180: be stored (not added) to PROD. We also avoid a loop for zeroing. */
! 181: v_limb = vp[0];
! 182: cy_dig = 0;
! 183: j = -usize;
! 184: do
! 185: {
! 186: umul_ppmm (prod_high, prod_low, up[j], v_limb);
! 187: add_ssaaaa (cy_dig, prodp[j], prod_high, prod_low, 0, cy_dig);
! 188: j++;
! 189: }
! 190: while (j < 0);
! 191:
! 192: prodp[j] = cy_dig;
! 193: prodp++;
! 194:
! 195: /* For each iteration in the outer loop, multiply one limb from
! 196: U with one limb from V, and add it to PROD. */
! 197: for (i = 1; i < vsize; i++)
! 198: {
! 199: v_limb = vp[i];
! 200: cy_dig = 0;
! 201: j = -usize;
! 202:
! 203: /* Inner loops. Simulate the carry flag by jumping between
! 204: these loops. The first is used when there was no carry
! 205: in the previois iteration; the second when there was carry. */
! 206:
! 207: do
! 208: {
! 209: umul_ppmm (prod_high, prod_low, up[j], v_limb);
! 210: add_ssaaaa (cy_dig, prod_low, prod_high, prod_low, 0, cy_dig);
! 211: c = prodp[j];
! 212: prod_low += c;
! 213: prodp[j] = prod_low;
! 214: if (prod_low < c)
! 215: goto cy_loop;
! 216: ncy_loop:
! 217: j++;
! 218: }
! 219: while (j < 0);
! 220:
! 221: prodp[j] = cy_dig;
! 222: prodp++;
! 223: continue;
! 224:
! 225: do
! 226: {
! 227: umul_ppmm (prod_high, prod_low, up[j], v_limb);
! 228: add_ssaaaa (cy_dig, prod_low, prod_high, prod_low, 0, cy_dig);
! 229: c = prodp[j];
! 230: prod_low += c + 1;
! 231: prodp[j] = prod_low;
! 232: if (prod_low > c)
! 233: goto ncy_loop;
! 234: cy_loop:
! 235: j++;
! 236: }
! 237: while (j < 0);
! 238:
! 239: cy_dig += 1;
! 240: prodp[j] = cy_dig;
! 241: prodp++;
! 242: }
! 243:
! 244: return usize + vsize - (cy_dig == 0);
! 245: }
! 246:
! 247: dump_abort (multiplier, multiplicand)
! 248: mpz_t multiplier, multiplicand;
! 249: {
! 250: fprintf (stderr, "ERROR\n");
! 251: fprintf (stderr, "multiplier = "); debug_mp (multiplier, -16);
! 252: fprintf (stderr, "multiplicand = "); debug_mp (multiplicand, -16);
! 253: abort();
! 254: }
! 255:
! 256: void
! 257: debug_mp (x, base)
! 258: mpz_t x;
! 259: {
! 260: mpz_out_str (stderr, base, x); fputc ('\n', stderr);
! 261: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to t-mul.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz / tests