From a8563d53e74e1d30d2cb1f27f5e22f9420c81c28 Mon Sep 17 00:00:00 2001
From: Leo Tenenbaum <pommicket@gmail.com>
Date: Tue, 12 May 2020 18:32:13 -0400
Subject: initial commit

---
 .gitignore |  6 +++++
 Makefile   |  3 +++
 README.md  |  9 +++++++
 main.c     | 82 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 100 insertions(+)
 create mode 100644 .gitignore
 create mode 100644 Makefile
 create mode 100644 README.md
 create mode 100644 main.c

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..50b2928
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+pi
+*.swp
+*~
+tags
+TAGS
+pi.txt
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..e7b76e8
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,3 @@
+CFLAGS=-O3 -lgmp -lm -Wall -Wextra -Wpedantic -pedantic -std=c89
+pi: main.c
+	$(CC) $(CFLAGS) main.c -o pi 
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..025534f
--- /dev/null
+++ b/README.md
@@ -0,0 +1,9 @@
+# pi
+
+Digits of pi, using the fact that
+
+```
+pi = pi + sin(pi)
+```
+
+Usage: `./pi` - generate 100 digits of pi; `./pi n` - generate n digits of pi
diff --git a/main.c b/main.c
new file mode 100644
index 0000000..9fb724e
--- /dev/null
+++ b/main.c
@@ -0,0 +1,82 @@
+#include <gmp.h>
+#include <stdio.h>
+#include <stdint.h>
+#include <math.h>
+#include <stdlib.h>
+
+typedef uint64_t u64;
+
+#define logb(a, b) ((log(a)) / (log(b)))
+
+/* approx -log10(pi^n/n!) */
+static double sin_precision_for_n(u64 n) {
+	return -(n * logb(3.2, 10) - 1/log(10) * n * (log(n) - 1));
+}
+
+/* n we should go up to for number of decimal digits */
+static u64 n_for_precision(u64 digits) {
+	double want = (double)digits;
+	double lo = 0.0;
+	double hi = want;
+	double epsilon = 1.0;
+	while (hi > lo + epsilon) {
+		double mid = (lo + hi) * 0.5;
+		double prec = sin_precision_for_n(mid);
+		if (prec < want) {
+			lo = mid;
+		} else {
+			hi = mid;
+		}
+	}
+	return (u64)(hi + 20.0); /* add 20 to be safe */
+}
+
+int main(int argc, char **argv) {
+	u64 original_digits, digits;
+	u64 sin_terms;
+	mpf_t pi, sin_pi, term, last_pi, desired_error;
+	u64 i;
+
+	original_digits = argc < 2 ? 100 : (u64)atol(argv[1]);
+	digits = original_digits < 10 ? 10 : original_digits;
+	sin_terms = n_for_precision(digits);
+
+	mpf_set_default_prec(digits < 100 ? 333 : (digits * 10) / 3);
+	mpf_init(pi);
+	mpf_init(sin_pi);
+	mpf_init(term);
+	mpf_init(desired_error);
+	mpf_init(last_pi);
+
+	mpf_set_ui(pi, 3);
+	mpf_set_ui(desired_error, 1);
+	mpf_div_ui(desired_error, desired_error, 10);
+	mpf_pow_ui(desired_error, desired_error, digits);
+
+	for (i = 0; ; ++i) {
+		unsigned j;
+		mpf_set(last_pi, pi);
+		mpf_set(sin_pi, pi);
+		mpf_set(term, pi);
+		mpf_add(sin_pi, sin_pi, term);
+		/* pi is now actually -pi^2 */
+		mpf_mul(pi, pi, pi);
+		mpf_neg(pi, pi);
+		for (j = 3; j < sin_terms; j += 2) {
+			mpf_div_ui(term, term, j * (j-1));
+			mpf_mul(term, term, pi);
+			mpf_add(sin_pi, sin_pi, term);
+		}
+		mpf_set(pi, sin_pi);
+		mpf_sub(last_pi, pi, last_pi);
+		/* last_pi is now delta */
+		mpf_abs(last_pi, last_pi);
+		gmp_fprintf(stderr, "Iteration: %9lu. Delta: %Fe\n",1+(unsigned long)i, last_pi);
+		if (mpf_cmp(last_pi, desired_error) < 0)
+			break;
+	}
+	gmp_printf("%.*Ff\n", original_digits, pi);
+
+
+	return 0;
+}
-- 
cgit v1.2.3