Check Machin-like formulas
You are encouraged to solve this task according to the task description, using any language you may know.
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:
You can store the equations in any convenient data structure, but for extra credit parse them from human-readable text input.
Note that to formally prove the formula correct you would also have to show that < right hand side < due to periodicity.
D
This uses the module of the Arithmetic Rational Task.
<lang d>import std.stdio, std.regex, std.conv, std.string, std.range,
arithmetic_rational;
struct Pair { int x; Rational r; }
Pair[][] parseEquations(in string text) /*pure nothrow*/ {
auto r = regex(r"\s*(?P<sign>[+-])?\s*(?:(?P<mul>\d+)\s*\*)?\s*" ~ r"arctan\((?P<num>\d+)/(?P<denom>\d+)\)"); Pair[][] machins; foreach (const line; text.splitLines) { Pair[] formula; foreach (part; line.split("=")[1].matchAll(r)) { immutable mul = part["mul"], num = part["num"], denom = part["denom"]; formula ~= Pair((part["sign"] == "-" ? -1 : 1) * (mul.empty ? 1 : mul.to!int), Rational(num.to!int, denom.empty ? 1 : denom.to!int)); } machins ~= formula; } return machins;
}
Rational tans(in Pair[] xs) pure nothrow {
static Rational tanEval(in int coef, in Rational f) pure nothrow { if (coef == 1) return f; if (coef < 0) return -tanEval(-coef, f); immutable a = tanEval(coef / 2, f), b = tanEval(coef - coef / 2, f); return (a + b) / (1 - a * b); }
if (xs.length == 1) return tanEval(xs[0].tupleof); immutable a = xs[0 .. $ / 2].tans, b = xs[$ / 2 .. $].tans; return (a + b) / (1 - a * b);
}
void main() {
immutable 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)";
const machins = equationText.parseEquations; foreach (const machin, const eqn; machins.zip(equationText.splitLines)) { immutable ans = machin.tans; writefln("%5s: %s", ans == 1 ? "OK" : "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)
EchoLisp
<lang scheme> (lib 'math) (lib 'match) (math-precision 1.e-10)
- formally derive (tan ..) expressions
- copied from Racket
- adapted and improved for performance
(define (reduce e)
- (set! rcount (1+ rcount)) ;; # of calls
(match e [(? number? a) a] [('+ (? number? a) (? number? b)) (+ a b)] [('- (? number? a) (? number? b)) (- a b)] [('- (? number? a)) (- a)] [('* (? number? a) (? number? b)) (* a b)] [('/ (? number? a) (? number? b)) (/ a b)] ; patch [( '+ a b) (reduce `(+ ,(reduce a) ,(reduce b)))] [( '- a b) (reduce `(- ,(reduce a) ,(reduce b)))] [( '- a) (reduce `(- ,(reduce a)))] [( '* a b) (reduce `(* ,(reduce a) ,(reduce b)))] [( '/ a b) (reduce `(/ ,(reduce a) ,(reduce b)))] [( 'tan ('arctan a)) (reduce a)] [( 'tan ( '- a)) (reduce `(- (tan ,a)))]
;; x 100 # calls reduction : derive (tan ,a) only once [( 'tan ( '+ a b)) (let ((alpha (reduce `(tan ,a))) (beta (reduce `(tan ,b)))) (reduce `(/ (+ ,alpha ,beta) (- 1 (* ,alpha ,beta)))))] [( 'tan ( '+ a b c ...)) (reduce `(tan (+ ,a (+ ,b ,@c))))] [( 'tan ( '- a b)) (let ((alpha (reduce `(tan ,a))) (beta (reduce `(tan ,b)))) (reduce `(/ (- ,alpha ,beta) (+ 1 (* ,alpha ,beta)))))]
;; add formula for (tan 2 (arctan a)) = 2 a / (1 - a^2)) [( 'tan ( '* 2 ('arctan a))) (reduce `(/ (* 2 ,a) (- 1 (* ,a ,a))))] [( 'tan ( '* 1 ('arctan a))) (reduce a)] ; added [( 'tan ( '* (? number? n) a)) (cond [(< n 0) (reduce `(- (tan (* ,(- n) ,a))))] [(= n 0) 0] [(= n 1) (reduce `(tan ,a))] [(even? n) (let ((alpha (reduce `(tan (* ,(/ n 2) ,a))))) ;; # calls reduction (reduce `(/ (* 2 ,alpha) (- 1 (* ,alpha ,alpha)))))] [else (reduce `(tan (+ ,a (* ,(- n 1) ,a))))])] ))
(define (task) (for ((f machins)) (if (~= 1 (reduce f)) (writeln '👍 f '⟾ 1 ) (writeln '❌ f '➽ (reduce f) ))))
</lang>
- Output:
<lang scheme>
(define machins
'((tan (+ (arctan 1/2) (arctan 1/3))) (tan (+ (* 2 (arctan 1/3)) (arctan 1/7))) (tan (- (* 4 (arctan 1/5)) (arctan 1/239))) (tan (+ (* 5 (arctan 1/7)) (* 2 (arctan 3/79)))) (tan (+ (* 5 (arctan 29/278)) (* 7 (arctan 3/79)))) (tan (+ (arctan 1/2) (arctan 1/5) (arctan 1/8))) (tan (+ (* 4 (arctan 1/5)) (* -1 (arctan 1/70)) (arctan 1/99))) (tan (+ (* 5 (arctan 1/7)) (* 4 (arctan 1/53)) (* 2 (arctan 1/4443)))) (tan (+ (* 6 (arctan 1/8)) (* 2 (arctan 1/57)) (arctan 1/239))) (tan (+ (* 8 (arctan 1/10)) (* -1 (arctan 1/239)) (* -4 (arctan 1/515)))) (tan (+ (* 12 (arctan 1/18)) (* 8 (arctan 1/57)) (* -5 (arctan 1/239)))) (tan (+ (* 16 (arctan 1/21)) (* 3 (arctan 1/239)) (* 4 (arctan 3/1042)))) (tan (+ (* 22 (arctan 1/28)) (* 2 (arctan 1/443)) (* -5 (arctan 1/1393)) (* -10 (arctan 1/11018)))) (tan (+ (* 22 (arctan 1/38)) (* 17 (arctan 7/601)) (* 10 (arctan 7/8149)))) (tan (+ (* 44 (arctan 1/57)) (* 7 (arctan 1/239)) (* -12 (arctan 1/682)) (* 24 (arctan 1/12943)))) (tan (+ (* 88 (arctan 1/172)) (* 51 (arctan 1/239)) (* 32 (arctan 1/682)) (* 44 (arctan 1/5357)) (* 68 (arctan 1/12943)))) (tan (+ (* 88 (arctan 1/172)) (* 51 (arctan 1/239)) (* 32 (arctan 1/682)) (* 44 (arctan 1/5357)) (* 68 (arctan 1/12944))))))
(task)
👍 (tan (+ (arctan 1/2) (arctan 1/3))) ⟾ 1 👍 (tan (+ (* 2 (arctan 1/3)) (arctan 1/7))) ⟾ 1 👍 (tan (- (* 4 (arctan 1/5)) (arctan 1/239))) ⟾ 1 👍 (tan (+ (* 5 (arctan 1/7)) (* 2 (arctan 3/79)))) ⟾ 1 👍 (tan (+ (* 5 (arctan 29/278)) (* 7 (arctan 3/79)))) ⟾ 1 👍 (tan (+ (arctan 1/2) (arctan 1/5) (arctan 1/8))) ⟾ 1 👍 (tan (+ (* 4 (arctan 1/5)) (* -1 (arctan 1/70)) (arctan 1/99))) ⟾ 1 👍 (tan (+ (* 5 (arctan 1/7)) (* 4 (arctan 1/53)) (* 2 (arctan 1/4443)))) ⟾ 1 👍 (tan (+ (* 6 (arctan 1/8)) (* 2 (arctan 1/57)) (arctan 1/239))) ⟾ 1 👍 (tan (+ (* 8 (arctan 1/10)) (* -1 (arctan 1/239)) (* -4 (arctan 1/515)))) ⟾ 1 👍 (tan (+ (* 12 (arctan 1/18)) (* 8 (arctan 1/57)) (* -5 (arctan 1/239)))) ⟾ 1 👍 (tan (+ (* 16 (arctan 1/21)) (* 3 (arctan 1/239)) (* 4 (arctan 3/1042)))) ⟾ 1 👍 (tan (+ (* 22 (arctan 1/28)) (* 2 (arctan 1/443)) (* -5 (arctan 1/1393)) (* -10 (arctan 1/11018)))) ⟾ 1 👍 (tan (+ (* 22 (arctan 1/38)) (* 17 (arctan 7/601)) (* 10 (arctan 7/8149)))) ⟾ 1 👍 (tan (+ (* 44 (arctan 1/57)) (* 7 (arctan 1/239)) (* -12 (arctan 1/682)) (* 24 (arctan 1/12943)))) ⟾ 1 👍 (tan (+ (* 88 (arctan 1/172)) (* 51 (arctan 1/239)) (* 32 (arctan 1/682))
(* 44 (arctan 1/5357)) (* 68 (arctan 1/12943)))) ⟾ 1
❌ (tan (+ (* 88 (arctan 1/172)) (* 51 (arctan 1/239)) (* 32 (arctan 1/682))
(* 44 (arctan 1/5357)) (* 68 (arctan 1/12944)))) ➽ 0.9999991882257442
</lang>
GAP
The formula is entered as a list of pairs [k, x], where each pair means , and all the terms in the list are summed. Like most other solutions, the program will only check that the tangent of the resulting sum is . For instance, Check([[5, 1/2], [5, 1/3]]);
returns also true
, though the result is .
<lang gap>TanPlus := function(a, b)
return (a + b) / (1 - a * b);
end;
TanTimes := function(n, a)
local x; x := 0; while n > 0 do if IsOddInt(n) then x := TanPlus(x, a); fi; a := TanPlus(a, a); n := QuoInt(n, 2); od; return x;
end;
Check := function(a)
local x, p; x := 0; for p in a do x := TanPlus(x, SignInt(p[1]) * TanTimes(AbsInt(p[1]), p[2])); od; return x = 1;
end;
ForAll([
[[1, 1/2], [1, 1/3]], [[2, 1/3], [1, 1/7]], [[4, 1/5], [-1, 1/239]], [[5, 1/7], [2, 3/79]], [[5, 29/278], [7, 3/79]], [[1, 1/2], [1, 1/5], [1, 1/8]], [[5, 1/7], [4, 1/53], [2, 1/4443]], [[6, 1/8], [2, 1/57], [1, 1/239]], [[8, 1/10], [-1, 1/239], [-4, 1/515]], [[12, 1/18], [8, 1/57], [-5, 1/239]], [[16, 1/21], [3, 1/239], [4, 3/1042]], [[22, 1/28], [2, 1/443], [-5, 1/1393], [-10, 1/11018]], [[22, 1/38], [17, 7/601], [10, 7/8149]], [[44, 1/57], [7, 1/239], [-12, 1/682], [24, 1/12943]], [[88, 1/172], [51, 1/239], [32, 1/682], [44, 1/5357], [68, 1/12943]]], Check);
Check([[88, 1/172], [51, 1/239], [32, 1/682], [44, 1/5357], [68, 1/12944]]);</lang>
Go
<lang go>package main
import (
"fmt" "math/big"
)
type mTerm struct {
a, n, d int64
}
var testCases = [][]mTerm{
{{1, 1, 2}, {1, 1, 3}}, {{2, 1, 3}, {1, 1, 7}}, {{4, 1, 5}, {-1, 1, 239}}, {{5, 1, 7}, {2, 3, 79}}, {{1, 1, 2}, {1, 1, 5}, {1, 1, 8}}, {{4, 1, 5}, {-1, 1, 70}, {1, 1, 99}}, {{5, 1, 7}, {4, 1, 53}, {2, 1, 4443}}, {{6, 1, 8}, {2, 1, 57}, {1, 1, 239}}, {{8, 1, 10}, {-1, 1, 239}, {-4, 1, 515}}, {{12, 1, 18}, {8, 1, 57}, {-5, 1, 239}}, {{16, 1, 21}, {3, 1, 239}, {4, 3, 1042}}, {{22, 1, 28}, {2, 1, 443}, {-5, 1, 1393}, {-10, 1, 11018}}, {{22, 1, 38}, {17, 7, 601}, {10, 7, 8149}}, {{44, 1, 57}, {7, 1, 239}, {-12, 1, 682}, {24, 1, 12943}}, {{88, 1, 172}, {51, 1, 239}, {32, 1, 682}, {44, 1, 5357}, {68, 1, 12943}}, {{88, 1, 172}, {51, 1, 239}, {32, 1, 682}, {44, 1, 5357}, {68, 1, 12944}},
}
func main() {
for _, m := range testCases { fmt.Printf("tan %v = %v\n", m, tans(m)) }
}
var one = big.NewRat(1, 1)
func tans(m []mTerm) *big.Rat {
if len(m) == 1 { return tanEval(m[0].a, big.NewRat(m[0].n, m[0].d)) } half := len(m) / 2 a := tans(m[:half]) b := tans(m[half:]) r := new(big.Rat) return r.Quo(new(big.Rat).Add(a, b), r.Sub(one, r.Mul(a, b)))
}
func tanEval(coef int64, f *big.Rat) *big.Rat {
if coef == 1 { return f } if coef < 0 { r := tanEval(-coef, f) return r.Neg(r) } ca := coef / 2 cb := coef - ca a := tanEval(ca, f) b := tanEval(cb, f) r := new(big.Rat) return r.Quo(new(big.Rat).Add(a, b), r.Sub(one, r.Mul(a, b)))
}</lang>
- Output:
Last line edited to show only most significant digits of fraction which is near, but not exactly equal to 1.
tan [{1 1 2} {1 1 3}] = 1/1 tan [{2 1 3} {1 1 7}] = 1/1 tan [{4 1 5} {-1 1 239}] = 1/1 tan [{5 1 7} {2 3 79}] = 1/1 tan [{1 1 2} {1 1 5} {1 1 8}] = 1/1 tan [{4 1 5} {-1 1 70} {1 1 99}] = 1/1 tan [{5 1 7} {4 1 53} {2 1 4443}] = 1/1 tan [{6 1 8} {2 1 57} {1 1 239}] = 1/1 tan [{8 1 10} {-1 1 239} {-4 1 515}] = 1/1 tan [{12 1 18} {8 1 57} {-5 1 239}] = 1/1 tan [{16 1 21} {3 1 239} {4 3 1042}] = 1/1 tan [{22 1 28} {2 1 443} {-5 1 1393} {-10 1 11018}] = 1/1 tan [{22 1 38} {17 7 601} {10 7 8149}] = 1/1 tan [{44 1 57} {7 1 239} {-12 1 682} {24 1 12943}] = 1/1 tan [{88 1 172} {51 1 239} {32 1 682} {44 1 5357} {68 1 12943}] = 1/1 tan [{88 1 172} {51 1 239} {32 1 682} {44 1 5357} {68 1 12944}] = 100928801... / 100928883...
Haskell
<lang haskell>import Data.Ratio import Data.List (foldl')
tanPlus :: Fractional a => a -> a -> a tanPlus a b = (a + b) / (1 - a * b)
tanEval :: (Integral a, Fractional b) => (a, b) -> b tanEval (0,_) = 0 tanEval (coef,f) | coef < 0 = -tanEval (-coef, f) | odd coef = tanPlus f $ tanEval (coef - 1, f) | otherwise = tanPlus a a where a = tanEval (coef `div` 2, f)
tans :: (Integral a, Fractional b) => [(a, b)] -> b tans = foldl' tanPlus 0 . map tanEval
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>
A crazier way to do the above, exploiting the built-in exponentiation algorithms: <lang haskell>import Data.Ratio
-- Private type. Do not use outside of the tans function newtype Tan a = Tan a deriving (Eq, Show) instance Fractional a => Num (Tan a) where
_ + _ = undefined Tan a * Tan b = Tan $ (a + b) / (1 - a * b) negate _ = undefined abs _ = undefined signum _ = undefined fromInteger 1 = Tan 0 -- identity for the (*) above fromInteger _ = undefined
instance Fractional a => Fractional (Tan a) where
fromRational _ = undefined recip (Tan f) = Tan (-f) -- inverse for the (*) above
tans :: (Integral a, Fractional b) => [(a, b)] -> b tans xs = x where
Tan x = product [Tan f ^^ coef | (coef,f) <- 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>
J
Solution:<lang j> machin =: 1r4p1 = [: +/ ({. * _3 o. %/@:}.)"1@:x:</lang> Example (test cases from task description):<lang j> R =: <@:(0&".);._2 ];._2 noun define 1 1 2 1 1 3
2 1 3 1 1 7
4 1 5 _1 1 239
5 1 7 2 3 79
5 29 278 7 3 79
1 1 2 1 1 5 1 1 8
4 1 5 _1 1 70 1 1 99
5 1 7 4 1 53 2 1 4443
6 1 8 2 1 57 1 1 239
8 1 10 _1 1 239 _4 1 515
12 1 18 8 1 57 _5 1 239
16 1 21 3 1 239 4 3 1042
22 1 28 2 1 443 _5 1 1393 _10 1 11018
22 1 38 17 7 601 10 7 8149
44 1 57 7 1 239 _12 1 682 24 1 12943
88 1 172 51 1 239 32 1 682 44 1 5357 68 1 12943
)
machin&> R 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1</lang> Example (counterexample):<lang j> counterExample=. 12944 (<_1;_1)} >{:R counterExample NB. Same as final test case with 12943 incremented to 12944 88 1 172 51 1 239 32 1 682 44 1 5357 68 1 12944 machin counterExample 0</lang> Notes: The function machin compares the results of each formula to π/4 (expressed as 1r4p1 in J's numeric notation). The first example above shows the results of these comparisons for each formula (with 1 for true and 0 for false). In J, arctan is expressed as _3 o. values and the function x: coerces values to exact representation; thereafter J will maintain exactness throughout its calculations, as long as it can.
Mathematica / Wolfram Language
<lang>Tan[ArcTan[1/2] + ArcTan[1/3]] == 1 Tan[2 ArcTan[1/3] + ArcTan[1/7]] == 1 Tan[4 ArcTan[1/5] - ArcTan[1/239]] == 1 Tan[5 ArcTan[1/7] + 2 ArcTan[3/79]] == 1 Tan[5 ArcTan[29/278] + 7 ArcTan[3/79]] == 1 Tan[ArcTan[1/2] + ArcTan[1/5] + ArcTan[1/8]] == 1 Tan[4 ArcTan[1/5] - ArcTan[1/70] + ArcTan[1/99]] == 1 Tan[5 ArcTan[1/7] + 4 ArcTan[1/53] + 2 ArcTan[1/4443]] == 1 Tan[6 ArcTan[1/8] + 2 ArcTan[1/57] + ArcTan[1/239]] == 1 Tan[8 ArcTan[1/10] - ArcTan[1/239] - 4 ArcTan[1/515]] == 1 Tan[12 ArcTan[1/18] + 8 ArcTan[1/57] - 5 ArcTan[1/239]] == 1 Tan[16 ArcTan[1/21] + 3 ArcTan[1/239] + 4 ArcTan[3/1042]] == 1 Tan[22 ArcTan[1/28] + 2 ArcTan[1/443] - 5 ArcTan[1/1393] -
10 ArcTan[1/11018]] == 1
Tan[22 ArcTan[1/38] + 17 ArcTan[7/601] + 10 ArcTan[7/8149]] == 1 Tan[44 ArcTan[1/57] + 7 ArcTan[1/239] - 12 ArcTan[1/682] +
24 ArcTan[1/12943]] == 1
Tan[88 ArcTan[1/172] + 51 ArcTan[1/239] + 32 ArcTan[1/682] +
44 ArcTan[1/5357] + 68 ArcTan[1/12943]] == 1
Tan[88 ArcTan[1/172] + 51 ArcTan[1/239] + 32 ArcTan[1/682] +
44 ArcTan[1/5357] + 68 ArcTan[1/12944]] == 1</lang>
- Output:
True True True True True True True True True True True True True True True True False
Maxima
<lang maxima>trigexpand:true$ is(tan(atan(1/2)+atan(1/3))=1); is(tan(2*atan(1/3)+atan(1/7))=1); is(tan(4*atan(1/5)-atan(1/239))=1); is(tan(5*atan(1/7)+2*atan(3/79))=1); is(tan(5*atan(29/278)+7*atan(3/79))=1); is(tan(atan(1/2)+atan(1/5)+atan(1/8))=1); is(tan(4*atan(1/5)-atan(1/70)+atan(1/99))=1); is(tan(5*atan(1/7)+4*atan(1/53)+2*atan(1/4443))=1); is(tan(6*atan(1/8)+2*atan(1/57)+atan(1/239))=1); is(tan(8*atan(1/10)-atan(1/239)-4*atan(1/515))=1); is(tan(12*atan(1/18)+8*atan(1/57)-5*atan(1/239))=1); is(tan(16*atan(1/21)+3*atan(1/239)+4*atan(3/1042))=1); is(tan(22*atan(1/28)+2*atan(1/443)-5*atan(1/1393)-10*atan(1/11018))=1); is(tan(22*atan(1/38)+17*atan(7/601)+10*atan(7/8149))=1); is(tan(44*atan(1/57)+7*atan(1/239)-12*atan(1/682)+24*atan(1/12943))=1); is(tan(88*atan(1/172)+51*atan(1/239)+32*atan(1/682)+44*atan(1/5357)+68*atan(1/12943))=1); is(tan(88*atan(1/172)+51*atan(1/239)+32*atan(1/682)+44*atan(1/5357)+68*atan(1/12944))=1);</lang>
- Output:
(%i2) (%o2) true (%i3) (%o3) true (%i4) (%o4) true (%i5) (%o5) true (%i6) (%o6) true (%i7) (%o7) true (%i8) (%o8) true (%i9) (%o9) true (%i10) (%o10) true (%i11) (%o11) true (%i12) (%o12) true (%i13) (%o13) true (%i14) (%o14) true (%i15) (%o15) true (%i16) (%o16) true (%i17) (%o17) true (%i18) (%o18) false
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
ooRexx
<lang oorexx>/*REXX ----------------------------------------------------------------
- 09.04.2014 Walter Pachl the REXX solution adapted for ooRexx
- which provides a function package rxMath
- --------------------------------------------------------------------*/
Numeric Digits 16 Numeric Fuzz 3; pi=rxCalcpi(); a.=
a.1 = 'pi/4 = rxCalcarctan(1/2,16,'R') + rxCalcarctan(1/3,16,'R')' a.2 = 'pi/4 = 2*rxCalcarctan(1/3,16,'R') + rxCalcarctan(1/7,16,'R')' a.3 = 'pi/4 = 4*rxCalcarctan(1/5,16,'R') - rxCalcarctan(1/239,16,'R')' a.4 = 'pi/4 = 5*rxCalcarctan(1/7,16,'R') + 2*rxCalcarctan(3/79,16,'R')' a.5 = 'pi/4 = 5*rxCalcarctan(29/278,16,'R') + 7*rxCalcarctan(3/79,16,'R')' a.6 = 'pi/4 = rxCalcarctan(1/2,16,'R') + rxCalcarctan(1/5,16,'R') + rxCalcarctan(1/8,16,'R')' a.7 = 'pi/4 = 4*rxCalcarctan(1/5,16,'R') - rxCalcarctan(1/70,16,'R') + rxCalcarctan(1/99,16,'R')' a.8 = 'pi/4 = 5*rxCalcarctan(1/7,16,'R') + 4*rxCalcarctan(1/53,16,'R') + 2*rxCalcarctan(1/4443,16,'R')' a.9 = 'pi/4 = 6*rxCalcarctan(1/8,16,'R') + 2*rxCalcarctan(1/57,16,'R') + rxCalcarctan(1/239,16,'R')'
a.10 = 'pi/4 = 8*rxCalcarctan(1/10,16,'R') - rxCalcarctan(1/239,16,'R') - 4*rxCalcarctan(1/515,16,'R')' a.11 = 'pi/4 = 12*rxCalcarctan(1/18,16,'R') + 8*rxCalcarctan(1/57,16,'R') - 5*rxCalcarctan(1/239,16,'R')' a.12 = 'pi/4 = 16*rxCalcarctan(1/21,16,'R') + 3*rxCalcarctan(1/239,16,'R') + 4*rxCalcarctan(3/1042,16,'R')' a.13 = 'pi/4 = 22*rxCalcarctan(1/28,16,'R') + 2*rxCalcarctan(1/443,16,'R') - 5*rxCalcarctan(1/1393,16,'R') - 10*rxCalcarctan(1/11018,16,'R')' a.14 = 'pi/4 = 22*rxCalcarctan(1/38,16,'R') + 17*rxCalcarctan(7/601,16,'R') + 10*rxCalcarctan(7/8149,16,'R')' a.15 = 'pi/4 = 44*rxCalcarctan(1/57,16,'R') + 7*rxCalcarctan(1/239,16,'R') - 12*rxCalcarctan(1/682,16,'R') + 24*rxCalcarctan(1/12943,16,'R')' a.16 = 'pi/4 = 88*rxCalcarctan(1/172,16,'R') + 51*rxCalcarctan(1/239,16,'R') + 32*rxCalcarctan(1/682,16,'R') + 44*rxCalcarctan(1/5357,16,'R') + 68*rxCalcarctan(1/12943,16,'R')' a.17 = 'pi/4 = 88*rxCalcarctan(1/172,16,'R') + 51*rxCalcarctan(1/239,16,'R') + 32*rxCalcarctan(1/682,16,'R') + 44*rxCalcarctan(1/5357,16,'R') + 68*rxCalcarctan(1/12944,16,'R')'
do j=1 while a.j\== /*evaluate each of the formulas. */ interpret 'answer=' "(" a.j ")" /*the heavy lifting.*/ say right(word('bad OK',answer+1),3)": " space(a.j,0) end /*j*/ /* [?] show OK | bad, formula. */
- requires rxmath library
</lang>
- Output:
OK: pi/4=rxCalcarctan(1/2,16,R)+rxCalcarctan(1/3,16,R) OK: pi/4=2*rxCalcarctan(1/3,16,R)+rxCalcarctan(1/7,16,R) OK: pi/4=4*rxCalcarctan(1/5,16,R)-rxCalcarctan(1/239,16,R) OK: pi/4=5*rxCalcarctan(1/7,16,R)+2*rxCalcarctan(3/79,16,R) OK: pi/4=5*rxCalcarctan(29/278,16,R)+7*rxCalcarctan(3/79,16,R) OK: pi/4=rxCalcarctan(1/2,16,R)+rxCalcarctan(1/5,16,R)+rxCalcarctan(1/8,16,R) OK: pi/4=4*rxCalcarctan(1/5,16,R)-rxCalcarctan(1/70,16,R)+rxCalcarctan(1/99,16,R) OK: pi/4=5*rxCalcarctan(1/7,16,R)+4*rxCalcarctan(1/53,16,R)+2*rxCalcarctan(1/4443,16,R) OK: pi/4=6*rxCalcarctan(1/8,16,R)+2*rxCalcarctan(1/57,16,R)+rxCalcarctan(1/239,16,R) OK: pi/4=8*rxCalcarctan(1/10,16,R)-rxCalcarctan(1/239,16,R)-4*rxCalcarctan(1/515,16,R) OK: pi/4=12*rxCalcarctan(1/18,16,R)+8*rxCalcarctan(1/57,16,R)-5*rxCalcarctan(1/239,16,R) OK: pi/4=16*rxCalcarctan(1/21,16,R)+3*rxCalcarctan(1/239,16,R)+4*rxCalcarctan(3/1042,16,R) OK: pi/4=22*rxCalcarctan(1/28,16,R)+2*rxCalcarctan(1/443,16,R)-5*rxCalcarctan(1/1393,16,R)-10*rxCalcarctan(1/11018,16,R) OK: pi/4=22*rxCalcarctan(1/38,16,R)+17*rxCalcarctan(7/601,16,R)+10*rxCalcarctan(7/8149,16,R) OK: pi/4=44*rxCalcarctan(1/57,16,R)+7*rxCalcarctan(1/239,16,R)-12*rxCalcarctan(1/682,16,R)+24*rxCalcarctan(1/12943,16,R) OK: pi/4=88*rxCalcarctan(1/172,16,R)+51*rxCalcarctan(1/239,16,R)+32*rxCalcarctan(1/682,16,R)+44*rxCalcarctan(1/5357,16,R)+68*rxCalcarctan(1/12943,16,R) bad: pi/4=88*rxCalcarctan(1/172,16,R)+51*rxCalcarctan(1/239,16,R)+32*rxCalcarctan(1/682,16,R)+44*rxCalcarctan(1/5357,16,R)+68*rxCalcarctan(1/12944,16,R)
PARI/GP
<lang parigp>tanEval(coef, f)={ if (coef <= 1, return(if(coef<1,-tanEval(-coef, f),f))); my(a=tanEval(coef\2, f), b=tanEval(coef-coef\2, f)); (a + b)/(1 - a*b) }; tans(xs)={ if (#xs == 1, return(tanEval(xs[1][1], xs[1][2]))); my(a=tans(xs[1..#xs\2]),b=tans(xs[#xs\2+1..#xs])); (a + b)/(1 - a*b) }; test(v)={ my(t=tans(v)); if(t==1,print("OK"),print("Error: "v)) }; test([[1,1/2],[1,1/3]]); test([[2,1/3],[1,1/7]]); test([[4,1/5],[-1,1/239]]); test([[5,1/7],[2,3/79]]); test([[5,29/278],[7,3/79]]); test([[1,1/2],[1,1/5],[1,1/8]]); test([[4,1/5],[-1,1/70],[1,1/99]]); test([[5,1/7],[4,1/53],[2,1/4443]]); test([[6,1/8],[2,1/57],[1,1/239]]); test([[8,1/10],[-1,1/239],[-4,1/515]]); test([[12,1/18],[8,1/57],[-5,1/239]]); test([[16,1/21],[3,1/239],[4,3/1042]]); test([[22,1/28],[2,1/443],[-5,1/1393],[-10,1/11018]]); test([[22,1/38],[17,7/601],[10,7/8149]]); test([[44,1/57],[7,1/239],[-12,1/682],[24,1/12943]]); test([[88,1/172],[51,1/239],[32,1/682],[44,1/5357],[68,1/12943]]); test([[88,1/172],[51,1/239],[32,1/682],[44,1/5357],[68,1/12944]]);</lang>
- Output:
OK OK OK OK OK OK OK OK OK OK OK OK OK OK OK OK Error: [[88, 1/172], [51, 1/239], [32, 1/682], [44, 1/5357], [68, 1/12944]]
Perl
<lang Perl>use Math::BigRat try=>"GMP";
sub taneval {
my($coef,$f) = @_; $f = Math::BigRat->new($f) unless ref($f); return 0 if $coef == 0; return $f if $coef == 1; return -taneval(-$coef, $f) if $coef < 0; my($a,$b) = ( taneval($coef>>1, $f), taneval($coef-($coef>>1),$f) ); ($a+$b)/(1-$a*$b);
}
sub tans {
my @xs=@_; return taneval(@{$xs[0]}) if scalar(@xs)==1; my($a,$b) = ( tans(@xs[0..($#xs>>1)]), tans(@xs[($#xs>>1)+1..$#xs]) ); ($a+$b)/(1-$a*$b);
}
sub test {
printf "%5s (%s)\n", (tans(@_)==1)?"OK":"Error", join(" ",map{"[@$_]"} @_);
}
test([1,'1/2'], [1,'1/3']); test([2,'1/3'], [1,'1/7']); test([4,'1/5'], [-1,'1/239']); test([5,'1/7'],[2,'3/79']); test([5,'29/278'],[7,'3/79']); test([1,'1/2'],[1,'1/5'],[1,'1/8']); test([4,'1/5'],[-1,'1/70'],[1,'1/99']); test([5,'1/7'],[4,'1/53'],[2,'1/4443']); test([6,'1/8'],[2,'1/57'],[1,'1/239']); test([8,'1/10'],[-1,'1/239'],[-4,'1/515']); test([12,'1/18'],[8,'1/57'],[-5,'1/239']); test([16,'1/21'],[3,'1/239'],[4,'3/1042']); test([22,'1/28'],[2,'1/443'],[-5,'1/1393'],[-10,'1/11018']); test([22,'1/38'],[17,'7/601'],[10,'7/8149']); test([44,'1/57'],[7,'1/239'],[-12,'1/682'],[24,'1/12943']); test([88,'1/172'],[51,'1/239'],[32,'1/682'],[44,'1/5357'],[68,'1/12943']); test([88,'1/172'],[51,'1/239'],[32,'1/682'],[44,'1/5357'],[68,'1/12944']);</lang>
- Output:
OK ([1 1/2] [1 1/3]) OK ([2 1/3] [1 1/7]) OK ([4 1/5] [-1 1/239]) OK ([5 1/7] [2 3/79]) OK ([5 29/278] [7 3/79]) OK ([1 1/2] [1 1/5] [1 1/8]) OK ([4 1/5] [-1 1/70] [1 1/99]) OK ([5 1/7] [4 1/53] [2 1/4443]) OK ([6 1/8] [2 1/57] [1 1/239]) OK ([8 1/10] [-1 1/239] [-4 1/515]) OK ([12 1/18] [8 1/57] [-5 1/239]) OK ([16 1/21] [3 1/239] [4 3/1042]) OK ([22 1/28] [2 1/443] [-5 1/1393] [-10 1/11018]) OK ([22 1/38] [17 7/601] [10 7/8149]) OK ([44 1/57] [7 1/239] [-12 1/682] [24 1/12943]) OK ([88 1/172] [51 1/239] [32 1/682] [44 1/5357] [68 1/12943]) Error ([88 1/172] [51 1/239] [32 1/682] [44 1/5357] [68 1/12944])
Perl 6
The task description requires exact computation. In this solution atan(a/b) produce floating points. Note that the Maxima solution is correct, since Maxima uses a symbolical representation of expressions.
<lang Perl6>use Test; plan *;
is tan(atan(1/2)+atan(1/3)), 1; is tan(2*atan(1/3)+atan(1/7)), 1; is tan(4*atan(1/5)-atan(1/239)), 1; is tan(5*atan(1/7)+2*atan(3/79)), 1; is tan(5*atan(29/278)+7*atan(3/79)), 1; is tan(atan(1/2)+atan(1/5)+atan(1/8)), 1; is tan(4*atan(1/5)-atan(1/70)+atan(1/99)), 1; is tan(5*atan(1/7)+4*atan(1/53)+2*atan(1/4443)), 1; is tan(6*atan(1/8)+2*atan(1/57)+atan(1/239)), 1; is tan(8*atan(1/10)-atan(1/239)-4*atan(1/515)), 1; is tan(12*atan(1/18)+8*atan(1/57)-5*atan(1/239)), 1; is tan(16*atan(1/21)+3*atan(1/239)+4*atan(3/1042)), 1; is tan(22*atan(1/28)+2*atan(1/443)-5*atan(1/1393)-10*atan(1/11018)), 1; is tan(22*atan(1/38)+17*atan(7/601)+10*atan(7/8149)), 1; is tan(44*atan(1/57)+7*atan(1/239)-12*atan(1/682)+24*atan(1/12943)), 1; is tan(88*atan(1/172)+51*atan(1/239)+32*atan(1/682)+44*atan(1/5357)+68*atan(1/12943)), 1; is tan(88*atan(1/172)+51*atan(1/239)+32*atan(1/682)+44*atan(1/5357)+68*atan(1/12944)), 1; </lang>
- Output:
ok 1 - ok 2 - ok 3 - ok 4 - ok 5 - ok 6 - ok 7 - ok 8 - ok 9 - ok 10 - ok 11 - ok 12 - ok 13 - ok 14 - ok 15 - ok 16 - not ok 17 - # got: '0.999999188225744' # expected: '1'
The text input is almost directly parsable as Perl6 code, so we can just write: <lang Perl 6>use Test; plan *;
for lines() {
is tan(.eval), /INCORRECT/ ?? none(1) !! 1 given .subst(/'pi/4 = '/, )\ .subst(/arctan/, 'atan', :g)\ .subst(/'[INCORRECT]'/, '# INCORRECT')
}</lang>
This program reads stdin.
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)
Note: the Kodos tool was used in developing the regular expression.
Racket
<lang racket>
- lang racket
(define (reduce e)
(match e [(? number? a) a] [(list '+ (? number? a) (? number? b)) (+ a b)] [(list '- (? number? a) (? number? b)) (- a b)] [(list '- (? number? a)) (- a)] [(list '* (? number? a) (? number? b)) (* a b)] [(list '/ (? number? a) (? number? b)) (/ a b)] [(list '+ a b) (reduce `(+ ,(reduce a) ,(reduce b)))] [(list '- a b) (reduce `(- ,(reduce a) ,(reduce b)))] [(list '- a) (reduce `(- ,(reduce a)))] [(list '* a b) (reduce `(* ,(reduce a) ,(reduce b)))] [(list '/ a b) (reduce `(/ ,(reduce a) ,(reduce b)))] [(list 'tan (list 'arctan a)) (reduce a)] [(list 'tan (list '- a)) (reduce `(- ,(reduce `(tan ,a))))] [(list 'tan (list '+ a b)) (reduce `(/ (+ (tan ,a) (tan ,b)) (- 1 (* (tan ,a) (tan ,b)))))] [(list 'tan (list '+ a b c ...)) (reduce `(tan (+ ,a (+ ,b ,@c))))] [(list 'tan (list '- a b)) (reduce `(/ (+ (tan ,a) (tan (- ,b))) (- 1 (* (tan ,a) (tan (- ,b))))))] [(list 'tan (list '* 1 a)) (reduce `(tan ,a))] [(list 'tan (list '* (? number? n) a)) (cond [(< n 0) (reduce `(- (tan (* ,(- n) ,a))))] [(= n 0) 0] [(even? n) (reduce `(tan (+ (* ,(/ n 2) ,a) (* ,(/ n 2) ,a))))] [else (reduce `(tan (+ ,a (* ,(- n 1) ,a))))])]))
(define correct-formulas
'((tan (+ (arctan 1/2) (arctan 1/3))) (tan (+ (* 2 (arctan 1/3)) (arctan 1/7))) (tan (- (* 4 (arctan 1/5)) (arctan 1/239))) (tan (+ (* 5 (arctan 1/7)) (* 2 (arctan 3/79)))) (tan (+ (* 5 (arctan 29/278)) (* 7 (arctan 3/79)))) (tan (+ (arctan 1/2) (arctan 1/5) (arctan 1/8))) (tan (+ (* 4 (arctan 1/5)) (* -1 (arctan 1/70)) (arctan 1/99))) (tan (+ (* 5 (arctan 1/7)) (* 4 (arctan 1/53)) (* 2 (arctan 1/4443)))) (tan (+ (* 6 (arctan 1/8)) (* 2 (arctan 1/57)) (arctan 1/239))) (tan (+ (* 8 (arctan 1/10)) (* -1 (arctan 1/239)) (* -4 (arctan 1/515)))) (tan (+ (* 12 (arctan 1/18)) (* 8 (arctan 1/57)) (* -5 (arctan 1/239)))) (tan (+ (* 16 (arctan 1/21)) (* 3 (arctan 1/239)) (* 4 (arctan 3/1042)))) (tan (+ (* 22 (arctan 1/28)) (* 2 (arctan 1/443)) (* -5 (arctan 1/1393)) (* -10 (arctan 1/11018)))) (tan (+ (* 22 (arctan 1/38)) (* 17 (arctan 7/601)) (* 10 (arctan 7/8149)))) (tan (+ (* 44 (arctan 1/57)) (* 7 (arctan 1/239)) (* -12 (arctan 1/682)) (* 24 (arctan 1/12943)))) (tan (+ (* 88 (arctan 1/172)) (* 51 (arctan 1/239)) (* 32 (arctan 1/682)) (* 44 (arctan 1/5357)) (* 68 (arctan 1/12943))))))
(define wrong-formula
'(tan (+ (* 88 (arctan 1/172)) (* 51 (arctan 1/239)) (* 32 (arctan 1/682)) (* 44 (arctan 1/5357)) (* 68 (arctan 1/12944)))))
(displayln "Do all correct formulas reduce to 1?") (for/and ([f correct-formulas]) (= 1 (reduce f)))
(displayln "The incorrect formula reduces to:") (reduce wrong-formula) </lang> Output: <lang racket> Do all correct formulas reduce to 1?
- t
The incorrect formula reduces to: 1009288018000944050967896710431587186456256928584351786643498522649995492271475761189348270710224618853590682465929080006511691833816436374107451368838065354726517908250456341991684635768915704374493675498637876700129004484434187627909285979251682006538817341793224963346197503893270875008524149334251672855130857035205217929335932890740051319216343365800342290782260673215928499123722781078448297609548233999010983373327601187505623621602789012550584784738082074783523787011976757247516095289966708782862528690942242793667539020699840402353522108223/1009288837315638583415701528780402795721935641614456853534313491853293025565940011104051964874275710024625850092154664245109626053906509780125743180758231049920425664246286578958307532545458843067352531217230461290763258378749459637420702619029075083089762088232401888676895047947363883809724322868121990870409574061477638203859217672620508200713073485398199091153535700094640095900731630771349477187594074169815106104524371099618096164871416282464532355211521113449237814080332335526420331468258917484010722587072087349909684004660371264507984339711 </lang>
REXX
Programming notes:
- REXX doesn't have many math functions, so a few of them are included here.
- The value of π should be at least as accurate as digs.
- If any formula needs any high math functions, they need to be included.
<lang rexx>/*REXX program evaluates some Machin-like formulas and verifies their veracity*/ parse arg digs .; if digs== then digs=100 /*use default for decimal digs?*/ numeric digits digs+10; numeric fuzz 3; pi=pi(); @.=
@.1 = 'pi/4 = atan(1/2) + atan(1/3)' @.2 = 'pi/4 = 2*atan(1/3) + atan(1/7)' @.3 = 'pi/4 = 4*atan(1/5) - atan(1/239)' @.4 = 'pi/4 = 5*atan(1/7) + 2*atan(3/79)' @.5 = 'pi/4 = 5*atan(29/278) + 7*atan(3/79)' @.6 = 'pi/4 = atan(1/2) + atan(1/5) + atan(1/8)' @.7 = 'pi/4 = 4*atan(1/5) - atan(1/70) + atan(1/99)' @.8 = 'pi/4 = 5*atan(1/7) + 4*atan(1/53) + 2*atan(1/4443)' @.9 = 'pi/4 = 6*atan(1/8) + 2*atan(1/57) + atan(1/239)'
@.10 = 'pi/4 = 8*atan(1/10) - atan(1/239) - 4*atan(1/515)' @.11 = 'pi/4 = 12*atan(1/18) + 8*atan(1/57) - 5*atan(1/239)' @.12 = 'pi/4 = 16*atan(1/21) + 3*atan(1/239) + 4*atan(3/1042)' @.13 = 'pi/4 = 22*atan(1/28) + 2*atan(1/443) - 5*atan(1/1393) - 10*atan(1/11018)' @.14 = 'pi/4 = 22*atan(1/38) + 17*atan(7/601) + 10*atan(7/8149)' @.15 = 'pi/4 = 44*atan(1/57) + 7*atan(1/239) - 12*atan(1/682) + 24*atan(1/12943)' @.16 = 'pi/4 = 88*atan(1/172) + 51*atan(1/239) + 32*atan(1/682) + 44*atan(1/5357) + 68*atan(1/12943)' @.17 = 'pi/4 = 88*atan(1/172) + 51*atan(1/239) + 32*atan(1/682) + 44*atan(1/5357) + 68*atan(1/12944)'
do j=1 while @.j\== /*evaluate each "Machin-like" formulas.*/ interpret 'answer=' "(" @.j ")" /*this is the heavy lifting.*/ say right(word('bad OK',answer+1),3)": " space(@.j,0) end /*j*/ /* [↑] show OK or bad, and the formula*/
exit /*stick a fork in it, we're all done. */ /*────broutines───────────────────────────────────────────────────────────────*/ pi: return 3.14159265358979323846264338327950288419716939937510582097494459 ||,
230781640628620899862803482534211706798214808651
AcosErr: call tellErr 'Acos(x), X must be in the range of -1 ──► +1, X='||x AsinErr: call tellErr 'Asin(x), X must be in the range of -1 ──► +1, X='||x tanErr: call tellErr 'tan(' || x") causes division by zero, X=" || x tellErr: say; say '*** error! ***'; say; say arg(1); say; exit 13
Acos: procedure; parse arg x; if x<-1 | x>1 then call AcosErr
return .5*pi()-Asin(x)
Asin: procedure expose $.; parse arg x 1 z 1 o 1 p; a=abs(x); aa=a*a
if a>1 then call AsinErr x /*X argument is out of valid range. */ if a>=sqrt(2)*.5 then return sign(x)*acos(sqrt(1-aa), '-ASIN') do j=2 by 2 until p=z; p=z; o=o*aa*(j-1)/j; z=z+o/(j+1); end return z /* [↑] compute until no more noise. */
Atan: procedure; parse arg x; if abs(x)=1 then return pi()/4*sign(x)
return Asin(x/sqrt(1+x**2))
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); i=; m.=9
numeric digits 9; numeric form; h=d+6; if x<0 then do; x=-x; i='i'; end parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g*.5'e'_%2 do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/ do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/ numeric digits d; return (g/1)i /*make complex if X < 0.*/</lang>
output
OK: pi/4=atan(1/2)+atan(1/3) OK: pi/4=2*atan(1/3)+atan(1/7) OK: pi/4=4*atan(1/5)-atan(1/239) OK: pi/4=5*atan(1/7)+2*atan(3/79) OK: pi/4=5*atan(29/278)+7*atan(3/79) OK: pi/4=atan(1/2)+atan(1/5)+atan(1/8) OK: pi/4=4*atan(1/5)-atan(1/70)+atan(1/99) OK: pi/4=5*atan(1/7)+4*atan(1/53)+2*atan(1/4443) OK: pi/4=6*atan(1/8)+2*atan(1/57)+atan(1/239) OK: pi/4=8*atan(1/10)-atan(1/239)-4*atan(1/515) OK: pi/4=12*atan(1/18)+8*atan(1/57)-5*atan(1/239) OK: pi/4=16*atan(1/21)+3*atan(1/239)+4*atan(3/1042) OK: pi/4=22*atan(1/28)+2*atan(1/443)-5*atan(1/1393)-10*atan(1/11018) OK: pi/4=22*atan(1/38)+17*atan(7/601)+10*atan(7/8149) OK: pi/4=44*atan(1/57)+7*atan(1/239)-12*atan(1/682)+24*atan(1/12943) OK: pi/4=88*atan(1/172)+51*atan(1/239)+32*atan(1/682)+44*atan(1/5357)+68*atan(1/12943) bad: pi/4=88*atan(1/172)+51*atan(1/239)+32*atan(1/682)+44*atan(1/5357)+68*atan(1/12944)
Seed7
<lang seed7>$ include "seed7_05.s7i";
include "bigint.s7i"; include "bigrat.s7i";
const type: mTerms is array array bigInteger;
const array mTerms: testCases is [] (
[] ([] ( 1_, 1_, 2_), [] ( 1_, 1_, 3_)), [] ([] ( 2_, 1_, 3_), [] ( 1_, 1_, 7_)), [] ([] ( 4_, 1_, 5_), [] (-1_, 1_, 239_)), [] ([] ( 5_, 1_, 7_), [] ( 2_, 3_, 79_)), [] ([] ( 1_, 1_, 2_), [] ( 1_, 1_, 5_), [] ( 1_, 1_, 8_)), [] ([] ( 4_, 1_, 5_), [] (-1_, 1_, 70_), [] ( 1_, 1_, 99_)), [] ([] ( 5_, 1_, 7_), [] ( 4_, 1_, 53_), [] ( 2_, 1_, 4443_)), [] ([] ( 6_, 1_, 8_), [] ( 2_, 1_, 57_), [] ( 1_, 1_, 239_)), [] ([] ( 8_, 1_, 10_), [] (-1_, 1_, 239_), [] ( -4_, 1_, 515_)), [] ([] (12_, 1_, 18_), [] ( 8_, 1_, 57_), [] ( -5_, 1_, 239_)), [] ([] (16_, 1_, 21_), [] ( 3_, 1_, 239_), [] ( 4_, 3_, 1042_)), [] ([] (22_, 1_, 28_), [] ( 2_, 1_, 443_), [] ( -5_, 1_, 1393_), [] (-10_, 1_, 11018_)), [] ([] (22_, 1_, 38_), [] (17_, 7_, 601_), [] ( 10_, 7_, 8149_)), [] ([] (44_, 1_, 57_), [] ( 7_, 1_, 239_), [] (-12_, 1_, 682_), [] ( 24_, 1_, 12943_)), [] ([] (88_, 1_, 172_), [] (51_, 1_, 239_), [] ( 32_, 1_, 682_), [] ( 44_, 1_, 5357_), [] (68_, 1_, 12943_)), [] ([] (88_, 1_, 172_), [] (51_, 1_, 239_), [] ( 32_, 1_, 682_), [] ( 44_, 1_, 5357_), [] (68_, 1_, 12944_)) );
const func bigRational: tanEval (in bigInteger: coef, in bigRational: f) is func
result var bigRational: tanEval is bigRational.value; local var bigRational: a is bigRational.value; var bigRational: b is bigRational.value; begin if coef = 1_ then tanEval := f; elsif coef < 0_ then tanEval := -tanEval(-coef, f); else a := tanEval(coef div 2_, f); b := tanEval(coef - coef div 2_, f); tanEval := (a + b) / (1_/1_ - a * b); end if; end func;
const func bigRational: tans (in mTerms: terms) is func
result var bigRational: tans is bigRational.value; local var bigRational: a is bigRational.value; var bigRational: b is bigRational.value; begin if length(terms) = 1 then tans := tanEval(terms[1][1], terms[1][2] / terms[1][3]); else a := tans(terms[.. length(terms) div 2]); b := tans(terms[succ(length(terms) div 2) ..]); tans := (a + b) / (1_/1_ - a * b); end if; end func;
const proc: main is func
local var integer: index is 0; var array bigInteger: term is 0 times 0_; begin for key index range testCases do write(tans(testCases[index]) = 1_/1_ <& ": pi/4 = "); for term range testCases[index] do write([0] ("+", "-")[ord(term[1] < 0_)] <& abs(term[1]) <& "*arctan(" <& term[2] <& "/" <& term[3] <& ")"); end for; writeln; end for; end func;</lang>
- Output:
TRUE: pi/4 = +1*arctan(1/2)+1*arctan(1/3) TRUE: pi/4 = +2*arctan(1/3)+1*arctan(1/7) TRUE: pi/4 = +4*arctan(1/5)-1*arctan(1/239) TRUE: pi/4 = +5*arctan(1/7)+2*arctan(3/79) TRUE: pi/4 = +1*arctan(1/2)+1*arctan(1/5)+1*arctan(1/8) TRUE: pi/4 = +4*arctan(1/5)-1*arctan(1/70)+1*arctan(1/99) TRUE: pi/4 = +5*arctan(1/7)+4*arctan(1/53)+2*arctan(1/4443) TRUE: pi/4 = +6*arctan(1/8)+2*arctan(1/57)+1*arctan(1/239) TRUE: pi/4 = +8*arctan(1/10)-1*arctan(1/239)-4*arctan(1/515) TRUE: pi/4 = +12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) TRUE: pi/4 = +16*arctan(1/21)+3*arctan(1/239)+4*arctan(3/1042) TRUE: pi/4 = +22*arctan(1/28)+2*arctan(1/443)-5*arctan(1/1393)-10*arctan(1/11018) TRUE: pi/4 = +22*arctan(1/38)+17*arctan(7/601)+10*arctan(7/8149) TRUE: pi/4 = +44*arctan(1/57)+7*arctan(1/239)-12*arctan(1/682)+24*arctan(1/12943) TRUE: pi/4 = +88*arctan(1/172)+51*arctan(1/239)+32*arctan(1/682)+44*arctan(1/5357)+68*arctan(1/12943) FALSE: pi/4 = +88*arctan(1/172)+51*arctan(1/239)+32*arctan(1/682)+44*arctan(1/5357)+68*arctan(1/12944)
Tcl
<lang tcl>package require Tcl 8.5
- Compute tan(atan(p)+atan(q)) using rationals
proc tadd {p q} {
lassign $p pp pq lassign $q qp qq set topp [expr {$pp*$qq + $qp*$pq}] set topq [expr {$pq*$qq}] set prodp [expr {$pp*$qp}] set prodq [expr {$pq*$qq}] set lowp [expr {$prodq - $prodp}] set resultp [set gcd1 [expr {$topp * $prodq}]] set resultq [set gcd2 [expr {$topq * $lowp}]] # Critical! Normalize using the GCD while {$gcd2 != 0} {
lassign [list $gcd2 [expr {$gcd1 % $gcd2}]] gcd1 gcd2
} list [expr {$resultp / abs($gcd1)}] [expr {$resultq / abs($gcd1)}]
} proc termTan {n a b} {
if {$n < 0} {
set n [expr {-$n}] set a [expr {-$a}]
} if {$n == 1} {
return [list $a $b]
} set k [expr {$n - [set m [expr {$n / 2}]]*2}] set t2 [termTan $m $a $b] set m2 [tadd $t2 $t2] if {$k == 0} {
return $m2
} return [tadd [termTan $k $a $b] $m2]
} proc machinTan {terms} {
set sum {0 1} foreach term $terms {
set sum [tadd $sum [termTan {*}$term]]
} return $sum
}
- Assumes that the formula is in the very specific form below!
proc parseFormula {formula} {
set RE {(-?\s*\d*\s*\*?)\s*arctan\s*\(\s*(-?\s*\d+)\s*/\s*(-?\s*\d+)\s*\)} set nospace {" " "" "*" ""} foreach {all n a b} [regexp -inline -all $RE $formula] {
if {![regexp {\d} $n]} {append n 1} lappend result [list [string map $nospace $n] [string map $nospace $a] [string map $nospace $b]]
} return $result
}
foreach formula {
"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)"
} {
if {[tcl::mathop::== {*}[machinTan [parseFormula $formula]]]} {
puts "Yes! '$formula' is true"
} else {
puts "No! '$formula' not true"
}
}</lang>
- Output:
Yes! 'pi/4 = arctan(1/2) + arctan(1/3)' is true Yes! 'pi/4 = 2*arctan(1/3) + arctan(1/7)' is true Yes! 'pi/4 = 4*arctan(1/5) - arctan(1/239)' is true Yes! 'pi/4 = 5*arctan(1/7) + 2*arctan(3/79)' is true Yes! 'pi/4 = 5*arctan(29/278) + 7*arctan(3/79)' is true Yes! 'pi/4 = arctan(1/2) + arctan(1/5) + arctan(1/8)' is true Yes! 'pi/4 = 4*arctan(1/5) - arctan(1/70) + arctan(1/99)' is true Yes! 'pi/4 = 5*arctan(1/7) + 4*arctan(1/53) + 2*arctan(1/4443)' is true Yes! 'pi/4 = 6*arctan(1/8) + 2*arctan(1/57) + arctan(1/239)' is true Yes! 'pi/4 = 8*arctan(1/10) - arctan(1/239) - 4*arctan(1/515)' is true Yes! 'pi/4 = 12*arctan(1/18) + 8*arctan(1/57) - 5*arctan(1/239)' is true Yes! 'pi/4 = 16*arctan(1/21) + 3*arctan(1/239) + 4*arctan(3/1042)' is true Yes! 'pi/4 = 22*arctan(1/28) + 2*arctan(1/443) - 5*arctan(1/1393) - 10*arctan(1/11018)' is true Yes! 'pi/4 = 22*arctan(1/38) + 17*arctan(7/601) + 10*arctan(7/8149)' is true Yes! 'pi/4 = 44*arctan(1/57) + 7*arctan(1/239) - 12*arctan(1/682) + 24*arctan(1/12943)' is true Yes! 'pi/4 = 88*arctan(1/172) + 51*arctan(1/239) + 32*arctan(1/682) + 44*arctan(1/5357) + 68*arctan(1/12943)' is true No! 'pi/4 = 88*arctan(1/172) + 51*arctan(1/239) + 32*arctan(1/682) + 44*arctan(1/5357) + 68*arctan(1/12944)' not true
XPL0
<lang XPL0>code ChOut=8, Text=12; \intrinsic routines int Number(18); \numbers from equations def LF=$0A; \ASCII line feed (end-of-line character)
func Parse(S); \Convert numbers in string S to binary in Number array char S; int I, Neg;
proc GetNum; \Get number from string S int N; [while S(0)<^0 ! S(0)>^9 do S:= S+1; N:= S(0)-^0; S:= S+1; while S(0)>=^0 & S(0)<=^9 do [N:= N*10 + S(0) - ^0; S:= S+1]; Number(I):= N; I:= I+1; ];
[while S(0)#^= do S:= S+1; \skip to "=" I:= 0; loop [Neg:= false; \assume positive term
loop [S:= S+1; \next char case S(0) of LF: [Number(I):= 0; return S+1]; \mark end of array ^-: Neg:= true; \term is negative ^a: [Number(I):= 1; I:= I+1; quit] \no coefficient so use 1 other if S(0)>=^0 & S(0)<=^9 then \if digit [S:= S-1; GetNum; quit]; \backup and get number ]; GetNum; \numerator if Neg then Number(I-1):= -Number(I-1); \tan(-a) = -tan(a) GetNum; \denominator ];
];
func GCD(U, V); \Return the greatest common divisor of U and V
int U, V;
int T;
[while V do \Euclid's method
[T:= U; U:= V; V:= rem(T/V)];
return abs(U); ];
proc Verify; \Verify that tangent of equation = 1 (i.e: E = F) int E, F, I, J;
proc Machin(A, B, C, D); int A, B, C, D; int Div; \tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a)*tan(b)) \tan(arctan(A/B) + arctan(C/D)) \ = (tan(arctan(A/B)) + tan(arctan(C/D))) / (1 - tan(arctan(A/B))*tan(arctan(C/D))) \ = (A/B + C/D) / (1 - A/B*C/D) \ = (A*D/B*D + B*C/B*D) / (B*D/B*D - A*C/B*D) \ = (A*D + B*C) / (B*D - A*C) [E:= A*D + B*C; F:= B*D - A*C; Div:= GCD(E, F); \keep integers from getting too big E:= E/Div; F:= F/Div; ];
[E:= 0; F:= 1; I:= 0; while Number(I) do
[for J:= 1 to Number(I) do Machin(E, F, Number(I+1), Number(I+2)); I:= I+3; ];
Text(0, if E=F then "Yes " else "No "); ];
char S, SS; int I;
[S:= "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)
"; \Python version of equations (thanks!)
for I:= 1 to 17 do
[SS:= S; \save start of string line S:= Parse(S); \returns start of next line Verify; \correct Machin equation? Yes or No repeat ChOut(0, SS(0)); SS:= SS+1 until SS(0)=LF; ChOut(0, LF); \show equation ];
]</lang>
- Output:
Yes pi/4 = arctan(1/2) + arctan(1/3) Yes pi/4 = 2*arctan(1/3) + arctan(1/7) Yes pi/4 = 4*arctan(1/5) - arctan(1/239) Yes pi/4 = 5*arctan(1/7) + 2*arctan(3/79) Yes pi/4 = 5*arctan(29/278) + 7*arctan(3/79) Yes pi/4 = arctan(1/2) + arctan(1/5) + arctan(1/8) Yes pi/4 = 4*arctan(1/5) - arctan(1/70) + arctan(1/99) Yes pi/4 = 5*arctan(1/7) + 4*arctan(1/53) + 2*arctan(1/4443) Yes pi/4 = 6*arctan(1/8) + 2*arctan(1/57) + arctan(1/239) Yes pi/4 = 8*arctan(1/10) - arctan(1/239) - 4*arctan(1/515) Yes pi/4 = 12*arctan(1/18) + 8*arctan(1/57) - 5*arctan(1/239) Yes pi/4 = 16*arctan(1/21) + 3*arctan(1/239) + 4*arctan(3/1042) Yes pi/4 = 22*arctan(1/28) + 2*arctan(1/443) - 5*arctan(1/1393) - 10*arctan(1/11018) Yes pi/4 = 22*arctan(1/38) + 17*arctan(7/601) + 10*arctan(7/8149) Yes pi/4 = 44*arctan(1/57) + 7*arctan(1/239) - 12*arctan(1/682) + 24*arctan(1/12943) Yes pi/4 = 88*arctan(1/172) + 51*arctan(1/239) + 32*arctan(1/682) + 44*arctan(1/5357) + 68*arctan(1/12943) No pi/4 = 88*arctan(1/172) + 51*arctan(1/239) + 32*arctan(1/682) + 44*arctan(1/5357) + 68*arctan(1/12944)