pi.go

package main

import (
	"fmt"
	"math/big"
	"runtime"
	"sync"
	"time"
)

// Calculate Pi using the Bailey-Borwein-Plouffe formula
// This formula allows direct computation of any hexadecimal digit of Pi
func calculatePiBBP(digits int) string {
	// Set precision high enough for the calculation
	precision := uint(digits*4 + 100)

	// Initialize result
	pi := new(big.Float).SetPrec(precision)
	pi.SetInt64(0)

	// Constants
	one := new(big.Float).SetPrec(precision)
	one.SetInt64(1)

	sixteen := new(big.Float).SetPrec(precision)
	sixteen.SetInt64(16)

	// Temporary variables
	term := new(big.Float).SetPrec(precision)
	sum := new(big.Float).SetPrec(precision)
	power := new(big.Float).SetPrec(precision)

	// Number of terms needed for the desired precision
	terms := digits * 2

	// Calculate Pi using the BBP formula
	// Pi = sum_{k=0}^infinity (1/16^k) * (4/(8k+1) - 2/(8k+4) - 1/(8k+5) - 1/(8k+6))
	for k := 0; k < terms; k++ {
		// Calculate 1/16^k
		power.SetInt64(1)
		for i := 0; i < k; i++ {
			power.Quo(power, sixteen)
		}

		// Calculate the sum for this k
		sum.SetInt64(0)

		// 4/(8k+1)
		term.SetInt64(4)
		term.Quo(term, new(big.Float).SetInt64(int64(8*k+1)))
		sum.Add(sum, term)

		// -2/(8k+4)
		term.SetInt64(2)
		term.Quo(term, new(big.Float).SetInt64(int64(8*k+4)))
		sum.Sub(sum, term)

		// -1/(8k+5)
		term.SetInt64(1)
		term.Quo(term, new(big.Float).SetInt64(int64(8*k+5)))
		sum.Sub(sum, term)

		// -1/(8k+6)
		term.SetInt64(1)
		term.Quo(term, new(big.Float).SetInt64(int64(8*k+6)))
		sum.Sub(sum, term)

		// Multiply by 1/16^k and add to pi
		sum.Mul(sum, power)
		pi.Add(pi, sum)
	}

	// Convert to string with the desired precision
	return pi.Text('f', digits)
}

// Parallel version of the BBP formula
func calculatePiBBPParallel(digits int) string {
	// Set precision high enough for the calculation
	precision := uint(digits*4 + 100)

	// Number of terms needed for the desired precision
	terms := digits * 2

	// Determine number of goroutines based on CPU cores
	numCPU := runtime.NumCPU()
	numWorkers := numCPU * 2 // Use 2x CPU cores for better utilization

	// Ensure we don't create more workers than terms
	if numWorkers > terms {
		numWorkers = terms
	}

	// Calculate terms per worker
	termsPerWorker := terms / numWorkers

	// Create a channel for results
	results := make(chan *big.Float, numWorkers)

	// Use a wait group to synchronize goroutines
	wg := sync.WaitGroup{}

	// Launch worker goroutines
	for w := 0; w < numWorkers; w++ {
		wg.Add(1)

		// Calculate the range for this worker
		startK := w * termsPerWorker
		endK := (w + 1) * termsPerWorker

		// Last worker takes any remaining terms
		if w == numWorkers-1 {
			endK = terms
		}

		// Launch the worker goroutine
		go func(start, end int) {
			defer wg.Done()

			// Initialize worker's partial sum
			partialSum := new(big.Float).SetPrec(precision)
			partialSum.SetInt64(0)

			// Constants
			sixteen := new(big.Float).SetPrec(precision)
			sixteen.SetInt64(16)

			// Temporary variables
			term := new(big.Float).SetPrec(precision)
			sum := new(big.Float).SetPrec(precision)
			power := new(big.Float).SetPrec(precision)

			// Calculate terms in this worker's range
			for k := start; k < end; k++ {
				// Calculate 1/16^k
				power.SetInt64(1)
				for i := 0; i < k; i++ {
					power.Quo(power, sixteen)
				}

				// Calculate the sum for this k
				sum.SetInt64(0)

				// 4/(8k+1)
				term.SetInt64(4)
				term.Quo(term, new(big.Float).SetInt64(int64(8*k+1)))
				sum.Add(sum, term)

				// -2/(8k+4)
				term.SetInt64(2)
				term.Quo(term, new(big.Float).SetInt64(int64(8*k+4)))
				sum.Sub(sum, term)

				// -1/(8k+5)
				term.SetInt64(1)
				term.Quo(term, new(big.Float).SetInt64(int64(8*k+5)))
				sum.Sub(sum, term)

				// -1/(8k+6)
				term.SetInt64(1)
				term.Quo(term, new(big.Float).SetInt64(int64(8*k+6)))
				sum.Sub(sum, term)

				// Multiply by 1/16^k and add to partial sum
				sum.Mul(sum, power)
				partialSum.Add(partialSum, sum)
			}

			// Send the partial sum to the results channel
			results <- partialSum
		}(startK, endK)
	}

	// Wait for all workers to complete
	go func() {
		wg.Wait()
		close(results)
	}()

	// Combine the results
	pi := new(big.Float).SetPrec(precision)
	pi.SetInt64(0)

	for partialSum := range results {
		pi.Add(pi, partialSum)
	}

	// Convert to string with the desired precision
	return pi.Text('f', digits)
}

