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