Annotation of OpenXM_contrib/gmp/tests/mpz/t-mul.c, Revision 1.1
1.1 ! ohara 1: /* Test mpz_cmp, mpz_cmp_ui, mpz_tdiv_qr, mpz_mul.
! 2:
! 3: Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002 Free Software
! 4: Foundation, Inc.
! 5:
! 6: This file is part of the GNU MP Library.
! 7:
! 8: The GNU MP Library is free software; you can redistribute it and/or modify
! 9: it under the terms of the GNU Lesser General Public License as published by
! 10: the Free Software Foundation; either version 2.1 of the License, or (at your
! 11: option) any later version.
! 12:
! 13: The GNU MP Library is distributed in the hope that it will be useful, but
! 14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 16: License for more details.
! 17:
! 18: You should have received a copy of the GNU Lesser General Public License
! 19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 21: MA 02111-1307, USA. */
! 22:
! 23: #include <stdio.h>
! 24: #include <stdlib.h>
! 25:
! 26: #include "gmp.h"
! 27: #include "gmp-impl.h"
! 28: #include "longlong.h"
! 29: #include "tests.h"
! 30:
! 31: void debug_mp _PROTO ((mpz_t, int));
! 32: static void base_mul _PROTO ((mp_ptr,mp_srcptr,mp_size_t,mp_srcptr,mp_size_t));
! 33: static void ref_mpz_mul _PROTO ((mpz_t, const mpz_t, const mpz_t));
! 34: void dump_abort _PROTO ((int, char *, mpz_t, mpz_t, mpz_t, mpz_t));
! 35:
! 36: int
! 37: main (int argc, char **argv)
! 38: {
! 39: mpz_t multiplier, multiplicand;
! 40: mpz_t product, ref_product;
! 41: mpz_t quotient, remainder;
! 42: mp_size_t multiplier_size, multiplicand_size;
! 43: int i;
! 44: int reps = 100;
! 45: gmp_randstate_ptr rands;
! 46: mpz_t bs;
! 47: unsigned long bsi, size_range;
! 48:
! 49: tests_start ();
! 50: rands = RANDS;
! 51:
! 52: mpz_init (bs);
! 53:
! 54: if (argc == 2)
! 55: reps = atoi (argv[1]);
! 56:
! 57: mpz_init (multiplier);
! 58: mpz_init (multiplicand);
! 59: mpz_init (product);
! 60: mpz_init (ref_product);
! 61: mpz_init (quotient);
! 62: mpz_init (remainder);
! 63:
! 64: for (i = 0; i < reps; i++)
! 65: {
! 66: mpz_urandomb (bs, rands, 32);
! 67: size_range = mpz_get_ui (bs) % 18 + 2;
! 68:
! 69: mpz_urandomb (bs, rands, size_range);
! 70: multiplier_size = mpz_get_ui (bs);
! 71: mpz_rrandomb (multiplier, rands, multiplier_size);
! 72:
! 73: mpz_urandomb (bs, rands, size_range);
! 74: multiplicand_size = mpz_get_ui (bs);
! 75: mpz_rrandomb (multiplicand, rands, multiplicand_size);
! 76:
! 77: mpz_urandomb (bs, rands, 2);
! 78: bsi = mpz_get_ui (bs);
! 79: if ((bsi & 1) != 0)
! 80: mpz_neg (multiplier, multiplier);
! 81: if ((bsi & 2) != 0)
! 82: mpz_neg (multiplicand, multiplicand);
! 83:
! 84: /* printf ("%ld %ld\n", SIZ (multiplier), SIZ (multiplicand)); */
! 85:
! 86: mpz_mul (product, multiplier, multiplicand);
! 87:
! 88: if (size_range <= 16) /* avoid calling ref_mpz_mul for huge operands */
! 89: {
! 90: ref_mpz_mul (ref_product, multiplier, multiplicand);
! 91: if (mpz_cmp (product, ref_product))
! 92: dump_abort (i, "incorrect plain product",
! 93: multiplier, multiplicand, product, ref_product);
! 94: }
! 95:
! 96: if (mpz_cmp_ui (multiplicand, 0) != 0)
! 97: {
! 98: mpz_tdiv_qr (quotient, remainder, product, multiplicand);
! 99: if (mpz_cmp_ui (remainder, 0) || mpz_cmp (quotient, multiplier))
! 100: {
! 101: debug_mp (quotient, -16);
! 102: debug_mp (remainder, -16);
! 103: dump_abort (i, "incorrect quotient or remainder",
! 104: multiplier, multiplicand, product, ref_product);
! 105: }
! 106: }
! 107:
! 108: /* Test squaring. */
! 109: if (size_range <= 16) /* avoid calling ref_mpz_mul for huge operands */
! 110: {
! 111: mpz_mul (product, multiplier, multiplier);
! 112: ref_mpz_mul (ref_product, multiplier, multiplier);
! 113: }
! 114: else
! 115: {
! 116: mpz_mul (product, multiplier, multiplier);
! 117: mpz_set (multiplicand, multiplier);
! 118: mpz_mul (ref_product, multiplier, multiplicand);
! 119: }
! 120:
! 121: if (mpz_cmp (product, ref_product))
! 122: dump_abort (i, "incorrect square product",
! 123: multiplier, multiplier, product, ref_product);
! 124: }
! 125:
! 126: mpz_clear (bs);
! 127: mpz_clear (multiplier);
! 128: mpz_clear (multiplicand);
! 129: mpz_clear (product);
! 130: mpz_clear (ref_product);
! 131: mpz_clear (quotient);
! 132: mpz_clear (remainder);
! 133:
! 134: tests_end ();
! 135: exit (0);
! 136: }
! 137:
! 138: static void
! 139: ref_mpz_mul (mpz_t w, const mpz_t u, const mpz_t v)
! 140: {
! 141: mp_size_t usize = u->_mp_size;
! 142: mp_size_t vsize = v->_mp_size;
! 143: mp_size_t wsize;
! 144: mp_size_t sign_product;
! 145: mp_ptr up, vp;
! 146: mp_ptr wp;
! 147: mp_ptr free_me = NULL;
! 148: size_t free_me_size;
! 149: TMP_DECL (marker);
! 150:
! 151: TMP_MARK (marker);
! 152: sign_product = usize ^ vsize;
! 153: usize = ABS (usize);
! 154: vsize = ABS (vsize);
! 155:
! 156: if (usize < vsize)
! 157: {
! 158: /* Swap U and V. */
! 159: {const __mpz_struct *t = u; u = v; v = t;}
! 160: {mp_size_t t = usize; usize = vsize; vsize = t;}
! 161: }
! 162:
! 163: if (vsize == 0)
! 164: {
! 165: SIZ (w) = 0;
! 166: return;
! 167: }
! 168:
! 169: up = u->_mp_d;
! 170: vp = v->_mp_d;
! 171: wp = w->_mp_d;
! 172:
! 173: /* Ensure W has space enough to store the result. */
! 174: wsize = usize + vsize;
! 175: if (w->_mp_alloc < wsize)
! 176: {
! 177: if (wp == up || wp == vp)
! 178: {
! 179: free_me = wp;
! 180: free_me_size = w->_mp_alloc;
! 181: }
! 182: else
! 183: (*__gmp_free_func) (wp, w->_mp_alloc * BYTES_PER_MP_LIMB);
! 184:
! 185: w->_mp_alloc = wsize;
! 186: wp = (mp_ptr) (*__gmp_allocate_func) (wsize * BYTES_PER_MP_LIMB);
! 187: w->_mp_d = wp;
! 188: }
! 189: else
! 190: {
! 191: /* Make U and V not overlap with W. */
! 192: if (wp == up)
! 193: {
! 194: /* W and U are identical. Allocate temporary space for U. */
! 195: up = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
! 196: /* Is V identical too? Keep it identical with U. */
! 197: if (wp == vp)
! 198: vp = up;
! 199: /* Copy to the temporary space. */
! 200: MPN_COPY (up, wp, usize);
! 201: }
! 202: else if (wp == vp)
! 203: {
! 204: /* W and V are identical. Allocate temporary space for V. */
! 205: vp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
! 206: /* Copy to the temporary space. */
! 207: MPN_COPY (vp, wp, vsize);
! 208: }
! 209: }
! 210:
! 211: base_mul (wp, up, usize, vp, vsize);
! 212: wsize = usize + vsize;
! 213: wsize -= wp[wsize - 1] == 0;
! 214: w->_mp_size = sign_product < 0 ? -wsize : wsize;
! 215: if (free_me != NULL)
! 216: (*__gmp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
! 217:
! 218: TMP_FREE (marker);
! 219: }
! 220:
! 221: static void
! 222: base_mul (mp_ptr wp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
! 223: {
! 224: mp_size_t i, j;
! 225: mp_limb_t prod_low, prod_high;
! 226: mp_limb_t cy_dig;
! 227: mp_limb_t v_limb;
! 228:
! 229: /* Multiply by the first limb in V separately, as the result can
! 230: be stored (not added) to PROD. We also avoid a loop for zeroing. */
! 231: v_limb = vp[0];
! 232: cy_dig = 0;
! 233: for (j = un; j > 0; j--)
! 234: {
! 235: mp_limb_t u_limb, w_limb;
! 236: u_limb = *up++;
! 237: umul_ppmm (prod_high, prod_low, u_limb, v_limb << GMP_NAIL_BITS);
! 238: add_ssaaaa (cy_dig, w_limb, prod_high, prod_low, 0, cy_dig << GMP_NAIL_BITS);
! 239: *wp++ = w_limb >> GMP_NAIL_BITS;
! 240: }
! 241:
! 242: *wp++ = cy_dig;
! 243: wp -= un;
! 244: up -= un;
! 245:
! 246: /* For each iteration in the outer loop, multiply one limb from
! 247: U with one limb from V, and add it to PROD. */
! 248: for (i = 1; i < vn; i++)
! 249: {
! 250: v_limb = vp[i];
! 251: cy_dig = 0;
! 252:
! 253: for (j = un; j > 0; j--)
! 254: {
! 255: mp_limb_t u_limb, w_limb;
! 256: u_limb = *up++;
! 257: umul_ppmm (prod_high, prod_low, u_limb, v_limb << GMP_NAIL_BITS);
! 258: w_limb = *wp;
! 259: add_ssaaaa (prod_high, prod_low, prod_high, prod_low, 0, w_limb << GMP_NAIL_BITS);
! 260: prod_low >>= GMP_NAIL_BITS;
! 261: prod_low += cy_dig;
! 262: #if GMP_NAIL_BITS == 0
! 263: cy_dig = prod_high + (prod_low < cy_dig);
! 264: #else
! 265: cy_dig = prod_high;
! 266: cy_dig += prod_low >> GMP_NUMB_BITS;
! 267: #endif
! 268: *wp++ = prod_low & GMP_NUMB_MASK;
! 269: }
! 270:
! 271: *wp++ = cy_dig;
! 272: wp -= un;
! 273: up -= un;
! 274: }
! 275: }
! 276:
! 277: void
! 278: dump_abort (int i, char *s,
! 279: mpz_t multiplier, mpz_t multiplicand, mpz_t product, mpz_t ref_product)
! 280: {
! 281: mpz_t diff;
! 282: fprintf (stderr, "ERROR: %s in test %d\n", s, i);
! 283: fprintf (stderr, "multiplier = "); debug_mp (multiplier, -16);
! 284: fprintf (stderr, "multiplicand = "); debug_mp (multiplicand, -16);
! 285: fprintf (stderr, " product = "); debug_mp (product, -16);
! 286: fprintf (stderr, "ref_product = "); debug_mp (ref_product, -16);
! 287: mpz_init (diff);
! 288: mpz_sub (diff, ref_product, product);
! 289: fprintf (stderr, "diff: "); debug_mp (diff, -16);
! 290: abort();
! 291: }
! 292:
! 293: void
! 294: debug_mp (mpz_t x, int base)
! 295: {
! 296: mpz_out_str (stderr, base, x); fputc ('\n', stderr);
! 297: }
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 / tests / mpz