1/* Test for sqrmod_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 12 33#endif 34 35#ifndef COUNT 36#define COUNT 3000 37#endif 38 39#define MAX_N (1L << SIZE_LOG) 40#define MIN_N 1 41 42/* 43 Reference function for squaring 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 sqrmod_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_sqrmod_bnm1 (mp_ptr rp, mp_size_t rn, mp_srcptr ap, mp_size_t an) 54{ 55 mp_limb_t cy; 56 57 ASSERT (0 < an && an <= rn); 58 59 refmpn_mul (rp, ap, an, ap, an); 60 an *= 2; 61 if (an > rn) { 62 cy = mpn_add (rp, rp, rn, rp + rn, an - rn); 63 /* If cy == 1, then the value of rp is at most B^rn - 2, so there can 64 * be no overflow when adding in the carry. */ 65 MPN_INCR_U (rp, rn, cy); 66 } 67} 68 69/* 70 Compare the result of the mpn_sqrmod_bnm1 function in the library 71 with the reference function above. 72*/ 73 74int 75main (int argc, char **argv) 76{ 77 mp_ptr ap, refp, pp, scratch; 78 int count = COUNT; 79 int test; 80 gmp_randstate_ptr rands; 81 TMP_DECL; 82 TMP_MARK; 83 84 if (argc > 1) 85 { 86 char *end; 87 count = strtol (argv[1], &end, 0); 88 if (*end || count <= 0) 89 { 90 fprintf (stderr, "Invalid test count: %s.\n", argv[1]); 91 return 1; 92 } 93 } 94 95 tests_start (); 96 rands = RANDS; 97 98 ASSERT_ALWAYS (mpn_sqrmod_bnm1_next_size (MAX_N) == MAX_N); 99 100 ap = TMP_ALLOC_LIMBS (MAX_N); 101 refp = TMP_ALLOC_LIMBS (MAX_N * 4); 102 pp = 1+TMP_ALLOC_LIMBS (MAX_N + 2); 103 scratch 104 = 1+TMP_ALLOC_LIMBS (mpn_sqrmod_bnm1_itch (MAX_N, MAX_N) + 2); 105 106 for (test = 0; test < count; test++) 107 { 108 unsigned size_min; 109 unsigned size_range; 110 mp_size_t an,rn,n; 111 mp_size_t itch; 112 mp_limb_t p_before, p_after, s_before, s_after; 113 114 for (size_min = 1; (1L << size_min) < MIN_N; size_min++) 115 ; 116 117 /* We generate an in the MIN_N <= n <= (1 << size_range). */ 118 size_range = size_min 119 + gmp_urandomm_ui (rands, SIZE_LOG + 1 - size_min); 120 121 n = MIN_N 122 + gmp_urandomm_ui (rands, (1L << size_range) + 1 - MIN_N); 123 n = mpn_sqrmod_bnm1_next_size (n); 124 125 if (n == 1) 126 an = 1; 127 else 128 an = ((n+1) >> 1) + gmp_urandomm_ui (rands, (n+1) >> 1); 129 130 mpn_random2 (ap, an); 131 132 /* Sometime trigger the borderline conditions 133 A = -1,0,+1 Mod(B^{n/2}+1). 134 This only makes sense if there is at least a split, i.e. n is even. */ 135 if ((test & 0x1f) == 1 && (n & 1) == 0) { 136 mp_size_t x; 137 MPN_COPY (ap, ap + (n >> 1), an - (n >> 1)); 138 MPN_ZERO (ap + an - (n >> 1) , n - an); 139 x = (n == an) ? 0 : gmp_urandomm_ui (rands, n - an); 140 ap[x] += gmp_urandomm_ui (rands, 3) - 1; 141 } 142 rn = MIN(n, 2*an); 143 mpn_random2 (pp-1, rn + 2); 144 p_before = pp[-1]; 145 p_after = pp[rn]; 146 147 itch = mpn_sqrmod_bnm1_itch (n, an); 148 ASSERT_ALWAYS (itch <= mpn_sqrmod_bnm1_itch (MAX_N, MAX_N)); 149 mpn_random2 (scratch-1, itch+2); 150 s_before = scratch[-1]; 151 s_after = scratch[itch]; 152 153 mpn_sqrmod_bnm1 ( pp, n, ap, an, scratch); 154 ref_sqrmod_bnm1 (refp, n, ap, an); 155 if (pp[-1] != p_before || pp[rn] != p_after 156 || scratch[-1] != s_before || scratch[itch] != s_after 157 || mpn_cmp (refp, pp, rn) != 0) 158 { 159 printf ("ERROR in test %d, an = %d, n = %d\n", 160 test, (int) an, (int) n); 161 if (pp[-1] != p_before) 162 { 163 printf ("before pp:"); mpn_dump (pp -1, 1); 164 printf ("keep: "); mpn_dump (&p_before, 1); 165 } 166 if (pp[rn] != p_after) 167 { 168 printf ("after pp:"); mpn_dump (pp + rn, 1); 169 printf ("keep: "); mpn_dump (&p_after, 1); 170 } 171 if (scratch[-1] != s_before) 172 { 173 printf ("before scratch:"); mpn_dump (scratch-1, 1); 174 printf ("keep: "); mpn_dump (&s_before, 1); 175 } 176 if (scratch[itch] != s_after) 177 { 178 printf ("after scratch:"); mpn_dump (scratch + itch, 1); 179 printf ("keep: "); mpn_dump (&s_after, 1); 180 } 181 mpn_dump (ap, an); 182 mpn_dump (pp, rn); 183 mpn_dump (refp, rn); 184 185 abort(); 186 } 187 } 188 TMP_FREE; 189 tests_end (); 190 return 0; 191} 192