// Calculate Pi using the Machin formula
// Pi = 16*arctan(1/5) - 4*arctan(1/239)
func calculatePiMachin(digits int) string {
	// Set precision high enough for the calculation
	precision := uint(digits*4 + 100)

	// Initialize result with high precision
	pi := new(big.Float).SetPrec(precision)
	pi.SetInt64(0)

	// Calculate 16*arctan(1/5)
	arctan1 := arctan(new(big.Float).SetInt64(5), precision, digits)
	arctan1.Mul(arctan1, new(big.Float).SetInt64(16))

	// Calculate 4*arctan(1/239)
	arctan2 := arctan(new(big.Float).SetInt64(239), precision, digits)
	arctan2.Mul(arctan2, new(big.Float).SetInt64(4))

	// Pi = 16*arctan(1/5) - 4*arctan(1/239)
	pi.Sub(arctan1, arctan2)

	// Convert to string with the desired precision
	return pi.Text('f', digits)
}

// Parallel version of the Machin formula
func calculatePiMachinParallel(digits int) string {
	// Set precision high enough for the calculation
	precision := uint(digits*4 + 100)

	// Initialize result with high precision
	pi := new(big.Float).SetPrec(precision)
	pi.SetInt64(0)

	// Use a wait group to synchronize goroutines
	wg := sync.WaitGroup{}

	// Create channels for results
	arctan1Chan := make(chan *big.Float, 1)
	arctan2Chan := make(chan *big.Float, 1)

	// Calculate 16*arctan(1/5) in a goroutine
	wg.Add(1)
	go func() {
		defer wg.Done()
		result := arctanParallel(new(big.Float).SetInt64(5), precision, digits)
		result.Mul(result, new(big.Float).SetInt64(16))
		arctan1Chan <- result
	}()

	// Calculate 4*arctan(1/239) in a goroutine
	wg.Add(1)
	go func() {
		defer wg.Done()
		result := arctanParallel(new(big.Float).SetInt64(239), precision, digits)
		result.Mul(result, new(big.Float).SetInt64(4))
		arctan2Chan <- result
	}()

	// Wait for both calculations to complete
	go func() {
		wg.Wait()
		close(arctan1Chan)
		close(arctan2Chan)
	}()

	// Get the results
	arctan1 := <-arctan1Chan
	arctan2 := <-arctan2Chan

	// Pi = 16*arctan(1/5) - 4*arctan(1/239)
	pi.Sub(arctan1, arctan2)

	// Convert to string with the desired precision
	return pi.Text('f', digits)
}

