check/exam_inifcns.cpp

   1 /** @file exam_inifcns.cpp
   2  *
   3  *  This test routine applies assorted tests on initially known higher level
   4  *  functions. */
   5
   6 /*
   7  *  GiNaC Copyright (C) 1999-2015 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 /* Assorted tests on other transcendental functions. */
  31 static unsigned inifcns_consist_trans()
  32 {
  33         using GiNaC::asin; using GiNaC::acos;
  34         using GiNaC::asinh; using GiNaC::acosh; using GiNaC::atanh;
  35
  36         unsigned result = 0;
  37         symbol x("x");
  38         ex chk;
  39
  40         chk = asin(1)-acos(0);
  41         if (!chk.is_zero()) {
  42                 clog << "asin(1)-acos(0) erroneously returned " << chk
  43                      << " instead of 0" << endl;
  44                 ++result;
  45         }
  46
  47         // arbitrary check of type sin(f(x)):
  48         chk = pow(sin(acos(x)),2) + pow(sin(asin(x)),2)
  49                 - (1+pow(x,2))*pow(sin(atan(x)),2);
  50         if (chk != 1-pow(x,2)) {
  51                 clog << "sin(acos(x))^2 + sin(asin(x))^2 - (1+x^2)*sin(atan(x))^2 "
  52                      << "erroneously returned " << chk << " instead of 1-x^2" << endl;
  53                 ++result;
  54         }
  55
  56         // arbitrary check of type cos(f(x)):
  57         chk = pow(cos(acos(x)),2) + pow(cos(asin(x)),2)
  58                 - (1+pow(x,2))*pow(cos(atan(x)),2);
  59         if (!chk.is_zero()) {
  60                 clog << "cos(acos(x))^2 + cos(asin(x))^2 - (1+x^2)*cos(atan(x))^2 "
  61                      << "erroneously returned " << chk << " instead of 0" << endl;
  62                 ++result;
  63         }
  64
  65         // arbitrary check of type tan(f(x)):
  66         chk = tan(acos(x))*tan(asin(x)) - tan(atan(x));
  67         if (chk != 1-x) {
  68                 clog << "tan(acos(x))*tan(asin(x)) - tan(atan(x)) "
  69                      << "erroneously returned " << chk << " instead of -x+1" << endl;
  70                 ++result;
  71         }
  72
  73         // arbitrary check of type sinh(f(x)):
  74         chk = -pow(sinh(acosh(x)),2).expand()*pow(sinh(atanh(x)),2)
  75                 - pow(sinh(asinh(x)),2);
  76         if (!chk.is_zero()) {
  77                 clog << "expand(-(sinh(acosh(x)))^2)*(sinh(atanh(x))^2) - sinh(asinh(x))^2 "
  78                      << "erroneously returned " << chk << " instead of 0" << endl;
  79                 ++result;
  80         }
  81
  82         // arbitrary check of type cosh(f(x)):
  83         chk = (pow(cosh(asinh(x)),2) - 2*pow(cosh(acosh(x)),2))
  84                 * pow(cosh(atanh(x)),2);
  85         if (chk != 1) {
  86                 clog << "(cosh(asinh(x))^2 - 2*cosh(acosh(x))^2) * cosh(atanh(x))^2 "
  87                      << "erroneously returned " << chk << " instead of 1" << endl;
  88                 ++result;
  89         }
  90
  91         // arbitrary check of type tanh(f(x)):
  92         chk = (pow(tanh(asinh(x)),-2) - pow(tanh(acosh(x)),2)).expand()
  93                 * pow(tanh(atanh(x)),2);
  94         if (chk != 2) {
  95                 clog << "expand(tanh(acosh(x))^2 - tanh(asinh(x))^(-2)) * tanh(atanh(x))^2 "
  96                      << "erroneously returned " << chk << " instead of 2" << endl;
  97                 ++result;
  98         }
  99
 100         // check consistency of log and eta phases:
 101         for (int r1=-1; r1<=1; ++r1) {
 102                 for (int i1=-1; i1<=1; ++i1) {
 103                         ex x1 = r1+I*i1;
 104                         if (x1.is_zero())
 105                                 continue;
 106                         for (int r2=-1; r2<=1; ++r2) {
 107                                 for (int i2=-1; i2<=1; ++i2) {
 108                                         ex x2 = r2+I*i2;
 109                                         if (x2.is_zero())
 110                                                 continue;
 111                                         if (abs(evalf(eta(x1,x2)-log(x1*x2)+log(x1)+log(x2)))>.1e-12) {
 112                                                 clog << "either eta(x,y), log(x), log(y) or log(x*y) is wrong"
 113                                                      << " at x==" << x1 << ", y==" << x2 << endl;
 114                                                 ++result;
 115                                         }
 116                                 }
 117                         }
 118                 }
 119         }
 120
 121         return result;
 122 }
 123
 124 /* Simple tests on the tgamma function.  We stuff in arguments where the results
 125  * exists in closed form and check if it's ok. */
 126 static unsigned inifcns_consist_gamma()
 127 {
 128         using GiNaC::tgamma;
 129         unsigned result = 0;
 130         ex e;
 131
 132         e = tgamma(1);
 133         for (int i=2; i<8; ++i)
 134                 e += tgamma(ex(i));
 135         if (e != numeric(874)) {
 136                 clog << "tgamma(1)+...+tgamma(7) erroneously returned "
 137                      << e << " instead of 874" << endl;
 138                 ++result;
 139         }
 140
 141         e = tgamma(1);
 142         for (int i=2; i<8; ++i)
 143                 e *= tgamma(ex(i));
 144         if (e != numeric(24883200)) {
 145                 clog << "tgamma(1)*...*tgamma(7) erroneously returned "
 146                      << e << " instead of 24883200" << endl;
 147                 ++result;
 148         }
 149
 150         e = tgamma(ex(numeric(5, 2)))*tgamma(ex(numeric(9, 2)))*64;
 151         if (e != 315*Pi) {
 152                 clog << "64*tgamma(5/2)*tgamma(9/2) erroneously returned "
 153                      << e << " instead of 315*Pi" << endl;
 154                 ++result;
 155         }
 156
 157         e = tgamma(ex(numeric(-13, 2)));
 158         for (int i=-13; i<7; i=i+2)
 159                 e += tgamma(ex(numeric(i, 2)));
 160         e = (e*tgamma(ex(numeric(15, 2)))*numeric(512));
 161         if (e != numeric(633935)*Pi) {
 162                 clog << "512*(tgamma(-13/2)+...+tgamma(5/2))*tgamma(15/2) erroneously returned "
 163                      << e << " instead of 633935*Pi" << endl;
 164                 ++result;
 165         }
 166
 167         return result;
 168 }
 169
 170 /* Simple tests on the Psi-function (aka polygamma-function).  We stuff in
 171    arguments where the result exists in closed form and check if it's ok. */
 172 static unsigned inifcns_consist_psi()
 173 {
 174         using GiNaC::log;
 175         using GiNaC::tgamma;
 176
 177         unsigned result = 0;
 178         symbol x;
 179         ex e, f;
 180
 181         // We check psi(1) and psi(1/2) implicitly by calculating the curious
 182         // little identity tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) == 2*log(2).
 183         e += (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1));
 184         e -= (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1,2));
 185         if (e!=2*log(2)) {
 186                 clog << "tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) erroneously returned "
 187                      << e << " instead of 2*log(2)" << endl;
 188                 ++result;
 189         }
 190
 191         return result;
 192 }
 193
 194 /* Simple tests on the Riemann Zeta function.  We stuff in arguments where the
 195  * result exists in closed form and check if it's ok.  Of course, this checks
 196  * the Bernoulli numbers as a side effect. */
 197 static unsigned inifcns_consist_zeta()
 198 {
 199         unsigned result = 0;
 200         ex e;
 201
 202         for (int i=0; i<13; i+=2)
 203                 e += zeta(i)/pow(Pi,i);
 204         if (e!=numeric(-204992279,638512875)) {
 205                 clog << "zeta(0) + zeta(2) + ... + zeta(12) erroneously returned "
 206                      << e << " instead of -204992279/638512875" << endl;
 207                 ++result;
 208         }
 209
 210         e = 0;
 211         for (int i=-1; i>-16; i--)
 212                 e += zeta(i);
 213         if (e!=numeric(487871,1633632)) {
 214                 clog << "zeta(-1) + zeta(-2) + ... + zeta(-15) erroneously returned "
 215                      << e << " instead of 487871/1633632" << endl;
 216                 ++result;
 217         }
 218
 219         return result;
 220 }
 221
 222 static unsigned inifcns_consist_abs()
 223 {
 224         unsigned result = 0;
 225         realsymbol a("a"), b("b"), x("x"), y("y");
 226         possymbol p("p");
 227         symbol z("z");
 228
 229         if (!abs(exp(x+I*y)).eval().is_equal(exp(x)))
 230                 ++result;
 231
 232         if (!abs(pow(p,a+I*b)).eval().is_equal(pow(p,a)))
 233                 ++result;
 234
 235         if (!abs(sqrt(p)).eval().is_equal(sqrt(p)))
 236                 ++result;
 237
 238         if (!abs(-sqrt(p)).eval().is_equal(sqrt(p)))
 239                 ++result;
 240
 241         // also checks that abs(p)=p
 242         if (!abs(pow(p,a+I*b)).eval().is_equal(pow(p,a)))
 243                 ++result;
 244
 245         if (!abs(pow(x+I*y,a)).eval().is_equal(pow(abs(x+I*y),a)))
 246                 ++result;
 247
 248         // it is not necessary a simplification if the following is really evaluated
 249         if (!abs(pow(x+I*y,a+I*b)).eval().is_equal(abs(pow(x+I*y,a+I*b))))
 250                 ++result;
 251
 252         // check expansion of abs
 253         if (!abs(-7*z*a*p).expand(expand_options::expand_transcendental).is_equal(7*abs(z)*abs(a)*p))
 254                 ++result;
 255
 256         if (!abs(z.conjugate()).eval().is_equal(abs(z)))
 257                 ++result;
 258
 259         if (!abs(step(z)).eval().is_equal(step(z)))
 260                 ++result;
 261
 262         if (!abs(p).info(info_flags::positive) || !abs(a).info(info_flags::real))
 263                 ++result;
 264
 265         if (abs(a).info(info_flags::positive) || !abs(a).info(info_flags::real))
 266                 ++result;
 267
 268         if (abs(z).info(info_flags::positive) || !abs(z).info(info_flags::real))
 269                 ++result;
 270
 271         return result;
 272 }
 273
 274 static unsigned inifcns_consist_exp()
 275 {
 276         unsigned result = 0;
 277         symbol a("a"), b("b");
 278
 279         if (!exp(a+b).expand(expand_options::expand_transcendental).is_equal(exp(a)*exp(b)))
 280                 ++result;
 281
 282         // shall not be expanded since the arg is not add
 283         if (!exp(pow(a+b,2)).expand(expand_options::expand_transcendental).is_equal(exp(pow(a+b,2))))
 284                 ++result;
 285
 286         // expand now
 287         if (!exp(pow(a+b,2)).expand(expand_options::expand_function_args | expand_options::expand_transcendental)
 288                 .is_equal(exp(a*a)*exp(b*b)*exp(2*a*b)))
 289                 ++result;
 290
 291         return result;
 292 }
 293
 294 static unsigned inifcns_consist_log()
 295 {
 296         using GiNaC::log;
 297         unsigned result = 0;
 298         symbol z("a"), w("b");
 299         realsymbol a("a"), b("b");
 300         possymbol p("p"), q("q");
 301
 302         // do not expand
 303         if (!log(z*w).expand(expand_options::expand_transcendental).is_equal(log(z*w)))
 304                 ++result;
 305
 306         // do not expand
 307         if (!log(a*b).expand(expand_options::expand_transcendental).is_equal(log(a*b)))
 308                 ++result;
 309
 310         // shall expand
 311         if (!log(p*q).expand(expand_options::expand_transcendental).is_equal(log(p) + log(q)))
 312                 ++result;
 313
 314         // a bit more complicated
 315         ex e1 = log(-7*p*pow(q,3)*a*pow(b,2)*z*w).expand(expand_options::expand_transcendental);
 316         ex e2 = log(7)+log(p)+log(pow(q,3))+log(-z*a*w*pow(b,2));
 317         if (!e1.is_equal(e2))
 318                 ++result;
 319
 320         // shall not do for non-real powers
 321         if (ex(log(pow(p,z))).is_equal(z*log(p)))
 322                 ++result;
 323
 324         // shall not do for non-positive basis
 325         if (ex(log(pow(a,b))).is_equal(b*log(a)))
 326                 ++result;
 327
 328         // infinite recursion log_series
 329         ex e(log(-p));
 330         ex ser = ex_to<pseries>(e.series(z, 1))
 331                 .convert_to_poly(/* no_order = */ true);
 332         if (!ser.is_equal(e)) {
 333                 clog << "series(" << e << ", " << z << "): wrong result" << endl;
 334                 ++result;
 335         }
 336
 337         return result;
 338 }
 339
 340 static unsigned inifcns_consist_various()
 341 {
 342         unsigned result = 0;
 343         symbol n;
 344
 345         if ( binomial(n, 0) != 1 ) {
 346                 clog << "ERROR: binomial(n,0) != 1" << endl;
 347                 ++result;
 348         }
 349
 350         return result;
 351 }
 352
 353 /* Several tests for derivatives */
 354 static unsigned inifcns_consist_derivatives()
 355 {
 356         unsigned result = 0;
 357         symbol z, w;
 358         realsymbol x;
 359         ex e, e1;
 360
 361         e=pow(x,z).conjugate().diff(x);
 362         e1=pow(x,z).conjugate()*z.conjugate()/x;
 363         if (! (e-e1).normal().is_zero() ) {
 364                 clog << "ERROR: pow(x,z).conjugate().diff(x) " << e << " != " << e1 << endl;
 365                 ++result;
 366         }
 367
 368         e=pow(w,z).conjugate().diff(w);
 369         e1=pow(w,z).conjugate()*z.conjugate()/w;
 370         if ( (e-e1).normal().is_zero() ) {
 371                 clog << "ERROR: pow(w,z).conjugate().diff(w) " << e << " = " << e1 << endl;
 372                 ++result;
 373         }
 374
 375         e=atanh(x).imag_part().diff(x);
 376         if (! e.is_zero() ) {
 377                 clog << "ERROR: atanh(x).imag_part().diff(x) " << e << " != 0" << endl;
 378                 ++result;
 379         }
 380
 381         e=atanh(w).imag_part().diff(w);
 382         if ( e.is_zero() ) {
 383                 clog << "ERROR: atanh(w).imag_part().diff(w) " << e << " = 0" << endl;
 384                 ++result;
 385         }
 386
 387         e=atanh(x).real_part().diff(x);
 388         e1=pow(1-x*x,-1);
 389         if (! (e-e1).normal().is_zero() ) {
 390                 clog << "ERROR: atanh(x).real_part().diff(x) " << e << " != " << e1 << endl;
 391                 ++result;
 392         }
 393
 394         e=atanh(w).real_part().diff(w);
 395         e1=pow(1-w*w,-1);
 396         if ( (e-e1).normal().is_zero() ) {
 397                 clog << "ERROR: atanh(w).real_part().diff(w) " << e << " = " << e1 << endl;
 398                 ++result;
 399         }
 400
 401         e=abs(log(z)).diff(z);
 402         e1=(conjugate(log(z))/z+log(z)/conjugate(z))/abs(log(z))/2;
 403         if (! (e-e1).normal().is_zero() ) {
 404                 clog << "ERROR: abs(log(z)).diff(z) " << e << " != " << e1 << endl;
 405                 ++result;
 406         }
 407
 408         e=Order(pow(x,4)).diff(x);
 409         e1=Order(pow(x,3));
 410         if (! (e-e1).normal().is_zero() ) {
 411                 clog << "ERROR: Order(pow(x,4)).diff(x) " << e << " != " << e1 << endl;
 412                 ++result;
 413         }
 414
 415         return result;
 416 }
 417
 418 unsigned exam_inifcns()
 419 {
 420         unsigned result = 0;
 421
 422         cout << "examining consistency of symbolic functions" << flush;
 423
 424         result += inifcns_consist_trans();  cout << '.' << flush;
 425         result += inifcns_consist_gamma();  cout << '.' << flush;
 426         result += inifcns_consist_psi();  cout << '.' << flush;
 427         result += inifcns_consist_zeta();  cout << '.' << flush;
 428         result += inifcns_consist_abs();  cout << '.' << flush;
 429         result += inifcns_consist_exp();  cout << '.' << flush;
 430         result += inifcns_consist_log();  cout << '.' << flush;
 431         result += inifcns_consist_various();  cout << '.' << flush;
 432         result += inifcns_consist_derivatives();  cout << '.' << flush;
 433
 434         return result;
 435 }
 436
 437 int main(int argc, char** argv)
 438 {
 439         return exam_inifcns();
 440 }