1 /*
2 Copyright (C) 2011 Sebastian Pancratz
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <gmp.h>
15 #include "flint.h"
16 #include "ulong_extras.h"
17
main(void)18 int main(void)
19 {
20 int i, result;
21 FLINT_TEST_INIT(state);
22
23 flint_printf("sqrtmod_primepow....");
24 fflush(stdout);
25
26
27
28 for (i = 0; i < 1000 * flint_test_multiplier(); i++) /* Test random squares mod a power of 2 */
29 {
30 mp_limb_t a, b, p, pow, pow2, pinv;
31 slong exp, num, i;
32 mp_limb_t * sqrt;
33 int btest;
34
35 p = 2;
36 exp = n_randint(state, FLINT_BITS - 1) + 1;
37 pow = n_pow(p, exp);
38 b = n_randtest(state) % pow;
39
40 pow2 = p;
41 while (FLINT_BIT_COUNT(p*pow2) <= 12)
42 pow2 *= p;
43
44 if ((b % (p*pow2)) == 0)
45 {
46 b += pow2;
47 b %= pow;
48 }
49
50 pinv = n_preinvert_limb(pow);
51 a = n_mulmod2_preinv(b, b, pow, pinv);
52
53 num = n_sqrtmod_primepow(&sqrt, a, p, exp);
54
55 btest = 0;
56 for (i = 0; i < num; i++)
57 {
58 if (a != n_mulmod2_preinv(sqrt[i], sqrt[i], pow, pinv))
59 break;
60 if (sqrt[i] == b)
61 btest = 1;
62 }
63
64 result = btest & (i == num);
65 if (!result)
66 {
67 flint_printf("FAIL:\n");
68 flint_printf("p = %wu\n", p);
69 flint_printf("exp = %wd\n", exp);
70 flint_printf("a = %wu\n", a);
71 flint_printf("b = %wu\n", b);
72 flint_printf("num = %wd\n", num);
73
74 if (!btest)
75 flint_printf("Square root not found.\n");
76 if (i != num)
77 flint_printf("%wu not a square root of %wu mod %wu\n", sqrt[i], a, pow);
78
79 abort();
80 }
81
82 flint_free(sqrt);
83 }
84
85 for (i = 0; i < 1000 * flint_test_multiplier(); i++) /* Test random squares mod other prime powers */
86 {
87 mp_limb_t a, b, p, pow, pow2, pinv;
88 slong exp, maxexp, num, i;
89 flint_bitcnt_t bits;
90 mp_limb_t * sqrt;
91 int btest;
92
93 bits = n_randint(state, 18) + 2;
94 p = n_randprime(state, bits, 0);
95 maxexp = FLINT_BITS/bits;
96 exp = n_randint(state, maxexp) + 1;
97 pow = n_pow(p, exp);
98 b = n_randtest(state) % pow;
99
100 if (bits <= FLINT_BITS/2)
101 {
102 pow2 = p;
103 while (FLINT_BIT_COUNT(p*pow2) <= 12)
104 pow2 *= p;
105
106 if ((b % (p*pow2)) == 0)
107 b += pow2;
108
109 b %= pow;
110 }
111
112 pinv = n_preinvert_limb(pow);
113 a = n_mulmod2_preinv(b, b, pow, pinv);
114
115 num = n_sqrtmod_primepow(&sqrt, a, p, exp);
116
117 btest = 0;
118 for (i = 0; i < num; i++)
119 {
120 if (a != n_mulmod2_preinv(sqrt[i], sqrt[i], pow, pinv))
121 break;
122 if (sqrt[i] == b)
123 btest = 1;
124 }
125
126 result = btest & (i == num);
127 if (!result)
128 {
129 flint_printf("FAIL:\n");
130 flint_printf("p = %wu\n", p);
131 flint_printf("exp = %wd\n", exp);
132 flint_printf("a = %wu\n", a);
133 flint_printf("b = %wu\n", b);
134 flint_printf("num = %wd\n", num);
135
136 if (!btest)
137 flint_printf("Square root not found.\n");
138 if (i != num)
139 flint_printf("%wu not a square root of %wu mod %wu\n", sqrt[i], a, pow);
140
141 abort();
142 }
143
144 flint_free(sqrt);
145 }
146
147 for (i = 0; i < 500 * flint_test_multiplier(); i++) /* Test random nonsquares */
148 {
149 mp_limb_t a, b, p, pow, pinv;
150 slong exp, maxexp;
151 flint_bitcnt_t bits;
152 mp_limb_t * sqrt;
153
154 bits = n_randint(state, 18) + 2;
155 p = n_randprime(state, bits, 0);
156 maxexp = 20/bits;
157 exp = n_randint(state, maxexp) + 1 + (p == 2);
158 pow = n_pow(p, exp);
159
160 pinv = n_preinvert_limb(pow);
161
162 a = n_randtest(state) % pow;
163 while (n_sqrtmod_primepow(&sqrt, a, p, exp))
164 {
165 if (n_mulmod2_preinv(sqrt[0], sqrt[0], pow, pinv) != a)
166 {
167 flint_printf("FAIL:\n");
168 flint_printf("%wu^2 is not %wu mod %wu\n", sqrt[0], a, pow);
169 abort();
170 }
171
172 flint_free(sqrt);
173 a = n_randtest(state) % pow;
174 }
175
176 for (b = 0; b < pow; b++)
177 {
178 if (n_mulmod2_preinv(b, b, pow, pinv) == a)
179 break;
180 }
181
182 result = (b == pow);
183 if (!result)
184 {
185 flint_printf("FAIL:\n");
186 flint_printf("p = %wu\n", p);
187 flint_printf("exp = %wd\n", exp);
188 flint_printf("a = %wu\n", a);
189 flint_printf("b = %wu\n", b);
190
191 abort();
192 }
193
194 flint_free(sqrt);
195 }
196
197 FLINT_TEST_CLEANUP(state);
198
199 flint_printf("PASS\n");
200 return 0;
201 }
202