// Calculate arctan(1/x) to the specified precision
func arctan(x *big.Float, precision uint, digits int) *big.Float {
	// Initialize result with high precision
	result := new(big.Float).SetPrec(precision)
	result.SetInt64(0)

	// Temporary variables
	term := new(big.Float).SetPrec(precision)
	term.SetInt64(1)
	term.Quo(term, x)

	xSquared := new(big.Float).SetPrec(precision)
	xSquared.Mul(x, x)

	// Number of terms needed for the desired precision
	terms := digits * 2

	// Calculate arctan(1/x) using the series: arctan(1/x) = 1/x - 1/(3*x^3) + 1/(5*x^5) - ...
	for k := 0; k < terms; k++ {
		// Add or subtract the term based on whether k is even or odd
		if k%2 == 0 {
			result.Add(result, term)
		} else {
			result.Sub(result, term)
		}

		// Calculate the next term: term = term / (x^2 * (2k+3)/(2k+1))
		term.Quo(term, xSquared)
		term.Mul(term, new(big.Float).SetInt64(int64(2*k + 1)))
		term.Quo(term, new(big.Float).SetInt64(int64(2*k + 3)))
	}

	return result
}

// Parallel version of arctan calculation
func arctanParallel(x *big.Float, precision uint, digits int) *big.Float {
	// Initialize result with high precision
	result := new(big.Float).SetPrec(precision)
	result.SetInt64(0)

	// Number of terms needed for the desired precision
	terms := digits * 2

	// Determine number of goroutines based on CPU cores
	numCPU := runtime.NumCPU()
	numWorkers := numCPU * 2 // Use 2x CPU cores for better utilization

	// Ensure we don't create more workers than terms
	if numWorkers > terms {
		numWorkers = terms
	}

	// Calculate terms per worker
	termsPerWorker := terms / numWorkers

	// Create a channel for results
	results := make(chan *big.Float, numWorkers)

	// Use a wait group to synchronize goroutines
	wg := sync.WaitGroup{}

	// Launch worker goroutines
	for w := 0; w < numWorkers; w++ {
		wg.Add(1)

		// Calculate the range for this worker
		startK := w * termsPerWorker
		endK := (w + 1) * termsPerWorker

		// Last worker takes any remaining terms
		if w == numWorkers-1 {
			endK = terms
		}

		// Launch the worker goroutine
		go func(start, end int) {
			defer wg.Done()

			// Initialize worker's partial sum
			partialSum := new(big.Float).SetPrec(precision)
			partialSum.SetInt64(0)

			// Calculate x^2 for the series
			xSquared := new(big.Float).SetPrec(precision)
			xSquared.Mul(x, x)

			// For each term in this worker's range
			for k := start; k < end; k++ {
				// Calculate the term
				term := new(big.Float).SetPrec(precision)
				term.SetInt64(1)
				term.Quo(term, x)

				// Apply x^(2k) to the denominator
				for i := 0; i < k; i++ {
					term.Quo(term, xSquared)
				}

				// Apply the factorial ratio (2k+1)/(2k+3)/(2k+5)/...
				for i := 0; i < k; i++ {
					term.Mul(term, new(big.Float).SetInt64(int64(2*i + 1)))
					term.Quo(term, new(big.Float).SetInt64(int64(2*i + 3)))
				}

				// Apply sign based on whether k is even or odd
				if k%2 == 1 {
					term.Neg(term)
				}

				// Add to partial sum
				partialSum.Add(partialSum, term)
			}

			// Send the partial sum to the results channel
			results <- partialSum
		}(startK, endK)
	}

	// Wait for all workers to complete
	go func() {
		wg.Wait()
		close(results)
	}()

	// Combine the results
	for partialSum := range results {
		result.Add(result, partialSum)
	}

	return result
}

