1 // Copyright (c) 2015 John Abbott
2
3 // This file is part of the source of CoCoALib, the CoCoA Library.
4
5 // CoCoALib is free software: you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation, either version 3 of the License, or
8 // (at your option) any later version.
9
10 // CoCoALib is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
14
15 // You should have received a copy of the GNU General Public License
16 // along with CoCoALib. If not, see <http://www.gnu.org/licenses/>.
17
18
19 #include "CoCoA/BuildInfo.H"
20 #include "CoCoA/GlobalManager.H"
21 #include "CoCoA/error.H"
22 #include "CoCoA/NumTheory-modular.H"
23 #include "CoCoA/NumTheory-factor.H"
24 #include "CoCoA/NumTheory-prime.H"
25
26
27 #include <algorithm>
28 using std::min;
29 #include <iostream>
30 using std::cerr;
31 using std::endl;
32 #include <vector>
33 using std::vector;
34
35
36 namespace CoCoA
37 {
38
39 // For test_radical:
40 // ad hoc impl -- no arg checks!
pwr(long base,int exp)41 long pwr(long base, int exp)
42 {
43 long ans = 1;
44 for (int i=0; i < exp; ++i)
45 ans *= base;
46 return ans;
47 }
48
49
50 //----------------------------------------------------------------------
51 // This test checks that radical gives the expected answer over a
52 // range of sall numbers with prime factors 2,3,5,7.
53 //----------------------------------------------------------------------
test_radical()54 void test_radical()
55 {
56 const long p1 = 2;
57 const long p2 = 3;
58 const long p3 = 5;
59 const long p4 = 7;
60 for (int e1 = 0; e1 < 4; ++e1)
61 for (int e2 = 0; e2 < 4; ++e2)
62 for (int e3 = 0; e3 < 4; ++e3)
63 for (int e4 = 0; e4 < 4; ++e4)
64 {
65 const long N = pwr(p1,e1) * pwr(p2,e2) * pwr(p3,e3) * pwr(p4,e4);
66 const long radN = pwr(p1,min(e1,1)) * pwr(p2,min(e2,1)) * pwr(p3,min(e3,1)) * pwr(p4,min(e4,1));
67 CoCoA_ASSERT_ALWAYS(radical(N) == radN);
68 CoCoA_ASSERT_ALWAYS(radical(-N) == radN);
69 CoCoA_ASSERT_ALWAYS(radical(radical(N)) == radN);
70 }
71 }
72
73
74 //----------------------------------------------------------------------
75 // This test checks that IsSqFree gives the expected answer over the
76 // range 1 to 100.
77 //----------------------------------------------------------------------
test_IsSqFree()78 void test_IsSqFree()
79 {
80 for (int i=1; i < 100; ++i)
81 {
82 const factorization<long> FacInfo = factor(i);
83 const vector<long>& mult = FacInfo.myMultiplicities();
84 int MaxExp = 0;
85 for (int j=0; j < len(mult); ++j)
86 if (mult[j] > MaxExp) MaxExp = mult[j];
87 const bool ans1 = IsSqFree(i);
88 const bool3 ans2 = IsSqFree(BigInt(i));
89 const bool ans3 = IsSqFree(-i);
90 const bool3 ans4 = IsSqFree(BigInt(-i));
91 CoCoA_ASSERT_ALWAYS(ans1 == (MaxExp <= 1));
92 CoCoA_ASSERT_ALWAYS(ans1 == IsTrue3(ans2) && ans1 == !IsFalse3(ans2));
93 CoCoA_ASSERT_ALWAYS(ans1 == ans3);
94 CoCoA_ASSERT_ALWAYS(ans1 == IsTrue3(ans4) && ans1 == !IsFalse3(ans4));
95 }
96
97 }
98
99
100 //-------------------------------------------------------
101 // Test the various "prime sequence" implementations
102 //-------------------------------------------------------
test_PrimeSeq()103 void test_PrimeSeq()
104 {
105 // FastFinitePrimeSeq is really just table-lookup.
106 // Check the sequence is coherent: starts with 2, and each
107 // successive entry is NextPrime of previous entry.
108 FastFinitePrimeSeq seq3;
109 long p3 = *seq3;
110 CoCoA_ASSERT_ALWAYS(p3 == 2);
111 while (!IsEnded(++seq3))
112 {
113 CoCoA_ASSERT_ALWAYS(*seq3 == NextPrime(p3));
114 p3 = *seq3;
115 }
116
117
118 // FastMostlyPrimeSeq is same as FastFinitePrimeSeq up to index 82024
119 // then it generates some composite numbers (but with no factor < 29)
120 FastMostlyPrimeSeq seq0;
121 for (int i=0; i < 100000; ++i)
122 {
123 const long n = *seq0;
124 if (n < 25) { CoCoA_ASSERT_ALWAYS(IsPrime(n)); ++seq0; continue; }
125 CoCoA_ASSERT_ALWAYS(n%2 != 0);
126 CoCoA_ASSERT_ALWAYS(n%3 != 0);
127 CoCoA_ASSERT_ALWAYS(n%5 != 0);
128 CoCoA_ASSERT_ALWAYS(n%7 != 0);
129 CoCoA_ASSERT_ALWAYS(n%11 != 0);
130 CoCoA_ASSERT_ALWAYS(n%13 != 0);
131 CoCoA_ASSERT_ALWAYS(n%17 != 0);
132 CoCoA_ASSERT_ALWAYS(n%19 != 0);
133 CoCoA_ASSERT_ALWAYS(n%23 != 0);
134 ++seq0;
135 }
136
137
138 // PrimeSeq generates the sequence of primes starting from 2.
139 // Check the sequence is coherent: starts with 2, and each
140 // successive entry is NextPrime of previous entry.
141 // Internally it "changes gear" at around 82024; check for
142 // correct behaviour around this point.
143 PrimeSeq seq1;
144 CoCoA_ASSERT_ALWAYS(*seq1 == 2);
145 for (int i=0; i < 82000; ++i)
146 ++seq1;
147 long p1 = *seq1;
148 for (int i=0; i < 100; ++i)
149 {
150 ++seq1;
151 CoCoA_ASSERT_ALWAYS(*seq1 == NextPrime(p1));
152 p1 = *seq1;
153 }
154
155 // PrimeSeqForCRT generates "large" small primes.
156 // It uses a table for the first 50000ish entries then
157 // generates on the fly. We check for a smooth transition.
158 PrimeSeqForCRT seq2;
159 for (int i=0; i < 50000; ++i)
160 ++seq2;
161 long p2 = *seq2;
162 for (int i=0; i < 100; ++i)
163 {
164 ++seq2;
165 CoCoA_ASSERT_ALWAYS(*seq2 == NextPrime(p2));
166 p2 = *seq2;
167 }
168 }
169
170
171 //----------------------------------------------------------------------
172 // Next test checks the function KroneckerSymbol
173 // (for testing whether a value is a quadratic residue).
174 // First 2 auxiliary fns.
175 //----------------------------------------------------------------------
176
177 // KroneckerSymbol: Check that sum of all values is 0 if p is odd
Kronecker_CheckPrime1(int p)178 void Kronecker_CheckPrime1(int p)
179 {
180 if (p == 2) return;
181 long sum = 0;
182 for (long r=0; r < p; ++r)
183 sum += KroneckerSymbol(r, p);
184 CoCoA_ASSERT_ALWAYS(sum == 0);
185
186 for (long r=0; r < p; ++r)
187 sum += KroneckerSymbol(r, BigInt(p));
188 CoCoA_ASSERT_ALWAYS(sum == 0);
189
190 for (long r=0; r < p; ++r)
191 sum += KroneckerSymbol(BigInt(r), p);
192 CoCoA_ASSERT_ALWAYS(sum == 0);
193
194 for (long r=0; r < p; ++r)
195 sum += KroneckerSymbol(BigInt(r), BigInt(p));
196 CoCoA_ASSERT_ALWAYS(sum == 0);
197 }
198
199
200 // KroneckerSymbol: Check that 0 gives 0, and that all squares give 1
Kronecker_CheckPrime2(int p)201 void Kronecker_CheckPrime2(int p)
202 {
203 CoCoA_ASSERT_ALWAYS(KroneckerSymbol(0, p) == 0);
204 for (long r=1; r < p; ++r)
205 CoCoA_ASSERT_ALWAYS(KroneckerSymbol(r*r, p) == 1);
206
207 CoCoA_ASSERT_ALWAYS(KroneckerSymbol(0, BigInt(p)) == 0);
208 for (long r=1; r < p; ++r)
209 CoCoA_ASSERT_ALWAYS(KroneckerSymbol(r*r, BigInt(p)) == 1);
210
211 CoCoA_ASSERT_ALWAYS(KroneckerSymbol(BigInt(0), p) == 0);
212 for (long r=1; r < p; ++r)
213 CoCoA_ASSERT_ALWAYS(KroneckerSymbol(BigInt(r*r), p) == 1);
214
215 CoCoA_ASSERT_ALWAYS(KroneckerSymbol(BigInt(0), BigInt(p)) == 0);
216 for (long r=1; r < p; ++r)
217 CoCoA_ASSERT_ALWAYS(KroneckerSymbol(BigInt(r*r), BigInt(p)) == 1);
218 }
219
220
test_kronecker()221 void test_kronecker()
222 {
223 Kronecker_CheckPrime1(2);
224 Kronecker_CheckPrime2(2);
225
226 Kronecker_CheckPrime1(3);
227 Kronecker_CheckPrime2(3);
228
229 Kronecker_CheckPrime1(5);
230 Kronecker_CheckPrime2(5);
231
232 Kronecker_CheckPrime1(101);
233 Kronecker_CheckPrime2(101);
234
235 Kronecker_CheckPrime1(32003);
236 Kronecker_CheckPrime2(32003);
237 }
238
239
program()240 void program()
241 {
242 GlobalManager CoCoAFoundations;
243
244 test_radical();
245 test_IsSqFree();
246 test_PrimeSeq();
247 test_kronecker();
248 }
249
250 } // end of namespace CoCoA
251
252
253 // Use main() to handle any uncaught exceptions and warn the user about them.
main()254 int main()
255 {
256 try
257 {
258 CoCoA::program();
259 return 0;
260 }
261 catch (const CoCoA::ErrorInfo& err)
262 {
263 cerr << "***ERROR*** UNCAUGHT CoCoA Error";
264 ANNOUNCE(cerr, err);
265 }
266 catch (const std::exception& exc)
267 {
268 cerr << "***ERROR*** UNCAUGHT std::exception: " << exc.what() << endl;
269 }
270 catch(...)
271 {
272 cerr << "***ERROR*** UNCAUGHT UNKNOWN EXCEPTION" << endl;
273 }
274
275 CoCoA::BuildInfo::PrintAll(cerr);
276 return 1;
277 }
278