1 /* Test for mulmod_bknp1 function.
2
3 Contributed to the GNU project by Marco Bodrato.
4
5 Copyright 2009, 2020-2022 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
30 #if MOD_BKNP1_USE11
31 #define USE11 11,
32 #else
33 #define USE11
34 #endif
35
36 #if GMP_NUMB_BITS % 32 == 0
37 #define MAX_K 17
38 #define SUPPORTED_K {3, 5, 7, 13, USE11 MAX_K}
39 #else
40 #if GMP_NUMB_BITS % 16 == 0
41 #define MAX_K 13
42 #define SUPPORTED_K {3, 5, 7, USE11 MAX_K}
43 #else
44 #if GMP_NUMB_BITS % 8 == 0
45 #define MAX_K 7
46 #define SUPPORTED_K {3, USE11 MAX_K}
47 #else
48 #define SUPPORTED_K {USE11} /* Supported ? */
49 #endif /* GMP_NUMB_BITS % 8 == 0 */
50 #endif /* GMP_NUMB_BITS % 16 == 0 */
51 #endif /* GMP_NUMB_BITS % 32 == 0 */
52
53 #if MOD_BKNP1_ONLY3
54 #undef SUPPORTED_K
55 #undef MAX_K
56 #define MAX_K 3
57 #define SUPPORTED_K {3}
58 #endif
59
60 /* Sizes are up to MAX_K * 2^SIZE_LOG limbs */
61 #ifndef SIZE_LOG
62 #define SIZE_LOG 7
63 #endif
64
65 #ifndef COUNT
66 #define COUNT 5000
67 #endif
68
69 #define MAX_N (MAX_K << SIZE_LOG)
70 #define MIN_N 1
71
72 /*
73 Reference function for multiplication modulo B^{k*rn}+1.
74 */
75
76 static void
77 ref_sqrmod_bnp1 (mp_ptr rp, mp_srcptr ap, mp_size_t rn)
78 {
79 mp_limb_t cy;
80
81 mpn_sqr (rp, ap, rn + 1);
82 cy = rp[2 * rn];
83 MPN_INCR_U (rp, 2 * rn + 1, rp[2 * rn]);
84 cy = rp[2 * rn] - cy + mpn_sub_n (rp, rp, rp + rn, rn);
85 rp[rn] = 0;
86 MPN_INCR_U (rp, rn + 1, cy);
87 }
88
89 /*
90 Compare the result of the mpn_mulmod_bnp1 function in the library
91 with the reference function above.
92 */
93 unsigned supported_k[] = SUPPORTED_K;
94
95 int
96 main (int argc, char **argv)
97 {
98 mp_ptr ap, refp, pp, scratch;
99 int count = COUNT;
100 int test;
101 gmp_randstate_ptr rands;
102 TMP_DECL;
103 TMP_MARK;
104
105 TESTS_REPS (count, argv, argc);
106
107 tests_start ();
108 rands = RANDS;
109
110 ap = TMP_ALLOC_LIMBS (MAX_N + 1);
111 refp = TMP_ALLOC_LIMBS (MAX_N * 2 + 2);
112 pp = 1 + TMP_ALLOC_LIMBS (MAX_N + 3);
113 scratch
114 = 1 + TMP_ALLOC_LIMBS (mpn_mulmod_bknp1_itch (MAX_N) + 2);
115
116 for (test = 0; test < count; test++)
117 {
118 unsigned size_min;
119 unsigned size_range;
120 unsigned k;
121 mp_size_t rn, n;
122 mp_size_t itch;
123 mp_limb_t p_before, p_after, s_before, s_after;
124
125 for (size_min = 1; (1L << size_min) < MIN_N; size_min++)
126 ;
127
128 /* We generate rn in the MIN_N <= n <= (1 << size_range). */
129 size_range = size_min
130 + gmp_urandomm_ui (rands, SIZE_LOG + 1 - size_min);
131
132 k = supported_k[test % numberof (supported_k)];
133 if (test < numberof (supported_k))
134 {
135 n = 1;
136 rn = k;
137 ap [rn] = 0;
138 mp_limb_t x = GMP_NUMB_MAX / k + 1;
139 ap [0] = x;
140 for (int i = 1; i < k; i += 2)
141 {
142 ap [i] = - x;
143 ap [i + 1] = x - 1;
144 }
145 }
146 else
147 {
148 n = MIN_N
149 + gmp_urandomm_ui (rands, (1L << size_range) + 1 - MIN_N);
150 rn = k * n;
151 if ((GMP_NUMB_MAX % k != 0) && (rn % 3 == 0))
152 n = rn / (k = 3);
153
154 mpn_random2 (ap, rn + 1);
155
156 ap [rn] &= 1;
157 }
158
159 mpn_random2 (pp-1, rn + 3);
160 p_before = pp[-1];
161 p_after = pp[rn + 1];
162
163 itch = mpn_sqrmod_bknp1_itch (rn);
164 ASSERT_ALWAYS (itch <= mpn_mulmod_bknp1_itch (MAX_N));
165 mpn_random2 (scratch - 1, itch + 2);
166 s_before = scratch[-1];
167 s_after = scratch[itch];
168
169 mpn_sqrmod_bknp1 ( pp, ap, n, k, scratch);
170 ref_sqrmod_bnp1 (refp, ap, rn);
171 if (pp[-1] != p_before || pp[rn + 1] != p_after
172 || scratch[-1] != s_before || scratch[itch] != s_after
173 || mpn_cmp (refp, pp, rn + 1) != 0)
174 {
175 printf ("ERROR in test %d(sqr), rn = %d, n = %d, k = %d\n",
176 test, (int) rn, (int) n, (int) k);
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 + 1] != p_after)
183 {
184 printf ("after pp:"); mpn_dump (pp + rn + 1, 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, rn + 1);
198 mpn_dump (pp, rn + 1);
199 mpn_dump (refp, rn + 1);
200
201 abort();
202 }
203
204 mpn_random2 (pp-1, rn + 3);
205 p_before = pp[-1];
206 p_after = pp[rn + 1];
207
208 itch = mpn_mulmod_bknp1_itch (rn);
209 ASSERT_ALWAYS (itch <= mpn_mulmod_bknp1_itch (MAX_N));
210 mpn_random2 (scratch - 1, itch + 2);
211 s_before = scratch[-1];
212 s_after = scratch[itch];
213
214 mpn_mulmod_bknp1 ( pp, ap, ap, n, k, scratch);
215 if (pp[-1] != p_before || pp[rn + 1] != p_after
216 || scratch[-1] != s_before || scratch[itch] != s_after
217 || mpn_cmp (refp, pp, rn + 1) != 0)
218 {
219 printf ("ERROR in test %d(mul), rn = %d, n = %d, k = %d\n",
220 test, (int) rn, (int) n, (int) k);
221 if (pp[-1] != p_before)
222 {
223 printf ("before pp:"); mpn_dump (pp - 1, 1);
224 printf ("keep: "); mpn_dump (&p_before, 1);
225 }
226 if (pp[rn + 1] != p_after)
227 {
228 printf ("after pp:"); mpn_dump (pp + rn + 1, 1);
229 printf ("keep: "); mpn_dump (&p_after, 1);
230 }
231 if (scratch[-1] != s_before)
232 {
233 printf ("before scratch:"); mpn_dump (scratch - 1, 1);
234 printf ("keep: "); mpn_dump (&s_before, 1);
235 }
236 if (scratch[itch] != s_after)
237 {
238 printf ("after scratch:"); mpn_dump (scratch + itch, 1);
239 printf ("keep: "); mpn_dump (&s_after, 1);
240 }
241 mpn_dump (ap, rn + 1);
242 mpn_dump (pp, rn + 1);
243 mpn_dump (refp, rn + 1);
244
245 abort();
246 }
247 }
248 TMP_FREE;
249 tests_end ();
250 return 0;
251 }