check/time_lw_N.cpp

   1 /** @file time_lw_N.cpp
   2  *
   3  *  Test N from the paper "Comparison of Polynomial-Oriented CAS" by Robert H.
   4  *  Lewis and Michael Wester (also known as the smaller version of the first
   5  *  Fermat-test). */
   6
   7 /*
   8  *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
   9  *
  10  *  This program is free software; you can redistribute it and/or modify
  11  *  it under the terms of the GNU General Public License as published by
  12  *  the Free Software Foundation; either version 2 of the License, or
  13  *  (at your option) any later version.
  14  *
  15  *  This program is distributed in the hope that it will be useful,
  16  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  *  GNU General Public License for more details.
  19  *
  20  *  You should have received a copy of the GNU General Public License
  21  *  along with this program; if not, write to the Free Software
  22  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  23  */
  24
  25 #include <iostream>
  26 #include <vector>
  27 #include "ginac.h"
  28 #include "timer.h"
  29 using namespace std;
  30 using namespace GiNaC;
  31
  32 static const bool do_test = false;  // set to true in order to run this beast
  33
  34 static unsigned test()
  35 {
  36         symbol p11("p11"), p12("p12"), p21("p21"), p22("p22");
  37         symbol a12("a12"), a21("a21"), a22("a22");
  38         symbol n11("n11"), n22("n22");
  39         symbol g("g");
  40         symbol q1("q1"), q2("q2"), q3("q3"), q4("q4");
  41
  42         ex ss1 = ex("(4*g*a22^3-g*a12*a21*a22^2-4*n22*a21*a22^2+4*n11*a21*a22^2+7*g*a12*a22^2+4*g^2*a22^2-4*n11*n22*a22^2+4*n11*a22^2+n22*a12*a21^2*a22-n11*a12*a21^2*a22-4*g*a21^2*a22-g*a12^2*a21*a22-5*g^2*a12*a21*a22+5*n11*n22*a12*a21*a22-7*n22*a12*a21*a22+2*n11*a12*a21*a22-4*g*a21*a22+3*g*a12^2*a22+3*g^2*a12*a22-3*n11*n22*a12*a22+3*n11*a12*a22+g*a12*a21^3+g^2*a12^2*a21^2-n11*n22*a12^2*a21^2+n22*a12^2*a21^2-2*g*a12*a21^2-3*g^2*a12^2*a21+3*n11*n22*a12^2*a21-3*n22*a12^2*a21-3*g*a12*a21)/(3*g*a12*a21*a22^2-3*n22*a21*a22^2+g*a12*a22^2-n22*a22^2+3*n11*a12*a21^2*a22-3*g*a21^2*a22+5*g*a12^2*a21*a22-5*n22*a12*a21*a22+4*n11*a12*a21*a22-4*g*a21*a22+g*a12^2*a22-n22*a12*a22+n11*a12*a22-g*a22+2*n11*a12^2*a21^2-2*g*a12*a21^2+2*g*a12^3*a21-2*n22*a12^2*a21+2*n11*a12^2*a21-2*g*a12*a21)",lst(g,a12,a21,a22,n11,n22));
  43
  44         ex ss2 = ex("(4*g*a12*a22^2-4*n22*a22^2+4*a22^2-g*a12^2*a21*a22+n22*a12*a21*a22+4*n11*a12*a21*a22-5*a12*a21*a22-4*g*a21*a22+3*g*a12^2*a22-3*n22*a12*a22+3*a12*a22-n11*a12^2*a21^2+a12^2*a21^2+g*a12*a21^2+3*n11*a12^2*a21-3*a12^2*a21-3*g*a12*a21)/(2*g*a12*a22^2-2*n22*a22^2+g*a12^2*a21*a22-n22*a12*a21*a22+2*n11*a12*a21*a22-2*g*a21*a22+2*g*a12^2*a22-2*n22*a12*a22+2*n11*a12*a22-2*g*a22+n11*a12^2*a21^2-g*a12*a21^2+g*a12^3*a21-n22*a12^2*a21+n11*a12^2*a21-g*a12*a21)",lst(g,a12,a21,a22,n11,n22));
  45
  46         ex ss3 = ex("(4*p21*a22^3-p21*a12*a21*a22^2-4*p22*a21*a22^2+4*p11*a21*a22^2+7*p21*a12*a22^2-4*p11*p22*a22^2+4*p12*p21*a22^2+4*p11*a22^2+p22*a12*a21^2*a22-p11*a12*a21^2*a22-4*p12*a21^2*a22-p21*a12^2*a21*a22+5*p11*p22*a12*a21*a22-7*p22*a12*a21*a22-5*p12*p21*a12*a21*a22+2*p11*a12*a21*a22-4*p12*a21*a22+3*p21*a12^2*a22-3*p11*p22*a12*a22+3*p12*p21*a12*a22+3*p11*a12*a22+p12*a12*a21^3-p11*p22*a12^2*a21^2+p22*a12^2*a21^2+p12*p21*a12^2*a21^2-2*p12*a12*a21^2+3*p11*p22*a12^2*a21-3*p22*a12^2*a21-3*p12*p21*a12^2*a21-3*p12*a12*a21)/(3*p21*a12*a21*a22^2-3*p22*a21*a22^2+p21*a12*a22^2-p22*a22^2+3*p11*a12*a21^2*a22-3*p12*a21^2*a22+5*p21*a12^2*a21*a22-5*p22*a12*a21*a22+4*p11*a12*a21*a22-4*p12*a21*a22+p21*a12^2*a22-p22*a12*a22+p11*a12*a22-p12*a22+2*p11*a12^2*a21^2-2*p12*a12*a21^2+2*p21*a12^3*a21-2*p22*a12^2*a21+2*p11*a12^2*a21-2*p12*a12*a21)",lst(a12,a21,a22,p11,p12,p21,p22));
  47
  48         ex ss4 = ex("(4*p21*a12*a22^2-4*p22*a22^2+4*a22^2-p21*a12^2*a21*a22+p22*a12*a21*a22+4*p11*a12*a21*a22-5*a12*a21*a22-4*p12*a21*a22+3*p21*a12^2*a22-3*p22*a12*a22+3*a12*a22-p11*a12^2*a21^2+a12^2*a21^2+p12*a12*a21^2+3*p11*a12^2*a21-3*a12^2*a21-3*p12*a12*a21)/(2*p21*a12*a22^2-2*p22*a22^2+p21*a12^2*a21*a22-p22*a12*a21*a22+2*p11*a12*a21*a22-2*p12*a21*a22+2*p21*a12^2*a22-2*p22*a12*a22+2*p11*a12*a22-2*p12*a22+p11*a12^2*a21^2-p12*a12*a21^2+p21*a12^3*a21-p22*a12^2*a21+p11*a12^2*a21-p12*a12*a21)",lst(p11,p12,p21,p22,a12,a21,a22));
  49
  50         ex res1 = ex("p11*p22*q1^2*q4^2-p12*p21*q1^2*q4^2-n22*p11*p22*q1*q4^2+2*n11*p11*p22*q1*q4^2-p11*p22*q1*q4^2+n22*p12*p21*q1*q4^2-2*n11*p12*p21*q1*q4^2+p12*p21*q1*q4^2+2*g^2*p11*p22*q4^2-2*n11*n22*p11*p22*q4^2+2*n22*p11*p22*q4^2+2*n11*p11*p22*q4^2-2*p11*p22*q4^2-2*g^2*p12*p21*q4^2+2*n11*n22*p12*p21*q4^2-2*n22*p12*p21*q4^2-2*n11*p12*p21*q4^2+2*p12*p21*q4^2-n11*p22*q1*q2*q3*q4+g*p21*q1*q2*q3*q4+g*p12*q1*q2*q3*q4-n22*p11*q1*q2*q3*q4-g^2*p22*q2*q3*q4+n11*n22*p22*q2*q3*q4-n11*p22*q2*q3*q4-2*g*p21*q2*q3*q4+g*p12*q2*q3*q4+2*g^2*p11*q2*q3*q4-2*n11*n22*p11*q2*q3*q4+2*n22*p11*q2*q3*q4-n11*p22*q1*q3*q4+p22*q1*q3*q4+g*p21*q1*q3*q4-2*g*p12*q1*q3*q4+2*n22*p11*q1*q3*q4-2*p11*q1*q3*q4-g^2*p22*q3*q4+n11*n22*p22*q3*q4-n22*p22*q3*q4-n11*p22*q3*q4+p22*q3*q4-4*g^2*p11*q3*q4+4*n11*n22*p11*q3*q4-4*n22*p11*q3*q4-4*n11*p11*q3*q4+4*p11*q3*q4+n22*p11*p22*q1*q2*q4-2*n11*p11*p22*q1*q2*q4+2*n11*p22*q1*q2*q4-n22*p12*p21*q1*q2*q4+2*n11*p12*p21*q1*q2*q4+g*p21*q1*q2*q4-2*g*p12*q1*q2*q4-n22*p11*q1*q2*q4-4*g^2*p11*p22*q2*q4+4*n11*n22*p11*p22*q2*q4-2*n22*p11*p22*q2*q4-2*n11*p11*p22*q2*q4+2*g^2*p22*q2*q4-2*n11*n22*p22*q2*q4+2*n11*p22*q2*q4+4*g^2*p12*p21*q2*q4-4*n11*n22*p12*p21*q2*q4+2*n22*p12*p21*q2*q4+2*n11*p12*p21*q2*q4-2*g*p21*q2*q4-2*g*p12*q2*q4+2*g^2*p11*q2*q4-2*n11*n22*p11*q2*q4+2*n22*p11*q2*q4-p11*p22*q1^2*q4-p22*q1^2*q4+p12*p21*q1^2*q4+2*p11*q1^2*q4-n22*p11*p22*q1*q4-4*n11*p11*p22*q1*q4+5*p11*p22*q1*q4+n22*p22*q1*q4-p22*q1*q4+n22*p12*p21*q1*q4+4*n11*p12*p21*q1*q4-5*p12*p21*q1*q4+g*p21*q1*q4+4*g*p12*q1*q4+4*n11*p11*q1*q4-4*p11*q1*q4-g^2*q2^2*q3^2+n11*n22*q2^2*q3^2+g^2*q2*q3^2-n11*n22*q2*q3^2-n22*q2*q3^2+2*n11*q2*q3^2+2*g^2*q3^2-2*n11*n22*q3^2+2*n22*q3^2+2*n11*q3^2-2*q3^2+g^2*p22*q2^2*q3-n11*n22*p22*q2^2*q3-2*g^2*p11*q2^2*q3+2*n11*n22*p11*q2^2*q3+g^2*q2^2*q3-n11*n22*q2^2*q3+2*n11*p22*q1*q2*q3-2*g*p21*q1*q2*q3+g*p12*q1*q2*q3-n22*p11*q1*q2*q3+n22*q1*q2*q3-2*n11*q1*q2*q3+g^2*p22*q2*q3-n11*n22*p22*q2*q3+n22*p22*q2*q3+4*g*p21*q2*q3+g*p12*q2*q3+4*g^2*p11*q2*q3-4*n11*n22*p11*q2*q3+4*n11*p11*q2*q3-5*g^2*q2*q3+5*n11*n22*q2*q3-n22*q2*q3-4*n11*q2*q3+2*n11*p22*q1*q3-2*p22*q1*q3-2*g*p21*q1*q3-2*g*p12*q1*q3+2*n22*p11*q1*q3-2*p11*q1*q3-2*n22*q1*q3-2*n11*q1*q3+4*q1*q3+2*g^2*p11*p22*q2^2-2*n11*n22*p11*p22*q2^2-2*g^2*p22*q2^2+2*n11*n22*p22*q2^2-2*g^2*p12*p21*q2^2+2*n11*n22*p12*p21*q2^2-2*g^2*p11*q2^2+2*n11*n22*p11*q2^2+2*g^2*q2^2-2*n11*n22*q2^2+n22*p11*p22*q1*q2+4*n11*p11*p22*q1*q2-n22*p22*q1*q2-4*n11*p22*q1*q2-n22*p12*p21*q1*q2-4*n11*p12*p21*q1*q2-n22*p11*q1*q2-4*n11*p11*q1*q2+n22*q1*q2+4*n11*q1*q2-2*p11*p22*q1^2+2*p22*q1^2+2*p12*p21*q1^2+2*p11*q1^2-2*q1^2",lst(p11,p12,p21,p22,n11,n22,g,q1,q2,q3,q4));
  51         ex result = res1.subs(lst(q1==ss1, q2==ss2, q3==ss3, q4==ss4));
  52         ex normalresult = normal(result);
  53         if (!normalresult.is_zero()) {
  54                 clog << "Normalization should have returned 0." << endl;
  55                 return 1;
  56         }
  57         return 0;
  58 }
  59
  60 unsigned time_lw_N()
  61 {
  62         unsigned result = 0;
  63         unsigned count = 0;
  64         timer tag_heuer;
  65         double time = .0;
  66
  67         cout << "timing Lewis-Wester test N (poly at rational fcns)" << flush;
  68
  69         if (do_test) {
  70                 tag_heuer.start();
  71                 // correct for very small times:
  72                 do {
  73                         result = test();
  74                         ++count;
  75                 } while ((time=tag_heuer.read())<0.1 && !result);
  76                 cout << '.' << flush;
  77                 cout << time/count << 's' << endl;
  78         } else {
  79                 cout << " disabled" << endl;
  80         }
  81
  82         return result;
  83 }
  84
  85 extern void randomify_symbol_serials();
  86
  87 int main(int argc, char** argv)
  88 {
  89         randomify_symbol_serials();
  90         cout << setprecision(2) << showpoint;
  91         return time_lw_N();
  92 }