Annotation of OpenXM_contrib/gmp/mpz/tests/t-mul.c, Revision 1.1.1.2
1.1 maekawa 1: /* Test mpz_add, mpz_cmp, mpz_cmp_ui, mpz_divmod, mpz_mul.
2:
1.1.1.2 ! maekawa 3: Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation, Inc.
1.1 maekawa 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
1.1.1.2 ! maekawa 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 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
1.1.1.2 ! maekawa 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 15: License for more details.
16:
1.1.1.2 ! maekawa 17: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 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
1.1.1.2 ! maekawa 33: #define SIZE 2 * TOOM3_MUL_THRESHOLD
1.1 maekawa 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;
1.1.1.2 ! maekawa 45: int reps = 2000;
1.1 maekawa 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: {
1.1.1.2 ! maekawa 59: multiplier_size = urandom () % 2 * SIZE - SIZE;
! 60: multiplicand_size = urandom () % 2 * SIZE - SIZE;
1.1 maekawa 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))
1.1.1.2 ! maekawa 71: dump_abort ("incorrect plain product",
! 72: multiplier, multiplicand, product, ref_product);
1.1 maekawa 73:
74: if (mpz_cmp_ui (multiplicand, 0) != 0)
75: if (mpz_cmp_ui (remainder, 0) || mpz_cmp (quotient, multiplier))
1.1.1.2 ! maekawa 76: {
! 77: debug_mp (quotient, -16);
! 78: debug_mp (remainder, -16);
! 79: dump_abort ("incorrect quotient or remainder",
! 80: multiplier, multiplicand, product, ref_product);
! 81: }
! 82:
! 83: /* Test squaring. */
! 84: mpz_mul (product, multiplier, multiplier);
! 85: mpz_refmul (ref_product, multiplier, multiplier);
! 86:
! 87: if (mpz_cmp (product, ref_product))
! 88: dump_abort ("incorrect square product",
! 89: multiplier, multiplier, product, ref_product);
1.1 maekawa 90: }
91:
92: exit (0);
93: }
94:
95: void
96: mpz_refmul (w, u, v)
97: mpz_t w;
98: const mpz_t u;
99: const mpz_t v;
100: {
101: mp_size_t usize = u->_mp_size;
102: mp_size_t vsize = v->_mp_size;
103: mp_size_t wsize;
104: mp_size_t sign_product;
105: mp_ptr up, vp;
106: mp_ptr wp;
107: mp_ptr free_me = NULL;
108: size_t free_me_size;
109: TMP_DECL (marker);
110:
111: TMP_MARK (marker);
112: sign_product = usize ^ vsize;
113: usize = ABS (usize);
114: vsize = ABS (vsize);
115:
116: if (usize < vsize)
117: {
118: /* Swap U and V. */
119: {const __mpz_struct *t = u; u = v; v = t;}
120: {mp_size_t t = usize; usize = vsize; vsize = t;}
121: }
122:
123: up = u->_mp_d;
124: vp = v->_mp_d;
125: wp = w->_mp_d;
126:
127: /* Ensure W has space enough to store the result. */
128: wsize = usize + vsize;
129: if (w->_mp_alloc < wsize)
130: {
131: if (wp == up || wp == vp)
132: {
133: free_me = wp;
134: free_me_size = w->_mp_alloc;
135: }
136: else
137: (*_mp_free_func) (wp, w->_mp_alloc * BYTES_PER_MP_LIMB);
138:
139: w->_mp_alloc = wsize;
140: wp = (mp_ptr) (*_mp_allocate_func) (wsize * BYTES_PER_MP_LIMB);
141: w->_mp_d = wp;
142: }
143: else
144: {
145: /* Make U and V not overlap with W. */
146: if (wp == up)
147: {
148: /* W and U are identical. Allocate temporary space for U. */
149: up = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
150: /* Is V identical too? Keep it identical with U. */
151: if (wp == vp)
152: vp = up;
153: /* Copy to the temporary space. */
154: MPN_COPY (up, wp, usize);
155: }
156: else if (wp == vp)
157: {
158: /* W and V are identical. Allocate temporary space for V. */
159: vp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
160: /* Copy to the temporary space. */
161: MPN_COPY (vp, wp, vsize);
162: }
163: }
164:
165: wsize = _mpn_mul_classic (wp, up, usize, vp, vsize);
166: w->_mp_size = sign_product < 0 ? -wsize : wsize;
167: if (free_me != NULL)
168: (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
169:
170: TMP_FREE (marker);
171: }
172:
173: mp_size_t
174: _mpn_mul_classic (prodp, up, usize, vp, vsize)
175: mp_ptr prodp;
176: mp_srcptr up;
177: mp_size_t usize;
178: mp_srcptr vp;
179: mp_size_t vsize;
180: {
181: mp_size_t i, j;
182: mp_limb_t prod_low, prod_high;
183: mp_limb_t cy_dig;
184: mp_limb_t v_limb, c;
185:
186: if (vsize == 0)
187: return 0;
188:
189: /* Offset UP and PRODP so that the inner loop can be faster. */
190: up += usize;
191: prodp += usize;
192:
193: /* Multiply by the first limb in V separately, as the result can
194: be stored (not added) to PROD. We also avoid a loop for zeroing. */
195: v_limb = vp[0];
196: cy_dig = 0;
197: j = -usize;
198: do
199: {
200: umul_ppmm (prod_high, prod_low, up[j], v_limb);
201: add_ssaaaa (cy_dig, prodp[j], prod_high, prod_low, 0, cy_dig);
202: j++;
203: }
204: while (j < 0);
205:
206: prodp[j] = cy_dig;
207: prodp++;
208:
209: /* For each iteration in the outer loop, multiply one limb from
210: U with one limb from V, and add it to PROD. */
211: for (i = 1; i < vsize; i++)
212: {
213: v_limb = vp[i];
214: cy_dig = 0;
215: j = -usize;
216:
217: /* Inner loops. Simulate the carry flag by jumping between
218: these loops. The first is used when there was no carry
219: in the previois iteration; the second when there was carry. */
220:
221: do
222: {
223: umul_ppmm (prod_high, prod_low, up[j], v_limb);
224: add_ssaaaa (cy_dig, prod_low, prod_high, prod_low, 0, cy_dig);
225: c = prodp[j];
226: prod_low += c;
227: prodp[j] = prod_low;
228: if (prod_low < c)
229: goto cy_loop;
230: ncy_loop:
231: j++;
232: }
233: while (j < 0);
234:
235: prodp[j] = cy_dig;
236: prodp++;
237: continue;
238:
239: do
240: {
241: umul_ppmm (prod_high, prod_low, up[j], v_limb);
242: add_ssaaaa (cy_dig, prod_low, prod_high, prod_low, 0, cy_dig);
243: c = prodp[j];
244: prod_low += c + 1;
245: prodp[j] = prod_low;
246: if (prod_low > c)
247: goto ncy_loop;
248: cy_loop:
249: j++;
250: }
251: while (j < 0);
252:
253: cy_dig += 1;
254: prodp[j] = cy_dig;
255: prodp++;
256: }
257:
258: return usize + vsize - (cy_dig == 0);
259: }
260:
1.1.1.2 ! maekawa 261: dump_abort (s, multiplier, multiplicand, product, ref_product)
! 262: char *s;
! 263: mpz_t multiplier, multiplicand, product, ref_product;
1.1 maekawa 264: {
1.1.1.2 ! maekawa 265: fprintf (stderr, "ERROR: %s\n", s);
1.1 maekawa 266: fprintf (stderr, "multiplier = "); debug_mp (multiplier, -16);
267: fprintf (stderr, "multiplicand = "); debug_mp (multiplicand, -16);
1.1.1.2 ! maekawa 268: fprintf (stderr, "product = "); debug_mp (product, -16);
! 269: fprintf (stderr, "ref_product = "); debug_mp (ref_product, -16);
1.1 maekawa 270: abort();
271: }
272:
273: void
274: debug_mp (x, base)
275: mpz_t x;
276: {
277: mpz_out_str (stderr, base, x); fputc ('\n', stderr);
278: }
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