number-theory/square_root_mod_tests.cpp at master · Mirraz/number-theory · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#include <assert.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <sys/time.h>
#include "square_root_mod.h"
#include "factorize.h"
#include "mul_mod.h"

uint_fast32_t my_min_nonresidue(uint_fast32_t primes[], size_t primes_count, uint_fast32_t p) {
	assert(p > 2);
	typedef MulMod<uint_fast32_t,((uint_fast32_t)1 << 31), uint_fast64_t> mul_mod_type;
	for (size_t i=0; i < primes_count; ++i) {
		uint_fast32_t nr = primes[i];
		if (mul_mod_type::pow_mod(p, nr, (p-1)/2) == p-1) return nr;
	}
	assert(0);
	return 0;
}

void test_least_nonresidue() {
	typedef uint_fast32_t num_type;
	typedef SquareRootMod<num_type, 32, uint_fast64_t> srm_type;
	typedef PrimesArray<num_type> primes_array_type;
	num_type primes[65536];
	size_t primes_count = primes_array_type::fill_primes(primes, 65536, UINT32_MAX);
	assert(primes_count == 65536);

	for (size_t idx=1; idx<primes_count; ++idx) {
		num_type p = primes[idx];
		num_type nr = srm_type::least_nonresidue(primes, primes_count, p);
		assert(nr != 0);
		num_type my_nr = my_min_nonresidue(primes, primes_count, p);
		assert(nr == my_nr);
	}
}

void test_square_root_mod_algo_01() {
	typedef uint_fast32_t num_type;
	typedef SquareRootMod<num_type, 32, uint_fast64_t> srm_type;
	typedef PrimesArray<num_type> primes_array_type;
	num_type primes[1024];
	size_t primes_count = primes_array_type::fill_primes(primes, 1024, UINT32_MAX);
	assert(primes_count == 1024);

	for (size_t idx=1; idx<primes_count; ++idx) {
		num_type p = primes[idx];
		num_type nr = srm_type::least_nonresidue(primes, primes_count, p);
		assert(nr != 0);
		num_type r_count = 0, nr_count = 0;
		for (num_type a=1; a<p; ++a) {
			num_type r = srm_type::square_root_mod_algo_01(p, nr, a);
			int legendre_symbol = srm_type::legendre_symbol(p, a);
			assert((legendre_symbol != 1) == (r == 0));
			if (r == 0) {
				++nr_count;
			} else {
				assert(srm_type::mul_mod_type::square_mod(p, r) == a);
				++r_count;
			}
		}
		assert(r_count == nr_count);
		assert(r_count == (p-1)/2);
	}
}

void test_tonelli_shanks_algo() {
	typedef uint_fast32_t num_type;
	typedef SquareRootMod<num_type, 32, uint_fast64_t> srm_type;
	typedef PrimesArray<num_type> primes_array_type;
	num_type primes[1024];
	size_t primes_count = primes_array_type::fill_primes(primes, 1024, UINT32_MAX);
	assert(primes_count == 1024);

	for (size_t idx=1; idx<primes_count; ++idx) {
		num_type p = primes[idx];
		srm_type square_root_mod(p, primes, primes_count);
		num_type r_count = 0, nr_count = 0;
		for (num_type a=1; a<p; ++a) {
			if (square_root_mod.legendre_symbol(a) != 1) {++nr_count; continue;}
			num_type r = square_root_mod.tonelli_shanks_algo(a);
			assert(srm_type::mul_mod_type::square_mod(p, r) == a);
			++r_count;
		}
		assert(r_count == nr_count);
		assert(r_count == (p-1)/2);
	}
}

void test_tonelli_shanks_algo_02_rand() {
	struct timeval tv;
	gettimeofday(&tv, NULL);
	unsigned int seed = (unsigned int)tv.tv_sec * 1000000 + tv.tv_usec;
	fprintf(stderr, "seed = %u\n", seed);
	srand(seed);

	typedef uint_fast32_t num_type;
	typedef SquareRootMod<num_type, 32, uint_fast64_t> srm_type;
	typedef PrimesArray<num_type> primes_array_type;
	num_type primes[1024];
	size_t primes_count = primes_array_type::fill_primes(primes, 1024, UINT32_MAX);
	assert(primes_count == 1024);
	PrimeChecker<num_type> prime_checker(primes_array_type(primes, primes_count));

	for (num_type p = UINT32_MAX-1024*4;; p+=2) {
		if (prime_checker.is_prime(p)) {
			srm_type square_root_mod(p, primes, primes_count);
			for (uint_fast16_t i=0; i<1024*4; ++i) {
				num_type a = rand() % p;
				if (square_root_mod.legendre_symbol(a) != 1) continue;
				num_type r = square_root_mod.tonelli_shanks_algo(a);
				assert(srm_type::mul_mod_type::square_mod(p, r) == a);
			}
		}
		if (p == UINT32_MAX) break;
	}
}

void tests_suite() {
	//test_least_nonresidue();
	//test_square_root_mod_algo_01();
	test_tonelli_shanks_algo();
	test_tonelli_shanks_algo_02_rand();
}

int main() {
	tests_suite();
	return 0;
}