1 /* Test mpz_add, mpz_cmp, mpz_cmp_ui, mpz_divmod, mpz_mul.
3 Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
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.
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.
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. */
29 mp_size_t _mpn_mul_classic ();
40 mpz_t multiplier, multiplicand;
41 mpz_t product, ref_product;
42 mpz_t quotient, remainder;
43 mp_size_t multiplier_size, multiplicand_size;
48 reps = atoi (argv[1]);
50 mpz_init (multiplier);
51 mpz_init (multiplicand);
53 mpz_init (ref_product);
57 for (i = 0; i < reps; i++)
59 multiplier_size = urandom () % SIZE - SIZE/2;
60 multiplicand_size = urandom () % SIZE - SIZE/2;
62 mpz_random2 (multiplier, multiplier_size);
63 mpz_random2 (multiplicand, multiplicand_size);
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);
70 if (mpz_cmp (product, ref_product))
71 dump_abort (multiplier, multiplicand);
73 if (mpz_cmp_ui (multiplicand, 0) != 0)
74 if (mpz_cmp_ui (remainder, 0) || mpz_cmp (quotient, multiplier))
75 dump_abort (multiplier, multiplicand);
87 mp_size_t usize = u->_mp_size;
88 mp_size_t vsize = v->_mp_size;
90 mp_size_t sign_product;
93 mp_ptr free_me = NULL;
98 sign_product = usize ^ vsize;
105 {const __mpz_struct *t = u; u = v; v = t;}
106 {mp_size_t t = usize; usize = vsize; vsize = t;}
113 /* Ensure W has space enough to store the result. */
114 wsize = usize + vsize;
115 if (w->_mp_alloc < wsize)
117 if (wp == up || wp == vp)
120 free_me_size = w->_mp_alloc;
123 (*_mp_free_func) (wp, w->_mp_alloc * BYTES_PER_MP_LIMB);
125 w->_mp_alloc = wsize;
126 wp = (mp_ptr) (*_mp_allocate_func) (wsize * BYTES_PER_MP_LIMB);
131 /* Make U and V not overlap with W. */
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. */
139 /* Copy to the temporary space. */
140 MPN_COPY (up, wp, usize);
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);
151 wsize = _mpn_mul_classic (wp, up, usize, vp, vsize);
152 w->_mp_size = sign_product < 0 ? -wsize : wsize;
154 (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
160 _mpn_mul_classic (prodp, up, usize, vp, vsize)
168 mp_limb_t prod_low, prod_high;
175 /* Offset UP and PRODP so that the inner loop can be faster. */
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. */
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);
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++)
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. */
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);
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);
244 return usize + vsize - (cy_dig == 0);
247 dump_abort (multiplier, multiplicand)
248 mpz_t multiplier, multiplicand;
250 fprintf (stderr, "ERROR\n");
251 fprintf (stderr, "multiplier = "); debug_mp (multiplier, -16);
252 fprintf (stderr, "multiplicand = "); debug_mp (multiplicand, -16);
260 mpz_out_str (stderr, base, x); fputc ('\n', stderr);