Check Machin-like formulas: Difference between revisions
m (simplified output) |
(updated output) |
||
Line 123: | Line 123: | ||
<pre> |
<pre> |
||
Testing: pi/4 = 1*atan(1/2) + 1*atan(1/3) |
Testing: pi/4 = 1*atan(1/2) + 1*atan(1/3) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 2*atan(1/3) + 1*atan(1/7) |
Testing: pi/4 = 2*atan(1/3) + 1*atan(1/7) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 4*atan(1/5) + -1*atan(1/239) |
Testing: pi/4 = 4*atan(1/5) + -1*atan(1/239) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 5*atan(1/7) + 2*atan(3/79) |
Testing: pi/4 = 5*atan(1/7) + 2*atan(3/79) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 5*atan(29/278) + 7*atan(3/79) |
Testing: pi/4 = 5*atan(29/278) + 7*atan(3/79) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 1*atan(1/2) + 1*atan(1/5) + 1*atan(1/8) |
Testing: pi/4 = 1*atan(1/2) + 1*atan(1/5) + 1*atan(1/8) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 4*atan(1/5) + -1*atan(1/70) + 1*atan(1/99) |
Testing: pi/4 = 4*atan(1/5) + -1*atan(1/70) + 1*atan(1/99) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 5*atan(1/7) + 4*atan(1/53) + 2*atan(1/4443) |
Testing: pi/4 = 5*atan(1/7) + 4*atan(1/53) + 2*atan(1/4443) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 6*atan(1/8) + 2*atan(1/57) + 1*atan(1/239) |
Testing: pi/4 = 6*atan(1/8) + 2*atan(1/57) + 1*atan(1/239) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 8*atan(1/10) + -1*atan(1/239) + -4*atan(1/515) |
Testing: pi/4 = 8*atan(1/10) + -1*atan(1/239) + -4*atan(1/515) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 12*atan(1/18) + 8*atan(1/57) + -5*atan(1/239) |
Testing: pi/4 = 12*atan(1/18) + 8*atan(1/57) + -5*atan(1/239) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 16*atan(1/21) + 3*atan(1/239) + 4*atan(3/1042) |
Testing: pi/4 = 16*atan(1/21) + 3*atan(1/239) + 4*atan(3/1042) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 22*atan(1/28) + 2*atan(1/443) + -5*atan(1/1393) + -10*atan(1/11018) |
Testing: pi/4 = 22*atan(1/28) + 2*atan(1/443) + -5*atan(1/1393) + -10*atan(1/11018) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 22*atan(1/38) + 17*atan(7/601) + 10*atan(7/8149) |
Testing: pi/4 = 22*atan(1/38) + 17*atan(7/601) + 10*atan(7/8149) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 44*atan(1/57) + 7*atan(1/239) + -12*atan(1/682) + 24*atan(1/12943) |
Testing: pi/4 = 44*atan(1/57) + 7*atan(1/239) + -12*atan(1/682) + 24*atan(1/12943) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 88*atan(1/172) + 51*atan(1/239) + 32*atan(1/682) + 44*atan(1/5357) + 68*atan(1/12943) |
Testing: pi/4 = 88*atan(1/172) + 51*atan(1/239) + 32*atan(1/682) + 44*atan(1/5357) + 68*atan(1/12943) |
||
tan(RHS) |
tan(RHS) is one |
||
Testing: pi/4 = 88*atan(1/172) + 51*atan(1/239) + 32*atan(1/682) + 44*atan(1/5357) + 68*atan(1/12944) |
Testing: pi/4 = 88*atan(1/172) + 51*atan(1/239) + 32*atan(1/682) + 44*atan(1/5357) + 68*atan(1/12944) |
||
tan(RHS) is not one |
|||
tan(RHS) = 1009288018000944050967896710431587186456256928584351786643498522649995492271475761189348270710224618853590682465929080006511691833816436374107451368838065354726517908250456341991684635768915704374493675498637876700129004484434187627909285979251682006538817341793224963346197503893270875008524149334251672855130857035205217929335932890740051319216343365800342290782260673215928499123722781078448297609548233999010983373327601187505623621602789012550584784738082074783523787011976757247516095289966708782862528690942242793667539020699840402353522108223/1009288837315638583415701528780402795721935641614456853534313491853293025565940011104051964874275710024625850092154664245109626053906509780125743180758231049920425664246286578958307532545458843067352531217230461290763258378749459637420702619029075083089762088232401888676895047947363883809724322868121990870409574061477638203859217672620508200713073485398199091153535700094640095900731630771349477187594074169815106104524371099618096164871416282464532355211521113449237814080332335526420331468258917484010722587072087349909684004660371264507984339711 |
|||
</pre> |
</pre> |
||
Revision as of 22:04, 8 November 2012
Machin-like formulas are useful for efficiently computing numerical approximations to Pi. Verify the following Machin-like formulas are correct by calculating the value of tan(right hand side) for each equation using exact arithmetic and showing they equal 1:
and confirm that the following formula is incorrect by showing tan(right hand side) is not 1:
These identities are useful in calculating the values:
Haskell
<lang haskell>import Data.Ratio
tanEval (1,f) = f tanEval (coef,f) | coef < 0 = -tanEval (-coef, f) | otherwise = (a + b) / (1 - a * b) where a = tanEval (ca, f) b = tanEval (cb, f) ca = coef `div` 2 cb = coef - ca
tans [x] = tanEval x tans xs = (a + b) / (1 - a * b) where a = tans aa b = tans bb (aa, bb) = splitAt ((length xs)`div`2) xs
machins = [ [(1, 1%2), (1, 1%3)], [(2, 1%3), (1, 1%7)], [(12, 1%18), (8, 1%57), (-5, 1%239)], [(88, 1%172), (51, 1%239), (32 , 1%682), (44, 1%5357), (68, 1%12943)]]
not_machin = [(88, 1%172), (51, 1%239), (32 , 1%682), (44, 1%5357), (68, 1%12944)]
main = do putStrLn "Machins:" mapM_ (\x->putStrLn $ show(tans x) ++ " <-- " ++ show(x)) machins
putStr "\nnot Machin: "; print not_machin print (tans not_machin)</lang>
OCaml
<lang ocaml>open Num;; (* use exact rationals for results *)
let tadd p q = (p +/ q) // ((Int 1) -/ (p */ q)) in
(* tan(n*arctan(a/b)) *) let rec tan_expr (n,a,b) =
if n = 1 then (Int a)//(Int b) else if n = -1 then (Int (-a))//(Int b) else let m = n/2 in let tm = tan_expr (m,a,b) in let m2 = tadd tm tm and k = n-m-m in if k = 0 then m2 else tadd (tan_expr (k,a,b)) m2 in
let verify (k, tlist) =
Printf.printf "Testing: pi/%d = " k; let t_str = List.map (fun (x,y,z) -> Printf.sprintf "%d*atan(%d/%d)" x y z) tlist in print_endline (String.concat " + " t_str); let ans_terms = List.map tan_expr tlist in let answer = List.fold_left tadd (Int 0) ans_terms in Printf.printf " tan(RHS) is %s\n" (if answer = (Int 1) then "one" else "not one") in
(* example: prog 4 5 29 278 7 3 79 represents pi/4 = 5*atan(29/278) + 7*atan(3/79) *) let args = Sys.argv in let nargs = Array.length args in let v k = int_of_string args.(k) in let rec triples n =
if n+2 > nargs-1 then [] else (v n, v (n+1), v (n+2)) :: triples (n+3) in
if nargs > 4 then let dat = (v 1, triples 2) in verify dat else List.iter verify [
(4,[(1,1,2);(1,1,3)]); (4,[(2,1,3);(1,1,7)]); (4,[(4,1,5);(-1,1,239)]); (4,[(5,1,7);(2,3,79)]); (4,[(5,29,278);(7,3,79)]); (4,[(1,1,2);(1,1,5);(1,1,8)]); (4,[(4,1,5);(-1,1,70);(1,1,99)]); (4,[(5,1,7);(4,1,53);(2,1,4443)]); (4,[(6,1,8);(2,1,57);(1,1,239)]); (4,[(8,1,10);(-1,1,239);(-4,1,515)]); (4,[(12,1,18);(8,1,57);(-5,1,239)]); (4,[(16,1,21);(3,1,239);(4,3,1042)]); (4,[(22,1,28);(2,1,443);(-5,1,1393);(-10,1,11018)]); (4,[(22,1,38);(17,7,601);(10,7,8149)]); (4,[(44,1,57);(7,1,239);(-12,1,682);(24,1,12943)]); (4,[(88,1,172);(51,1,239);(32,1,682);(44,1,5357);(68,1,12943)]); (4,[(88,1,172);(51,1,239);(32,1,682);(44,1,5357);(68,1,12944)])
]</lang>
Compile with
ocamlopt -o verify_machin.opt nums.cmxa verify_machin.ml
or run with
ocaml nums.cma verify_machin.ml
- Output:
Testing: pi/4 = 1*atan(1/2) + 1*atan(1/3) tan(RHS) is one Testing: pi/4 = 2*atan(1/3) + 1*atan(1/7) tan(RHS) is one Testing: pi/4 = 4*atan(1/5) + -1*atan(1/239) tan(RHS) is one Testing: pi/4 = 5*atan(1/7) + 2*atan(3/79) tan(RHS) is one Testing: pi/4 = 5*atan(29/278) + 7*atan(3/79) tan(RHS) is one Testing: pi/4 = 1*atan(1/2) + 1*atan(1/5) + 1*atan(1/8) tan(RHS) is one Testing: pi/4 = 4*atan(1/5) + -1*atan(1/70) + 1*atan(1/99) tan(RHS) is one Testing: pi/4 = 5*atan(1/7) + 4*atan(1/53) + 2*atan(1/4443) tan(RHS) is one Testing: pi/4 = 6*atan(1/8) + 2*atan(1/57) + 1*atan(1/239) tan(RHS) is one Testing: pi/4 = 8*atan(1/10) + -1*atan(1/239) + -4*atan(1/515) tan(RHS) is one Testing: pi/4 = 12*atan(1/18) + 8*atan(1/57) + -5*atan(1/239) tan(RHS) is one Testing: pi/4 = 16*atan(1/21) + 3*atan(1/239) + 4*atan(3/1042) tan(RHS) is one Testing: pi/4 = 22*atan(1/28) + 2*atan(1/443) + -5*atan(1/1393) + -10*atan(1/11018) tan(RHS) is one Testing: pi/4 = 22*atan(1/38) + 17*atan(7/601) + 10*atan(7/8149) tan(RHS) is one Testing: pi/4 = 44*atan(1/57) + 7*atan(1/239) + -12*atan(1/682) + 24*atan(1/12943) tan(RHS) is one Testing: pi/4 = 88*atan(1/172) + 51*atan(1/239) + 32*atan(1/682) + 44*atan(1/5357) + 68*atan(1/12943) tan(RHS) is one Testing: pi/4 = 88*atan(1/172) + 51*atan(1/239) + 32*atan(1/682) + 44*atan(1/5357) + 68*atan(1/12944) tan(RHS) is not one
Python
This example parses the original equations to form an intermediate representation then does the checks.
Function tans and tanEval are translations of the Haskel functions of the same names.
<lang python>import re
from fractions import Fraction
from pprint import pprint as pp
equationtext = \
pi/4 = arctan(1/2) + arctan(1/3) pi/4 = 2*arctan(1/3) + arctan(1/7) pi/4 = 4*arctan(1/5) - arctan(1/239) pi/4 = 5*arctan(1/7) + 2*arctan(3/79) pi/4 = 5*arctan(29/278) + 7*arctan(3/79) pi/4 = arctan(1/2) + arctan(1/5) + arctan(1/8) pi/4 = 4*arctan(1/5) - arctan(1/70) + arctan(1/99) pi/4 = 5*arctan(1/7) + 4*arctan(1/53) + 2*arctan(1/4443) pi/4 = 6*arctan(1/8) + 2*arctan(1/57) + arctan(1/239) pi/4 = 8*arctan(1/10) - arctan(1/239) - 4*arctan(1/515) pi/4 = 12*arctan(1/18) + 8*arctan(1/57) - 5*arctan(1/239) pi/4 = 16*arctan(1/21) + 3*arctan(1/239) + 4*arctan(3/1042) pi/4 = 22*arctan(1/28) + 2*arctan(1/443) - 5*arctan(1/1393) - 10*arctan(1/11018) pi/4 = 22*arctan(1/38) + 17*arctan(7/601) + 10*arctan(7/8149) pi/4 = 44*arctan(1/57) + 7*arctan(1/239) - 12*arctan(1/682) + 24*arctan(1/12943) pi/4 = 88*arctan(1/172) + 51*arctan(1/239) + 32*arctan(1/682) + 44*arctan(1/5357) + 68*arctan(1/12943) pi/4 = 88*arctan(1/172) + 51*arctan(1/239) + 32*arctan(1/682) + 44*arctan(1/5357) + 68*arctan(1/12944)
def parse_eqn(equationtext=equationtext):
eqn_re = re.compile(r"""(?mx) (?P<lhs> ^ \s* pi/4 \s* = \s*)? # LHS of equation (?: # RHS \s* (?P<sign> [+-])? \s* (?: (?P<mult> \d+) \s* \*)? \s* arctan\( (?P<numer> \d+) / (?P<denom> \d+) )""")
found = eqn_re.findall(equationtext) machins, part = [], [] for lhs, sign, mult, numer, denom in eqn_re.findall(equationtext): if lhs and part: machins.append(part) part = [] part.append( ( (-1 if sign == '-' else 1) * ( int(mult) if mult else 1), Fraction(int(numer), (int(denom) if denom else 1)) ) ) machins.append(part) return machins
def tans(xs):
xslen = len(xs) if xslen == 1: return tanEval(*xs[0]) aa, bb = xs[:xslen//2], xs[xslen//2:] a, b = tans(aa), tans(bb) return (a + b) / (1 - a * b)
def tanEval(coef, f):
if coef == 1: return f if coef < 0: return -tanEval(-coef, f) ca = coef // 2 cb = coef - ca a, b = tanEval(ca, f), tanEval(cb, f) return (a + b) / (1 - a * b)
if __name__ == '__main__':
machins = parse_eqn() #pp(machins, width=160) for machin, eqn in zip(machins, equationtext.split('\n')): ans = tans(machin) print('%5s: %s' % ( ('OK' if ans == 1 else 'ERROR'), eqn))</lang>
- Output:
OK: pi/4 = arctan(1/2) + arctan(1/3) OK: pi/4 = 2*arctan(1/3) + arctan(1/7) OK: pi/4 = 4*arctan(1/5) - arctan(1/239) OK: pi/4 = 5*arctan(1/7) + 2*arctan(3/79) OK: pi/4 = 5*arctan(29/278) + 7*arctan(3/79) OK: pi/4 = arctan(1/2) + arctan(1/5) + arctan(1/8) OK: pi/4 = 4*arctan(1/5) - arctan(1/70) + arctan(1/99) OK: pi/4 = 5*arctan(1/7) + 4*arctan(1/53) + 2*arctan(1/4443) OK: pi/4 = 6*arctan(1/8) + 2*arctan(1/57) + arctan(1/239) OK: pi/4 = 8*arctan(1/10) - arctan(1/239) - 4*arctan(1/515) OK: pi/4 = 12*arctan(1/18) + 8*arctan(1/57) - 5*arctan(1/239) OK: pi/4 = 16*arctan(1/21) + 3*arctan(1/239) + 4*arctan(3/1042) OK: pi/4 = 22*arctan(1/28) + 2*arctan(1/443) - 5*arctan(1/1393) - 10*arctan(1/11018) OK: pi/4 = 22*arctan(1/38) + 17*arctan(7/601) + 10*arctan(7/8149) OK: pi/4 = 44*arctan(1/57) + 7*arctan(1/239) - 12*arctan(1/682) + 24*arctan(1/12943) OK: pi/4 = 88*arctan(1/172) + 51*arctan(1/239) + 32*arctan(1/682) + 44*arctan(1/5357) + 68*arctan(1/12943) ERROR: pi/4 = 88*arctan(1/172) + 51*arctan(1/239) + 32*arctan(1/682) + 44*arctan(1/5357) + 68*arctan(1/12944)