1 /* Test for mulmod_bnm1 function. 2 3 Contributed to the GNU project by Marco Bodrato. 4 5 Copyright 2009 Free Software Foundation, Inc. 6 7 This file is part of the GNU MP Library test suite. 8 9 The GNU MP Library test suite is free software; you can redistribute it 10 and/or modify it under the terms of the GNU General Public License as 11 published by the Free Software Foundation; either version 3 of the License, 12 or (at your option) any later version. 13 14 The GNU MP Library test suite is distributed in the hope that it will be 15 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 17 Public License for more details. 18 19 You should have received a copy of the GNU General Public License along with 20 the GNU MP Library test suite. If not, see http://www.gnu.org/licenses/. */ 21 22 23 #include <stdlib.h> 24 #include <stdio.h> 25 26 #include "gmp.h" 27 #include "gmp-impl.h" 28 #include "tests.h" 29 30 /* Sizes are up to 2^SIZE_LOG limbs */ 31 #ifndef SIZE_LOG 32 #define SIZE_LOG 11 33 #endif 34 35 #ifndef COUNT 36 #define COUNT 5000 37 #endif 38 39 #define MAX_N (1L << SIZE_LOG) 40 #define MIN_N 1 41 42 /* 43 Reference function for multiplication modulo B^rn-1. 44 45 The result is expected to be ZERO if and only if one of the operand 46 already is. Otherwise the class [0] Mod(B^rn-1) is represented by 47 B^rn-1. This should not be a problem if mulmod_bnm1 is used to 48 combine results and obtain a natural number when one knows in 49 advance that the final value is less than (B^rn-1). 50 */ 51 52 static void 53 ref_mulmod_bnm1 (mp_ptr rp, mp_size_t rn, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) 54 { 55 mp_limb_t cy; 56 57 ASSERT (0 < an && an <= rn); 58 ASSERT (0 < bn && bn <= rn); 59 60 if (an >= bn) 61 refmpn_mul (rp, ap, an, bp, bn); 62 else 63 refmpn_mul (rp, bp, bn, ap, an); 64 an += bn; 65 if (an > rn) { 66 cy = mpn_add (rp, rp, rn, rp + rn, an - rn); 67 /* If cy == 1, then the value of rp is at most B^rn - 2, so there can 68 * be no overflow when adding in the carry. */ 69 MPN_INCR_U (rp, rn, cy); 70 } 71 } 72 73 /* 74 Compare the result of the mpn_mulmod_bnm1 function in the library 75 with the reference function above. 76 */ 77 78 int 79 main (int argc, char **argv) 80 { 81 mp_ptr ap, bp, refp, pp, scratch; 82 int count = COUNT; 83 int test; 84 gmp_randstate_ptr rands; 85 TMP_DECL; 86 TMP_MARK; 87 88 if (argc > 1) 89 { 90 char *end; 91 count = strtol (argv[1], &end, 0); 92 if (*end || count <= 0) 93 { 94 fprintf (stderr, "Invalid test count: %s.\n", argv[1]); 95 return 1; 96 } 97 } 98 99 tests_start (); 100 rands = RANDS; 101 102 ASSERT_ALWAYS (mpn_mulmod_bnm1_next_size (MAX_N) == MAX_N); 103 104 ap = TMP_ALLOC_LIMBS (MAX_N); 105 bp = TMP_ALLOC_LIMBS (MAX_N); 106 refp = TMP_ALLOC_LIMBS (MAX_N * 4); 107 pp = 1+TMP_ALLOC_LIMBS (MAX_N + 2); 108 scratch 109 = 1+TMP_ALLOC_LIMBS (mpn_mulmod_bnm1_itch (MAX_N, MAX_N, MAX_N) + 2); 110 111 for (test = 0; test < count; test++) 112 { 113 unsigned size_min; 114 unsigned size_range; 115 mp_size_t an,bn,rn,n; 116 mp_size_t itch; 117 mp_limb_t p_before, p_after, s_before, s_after; 118 119 for (size_min = 1; (1L << size_min) < MIN_N; size_min++) 120 ; 121 122 /* We generate an in the MIN_N <= n <= (1 << size_range). */ 123 size_range = size_min 124 + gmp_urandomm_ui (rands, SIZE_LOG + 1 - size_min); 125 126 n = MIN_N 127 + gmp_urandomm_ui (rands, (1L << size_range) + 1 - MIN_N); 128 n = mpn_mulmod_bnm1_next_size (n); 129 130 if ( (test & 1) || n == 1) { 131 /* Half of the tests are done with the main scenario in mind: 132 both an and bn >= rn/2 */ 133 an = ((n+1) >> 1) + gmp_urandomm_ui (rands, (n+1) >> 1); 134 bn = ((n+1) >> 1) + gmp_urandomm_ui (rands, (n+1) >> 1); 135 } else { 136 /* Second half of the tests are done using mulmod to compute a 137 full product with n/2 < an+bn <= n. */ 138 an = 1 + gmp_urandomm_ui (rands, n - 1); 139 if (an >= n/2) 140 bn = 1 + gmp_urandomm_ui (rands, n - an); 141 else 142 bn = n/2 + 1 - an + gmp_urandomm_ui (rands, (n+1)/2); 143 } 144 145 /* Make sure an >= bn */ 146 if (an < bn) 147 MP_SIZE_T_SWAP (an, bn); 148 149 mpn_random2 (ap, an); 150 mpn_random2 (bp, bn); 151 152 /* Sometime trigger the borderline conditions 153 A = -1,0,+1 or B = -1,0,+1 or A*B == -1,0,1 Mod(B^{n/2}+1). 154 This only makes sense if there is at least a split, i.e. n is even. */ 155 if ((test & 0x1f) == 1 && (n & 1) == 0) { 156 mp_size_t x; 157 MPN_COPY (ap, ap + (n >> 1), an - (n >> 1)); 158 MPN_ZERO (ap + an - (n >> 1) , n - an); 159 MPN_COPY (bp, bp + (n >> 1), bn - (n >> 1)); 160 MPN_ZERO (bp + bn - (n >> 1) , n - bn); 161 x = (n == an) ? 0 : gmp_urandomm_ui (rands, n - an); 162 ap[x] += gmp_urandomm_ui (rands, 3) - 1; 163 x = (n >> 1) - x % (n >> 1); 164 bp[x] += gmp_urandomm_ui (rands, 3) - 1; 165 /* We don't propagate carry, this means that the desired condition 166 is not triggered all the times. A few times are enough anyway. */ 167 } 168 rn = MIN(n, an + bn); 169 mpn_random2 (pp-1, rn + 2); 170 p_before = pp[-1]; 171 p_after = pp[rn]; 172 173 itch = mpn_mulmod_bnm1_itch (n, an, bn); 174 ASSERT_ALWAYS (itch <= mpn_mulmod_bnm1_itch (MAX_N, MAX_N, MAX_N)); 175 mpn_random2 (scratch-1, itch+2); 176 s_before = scratch[-1]; 177 s_after = scratch[itch]; 178 179 mpn_mulmod_bnm1 ( pp, n, ap, an, bp, bn, scratch); 180 ref_mulmod_bnm1 (refp, n, ap, an, bp, bn); 181 if (pp[-1] != p_before || pp[rn] != p_after 182 || scratch[-1] != s_before || scratch[itch] != s_after 183 || mpn_cmp (refp, pp, rn) != 0) 184 { 185 printf ("ERROR in test %d, an = %d, bn = %d, n = %d\n", 186 test, (int) an, (int) bn, (int) n); 187 if (pp[-1] != p_before) 188 { 189 printf ("before pp:"); mpn_dump (pp -1, 1); 190 printf ("keep: "); mpn_dump (&p_before, 1); 191 } 192 if (pp[rn] != p_after) 193 { 194 printf ("after pp:"); mpn_dump (pp + rn, 1); 195 printf ("keep: "); mpn_dump (&p_after, 1); 196 } 197 if (scratch[-1] != s_before) 198 { 199 printf ("before scratch:"); mpn_dump (scratch-1, 1); 200 printf ("keep: "); mpn_dump (&s_before, 1); 201 } 202 if (scratch[itch] != s_after) 203 { 204 printf ("after scratch:"); mpn_dump (scratch + itch, 1); 205 printf ("keep: "); mpn_dump (&s_after, 1); 206 } 207 mpn_dump (ap, an); 208 mpn_dump (bp, bn); 209 mpn_dump (pp, rn); 210 mpn_dump (refp, rn); 211 212 abort(); 213 } 214 } 215 TMP_FREE; 216 tests_end (); 217 return 0; 218 } 219