summaryrefslogtreecommitdiff
path: root/euler381.c
diff options
context:
space:
mode:
Diffstat (limited to 'euler381.c')
-rw-r--r--euler381.c73
1 files changed, 73 insertions, 0 deletions
diff --git a/euler381.c b/euler381.c
new file mode 100644
index 0000000..204774b
--- /dev/null
+++ b/euler381.c
@@ -0,0 +1,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);
+}