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 #if MOD_BKNP1_USE11
30 #define USE11 11,
31 #else
32 #define USE11
33 #endif
34
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_mulmod_bnp1 (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t rn)
78 {
79 mp_limb_t cy;
80
81 mpn_mul_n (rp, ap, bp, 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, bp, 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 bp = TMP_ALLOC_LIMBS (MAX_N + 1);
112 refp = TMP_ALLOC_LIMBS (MAX_N * 2 + 2);
113 pp = 1 + TMP_ALLOC_LIMBS (MAX_N + 3);
114 scratch
115 = 1 + TMP_ALLOC_LIMBS (mpn_mulmod_bknp1_itch (MAX_N) + 2);
116
117 for (test = 0; test < count; test++)
118 {
119 unsigned size_min;
120 unsigned size_range;
121 unsigned k;
122 mp_size_t rn, n;
123 mp_size_t itch;
124 mp_limb_t p_before, p_after, s_before, s_after;
125
126 for (size_min = 1; (1L << size_min) < MIN_N; size_min++)
127 ;
128
129 /* We generate rn in the MIN_N <= n <= (1 << size_range). */
130 size_range = size_min
131 + gmp_urandomm_ui (rands, SIZE_LOG + 1 - size_min);
132
133 k = supported_k[test % numberof (supported_k)];
134 n = MIN_N
135 + gmp_urandomm_ui (rands, (1L << size_range) + 1 - MIN_N);
136 rn = k * n;
137 if ((GMP_NUMB_MAX % k != 0) && (rn % 3 == 0))
138 n = rn / (k = 3);
139
140 if (test == 0)
141 {
142 mpn_random2 (ap, n);
143 mpn_add_1 (ap + n, ap, n, 1); /* {ap,an} = -1 mod B+1 */
144 MPN_ZERO (ap + 2 * n, rn - 2 * n + 1);
145 }
146 else
147 mpn_random2 (ap, rn + 1);
148 mpn_random2 (bp, rn + 1);
149
150 bp [rn] &= 1;
151 ap [rn] &= 1;
152
153 mpn_random2 (pp-1, rn + 3);
154 p_before = pp[-1];
155 p_after = pp[rn + 1];
156
157 itch = mpn_mulmod_bknp1_itch (rn);
158 ASSERT_ALWAYS (itch <= mpn_mulmod_bknp1_itch (MAX_N));
159 mpn_random2 (scratch - 1, itch + 2);
160 s_before = scratch[-1];
161 s_after = scratch[itch];
162
163 mpn_mulmod_bknp1 ( pp, ap, bp, n, k, scratch);
164 ref_mulmod_bnp1 (refp, ap, bp, rn);
165 if (pp[-1] != p_before || pp[rn + 1] != p_after
166 || scratch[-1] != s_before || scratch[itch] != s_after
167 || mpn_cmp (refp, pp, rn + 1) != 0)
168 {
169 printf ("ERROR in test %d, rn = %d, n = %d, k = %d\n",
170 test, (int) rn, (int) n, (int) k);
171 if (pp[-1] != p_before)
172 {
173 printf ("before pp:"); mpn_dump (pp - 1, 1);
174 printf ("keep: "); mpn_dump (&p_before, 1);
175 }
176 if (pp[rn + 1] != p_after)
177 {
178 printf ("after pp:"); mpn_dump (pp + rn + 1, 1);
179 printf ("keep: "); mpn_dump (&p_after, 1);
180 }
181 if (scratch[-1] != s_before)
182 {
183 printf ("before scratch:"); mpn_dump (scratch - 1, 1);
184 printf ("keep: "); mpn_dump (&s_before, 1);
185 }
186 if (scratch[itch] != s_after)
187 {
188 printf ("after scratch:"); mpn_dump (scratch + itch, 1);
189 printf ("keep: "); mpn_dump (&s_after, 1);
190 }
191 mpn_dump (ap, rn + 1);
192 mpn_dump (bp, rn + 1);
193 mpn_dump (pp, rn + 1);
194 mpn_dump (refp, rn + 1);
195
196 abort();
197 }
198 }
199 TMP_FREE;
200 tests_end ();
201 return 0;
202 }