check/powerlaws.cpp

   1 // check/powerlaws.cpp
   2
   3 /* Tests for power laws.  You shouldn't try to draw much inspiration from
   4  * this code, it is a sanity check rather deeply rooted in GiNaC's classes. */
   5
   6 #include <ginac/ginac.h>
   7
   8 static unsigned powerlaws1(void)
   9 {
  10     // (x^a)^b = x^(a*b)
  11
  12     symbol x("x");
  13     symbol a("a");
  14     symbol b("b");
  15
  16     ex e1=power(power(x,a),b);
  17     if (!(is_ex_exactly_of_type(e1,power) &&
  18           is_ex_exactly_of_type(e1.op(0),power) &&
  19           is_ex_exactly_of_type(e1.op(0).op(0),symbol) &&
  20           is_ex_exactly_of_type(e1.op(0).op(1),symbol) &&
  21           is_ex_exactly_of_type(e1.op(1),symbol) &&
  22           e1.is_equal(power(power(x,a),b)) )) {
  23         clog << "(x^a)^b, x,a,b symbolic wrong" << endl;
  24         clog << "returned: " << e1 << endl;
  25         return 1;
  26     }
  27
  28     ex e2=e1.subs(a==1);
  29     if (!(is_ex_exactly_of_type(e2,power) &&
  30           is_ex_exactly_of_type(e2.op(0),symbol) &&
  31           is_ex_exactly_of_type(e2.op(1),symbol) &&
  32           e2.is_equal(power(x,b)) )) {
  33         clog << "(x^a)^b, x,b symbolic, a==1 wrong" << endl;
  34         clog << "returned: " << e2 << endl;
  35         return 1;
  36     }
  37
  38     ex e3=e1.subs(a==-1);
  39     if (!(is_ex_exactly_of_type(e3,power) &&
  40           is_ex_exactly_of_type(e3.op(0),power) &&
  41           is_ex_exactly_of_type(e3.op(0).op(0),symbol) &&
  42           is_ex_exactly_of_type(e3.op(0).op(1),numeric) &&
  43           is_ex_exactly_of_type(e3.op(1),symbol) &&
  44           e3.is_equal(power(power(x,-1),b)) )) {
  45         clog << "(x^a)^b, x,b symbolic, a==-1 wrong" << endl;
  46         clog << "returned: " << e3 << endl;
  47         return 1;
  48     }
  49
  50     ex e4=e1.subs(lst(a==-1,b==2.5));
  51     if (!(is_ex_exactly_of_type(e4,power) &&
  52           is_ex_exactly_of_type(e4.op(0),power) &&
  53           is_ex_exactly_of_type(e4.op(0).op(0),symbol) &&
  54           is_ex_exactly_of_type(e4.op(0).op(1),numeric) &&
  55           is_ex_exactly_of_type(e4.op(1),numeric) &&
  56           e4.is_equal(power(power(x,-1),2.5)) )) {
  57         clog << "(x^a)^b, x symbolic, a==-1, b==2.5 wrong" << endl;
  58         clog << "returned: " << e4 << endl;
  59         return 1;
  60     }
  61
  62     ex e5=e1.subs(lst(a==-0.9,b==2.5));
  63     if (!(is_ex_exactly_of_type(e5,power) &&
  64           is_ex_exactly_of_type(e5.op(0),symbol) &&
  65           is_ex_exactly_of_type(e5.op(1),numeric) &&
  66           e5.is_equal(power(x,numeric(-0.9)*numeric(2.5))) )) {
  67         clog << "(x^a)^b, x symbolic, a==-0.9, b==2.5 wrong" << endl;
  68         clog << "returned: " << e5 << endl;
  69         return 1;
  70     }
  71
  72     ex e6=e1.subs(lst(a==numeric(3)+numeric(5.3)*I,b==-5));
  73     if (!(is_ex_exactly_of_type(e6,power) &&
  74           is_ex_exactly_of_type(e6.op(0),symbol) &&
  75           is_ex_exactly_of_type(e6.op(1),numeric) &&
  76           e6.is_equal(power(x,numeric(-15)+numeric(5.3)*numeric(-5)*I)) )) {
  77         clog << "(x^a)^b, x symbolic, a==3+5.3*I, b==-5 wrong" << endl;
  78         clog << "returned: " << e6 << endl;
  79         return 1;
  80     }
  81     return 0;
  82 }
  83
  84 static unsigned powerlaws2(void)
  85 {
  86     // (a*x)^b = a^b * x^b
  87
  88     symbol x("x");
  89     symbol a("a");
  90     symbol b("b");
  91
  92     ex e1=power(a*x,b);
  93     if (!(is_ex_exactly_of_type(e1,power) &&
  94           is_ex_exactly_of_type(e1.op(0),mul) &&
  95           (e1.op(0).nops()==2) &&
  96           is_ex_exactly_of_type(e1.op(0).op(0),symbol) &&
  97           is_ex_exactly_of_type(e1.op(0).op(1),symbol) &&
  98           is_ex_exactly_of_type(e1.op(1),symbol) &&
  99           e1.is_equal(power(a*x,b)) )) {
 100         clog << "(a*x)^b, x,a,b symbolic wrong" << endl;
 101         clog << "returned: " << e1 << endl;
 102         return 1;
 103     }
 104
 105     ex e2=e1.subs(a==3);
 106     if (!(is_ex_exactly_of_type(e2,power) &&
 107           is_ex_exactly_of_type(e2.op(0),mul) &&
 108           (e2.op(0).nops()==2) &&
 109           is_ex_exactly_of_type(e2.op(0).op(0),symbol) &&
 110           is_ex_exactly_of_type(e2.op(0).op(1),numeric) &&
 111           is_ex_exactly_of_type(e2.op(1),symbol) &&
 112           e2.is_equal(power(3*x,b)) )) {
 113         clog << "(a*x)^b, x,b symbolic, a==3 wrong" << endl;
 114         clog << "returned: " << e2 << endl;
 115         return 1;
 116     }
 117
 118     ex e3=e1.subs(b==-3);
 119     if (!(is_ex_exactly_of_type(e3,mul) &&
 120           (e3.nops()==2) &&
 121           is_ex_exactly_of_type(e3.op(0),power) &&
 122           is_ex_exactly_of_type(e3.op(1),power) &&
 123           e3.is_equal(power(a,-3)*power(x,-3)) )) {
 124         clog << "(a*x)^b, x,a symbolic, b==-3 wrong" << endl;
 125         clog << "returned: " << e3 << endl;
 126         return 1;
 127     }
 128
 129     ex e4=e1.subs(b==4.5);
 130     if (!(is_ex_exactly_of_type(e4,power) &&
 131           is_ex_exactly_of_type(e4.op(0),mul) &&
 132           (e4.op(0).nops()==2) &&
 133           is_ex_exactly_of_type(e4.op(0).op(0),symbol) &&
 134           is_ex_exactly_of_type(e4.op(0).op(1),symbol) &&
 135           is_ex_exactly_of_type(e4.op(1),numeric) &&
 136           e4.is_equal(power(a*x,4.5)) )) {
 137         clog << "(a*x)^b, x,a symbolic, b==4.5 wrong" << endl;
 138         clog << "returned: " << e4 << endl;
 139         return 1;
 140     }
 141
 142     ex e5=e1.subs(lst(a==3.2,b==3+numeric(5)*I));
 143     if (!(is_ex_exactly_of_type(e5,mul) &&
 144           (e5.nops()==2) &&
 145           is_ex_exactly_of_type(e5.op(0),power) &&
 146           is_ex_exactly_of_type(e5.op(1),numeric) &&
 147           e5.is_equal(power(x,3+numeric(5)*I)*
 148                       power(numeric(3.2),3+numeric(5)*I)) )) {
 149         clog << "(a*x)^b, x symbolic, a==3.2, b==3+5*I wrong" << endl;
 150         clog << "returned: " << e5 << endl;
 151         return 1;
 152     }
 153
 154     ex e6=e1.subs(lst(a==-3.2,b==3+numeric(5)*I));
 155     if (!(is_ex_exactly_of_type(e6,mul) &&
 156           (e6.nops()==2) &&
 157           is_ex_exactly_of_type(e6.op(0),power) &&
 158           is_ex_exactly_of_type(e6.op(1),numeric) &&
 159           e6.is_equal(power(-x,3+numeric(5)*I)*
 160                       power(numeric(3.2),3+numeric(5)*I)) )) {
 161         clog << "(a*x)^b, x symbolic, a==-3.2, b==3+5*I wrong" << endl;
 162         clog << "returned: " << e6 << endl;
 163         return 1;
 164     }
 165
 166     ex e7=e1.subs(lst(a==3+numeric(5)*I,b==3.2));
 167     if (!(is_ex_exactly_of_type(e7,power) &&
 168           is_ex_exactly_of_type(e7.op(0),mul) &&
 169           (e7.op(0).nops()==2) &&
 170           is_ex_exactly_of_type(e7.op(0).op(0),symbol) &&
 171           is_ex_exactly_of_type(e7.op(0).op(1),numeric) &&
 172           is_ex_exactly_of_type(e7.op(1),numeric) &&
 173           e7.is_equal(power((3+numeric(5)*I)*x,3.2)) )) {
 174         clog << "(a*x)^b, x symbolic, a==3+5*I, b==3.2 wrong" << endl;
 175         clog << "returned: " << e7 << endl;
 176         return 1;
 177     }
 178
 179     return 0;
 180 }
 181
 182 static unsigned powerlaws3(void)
 183 {
 184     // numeric evaluation
 185
 186     ex e1=power(numeric(4),numeric(1)/numeric(2));
 187     if (e1 != 2) {
 188         clog << "4^(1/2) wrongly returned " << e1 << endl;
 189         return 1;
 190     }
 191
 192     ex e2=power(numeric(27),numeric(2)/numeric(3));
 193     if (e2 != 9) {
 194         clog << "27^(2/3) wrongly returned " << e2 << endl;
 195         return 1;
 196     }
 197
 198     ex e3=power(numeric(5),numeric(1)/numeric(2));
 199     if (!(is_ex_exactly_of_type(e3,power) &&
 200           e3.op(0).is_equal(numeric(5)) &&
 201           e3.op(1).is_equal(numeric(1)/numeric(2)) )) {
 202         clog << "5^(1/2) wrongly returned " << e3 << endl;
 203         return 1;
 204     }
 205
 206     ex e4=power(numeric(5),evalf(numeric(1)/numeric(2)));
 207     if (!(is_ex_exactly_of_type(e4,numeric))) {
 208         clog << "5^(0.5) wrongly returned " << e4 << endl;
 209         return 1;
 210     }
 211
 212     ex e5=power(evalf(numeric(5)),numeric(1)/numeric(2));
 213     if (!(is_ex_exactly_of_type(e5,numeric))) {
 214         clog << "5.0^(1/2) wrongly returned " << e5 << endl;
 215         return 1;
 216     }
 217
 218     return 0;
 219 }
 220
 221 static unsigned powerlaws4(void)
 222 {
 223     // test for mul::eval()
 224
 225     symbol a("a");
 226     symbol b("b");
 227     symbol c("c");
 228
 229     ex f1=power(a*b,ex(1)/ex(2));
 230     ex f2=power(a*b,ex(3)/ex(2));
 231     ex f3=c;
 232
 233     exvector v;
 234     v.push_back(f1);
 235     v.push_back(f2);
 236     v.push_back(f3);
 237     ex e1=mul(v);
 238     if (e1!=a*a*b*b*c) {
 239         clog << "(a*b)^(1/2)*(a*b)^(3/2)*c wrongly returned " << e1 << endl;
 240         return 1;
 241     }
 242     return 0;
 243 }
 244
 245 unsigned powerlaws(void)
 246 {
 247     unsigned result = 0;
 248
 249     cout << "checking power laws..." << flush;
 250     clog << "---------power laws:" << endl;
 251
 252     result += powerlaws1();
 253     result += powerlaws2();
 254     result += powerlaws3();
 255     result += powerlaws4();
 256
 257     if (!result) {
 258         cout << " passed ";
 259         clog << "(no output)" << endl;
 260     } else {
 261         cout << " failed ";
 262     }
 263     return result;
 264 }