1 /* Test mpz_cmp, mpz_mul. 2 3 Copyright 1991, 1993, 1994, 1996, 1997, 2000-2004 Free Software Foundation, 4 Inc. 5 6 This file is part of the GNU MP Library test suite. 7 8 The GNU MP Library test suite is free software; you can redistribute it 9 and/or modify it under the terms of the GNU General Public License as 10 published by the Free Software Foundation; either version 3 of the License, 11 or (at your option) any later version. 12 13 The GNU MP Library test suite is distributed in the hope that it will be 14 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 16 Public License for more details. 17 18 You should have received a copy of the GNU General Public License along with 19 the 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 28 void debug_mp (mpz_t); 29 static void refmpz_mul (mpz_t, const mpz_t, const mpz_t); 30 void dump_abort (int, const char *, mpz_t, mpz_t, mpz_t, mpz_t); 31 32 #define FFT_MIN_BITSIZE 100000 33 34 char *extra_fft; 35 36 void 37 one (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 63 int 64 main (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 141 static void 142 refmpz_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 182 void 183 dump_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 202 void 203 debug_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
Indexes created Sat Feb 28 05:31:39 UTC 2026