Check Machin-like formulas: Difference between revisions
m (→{{header|J}}: Use more industry-standard terms in notation (thanks User:TobyK)) |
(+ D entry) |
||
Line 30: | Line 30: | ||
You can store the equations in any convenient data structure, but for extra credit parse them from human-readable [[Check_Machin-like_formulas/text_equations|text input]]. |
You can store the equations in any convenient data structure, but for extra credit parse them from human-readable [[Check_Machin-like_formulas/text_equations|text input]]. |
||
=={{header|D}}== |
|||
This uses the module of the Arithmetic Rational Task. |
|||
{{trans|Python}} |
|||
<lang d>import std.stdio, std.regex, std.conv, std.string, std.typecons, |
|||
std.array, std.range; |
|||
import arithmetic_rational: Rational; |
|||
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)"; |
|||
alias Pair = Tuple!(int, Rational); |
|||
Pair[][] parseEquations(in string text) { |
|||
auto r = regex(r"\s*(?P<sign>[+-])?\s*(?:(?P<mul>\d+)\s*\*)?\s*" ~ |
|||
r"arctan\((?P<num>\d+)/(?P<denom>\d+)\)", "g"); |
|||
Pair[][] machins; |
|||
foreach (const line; text.splitLines()) { |
|||
Pair[] formula; |
|||
foreach (part; std.string.split(line, "=")[1].match(r)) { |
|||
immutable string mul = part["mul"]; |
|||
immutable string num = part["num"]; |
|||
immutable string denom = part["denom"]; |
|||
formula ~= Pair((part["sign"] == "-" ? -1 : 1) * |
|||
(mul.empty ? 1 : to!int(mul)), |
|||
Rational(to!int(num), |
|||
denom.empty ? 1 : to!int(denom))); |
|||
} |
|||
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); |
|||
/*const*/ auto a = tanEval(coef / 2, f); |
|||
/*const*/ auto b = tanEval(coef - coef / 2, f); |
|||
return (a + b) / (1 - a * b); |
|||
} |
|||
if (xs.length == 1) |
|||
return tanEval(xs[0].tupleof); |
|||
/*const*/ auto a = tans(xs[0 .. xs.length / 2]); |
|||
/*const*/ auto b = tans(xs[xs.length / 2 .. $]); |
|||
return (a + b) / (1 - a * b); |
|||
} |
|||
void main() { |
|||
/*const*/ auto machins = parseEquations(equationText); |
|||
foreach (machin, eqn; zip(machins, equationText.splitLines())) { |
|||
/*const*/ auto ans = tans(machin); |
|||
writefln("%5s: %s", ans == 1 ? "OK" : "ERROR", eqn); |
|||
} |
|||
}</lang> |
|||
{{out}} |
|||
<pre> 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)</pre> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
Revision as of 17:42, 20 January 2013
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.
D
This uses the module of the Arithmetic Rational Task.
<lang d>import std.stdio, std.regex, std.conv, std.string, std.typecons,
std.array, std.range;
import arithmetic_rational: Rational;
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)";
alias Pair = Tuple!(int, Rational);
Pair[][] parseEquations(in string text) {
auto r = regex(r"\s*(?P<sign>[+-])?\s*(?:(?P<mul>\d+)\s*\*)?\s*" ~ r"arctan\((?P<num>\d+)/(?P<denom>\d+)\)", "g"); Pair[][] machins; foreach (const line; text.splitLines()) { Pair[] formula; foreach (part; std.string.split(line, "=")[1].match(r)) { immutable string mul = part["mul"]; immutable string num = part["num"]; immutable string denom = part["denom"]; formula ~= Pair((part["sign"] == "-" ? -1 : 1) * (mul.empty ? 1 : to!int(mul)), Rational(to!int(num), denom.empty ? 1 : to!int(denom))); } 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); /*const*/ auto a = tanEval(coef / 2, f); /*const*/ auto b = tanEval(coef - coef / 2, f); return (a + b) / (1 - a * b); }
if (xs.length == 1) return tanEval(xs[0].tupleof); /*const*/ auto a = tans(xs[0 .. xs.length / 2]); /*const*/ auto b = tans(xs[xs.length / 2 .. $]); return (a + b) / (1 - a * b);
}
void main() {
/*const*/ auto machins = parseEquations(equationText); foreach (machin, eqn; zip(machins, equationText.splitLines())) { /*const*/ auto ans = tans(machin); 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)
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
<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
Perl 6
<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'
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.
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)