Here is a bit of my code. The rest is here: https://github.com/Overboard-code/Pi-Pourri/blob/main/pi-pourri.py You call it with the number of desired digits, like 1_000_000 and the denominators, multipliers and operators. For your example use:
PiMachin(1_000_000,"John Machin 1706",[5,239],[4,1],[1,-1])
Here is the class. It needed some setup to work outside of my program.
from datetime import timedelta
from functools import partial
import sys,time,multiprocessing,logging,os,argparse
try:
# https://stackoverflow.com/questions/384076/how-can-i-color-python-logging-output
import Colorer
except ImportError:
pass
try:
from gmpy2 import mpz,isqrt,mpfr,atan2,sqrt,get_context,const_pi # Gumpy2 mpz large ints are ten times faster than python large int
except ImportError:
raise ImportError('This program requires gmpy2, please insatll. exiting....')
# CONSTANTS
LOG_LEVELS = ["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"]
DEFAULT_LOG_LEVEL = "INFO"
LOG2_10 = 3.32192809488736
class PiMachin:
def __init__(self,ndigits,name,denoms,mults,operators):
""" Initialization
:param int digits: digits of PI computation
:param string name: name for credit on formula
:param list denoms a lis of ints for the denomiators
:param list mults: a list of ints as multipliyers for the Machin formula
:param list operators: a list of 1 or -1 to cause addition or subtraction
"""
self.name = name # Not used. Planned to use in logging
self.ndigits = ndigits
self.denoms = denoms
self.mults = mults
self.operators = operators
self.xdigits = len(denoms)+7 # Extra digits to reduce trailing error More factors means more error
self.start_time = 0
def ArctanDenom(self,d):
cdigits = self.ndigits+self.xdigits
get_context().precision=int(cdigits * LOG2_10)
# Calculates arctan(1/d) = 1/d - 1/(3*d^3) + 1/(5*d^5) - 1/(7*d^7) + ...
logging.debug('arctan(1/%d) Started ',d )
total = mpfr(0)
arc_start_time = time.time() # Start the clock for this arctan calulation
total = mpfr(atan2(mpfr(1),mpfr(d)))
logging.debug('arctan(1/{}) Done! {:.2f} seconds.'
.format(int(d),time.time() - arc_start_time))
return total,int(1) # I used to calulate arctan by hand. Now I just use atan2() so just one iteration here
#
def compute(self):
self.start_time = time.time() # Start the clock for total time
ndigits = self.ndigits
cdigits = self.ndigits+self.xdigits
logging.info("Starting Machin-Like formula to {:,} decimal places"
.format(ndigits) )
get_context().precision=int(cdigits * LOG2_10)
logging.debug("Starting %d Pool threads to calculate arctan values.",
len(self.denoms) )
p = multiprocessing.Pool(processes=(len(self.denoms))) # get some threads for our pool
results=p.map(self.ArctanDenom, self.denoms) # one thread per arctan(1/xxxx)
p.close()
p.join() # wait for them to finish
# Now we have the arctan calculations from the pool threads in results[]
# Apply chosen Formula to the results and calculate pi using mults and signs
logging.debug ("Now multiplying and summing all arctan results")
arctanSum = pi = mpfr(0)
iters = 0
for i, result in enumerate(results):
iters += result[1] # Keep track of this thread's iterations for later
arctanSum += mpfr(mpfr(self.mults[i])*mpfr(result[0])*mpfr(self.operators[i])) # Add or subtract the product from the accumulated arctans
pi = mpfr(4) * arctanSum # change pi/4 = x to pi = 4 * x
# We calculated extra digits to compensate for roundoff error.
# Chop off the extra digits with format() using D to basically trunc() to ndigits
return "{0:.{1}Df}".format(pi,self.ndigits),iters,time.time()-self.start_time
ndigits= 1_000_000
obj = PiMachin(ndigits,"John Machin 1706",[5,239],[4,1],[1,-1])
# Calculate Pi using selected formula
pi,iters,time_to_calc = obj.compute()
print(f"First 200 of {ndigits:,}\n{pi[:200]}\nTook {time_to_calc:.2f} seconds")
