Annotation of OpenXM_contrib/gmp/tests/mpz/t-mul.c, Revision 1.1.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