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