Check Machin-like formulas

Revision as of 00:23, 8 November 2012 by rosettacode>TobyK (Use exact arithmetic to verify Machin-like formulas)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

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:

Check Machin-like formulas is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
 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)

and confirm that the following formula is incorrect by showing tan(right hand side) is not 1:

 pi/4 = 88*arctan(1/172) + 51*arctan(1/239) + 32*arctan(1/682) + 44*arctan(1/5357) + 68*arctan(1/12944)

These identities are useful in calculating the values:

 tan(a + b) = (tan(a) + tan(b))/(1 - tan(a)*tan(b))
 tan(arctan(a/b)) = a/b
 tan(-a) = -tan(a)

OCaml

<lang ocaml>open Num;; (* use exact rationals for results *)

let tadd p q = (p +/ q) // ((Int 1) -/ (p */ q)) in

let rec tsum = function

| [] -> failwith "empty tsum"
| [a] -> a
| [a;b] -> tadd a b
| h::t -> tadd h (tsum t) 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 = tsum ans_terms in
 Printf.printf "  tan(RHS) = %s\n" (string_of_num answer);
 let pi = 2.0*.acos(0.0) in
 Printf.printf "  tan(LHS) ~= %.15f\n" (tan (pi/.float k)) 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)])

]</lang>

Output:
Testing: pi/4 = 1*atan(1/2) + 1*atan(1/3)
  tan(RHS) = 1
Testing: pi/4 = 2*atan(1/3) + 1*atan(1/7)
  tan(RHS) = 1
Testing: pi/4 = 4*atan(1/5) + -1*atan(1/239)
  tan(RHS) = 1
Testing: pi/4 = 5*atan(1/7) + 2*atan(3/79)
  tan(RHS) = 1
Testing: pi/4 = 5*atan(29/278) + 7*atan(3/79)
  tan(RHS) = 1
Testing: pi/4 = 1*atan(1/2) + 1*atan(1/5) + 1*atan(1/8)
  tan(RHS) = 1
Testing: pi/4 = 4*atan(1/5) + -1*atan(1/70) + 1*atan(1/99)
  tan(RHS) = 1
Testing: pi/4 = 5*atan(1/7) + 4*atan(1/53) + 2*atan(1/4443)
  tan(RHS) = 1
Testing: pi/4 = 6*atan(1/8) + 2*atan(1/57) + 1*atan(1/239)
  tan(RHS) = 1
Testing: pi/4 = 8*atan(1/10) + -1*atan(1/239) + -4*atan(1/515)
  tan(RHS) = 1
Testing: pi/4 = 12*atan(1/18) + 8*atan(1/57) + -5*atan(1/239)
  tan(RHS) = 1
Testing: pi/4 = 16*atan(1/21) + 3*atan(1/239) + 4*atan(3/1042)
  tan(RHS) = 1
Testing: pi/4 = 22*atan(1/28) + 2*atan(1/443) + -5*atan(1/1393) + -10*atan(1/11018)
  tan(RHS) = 1
Testing: pi/4 = 22*atan(1/38) + 17*atan(7/601) + 10*atan(7/8149)
  tan(RHS) = 1
Testing: pi/4 = 44*atan(1/57) + 7*atan(1/239) + -12*atan(1/682) + 24*atan(1/12943)
  tan(RHS) = 1
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) = 1
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) = 1009288018000944050967896710431587186456256928584351786643498522649995492271475761189348270710224618853590682465929080006511691833816436374107451368838065354726517908250456341991684635768915704374493675498637876700129004484434187627909285979251682006538817341793224963346197503893270875008524149334251672855130857035205217929335932890740051319216343365800342290782260673215928499123722781078448297609548233999010983373327601187505623621602789012550584784738082074783523787011976757247516095289966708782862528690942242793667539020699840402353522108223/1009288837315638583415701528780402795721935641614456853534313491853293025565940011104051964874275710024625850092154664245109626053906509780125743180758231049920425664246286578958307532545458843067352531217230461290763258378749459637420702619029075083089762088232401888676895047947363883809724322868121990870409574061477638203859217672620508200713073485398199091153535700094640095900731630771349477187594074169815106104524371099618096164871416282464532355211521113449237814080332335526420331468258917484010722587072087349909684004660371264507984339711