main.go raw

   1  package main
   2  
   3  import (
   4  	"fmt"
   5  	"math/big"
   6  	"os"
   7  	"time"
   8  )
   9  
  10  const (
  11  	// Precision in bits. ~3.32 bits per decimal digit,
  12  	// so 330,000 bits ≈ 100,000 decimal digits.
  13  	precision = 330_000
  14  	// Iterations. Quadratic convergence: each step doubles correct digits.
  15  	// 2^17 = 131072, comfortably exceeds 100k digits.
  16  	iterations = 17
  17  )
  18  
  19  func newFloat() *big.Float { return new(big.Float).SetPrec(precision) }
  20  
  21  func seed() (a, b, t, p *big.Float) {
  22  	a = newFloat().SetInt64(1)
  23  
  24  	// b = 1/√2
  25  	two := newFloat().SetInt64(2)
  26  	b = newFloat().Sqrt(two)
  27  	b.Quo(newFloat().SetInt64(1), b)
  28  
  29  	t = newFloat().SetFloat64(0.25)
  30  	p = newFloat().SetInt64(1)
  31  	return
  32  }
  33  
  34  func step(a, b, t, p *big.Float) (*big.Float, *big.Float, *big.Float, *big.Float) {
  35  	// aₙ₊₁ = (aₙ + bₙ) / 2
  36  	aNext := newFloat().Add(a, b)
  37  	aNext.Quo(aNext, newFloat().SetInt64(2))
  38  
  39  	// bₙ₊₁ = √(aₙ · bₙ)
  40  	ab := newFloat().Mul(a, b)
  41  	bNext := newFloat().Sqrt(ab)
  42  
  43  	// tₙ₊₁ = tₙ − pₙ(aₙ − aₙ₊₁)²
  44  	diff := newFloat().Sub(a, aNext)
  45  	diff.Mul(diff, diff)
  46  	diff.Mul(p, diff)
  47  	tNext := newFloat().Sub(t, diff)
  48  
  49  	// pₙ₊₁ = 2pₙ
  50  	pNext := newFloat().Mul(p, newFloat().SetInt64(2))
  51  
  52  	return aNext, bNext, tNext, pNext
  53  }
  54  
  55  func computePi(a, b, t *big.Float) *big.Float {
  56  	// π ≈ (a + b)² / (4t)
  57  	sum := newFloat().Add(a, b)
  58  	sum.Mul(sum, sum)
  59  	fourT := newFloat().Mul(newFloat().SetInt64(4), t)
  60  	return newFloat().Quo(sum, fourT)
  61  }
  62  
  63  func main() {
  64  	start := time.Now()
  65  
  66  	a, b, t, p := seed()
  67  	for i := range iterations {
  68  		a, b, t, p = step(a, b, t, p)
  69  		_ = p
  70  		fmt.Fprintf(os.Stderr, "iteration %d/%d\n", i+1, iterations)
  71  	}
  72  
  73  	pi := computePi(a, b, t)
  74  	elapsed := time.Since(start)
  75  
  76  	// Format to decimal string. precision / 3.32 ≈ decimal digits.
  77  	var prec float64 = precision
  78  	digits := int(prec / 3.321928)
  79  	result := pi.Text('f', digits)
  80  
  81  	fmt.Fprintf(os.Stderr, "computed %d decimal digits in %v\n", digits, elapsed)
  82  
  83  	if err := os.WriteFile("pimax.txt", []byte(result+"\n"), 0644); err != nil {
  84  		fmt.Fprintf(os.Stderr, "error writing pimax.txt: %v\n", err)
  85  		os.Exit(1)
  86  	}
  87  	fmt.Fprintf(os.Stderr, "written to pimax.txt (%d bytes)\n", len(result)+1)
  88  }
  89