# Improve running speed of FullSimplify function

How can I make this code faster? The checking with FullSimplify takes a lot of time.

ClearAll["Global*"]
myPrint[args__, {style__}] := Print[Row[{args}, BaseStyle -> {style}]]
f1[a1_, a2_, a3_, a4_, a5_, a6_, b1_, b2_, b3_, b4_, b5_, b6_,
d_] := -(((a5 + a2 d - a5 d) (b4 + b1 d - b4 d) - (a4 + a1 d -
a4 d) (b5 + b2 d - b5 d))/(-(a6 + a3 d - a6 d) (b5 + b2 d -
b5 d) + (a5 + a2 d - a5 d) (b6 + b3 d - b6 d)));
f2[a1_, a2_, a3_, a4_, a5_, a6_, b1_, b2_, b3_, b4_, b5_, b6_,
d_] := (a4 b6 - a6 (-1 + d) (b4 (-1 + d) - b1 d) +
d (a4 b3 - a3 b4 + a1 b6 - 2 a4 b6 +
a3 (-b1 + b4) d + (a1 - a4) (b3 - b6) d))/(a6 (-1 +
d) (b5 (-1 + d) - b2 d) - a5 (-1 + d) (b6 (-1 + d) - b3 d) +
d (a2 b6 (-1 + d) - a2 b3 d + a3 (b5 + b2 d - b5 d)));
tup1 = Tuples[{{0, 1}, {-1, 0, 1}, {-1, 0, 1}, {0, 1}, {-1, 0,
1}, {-1, 0, 1}, {0, 1}, {-1, 0, 1}, {-1, 0, 1}, {0, 1}, {-1, 0,
1}, {-1, 0, 1}}];
myPrint["The total number of cases: ",
Length[tup1], {FontSize -> 25, FontWeight -> Bold,
Background -> LightGreen}]
tupn = Pick[tup1,
Quiet[FullSimplify[f2[##, d] == 1] ||
FullSimplify[f2[##, d] == f1[##, d]] ||
FullSimplify[f2[##, d] == (1 - f1[##, d]) ] ||
FullSimplify[f2[##, d] == (f1[##, d] - 1) ] & @@@ tup1]];
myPrint["The number of removed cases: ",
Length[tup1] - Length[tupn], {FontSize -> 25, FontWeight -> Bold,
Background -> LightGreen}]

• Have you tried Simplify? It's probably faster and sufficiently robust on rational functions. Also, you could be simplifying f2 four times and f1 three times, which seems likely to be wasteful (it might not be, but it probably is). Jan 13, 2021 at 14:46
• I tried but it's still very slow. How can I reduce the number of simplifying? Jan 13, 2021 at 14:58
• With Simplify only two cases remained out of 104 976 so I think Simplify doesn't work. Jan 13, 2021 at 15:04
• How many cases are you supposed to get? Jan 13, 2021 at 15:24
• @MichaelE2 I don't know exactly but I expect it would be more than several hundreds. Jan 13, 2021 at 15:25

Instead of filtering symbolically, we can filter numerically to remove a large number of undesired tuples. Since the numerators and denominators are quadratic polynomials, the rational functions are determined by their values at five distinct points. To keep things in machine integers, we write a rational function as an ordered pair and code some of the algebra by hand. (Has someone already written a package that implements rational arithmetic in terms of ordered pairs?)

f12[a1_, a2_, a3_, a4_, a5_, a6_, b1_, b2_, b3_, b4_, b5_, b6_,
d_] := {-((a5 + a2 d - a5 d) (b4 + b1 d - b4 d) - (a4 + a1 d -
a4 d) (b5 + b2 d - b5 d)),  (* numerator *)
(-(a6 + a3 d - a6 d) (b5 + b2 d - b5 d) + (a5 + a2 d - a5 d) (b6 +
b3 d - b6 d))};              (* denominator *)

f22[a1_, a2_, a3_, a4_, a5_, a6_, b1_, b2_, b3_, b4_, b5_, b6_,
d_] := {(a4 b6 - a6 (-1 + d) (b4 (-1 + d) - b1 d) +
d (a4 b3 - a3 b4 + a1 b6 - 2 a4 b6 + a3 (-b1 + b4) d +
(a1 - a4) (b3 - b6) d)),        (* numerator *)
(a6 (-1 + d) (b5 (-1 + d) - b2 d) - a5 (-1 + d) (b6 (-1 + d) - b3 d) +
d (a2 b6 (-1 + d) - a2 b3 d +
a3 (b5 + b2 d - b5 d)))};       (* denominator *)

cf1 = Compile[{{a1, _Integer, 1}, {a2, _Integer, 1}, {a3, _Integer,
1}, {a4, _Integer, 1}, {a5, _Integer, 1}, {a6, _Integer,
1}, {b1, _Integer, 1}, {b2, _Integer, 1}, {b3, _Integer,
1}, {b4, _Integer, 1}, {b5, _Integer, 1}, {b6, _Integer,
1}, {d, _Integer}},
Evaluate@f12[a1, a2, a3, a4, a5, a6, b1, b2, b3, b4, b5, b6, d]];
cf2 = Compile[{{a1, _Integer, 1}, {a2, _Integer, 1}, {a3, _Integer,
1}, {a4, _Integer, 1}, {a5, _Integer, 1}, {a6, _Integer,
1}, {b1, _Integer, 1}, {b2, _Integer, 1}, {b3, _Integer,
1}, {b4, _Integer, 1}, {b5, _Integer, 1}, {b6, _Integer,
1}, {d, _Integer}},
Evaluate@f22[a1, a2, a3, a4, a5, a6, b1, b2, b3, b4, b5, b6, d]];
test = Compile[{{pair1, _Integer, 1}, {pair2, _Integer, 1}},
If[Last@pair1 == 0 || Last@pair2 == 0,
1,
With[{
a1 = pair1[[1]], (* a1/a2 = f1[##, d] *)
a2 = pair1[[2]],
b1 = pair2[[1]], (* b1/b2 = f2[##, d] *)
b2 = pair2[[2]]
},
(b1 - b2) *                (* b1/b2=f2[##,d]=1 *)
(a1 b2 - a2 b1) *         (* b1/b2=f2[##,d]=f1[##,d]=a1/a2 *)
(a1 b2 + a2 b1 - a2 b2) * (* b1/b2=f2[##,d]=1-f1[##,d]=1-a1/a2 *)
(a1 b2 - a2 b1 - a2 b2)   (* b1/b2=f2[##,d]=f1[##,d]-1=a1/a2-1 *)
]],
RuntimeAttributes -> {Listable}, Parallelization -> True];

res = Fold[Function[{tup, d},
Pick[tup,
test[
cf1[##, d] & @@ Transpose[tup] // Transpose,
cf2[##, d] & @@ Transpose[tup] // Transpose],
0
]
],
tup1,
Range[3, 7]
]; // AbsoluteTiming

(*  {0.068052, Null}  *)

Length@res

(*  19544  *)


Check symbolically:

res = Pick[res,
Quiet[FullSimplify[f2[##, d] == 1] ||
FullSimplify[f2[##, d] == f1[##, d]] ||
FullSimplify[f2[##, d] == (1 - f1[##, d])] ||
FullSimplify[f2[##, d] == (f1[##, d] - 1)] & @@@
res]]; // AbsoluteTiming
Length@res

(*  {10.7007, Null}  *)
(*  19544  *)


Simplify is usually as robust as FullSimplify on rational functions. (I think "usually" can be replaced by "always," in fact.)

res = Pick[res,
Quiet[Simplify[f2[##, d] == 1] ||
Simplify[f2[##, d] == f1[##, d]] ||
Simplify[f2[##, d] == (1 - f1[##, d])] ||
Simplify[f2[##, d] == (f1[##, d] - 1)] & @@@
res]]; // AbsoluteTiming

(*  {8.97561, Null}  *)


Nesting Simplify sometimes speeds things up. There's no real way to be sure without testing. If the first Simplify simplifies a lot, then it usually will help. Simplify also caches some results, which can trip up the one who would predict which way is faster.

res = Pick[res,
Quiet[
With[{sf1 = Simplify[f1[##, d]], sf2 = Simplify[f2[##, d]]},
Simplify[sf2 == 1] || Simplify[sf2 == sf1] ||
Simplify[sf2 == (1 - sf1)] || Simplify[sf2 == (sf1 - 1)]
] & @@@ res]]; // AbsoluteTiming

(*  {5.78182, Null}  *)

res = Pick[res,
Quiet[With[{sf2 = Simplify[f2[##, d]]},
Simplify[sf2 == 1] || (* won't go on to simplify f1[] if True *)
With[{sf1 = Simplify[f1[##, d]]},
Simplify[sf2 == sf1] || Simplify[sf2 == (1 - sf1)] ||
Simplify[sf2 == (sf1 - 1)]
]] & @@@ res]]; // AbsoluteTiming

(*  {5.05181, Null}  *)

• What does 1 in {a1, _Integer, 1} mean? Also could you explain a bit pair1, pair2? Jan 15, 2021 at 7:19
• @anhnha In Compile, {var, type, arraydepth} declares the variable var to be an array of type type and depth arraydepth. A list has depth 1, a matrix depth 2, and so on. -- A pair1 should be a list of length two of the form {num, den}, where num/den = f1[..] and similarly for pair2 and f2[]. Thus pair1[[1]] = num and pair[[2]] = den. I added some comment I hope will clarify. Also note that X*Y=0 is equivalent to X=0 or Y=0. So the product is zero if any of the factors are zero. Jan 15, 2021 at 15:31
• Thanks. Why do you choose 1 as arraydepth for a1 while a1 is just a number? Jan 16, 2021 at 6:41
• @anhnha I thought you meant the code for any of the Compile functions. I misunderstood. The functions cf1 and cf2 are called cf1[##, d]&@@Transpose[tup], which means a1 is a list of integers. Since Transpose[tup] is a list of twelve lists, one list for each of a1,..., b6. Each of the twelve lists becomes an argument to cf1. This works because Plus and Times can operate on lists (and it helps here because they are extremely efficient on lists). Jan 16, 2021 at 6:58
• Thanks. The way that functions cf1 and cf2 are called makes sense that their inputs are lists. However, I'm confused this part too. cf1 = Compile[{{a1, _Integer, 1}, {a2, _Integer, 1}, ...] In this function you're declare that individual variable a1, a2` are lists? Jan 18, 2021 at 7:37