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