1
2 /* Test for mulmod_bnm1 function.
3
4 Contributed to the GNU project by Marco Bodrato.
5
6 Copyright 2009, 2020 Free Software Foundation, Inc.
7
8 This file is part of the GNU MP Library test suite.
9
10 The GNU MP Library test suite is free software; you can redistribute it
11 and/or modify it under the terms of the GNU General Public License as
12 published by the Free Software Foundation; either version 3 of the License,
13 or (at your option) any later version.
14
15 The GNU MP Library test suite is distributed in the hope that it will be
16 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
18 Public License for more details.
19
20 You should have received a copy of the GNU General Public License along with
21 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
22
23
24 #include <stdlib.h>
25 #include <stdio.h>
26
27 #include "gmp-impl.h"
28 #include "tests.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
52 static void
53 ref_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
78 int
79 main (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 TESTS_REPS (count, argv, argc);
89
90 tests_start ();
91 rands = RANDS;
92
93 ASSERT_ALWAYS (mpn_mulmod_bnm1_next_size (MAX_N) == MAX_N);
94
95 ap = TMP_ALLOC_LIMBS (MAX_N);
96 bp = TMP_ALLOC_LIMBS (MAX_N);
97 refp = TMP_ALLOC_LIMBS (MAX_N * 4);
98 pp = 1+TMP_ALLOC_LIMBS (MAX_N + 2);
99 scratch
100 = 1+TMP_ALLOC_LIMBS (mpn_mulmod_bnm1_itch (MAX_N, MAX_N, MAX_N) + 2);
101
102 for (test = 0; test < count; test++)
103 {
104 unsigned size_min;
105 unsigned size_range;
106 mp_size_t an,bn,rn,n;
107 mp_size_t itch;
108 mp_limb_t p_before, p_after, s_before, s_after;
109
110 for (size_min = 1; (1L << size_min) < MIN_N; size_min++)
111 ;
112
113 /* We generate an in the MIN_N <= n <= (1 << size_range). */
114 size_range = size_min
115 + gmp_urandomm_ui (rands, SIZE_LOG + 1 - size_min);
116
117 n = MIN_N
118 + gmp_urandomm_ui (rands, (1L << size_range) + 1 - MIN_N);
119 n = mpn_mulmod_bnm1_next_size (n);
120
121 if ( (test & 1) || n == 1) {
122 /* Half of the tests are done with the main scenario in mind:
123 both an and bn >= rn/2 */
124 an = ((n+1) >> 1) + gmp_urandomm_ui (rands, (n+1) >> 1);
125 bn = ((n+1) >> 1) + gmp_urandomm_ui (rands, (n+1) >> 1);
126 } else {
127 /* Second half of the tests are done using mulmod to compute a
128 full product with n/2 < an+bn <= n. */
129 an = 1 + gmp_urandomm_ui (rands, n - 1);
130 if (an >= n/2)
131 bn = 1 + gmp_urandomm_ui (rands, n - an);
132 else
133 bn = n/2 + 1 - an + gmp_urandomm_ui (rands, (n+1)/2);
134 }
135
136 /* Make sure an >= bn */
137 if (an < bn)
138 MP_SIZE_T_SWAP (an, bn);
139
140 mpn_random2 (ap, an);
141 mpn_random2 (bp, bn);
142
143 /* Sometime trigger the borderline conditions
144 A = -1,0,+1 or B = -1,0,+1 or A*B == -1,0,1 Mod(B^{n/2}+1).
145 This only makes sense if there is at least a split, i.e. n is even. */
146 if ((test & 0x1f) == 1 && (n & 1) == 0) {
147 mp_size_t x;
148 MPN_COPY (ap, ap + (n >> 1), an - (n >> 1));
149 MPN_ZERO (ap + an - (n >> 1) , n - an);
150 MPN_COPY (bp, bp + (n >> 1), bn - (n >> 1));
151 MPN_ZERO (bp + bn - (n >> 1) , n - bn);
152 x = 0;
153 /* x = (n == an) ? 0 : gmp_urandomm_ui (rands, n - an); */
154 ap[x] += gmp_urandomm_ui (rands, 3) - 1;
155 /* x = (n >> 1) - x % (n >> 1); */
156 bp[x] += gmp_urandomm_ui (rands, 3) - 1;
157 /* We don't propagate carry, this means that the desired condition
158 is not triggered all the times. A few times are enough anyway. */
159 }
160 rn = MIN(n, an + bn);
161 mpn_random2 (pp-1, rn + 2);
162 p_before = pp[-1];
163 p_after = pp[rn];
164
165 itch = mpn_mulmod_bnm1_itch (n, an, bn);
166 ASSERT_ALWAYS (itch <= mpn_mulmod_bnm1_itch (MAX_N, MAX_N, MAX_N));
167 mpn_random2 (scratch-1, itch+2);
168 s_before = scratch[-1];
169 s_after = scratch[itch];
170
171 mpn_mulmod_bnm1 ( pp, n, ap, an, bp, bn, scratch);
172 ref_mulmod_bnm1 (refp, n, ap, an, bp, bn);
173 if (pp[-1] != p_before || pp[rn] != p_after
174 || scratch[-1] != s_before || scratch[itch] != s_after
175 || mpn_cmp (refp, pp, rn) != 0)
176 {
177 printf ("ERROR in test %d, an = %d, bn = %d, n = %d\n",
178 test, (int) an, (int) bn, (int) n);
179 if (pp[-1] != p_before)
180 {
181 printf ("before pp:"); mpn_dump (pp -1, 1);
182 printf ("keep: "); mpn_dump (&p_before, 1);
183 }
184 if (pp[rn] != p_after)
185 {
186 printf ("after pp:"); mpn_dump (pp + rn, 1);
187 printf ("keep: "); mpn_dump (&p_after, 1);
188 }
189 if (scratch[-1] != s_before)
190 {
191 printf ("before scratch:"); mpn_dump (scratch-1, 1);
192 printf ("keep: "); mpn_dump (&s_before, 1);
193 }
194 if (scratch[itch] != s_after)
195 {
196 printf ("after scratch:"); mpn_dump (scratch + itch, 1);
197 printf ("keep: "); mpn_dump (&s_after, 1);
198 }
199 mpn_dump (ap, an);
200 mpn_dump (bp, bn);
201 mpn_dump (pp, rn);
202 mpn_dump (refp, rn);
203
204 abort();
205 }
206 }
207 TMP_FREE;
208 tests_end ();
209 return 0;
210 }