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);
+}