check/exam_powerlaws.cpp

   1 /** @file exam_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 /*
   7  *  GiNaC Copyright (C) 1999-2019 Johannes Gutenberg University Mainz, Germany
   8  *
   9  *  This program is free software; you can redistribute it and/or modify
  10  *  it under the terms of the GNU General Public License as published by
  11  *  the Free Software Foundation; either version 2 of the License, or
  12  *  (at your option) any later version.
  13  *
  14  *  This program is distributed in the hope that it will be useful,
  15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17  *  GNU General Public License for more details.
  18  *
  19  *  You should have received a copy of the GNU General Public License
  20  *  along with this program; if not, write to the Free Software
  21  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  22  */
  23
  24 #include "ginac.h"
  25 using namespace GiNaC;
  26
  27 #include <iostream>
  28 using namespace std;
  29
  30 static unsigned exam_powerlaws1()
  31 {
  32         // (x^a)^b = x^(a*b)
  33
  34         symbol x("x");
  35         symbol a("a");
  36         symbol b("b");
  37
  38         ex e1 = power(power(x,a), b);
  39         if (!(is_exactly_a<power>(e1) &&
  40               is_exactly_a<power>(e1.op(0)) &&
  41               is_exactly_a<symbol>(e1.op(0).op(0)) &&
  42               is_exactly_a<symbol>(e1.op(0).op(1)) &&
  43               is_exactly_a<symbol>(e1.op(1)) &&
  44               e1.is_equal(power(power(x,a),b)) )) {
  45                 clog << "(x^a)^b, x,a,b symbolic wrong" << endl;
  46                 clog << "returned: " << e1 << endl;
  47                 return 1;
  48         }
  49
  50         ex e2 = e1.subs(a==1);
  51         if (!(is_exactly_a<power>(e2) &&
  52               is_exactly_a<symbol>(e2.op(0)) &&
  53               is_exactly_a<symbol>(e2.op(1)) &&
  54               e2.is_equal(power(x,b)) )) {
  55                 clog << "(x^a)^b, x,b symbolic, a==1 wrong" << endl;
  56                 clog << "returned: " << e2 << endl;
  57                 return 1;
  58         }
  59
  60         ex e3 = e1.subs(a==-1);
  61         if (!(is_exactly_a<power>(e3) &&
  62               is_exactly_a<power>(e3.op(0)) &&
  63               is_exactly_a<symbol>(e3.op(0).op(0)) &&
  64               is_exactly_a<numeric>(e3.op(0).op(1)) &&
  65               is_exactly_a<symbol>(e3.op(1)) &&
  66               e3.is_equal(power(power(x,-1),b)) )) {
  67                 clog << "(x^a)^b, x,b symbolic, a==-1 wrong" << endl;
  68                 clog << "returned: " << e3 << endl;
  69                 return 1;
  70         }
  71
  72         ex e4 = e1.subs(lst{a==-1, b==-2.5});
  73         if (!(is_exactly_a<power>(e4) &&
  74               is_exactly_a<power>(e4.op(0)) &&
  75               is_exactly_a<symbol>(e4.op(0).op(0)) &&
  76               is_exactly_a<numeric>(e4.op(0).op(1)) &&
  77               is_exactly_a<numeric>(e4.op(1)) &&
  78               e4.is_equal(power(power(x,-1),-2.5)) )) {
  79                 clog << "(x^a)^b, x symbolic, a==-1, b==-2.5 wrong" << endl;
  80                 clog << "returned: " << e4 << endl;
  81                 return 1;
  82         }
  83
  84         ex e5 = e1.subs(lst{a==-0.9, b==2.5});
  85         if (!(is_exactly_a<power>(e5) &&
  86               is_exactly_a<symbol>(e5.op(0)) &&
  87               is_exactly_a<numeric>(e5.op(1)) &&
  88               e5.is_equal(power(x,numeric(-0.9)*numeric(2.5))) )) {
  89                 clog << "(x^a)^b, x symbolic, a==-0.9, b==2.5 wrong" << endl;
  90                 clog << "returned: " << e5 << endl;
  91                 return 1;
  92         }
  93
  94         ex e6 = e1.subs(lst{a==numeric(3)+numeric(5.3)*I, b==-5});
  95         if (!(is_exactly_a<power>(e6) &&
  96               is_exactly_a<symbol>(e6.op(0)) &&
  97               is_exactly_a<numeric>(e6.op(1)) &&
  98               e6.is_equal(power(x,numeric(-15)+numeric(5.3)*numeric(-5)*I)) )) {
  99                 clog << "(x^a)^b, x symbolic, a==3+5.3*I, b==-5 wrong" << endl;
 100                 clog << "returned: " << e6 << endl;
 101                 return 1;
 102         }
 103
 104         return 0;
 105 }
 106
 107 static unsigned exam_powerlaws2()
 108 {
 109         // (a*x)^b = a^b * x^b
 110
 111         symbol x("x");
 112         symbol a("a");
 113         symbol b("b");
 114
 115         ex e1 = power(a*x,b);
 116         if (!(is_exactly_a<power>(e1) &&
 117               is_exactly_a<mul>(e1.op(0)) &&
 118               (e1.op(0).nops()==2) &&
 119               is_exactly_a<symbol>(e1.op(0).op(0)) &&
 120               is_exactly_a<symbol>(e1.op(0).op(1)) &&
 121               is_exactly_a<symbol>(e1.op(1)) &&
 122               e1.is_equal(power(a*x,b)) )) {
 123                 clog << "(a*x)^b, x,a,b symbolic wrong" << endl;
 124                 clog << "returned: " << e1 << endl;
 125                 return 1;
 126         }
 127
 128         ex e2 = e1.subs(a==3);
 129         if (!(is_exactly_a<power>(e2) &&
 130               is_exactly_a<mul>(e2.op(0)) &&
 131               (e2.op(0).nops()==2) &&
 132               is_exactly_a<symbol>(e2.op(0).op(0)) &&
 133               is_exactly_a<numeric>(e2.op(0).op(1)) &&
 134               is_exactly_a<symbol>(e2.op(1)) &&
 135               e2.is_equal(power(3*x,b)) )) {
 136                 clog << "(a*x)^b, x,b symbolic, a==3 wrong" << endl;
 137                 clog << "returned: " << e2 << endl;
 138                 return 1;
 139         }
 140
 141         ex e3 = e1.subs(b==-3);
 142         if (!(is_exactly_a<mul>(e3) &&
 143               (e3.nops()==2) &&
 144               is_exactly_a<power>(e3.op(0)) &&
 145               is_exactly_a<power>(e3.op(1)) &&
 146               e3.is_equal(power(a,-3)*power(x,-3)) )) {
 147                 clog << "(a*x)^b, x,a symbolic, b==-3 wrong" << endl;
 148                 clog << "returned: " << e3 << endl;
 149                 return 1;
 150         }
 151
 152         ex e4 = e1.subs(b==4.5);
 153         if (!(is_exactly_a<power>(e4) &&
 154               is_exactly_a<mul>(e4.op(0)) &&
 155               (e4.op(0).nops()==2) &&
 156               is_exactly_a<symbol>(e4.op(0).op(0)) &&
 157               is_exactly_a<symbol>(e4.op(0).op(1)) &&
 158               is_exactly_a<numeric>(e4.op(1)) &&
 159               e4.is_equal(power(a*x,4.5)) )) {
 160                 clog << "(a*x)^b, x,a symbolic, b==4.5 wrong" << endl;
 161                 clog << "returned: " << e4 << endl;
 162                 return 1;
 163         }
 164
 165         ex e5 = e1.subs(lst{a==3.2, b==3+numeric(5)*I});
 166         if (!(is_exactly_a<mul>(e5) &&
 167               (e5.nops()==2) &&
 168               is_exactly_a<power>(e5.op(0)) &&
 169               is_exactly_a<numeric>(e5.op(1)) &&
 170               e5.is_equal(power(x,3+numeric(5)*I)*
 171                                           power(numeric(3.2),3+numeric(5)*I)) )) {
 172                 clog << "(a*x)^b, x symbolic, a==3.2, b==3+5*I wrong" << endl;
 173                 clog << "returned: " << e5 << endl;
 174                 return 1;
 175         }
 176
 177         ex e6 = e1.subs(lst{a==-3.2, b==3+numeric(5)*I});
 178         if (!(is_exactly_a<mul>(e6) &&
 179               (e6.nops()==2) &&
 180               is_exactly_a<power>(e6.op(0)) &&
 181               is_exactly_a<numeric>(e6.op(1)) &&
 182               e6.is_equal(power(-x,3+numeric(5)*I)*
 183                                           power(numeric(3.2),3+numeric(5)*I)) )) {
 184                 clog << "(a*x)^b, x symbolic, a==-3.2, b==3+5*I wrong" << endl;
 185                 clog << "returned: " << e6 << endl;
 186                 return 1;
 187         }
 188
 189         ex e7 = e1.subs(lst{a==3+numeric(5)*I, b==3.2});
 190         if (!(is_exactly_a<power>(e7) &&
 191               is_exactly_a<mul>(e7.op(0)) &&
 192               (e7.op(0).nops()==2) &&
 193               is_exactly_a<symbol>(e7.op(0).op(0)) &&
 194               is_exactly_a<numeric>(e7.op(0).op(1)) &&
 195               is_exactly_a<numeric>(e7.op(1)) &&
 196               e7.is_equal(power((3+numeric(5)*I)*x,3.2)) )) {
 197                 clog << "(a*x)^b, x symbolic, a==3+5*I, b==3.2 wrong" << endl;
 198                 clog << "returned: " << e7 << endl;
 199                 return 1;
 200         }
 201
 202         return 0;
 203 }
 204
 205 static unsigned exam_powerlaws3()
 206 {
 207         // numeric evaluation
 208
 209         ex e1 = power(numeric(4),numeric(1,2));
 210         if (e1 != 2) {
 211                 clog << "4^(1/2) wrongly returned " << e1 << endl;
 212                 return 1;
 213         }
 214
 215         ex e2 = power(numeric(27),numeric(2,3));
 216         if (e2 != 9) {
 217                 clog << "27^(2/3) wrongly returned " << e2 << endl;
 218                 return 1;
 219         }
 220
 221         ex e3 = power(numeric(5),numeric(1,2));
 222         if (!(is_exactly_a<power>(e3) &&
 223               e3.op(0).is_equal(numeric(5)) &&
 224               e3.op(1).is_equal(numeric(1,2)))) {
 225                 clog << "5^(1/2) wrongly returned " << e3 << endl;
 226                 return 1;
 227         }
 228
 229         ex e4 = power(numeric(5),evalf(numeric(1,2)));
 230         if (!(is_exactly_a<numeric>(e4))) {
 231                 clog << "5^(0.5) wrongly returned " << e4 << endl;
 232                 return 1;
 233         }
 234
 235         ex e5 = power(evalf(numeric(5)),numeric(1,2));
 236         if (!(is_exactly_a<numeric>(e5))) {
 237                 clog << "5.0^(1/2) wrongly returned " << e5 << endl;
 238                 return 1;
 239         }
 240
 241         return 0;
 242 }
 243
 244 static unsigned exam_powerlaws4()
 245 {
 246         // test for mul::eval()
 247
 248         symbol a("a");
 249         symbol b("b");
 250         symbol c("c");
 251
 252         ex f1 = power(a*b,ex(1)/ex(2));
 253         ex f2 = power(a*b,ex(3)/ex(2));
 254         ex f3 = c;
 255
 256         exvector v;
 257         v.push_back(f1);
 258         v.push_back(f2);
 259         v.push_back(f3);
 260         ex e1 = mul(v);
 261         if (e1!=a*a*b*b*c) {
 262                 clog << "(a*b)^(1/2)*(a*b)^(3/2)*c wrongly returned " << e1 << endl;
 263                 return 1;
 264         }
 265
 266         return 0;
 267 }
 268
 269 static unsigned exam_powerlaws5()
 270 {
 271         // cabinet of slightly pathological cases
 272
 273         symbol a("a");
 274
 275         ex e1 = pow(1,a);
 276         if (e1 != 1) {
 277                 clog << "1^a wrongly returned " << e1 << endl;
 278                 return 1;
 279         }
 280
 281         ex e2 = pow(0,a);
 282         if (!(is_exactly_a<power>(e2))) {
 283                 clog << "0^a was evaluated to " << e2
 284                      << " though nothing is known about a." << endl;
 285                 return 1;
 286         }
 287
 288         return 0;
 289 }
 290
 291 static unsigned exam_powerlaws6()
 292 {
 293         // check expansion rules for positive symbols
 294
 295         symbol a("a");
 296         symbol b("b");
 297         symbol c("c");
 298         realsymbol x("x");
 299         realsymbol y("y");
 300         possymbol p("p");
 301         possymbol q("q");
 302         numeric half=numeric(1,2);
 303
 304         ex e1 = pow(5*pow(3*a*b*x*y*p*q,2),7*half*c).expand();
 305         ex e2 = pow(p,7*c)*pow(q,7*c)*pow(pow(a*b*x*y,2),numeric(7,2)*c)*pow(45,numeric(7,2)*c);
 306         if (!e1.is_equal(e2)) {
 307                 clog << "Could not expand exponents with positive bases in " << e1 << endl;
 308                 return 1;
 309         }
 310
 311         ex e3 = pow(-pow(-a*x*p,3)*pow(b*y*p,3),half*c).expand().normal();
 312         ex e4 = pow(p,3*c)*pow(pow(a*b*x*y,3),half*c);
 313
 314         if (!e3.is_equal(e4)) {
 315                 clog << "Could not expand exponents with positive bases in " << e3 << endl;
 316                 return 1;
 317         }
 318
 319         return 0;
 320 }
 321
 322 unsigned exam_powerlaws()
 323 {
 324         unsigned result = 0;
 325
 326         cout << "examining power laws" << flush;
 327
 328         result += exam_powerlaws1();  cout << '.' << flush;
 329         result += exam_powerlaws2();  cout << '.' << flush;
 330         result += exam_powerlaws3();  cout << '.' << flush;
 331         result += exam_powerlaws4();  cout << '.' << flush;
 332         result += exam_powerlaws5();  cout << '.' << flush;
 333         result += exam_powerlaws6();  cout << '.' << flush;
 334
 335         return result;
 336 }
 337
 338 int main(int argc, char** argv)
 339 {
 340         return exam_powerlaws();
 341 }