// Known value of Pi to 1000 decimal places for verification
const knownPi = "3." +
	"1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679" +
	"8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196" +
	"4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273" +
	"7245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094" +
	"3305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912" +
	"9833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132" +
	"0005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235" +
	"4201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859" +
	"5024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303" +
	"5982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"

func main() {
	digits := 500

	// Get number of CPU cores for information
	numCPU := runtime.NumCPU()
	fmt.Printf("Using %d CPU cores\n", numCPU)

	// Calculate Pi using the BBP formula (sequential)
	start := time.Now()
	piBBP := calculatePiBBP(digits)
	durationBBP := time.Since(start)

	// Calculate Pi using the BBP formula (parallel)
	start = time.Now()
	piBBPParallel := calculatePiBBPParallel(digits)
	durationBBPParallel := time.Since(start)

	// Calculate Pi using the Machin formula (sequential)
	start = time.Now()
	piMachin := calculatePiMachin(digits)
	durationMachin := time.Since(start)

	// Calculate Pi using the Machin formula (parallel)
	start = time.Now()
	piMachinParallel := calculatePiMachinParallel(digits)
	durationMachinParallel := time.Since(start)

	// Choose the result to display
	var piStr string
	var algorithm string
	var duration time.Duration

	// Use the parallel BBP formula result as it's generally fastest and accurate
	piStr = piBBPParallel
	algorithm = "Parallel Bailey-Borwein-Plouffe"
	duration = durationBBPParallel

	fmt.Printf("Pi with %d decimal places (using %s algorithm):\n", digits, algorithm)
	fmt.Println(piStr)
	fmt.Printf("\nCalculation took %v\n", duration)

	// Verify against known value
	fmt.Println("\nVerification of first 100 digits:")
	fmt.Println("Known Pi:     ", knownPi[:102])
	fmt.Println("Computed Pi:  ", piStr[:102])

	if piStr[:102] == knownPi[:102] {
		fmt.Println("✓ First 100 digits match!")
	} else {
		fmt.Println("✗ Digits do not match!")

		// Find where the digits start to differ
		for i := 0; i < 102 && i < len(piStr) && i < len(knownPi); i++ {
			if piStr[i] != knownPi[i] {
				fmt.Printf("Digits differ at position %d: expected '%c', got '%c'\n",
					i, knownPi[i], piStr[i])
				break
			}
		}
	}

	// Show performance comparison
	fmt.Println("\nPerformance Comparison:")
	fmt.Printf("1. Sequential BBP:    %v\n", durationBBP)
	fmt.Printf("2. Parallel BBP:      %v (%.1fx speedup)\n",
		durationBBPParallel, float64(durationBBP.Nanoseconds())/float64(durationBBPParallel.Nanoseconds()))
	fmt.Printf("3. Sequential Machin: %v\n", durationMachin)
	fmt.Printf("4. Parallel Machin:   %v (%.1fx speedup)\n",
		durationMachinParallel, float64(durationMachin.Nanoseconds())/float64(durationMachinParallel.Nanoseconds()))

	// Verify that all methods produce the same result
	fmt.Println("\nVerifying that all methods produce the same result:")

	if piBBP[:102] == piBBPParallel[:102] {
		fmt.Println("✓ Sequential and parallel BBP match!")
	} else {
		fmt.Println("✗ Sequential and parallel BBP differ!")
	}

	if piMachin[:102] == piMachinParallel[:102] {
		fmt.Println("✓ Sequential and parallel Machin match!")
	} else {
		fmt.Println("✗ Sequential and parallel Machin differ!")
	}

	if piBBPParallel[:102] == piMachinParallel[:102] {
		fmt.Println("✓ BBP and Machin algorithms match!")
	} else {
		fmt.Println("✗ BBP and Machin algorithms differ!")
	}
}