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-2010 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_various()
 223 {
 224         unsigned result = 0;
 225         symbol n;
 226         ex e;
 227
 228         if ( binomial(n, 0) != 1 ) {
 229                 clog << "ERROR: binomial(n,0) != 1" << endl;
 230                 ++result;
 231         }
 232
 233         return result;
 234 }
 235
 236 unsigned exam_inifcns()
 237 {
 238         unsigned result = 0;
 239
 240         cout << "examining consistency of symbolic functions" << flush;
 241
 242         result += inifcns_consist_trans();  cout << '.' << flush;
 243         result += inifcns_consist_gamma();  cout << '.' << flush;
 244         result += inifcns_consist_psi();  cout << '.' << flush;
 245         result += inifcns_consist_zeta();  cout << '.' << flush;
 246         result += inifcns_consist_various();  cout << '.' << flush;
 247
 248         return result;
 249 }
 250
 251 int main(int argc, char** argv)
 252 {
 253         return exam_inifcns();
 254 }