1/* Test for mulmod_bnm1 function. 2 3 Contributed to the GNU project by Marco Bodrato. 4 5Copyright 2009 Free Software Foundation, Inc. 6 7This file is part of the GNU MP Library. 8 9The GNU MP Library is free software; you can redistribute it and/or modify 10it under the terms of the GNU Lesser General Public License as published by 11the Free Software Foundation; either version 3 of the License, or (at your 12option) any later version. 13 14The GNU MP Library is distributed in the hope that it will be useful, but 15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17License for more details. 18 19You should have received a copy of the GNU Lesser General Public License 20along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 21 22 23#include "gmp.h" 24#include "gmp-impl.h" 25#include "tests.h" 26 27#include <stdlib.h> 28#include <stdio.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 52static void 53ref_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 78int 79main (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