package main

import (
	"fmt"
	"math/big"
	"os"
	"time"
)

const (
	// Precision in bits. ~3.32 bits per decimal digit,
	// so 330,000 bits ≈ 100,000 decimal digits.
	precision = 330_000
	// Iterations. Quadratic convergence: each step doubles correct digits.
	// 2^17 = 131072, comfortably exceeds 100k digits.
	iterations = 17
)

func newFloat() *big.Float { return new(big.Float).SetPrec(precision) }

func seed() (a, b, t, p *big.Float) {
	a = newFloat().SetInt64(1)

	// b = 1/√2
	two := newFloat().SetInt64(2)
	b = newFloat().Sqrt(two)
	b.Quo(newFloat().SetInt64(1), b)

	t = newFloat().SetFloat64(0.25)
	p = newFloat().SetInt64(1)
	return
}

func step(a, b, t, p *big.Float) (*big.Float, *big.Float, *big.Float, *big.Float) {
	// aₙ₊₁ = (aₙ + bₙ) / 2
	aNext := newFloat().Add(a, b)
	aNext.Quo(aNext, newFloat().SetInt64(2))

	// bₙ₊₁ = √(aₙ · bₙ)
	ab := newFloat().Mul(a, b)
	bNext := newFloat().Sqrt(ab)

	// tₙ₊₁ = tₙ − pₙ(aₙ − aₙ₊₁)²
	diff := newFloat().Sub(a, aNext)
	diff.Mul(diff, diff)
	diff.Mul(p, diff)
	tNext := newFloat().Sub(t, diff)

	// pₙ₊₁ = 2pₙ
	pNext := newFloat().Mul(p, newFloat().SetInt64(2))

	return aNext, bNext, tNext, pNext
}

func computePi(a, b, t *big.Float) *big.Float {
	// π ≈ (a + b)² / (4t)
	sum := newFloat().Add(a, b)
	sum.Mul(sum, sum)
	fourT := newFloat().Mul(newFloat().SetInt64(4), t)
	return newFloat().Quo(sum, fourT)
}

func main() {
	start := time.Now()

	a, b, t, p := seed()
	for i := range iterations {
		a, b, t, p = step(a, b, t, p)
		_ = p
		fmt.Fprintf(os.Stderr, "iteration %d/%d\n", i+1, iterations)
	}

	pi := computePi(a, b, t)
	elapsed := time.Since(start)

	// Format to decimal string. precision / 3.32 ≈ decimal digits.
	var prec float64 = precision
	digits := int(prec / 3.321928)
	result := pi.Text('f', digits)

	fmt.Fprintf(os.Stderr, "computed %d decimal digits in %v\n", digits, elapsed)

	if err := os.WriteFile("pimax.txt", []byte(result+"\n"), 0644); err != nil {
		fmt.Fprintf(os.Stderr, "error writing pimax.txt: %v\n", err)
		os.Exit(1)
	}
	fmt.Fprintf(os.Stderr, "written to pimax.txt (%d bytes)\n", len(result)+1)
}