@@ -185,6 +185,12 @@ is_int v subst egr =
185185 Just x -> snd (properFraction x :: (Int , Double )) == 0
186186 Nothing -> False
187187
188+ is_positive :: Pattern Expr -> RewriteCondition (Maybe Double ) Expr
189+ is_positive v subst egr =
190+ case egr^. _class (unsafeGetSubst v subst). _data of
191+ Just x -> x > 0
192+ Nothing -> False
193+
188194is_sym :: Pattern Expr -> RewriteCondition (Maybe Double ) Expr
189195is_sym v subst egr =
190196 any ((\ case (Sym _) -> True ; _ -> False ) . unNode) (egr^. _class (unsafeGetSubst v subst). _nodes)
@@ -236,10 +242,19 @@ rewrites =
236242
237243 , " x" * (1 / " x" ) := 1 :| is_not_zero " x"
238244
245+ -- In principle, this with distributivity should capture the
246+ -- binomial theorem, but there is some trickiness.
247+ , powP " a" " n" := " a" * powP " a" (" n" - 1 ) :| is_int " n" :| is_positive " n"
248+ ]
249+
250+ -- How can the binomial theorem be represented?
251+ -- Is it really only available for one integer at a time?
252+ ++ [ powP (" a" + " b" ) (NonVariablePattern . Const $ fromIntegral n) := sum [(fromInteger $ n `choose` k) * powP " a" (fromInteger k) * powP " b" (fromInteger $ n - k) | k <- [0 .. n]] | n <- [2 .. 1000 ]] ++
253+
239254 -- It's a bit unclear to me how to determine that high powers
240255 -- can be reduced. Ideally something like:
241256 -- (cos x)^(2n) --> (1-sin^2 x)^n could happen.
242- , powP (cosP " x" ) 2 := 1 - powP (sinP " x" ) 2
257+ [ powP (cosP " x" ) 2 := 1 - powP (sinP " x" ) 2
243258 , powP (coshP " x" ) 2 := 1 + powP (sinhP " x" ) 2
244259 , powP (cnP " x" " k" ) 2 := 1 - powP (snP " x" " k" ) 2
245260 , powP (dnP " x" " k" ) 2 := (1 - powP (snP " x" " k" ) 2 ) / powP " k" 2
@@ -317,7 +332,13 @@ rewrites =
317332 -- Additional ad-hoc: because of negate representations?
318333 , " a" - (fromInteger (- 1 )* " b" ) := " a" + " b"
319334
320- ]
335+ ] where
336+ n `choose` k
337+ | k < 0 || k > n = 0
338+ | k == 0 || k == n = 1
339+ | k == 1 || k == n - 1 = n
340+ | 2 * k > n = n `choose` (n - k)
341+ | otherwise = (n - 1 ) `choose` (k - 1 ) * n `div` k
321342
322343rewrite :: Fix Expr -> Fix Expr
323344rewrite e = fst $ equalitySaturation e rewrites symCost
@@ -410,9 +431,19 @@ symTests = testGroup "Jacobi"
410431 , testCase " i6" $
411432 rewrite (_i (_ln " x" ) " x" ) @?= " x" * (_ln " x" + fromInteger (- 1 ))
412433
434+ -- Trig identities might be a stepping stone to the elliptic.
435+ , testCase " trig add thm: sin(a+b) = sin a * cos b + cos a * sin b" $
436+ rewrite (_sin(" a" + " b" )) @?= _sin " a" * _cos " b" + _cos " a" * _sin " b"
437+
413438 -- TODO: More elliptic function identities may be worthwhile.
414- , testCase " i7" $
415- rewrite (_pow (_dn " x" " k" ) 11 ) @?= _pow (1 - _pow " k" 2 * _pow (_sn " x" " k" ) 2 ) 5 * _dn " x" " k"
439+ , testCase " reduce (dn(x,k))^11 in terms of sn(x,k)" $
440+ rewrite (_pow (_dn " x" " k" ) 11 ) @?= _pow ((1 - _pow (_sn " x" " k" ) 2 ) / _pow " k" 2 ) 5 * _dn " x" " k" -- this should actually not be equal
441+
442+ , testCase " reduce (dn(x,k))^1001 in terms of sn(x,k)" $
443+ rewrite (_pow (_dn " x" " k" ) 1001 ) @?= _pow ((1 - _pow (_sn " x" " k" ) 2 ) / _pow " k" 2 ) 500 * _dn " x" " k"
444+
445+ , testCase " cubic binomial (a+b)^3 = a^3 + 3*a^2*b + 3*a*b^2 + b^3" $
446+ rewrite (_pow (" a" + " b" ) 3 ) @?= _pow " a" 3 + fromInteger 3 * _pow " a" 2 * " b" + fromInteger 3 * " a" * _pow " b" 2 + _pow " b" 3
416447
417448 -- TODO: Require ability to fine tune parameters
418449 -- , testCase "diff_power_harder" $
0 commit comments