summaryrefslogtreecommitdiff
path: root/euler381.c
blob: 204774bf16a1518b5f3929ffe901fdf1133bedfd (plain)
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
#include <stdio.h>

/* Wilson's theorem: (p-1)! = -1 (mod p)
 * (p-2)! = -1 * (p-1)^(-1) = 1
 * (p-3)! = -1 * (p-1)^(-1) * (p-2)^(-1) = (-2)^(-1)
 * (p-4)! = (-2)^(-1) * (-3)^(-1) = 6^(-1)
 * (p-5)! = (-24)^(-1)
 * S(p) = (-2)^(-1) + 6^(-1) + (-24)^(-1)
 */

/* a = q[0]*b + r[0]
 * b = q[1]*r[0] + r[1]
 * r[0] = q[2]*r[1] + r[2]
 * ...
 * r[n-2] = q[n]*r[n-1] + r[n] (r[n]=gcd(a,b)=1)
 */
int inv(int a, int p)
{
	/* if inv(p%a, a) = x, then
	 * x*(p%a) = qa + 1 (q=x*(p%a)/a)
	 * x*(p-(p/a)*a) = qa + 1
	 * (-x(p/a)-q)a = 1 (mod p)
	 * => inv(a, p) = -x(p/a)-q
	 */
	if (a == 1 || a == p-1)
		return a;
	int x = inv(p%a, a);
	long long q = (long long)x * (p%a) / a;
	long long res = (-(long long)x*(p/a)-q)%p;
	return (p+res)%p;
}

int S(int p)
{
	return (inv(p-2, p) + inv(6, p) + inv(p-24%p, p)) % p;
}

#define MAXNUM 100000000
int nprimes;
int primes[MAXNUM];

void getprimes()
{
	primes[0] = 2;
	nprimes = 1;
	for (int i = 3; i <= MAXNUM; i+=2) {
		int isprime = 1;
		for (int j = 0; j < nprimes; j++) {
			if (primes[j] * primes[j] > i)
				break;
			if (i % primes[j] == 0) {
				isprime = 0;
				break;
			}
		}
		if (isprime) {
			primes[nprimes] = i;
			nprimes ++;
		}
	}
}

int main()
{
	long long sum = 0;
	getprimes();

	/* primes[2] = 5 */
	for (int i = 2; i < nprimes; i++) {
		sum += S(primes[i]);
	}
	printf("%lld\n", sum);
}