1/* Test mpz_cmp, mpz_mul. 2 3Copyright 1991, 1993, 1994, 1996, 1997, 2000-2004 Free Software Foundation, 4Inc. 5 6This file is part of the GNU MP Library test suite. 7 8The GNU MP Library test suite is free software; you can redistribute it 9and/or modify it under the terms of the GNU General Public License as 10published by the Free Software Foundation; either version 3 of the License, 11or (at your option) any later version. 12 13The GNU MP Library test suite is distributed in the hope that it will be 14useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 16Public License for more details. 17 18You should have received a copy of the GNU General Public License along with 19the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ 20 21#include <stdio.h> 22#include <stdlib.h> 23 24#include "gmp-impl.h" 25#include "longlong.h" 26#include "tests.h" 27 28void debug_mp (mpz_t); 29static void refmpz_mul (mpz_t, const mpz_t, const mpz_t); 30void dump_abort (int, const char *, mpz_t, mpz_t, mpz_t, mpz_t); 31 32#define FFT_MIN_BITSIZE 100000 33 34char *extra_fft; 35 36void 37one (int i, mpz_t multiplicand, mpz_t multiplier) 38{ 39 mpz_t product, ref_product; 40 41 mpz_init (product); 42 mpz_init (ref_product); 43 44 /* Test plain multiplication comparing results against reference code. */ 45 mpz_mul (product, multiplier, multiplicand); 46 refmpz_mul (ref_product, multiplier, multiplicand); 47 if (mpz_cmp (product, ref_product)) 48 dump_abort (i, "incorrect plain product", 49 multiplier, multiplicand, product, ref_product); 50 51 /* Test squaring, comparing results against plain multiplication */ 52 mpz_mul (product, multiplier, multiplier); 53 mpz_set (multiplicand, multiplier); 54 mpz_mul (ref_product, multiplier, multiplicand); 55 if (mpz_cmp (product, ref_product)) 56 dump_abort (i, "incorrect square product", 57 multiplier, multiplier, product, ref_product); 58 59 mpz_clear (product); 60 mpz_clear (ref_product); 61} 62 63int 64main (int argc, char **argv) 65{ 66 mpz_t op1, op2; 67 int i; 68 int fft_max_2exp; 69 70 gmp_randstate_ptr rands; 71 mpz_t bs; 72 unsigned long bsi, size_range, fsize_range; 73 74 tests_start (); 75 rands = RANDS; 76 77 extra_fft = getenv ("GMP_CHECK_FFT"); 78 fft_max_2exp = 0; 79 if (extra_fft != NULL) 80 { 81 fft_max_2exp = atoi (extra_fft); 82 printf ("GMP_CHECK_FFT=%d (include this in bug reports)\n", fft_max_2exp); 83 } 84 85 if (fft_max_2exp <= 1) /* compat with old use of GMP_CHECK_FFT */ 86 fft_max_2exp = 22; /* default limit, good for any machine */ 87 88 mpz_init (bs); 89 mpz_init (op1); 90 mpz_init (op2); 91 92 fsize_range = 4 << 8; /* a fraction 1/256 of size_range */ 93 for (i = 0;; i++) 94 { 95 size_range = fsize_range >> 8; 96 fsize_range = fsize_range * 33 / 32; 97 98 if (size_range > fft_max_2exp) 99 break; 100 101 mpz_urandomb (bs, rands, size_range); 102 mpz_rrandomb (op1, rands, mpz_get_ui (bs)); 103 if (i & 1) 104 mpz_urandomb (bs, rands, size_range); 105 mpz_rrandomb (op2, rands, mpz_get_ui (bs)); 106 107 mpz_urandomb (bs, rands, 4); 108 bsi = mpz_get_ui (bs); 109 if ((bsi & 0x3) == 0) 110 mpz_neg (op1, op1); 111 if ((bsi & 0xC) == 0) 112 mpz_neg (op2, op2); 113 114 /* printf ("%d %d\n", SIZ (op1), SIZ (op2)); */ 115 one (i, op2, op1); 116 } 117 118 for (i = -50; i < 0; i++) 119 { 120 mpz_urandomb (bs, rands, 32); 121 size_range = mpz_get_ui (bs) % fft_max_2exp; 122 123 mpz_urandomb (bs, rands, size_range); 124 mpz_rrandomb (op1, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE); 125 mpz_urandomb (bs, rands, size_range); 126 mpz_rrandomb (op2, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE); 127 128 /* printf ("%d: %d %d\n", i, SIZ (op1), SIZ (op2)); */ 129 fflush (stdout); 130 one (-1, op2, op1); 131 } 132 133 mpz_clear (bs); 134 mpz_clear (op1); 135 mpz_clear (op2); 136 137 tests_end (); 138 exit (0); 139} 140 141static void 142refmpz_mul (mpz_t w, const mpz_t u, const mpz_t v) 143{ 144 mp_size_t usize = u->_mp_size; 145 mp_size_t vsize = v->_mp_size; 146 mp_size_t wsize; 147 mp_size_t sign_product; 148 mp_ptr up, vp; 149 mp_ptr wp; 150 mp_size_t talloc; 151 152 sign_product = usize ^ vsize; 153 usize = ABS (usize); 154 vsize = ABS (vsize); 155 156 if (usize == 0 || vsize == 0) 157 { 158 SIZ (w) = 0; 159 return; 160 } 161 162 talloc = usize + vsize; 163 164 up = u->_mp_d; 165 vp = v->_mp_d; 166 167 wp = __GMP_ALLOCATE_FUNC_LIMBS (talloc); 168 169 if (usize > vsize) 170 refmpn_mul (wp, up, usize, vp, vsize); 171 else 172 refmpn_mul (wp, vp, vsize, up, usize); 173 wsize = usize + vsize; 174 wsize -= wp[wsize - 1] == 0; 175 MPZ_REALLOC (w, wsize); 176 MPN_COPY (PTR(w), wp, wsize); 177 178 SIZ(w) = sign_product < 0 ? -wsize : wsize; 179 __GMP_FREE_FUNC_LIMBS (wp, talloc); 180} 181 182void 183dump_abort (int i, const char *s, 184 mpz_t op1, mpz_t op2, mpz_t product, mpz_t ref_product) 185{ 186 mp_size_t b, e; 187 fprintf (stderr, "ERROR: %s in test %d\n", s, i); 188 fprintf (stderr, "op1 = "); debug_mp (op1); 189 fprintf (stderr, "op2 = "); debug_mp (op2); 190 fprintf (stderr, " product = "); debug_mp (product); 191 fprintf (stderr, "ref_product = "); debug_mp (ref_product); 192 for (b = 0; b < ABSIZ(ref_product); b++) 193 if (PTR(ref_product)[b] != PTR(product)[b]) 194 break; 195 for (e = ABSIZ(ref_product) - 1; e >= 0; e--) 196 if (PTR(ref_product)[e] != PTR(product)[e]) 197 break; 198 printf ("ERRORS in %ld--%ld\n", b, e); 199 abort(); 200} 201 202void 203debug_mp (mpz_t x) 204{ 205 size_t siz = mpz_sizeinbase (x, 16); 206 207 if (siz > 65) 208 { 209 mpz_t q; 210 mpz_init (q); 211 mpz_tdiv_q_2exp (q, x, 4 * (mpz_sizeinbase (x, 16) - 25)); 212 gmp_fprintf (stderr, "%ZX...", q); 213 mpz_tdiv_r_2exp (q, x, 4 * 25); 214 gmp_fprintf (stderr, "%025ZX [%d]\n", q, (int) siz); 215 mpz_clear (q); 216 } 217 else 218 { 219 gmp_fprintf (stderr, "%ZX\n", x); 220 } 221} 222