0
$\begingroup$

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}]
Improve this question
$\endgroup$
7
  • $\begingroup$ 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). $\endgroup$
    – Michael E2
    Commented Jan 13, 2021 at 14:46
  • $\begingroup$ I tried but it's still very slow. How can I reduce the number of simplifying? $\endgroup$
    – emnha
    Commented Jan 13, 2021 at 14:58
  • $\begingroup$ With Simplify only two cases remained out of 104 976 so I think Simplify doesn't work. $\endgroup$
    – emnha
    Commented Jan 13, 2021 at 15:04
  • $\begingroup$ How many cases are you supposed to get? $\endgroup$
    – Michael E2
    Commented Jan 13, 2021 at 15:24
  • $\begingroup$ @MichaelE2 I don't know exactly but I expect it would be more than several hundreds. $\endgroup$
    – emnha
    Commented Jan 13, 2021 at 15:25

1 Answer 1

Reset to default
4
$\begingroup$

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  *)

Addendum: Symbolic speed-ups.

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}  *)
Improve this answer
$\endgroup$
16
  • $\begingroup$ What does 1 in {a1, _Integer, 1} mean? Also could you explain a bit pair1, pair2? $\endgroup$
    – emnha
    Commented Jan 15, 2021 at 7:19
  • $\begingroup$ @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. $\endgroup$
    – Michael E2
    Commented Jan 15, 2021 at 15:31
  • $\begingroup$ Thanks. Why do you choose 1 as arraydepth for a1 while a1 is just a number? $\endgroup$
    – emnha
    Commented Jan 16, 2021 at 6:41
  • $\begingroup$ @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). $\endgroup$
    – Michael E2
    Commented Jan 16, 2021 at 6:58
  • $\begingroup$ 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? $\endgroup$
    – emnha
    Commented Jan 18, 2021 at 7:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.