check/inifcns_consist.cpp

   1 /** @file inifcns_consist.cpp
   2  *
   3  *  This test routine applies assorted tests on initially known higher level
   4  *  functions. */
   5
   6 /*
   7  *  GiNaC Copyright (C) 1999 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  22  */
  23
  24 #include <ginac/ginac.h>
  25
  26 #ifndef NO_GINAC_NAMESPACE
  27 using namespace GiNaC;
  28 #endif // ndef NO_GINAC_NAMESPACE
  29
  30 /* Simple tests on the sine trigonometric function. */
  31 static unsigned inifcns_consist_sin(void)
  32 {
  33     unsigned result = 0;
  34     bool errorflag;
  35
  36     // sin(n*Pi) == 0?
  37     errorflag = false;
  38     for (int n=-10; n<=10; ++n) {
  39         if ( sin(n*Pi).eval() != numeric(0) ||
  40             !sin(n*Pi).eval().info(info_flags::integer))
  41             errorflag = true;
  42     }
  43     if (errorflag) {
  44         clog << "sin(n*Pi) with integer n does not always return exact 0"
  45              << endl;
  46         ++result;
  47     }
  48
  49     // sin((n+1/2)*Pi) == {+|-}1?
  50     errorflag = false;
  51     for (int n=-10; n<=10; ++n) {
  52         if (! sin((n+numeric(1,2))*Pi).eval().info(info_flags::integer) ||
  53             !(sin((n+numeric(1,2))*Pi).eval() == numeric(1) ||
  54               sin((n+numeric(1,2))*Pi).eval() == numeric(-1)))
  55             errorflag = true;
  56     }
  57     if (errorflag) {
  58         clog << "sin((n+1/2)*Pi) with integer n does not always return exact {+|-}1"
  59              << endl;
  60         ++result;
  61     }
  62
  63     return result;
  64 }
  65
  66 /* Simple tests on the cosine trigonometric function. */
  67 static unsigned inifcns_consist_cos(void)
  68 {
  69     unsigned result = 0;
  70     bool errorflag;
  71
  72     // cos((n+1/2)*Pi) == 0?
  73     errorflag = false;
  74     for (int n=-10; n<=10; ++n) {
  75         if ( cos((n+numeric(1,2))*Pi).eval() != numeric(0) ||
  76             !cos((n+numeric(1,2))*Pi).eval().info(info_flags::integer))
  77             errorflag = true;
  78     }
  79     if (errorflag) {
  80         clog << "cos((n+1/2)*Pi) with integer n does not always return exact 0"
  81              << endl;
  82         ++result;
  83     }
  84
  85     // cos(n*Pi) == 0?
  86     errorflag = false;
  87     for (int n=-10; n<=10; ++n) {
  88         if (! cos(n*Pi).eval().info(info_flags::integer) ||
  89             !(cos(n*Pi).eval() == numeric(1) ||
  90               cos(n*Pi).eval() == numeric(-1)))
  91             errorflag = true;
  92     }
  93     if (errorflag) {
  94         clog << "cos(n*Pi) with integer n does not always return exact {+|-}1"
  95              << endl;
  96         ++result;
  97     }
  98
  99     return result;
 100 }
 101
 102 /* Assorted tests on other transcendental functions. */
 103 static unsigned inifcns_consist_trans(void)
 104 {
 105     unsigned result = 0;
 106     symbol x("x");
 107     ex chk;
 108
 109     chk = asin(1)-acos(0);
 110     if (!chk.is_zero()) {
 111         clog << "asin(1)-acos(0) erroneously returned " << chk
 112              << " instead of 0" << endl;
 113         ++result;
 114     }
 115
 116     // arbitrary check of type sin(f(x)):
 117     chk = pow(sin(acos(x)),2) + pow(sin(asin(x)),2)
 118         - (1+pow(x,2))*pow(sin(atan(x)),2);
 119     if (chk != 1-pow(x,2)) {
 120         clog << "sin(acos(x))^2 + sin(asin(x))^2 - (1+x^2)*sin(atan(x))^2 "
 121              << "erroneously returned " << chk << " instead of 1-x^2" << endl;
 122         ++result;
 123     }
 124
 125     // arbitrary check of type cos(f(x)):
 126     chk = pow(cos(acos(x)),2) + pow(cos(asin(x)),2)
 127         - (1+pow(x,2))*pow(cos(atan(x)),2);
 128     if (!chk.is_zero()) {
 129         clog << "cos(acos(x))^2 + cos(asin(x))^2 - (1+x^2)*cos(atan(x))^2 "
 130              << "erroneously returned " << chk << " instead of 0" << endl;
 131         ++result;
 132     }
 133
 134     // arbitrary check of type tan(f(x)):
 135     chk = tan(acos(x))*tan(asin(x)) - tan(atan(x));
 136     if (chk != 1-x) {
 137         clog << "tan(acos(x))*tan(asin(x)) - tan(atan(x)) "
 138              << "erroneously returned " << chk << " instead of -x+1" << endl;
 139         ++result;
 140     }
 141
 142     // arbitrary check of type sinh(f(x)):
 143     chk = -pow(sinh(acosh(x)),2).expand()*pow(sinh(atanh(x)),2)
 144         - pow(sinh(asinh(x)),2);
 145     if (!chk.is_zero()) {
 146         clog << "expand(-(sinh(acosh(x)))^2)*(sinh(atanh(x))^2) - sinh(asinh(x))^2 "
 147              << "erroneously returned " << chk << " instead of 0" << endl;
 148         ++result;
 149     }
 150
 151     // arbitrary check of type cosh(f(x)):
 152     chk = (pow(cosh(asinh(x)),2) - 2*pow(cosh(acosh(x)),2))
 153         * pow(cosh(atanh(x)),2);
 154     if (chk != 1) {
 155         clog << "(cosh(asinh(x))^2 - 2*cosh(acosh(x))^2) * cosh(atanh(x))^2 "
 156              << "erroneously returned " << chk << " instead of 1" << endl;
 157         ++result;
 158     }
 159
 160     // arbitrary check of type tanh(f(x)):
 161     chk = (pow(tanh(asinh(x)),-2) - pow(tanh(acosh(x)),2)).expand()
 162         * pow(tanh(atanh(x)),2);
 163     if (chk != 2) {
 164         clog << "expand(tanh(acosh(x))^2 - tanh(asinh(x))^(-2)) * tanh(atanh(x))^2 "
 165              << "erroneously returned " << chk << " instead of 2" << endl;
 166         ++result;
 167     }
 168
 169     return result;
 170 }
 171
 172 /* Simple tests on the Gamma function.  We stuff in arguments where the results
 173  * exists in closed form and check if it's ok. */
 174 static unsigned inifcns_consist_gamma(void)
 175 {
 176     unsigned result = 0;
 177     ex e;
 178
 179     e = gamma(ex(1));
 180     for (int i=2; i<8; ++i)
 181         e += gamma(ex(i));
 182     if (e != numeric(874)) {
 183         clog << "gamma(1)+...+gamma(7) erroneously returned "
 184              << e << " instead of 874" << endl;
 185         ++result;
 186     }
 187
 188     e = gamma(ex(1));
 189     for (int i=2; i<8; ++i)
 190         e *= gamma(ex(i));
 191     if (e != numeric(24883200)) {
 192         clog << "gamma(1)*...*gamma(7) erroneously returned "
 193              << e << " instead of 24883200" << endl;
 194         ++result;
 195     }
 196
 197     e = gamma(ex(numeric(5, 2)))*gamma(ex(numeric(9, 2)))*64;
 198     if (e != 315*Pi) {
 199         clog << "64*gamma(5/2)*gamma(9/2) erroneously returned "
 200              << e << " instead of 315*Pi" << endl;
 201         ++result;
 202     }
 203
 204     e = gamma(ex(numeric(-13, 2)));
 205     for (int i=-13; i<7; i=i+2)
 206         e += gamma(ex(numeric(i, 2)));
 207     e = (e*gamma(ex(numeric(15, 2)))*numeric(512));
 208     if (e != numeric(633935)*Pi) {
 209         clog << "512*(gamma(-13/2)+...+gamma(5/2))*gamma(15/2) erroneously returned "
 210              << e << " instead of 633935*Pi" << endl;
 211         ++result;
 212     }
 213
 214     return result;
 215 }
 216
 217 /* Simple tests on the Psi-function (aka polygamma-function).  We stuff in
 218    arguments where the result exists in closed form and check if it's ok. */
 219 static unsigned inifcns_consist_psi(void)
 220 {
 221     unsigned result = 0;
 222     symbol x;
 223     ex e, f;
 224
 225     // We check psi(1) and psi(1/2) implicitly by calculating the curious
 226     // little identity gamma(1)'/gamma(1) - gamma(1/2)'/gamma(1/2) == 2*log(2).
 227     e += (gamma(x).diff(x)/gamma(x)).subs(x==numeric(1));
 228     e -= (gamma(x).diff(x)/gamma(x)).subs(x==numeric(1,2));
 229     if (e!=2*log(2)) {
 230         clog << "gamma(1)'/gamma(1) - gamma(1/2)'/gamma(1/2) erroneously returned "
 231              << e << " instead of 2*log(2)" << endl;
 232         ++result;
 233     }
 234
 235     return result;
 236 }
 237
 238 /* Simple tests on the Riemann Zeta function.  We stuff in arguments where the
 239  * result exists in closed form and check if it's ok.  Of course, this checks
 240  * the Bernoulli numbers as a side effect. */
 241 static unsigned inifcns_consist_zeta(void)
 242 {
 243     unsigned result = 0;
 244     ex e;
 245
 246     for (int i=0; i<13; i+=2)
 247         e += zeta(i)/pow(Pi,i);
 248     if (e!=numeric(-204992279,638512875)) {
 249         clog << "zeta(0) + zeta(2) + ... + zeta(12) erroneously returned "
 250              << e << " instead of -204992279/638512875" << endl;
 251         ++result;
 252     }
 253
 254     e = 0;
 255     for (int i=-1; i>-16; i--)
 256         e += zeta(i);
 257     if (e!=numeric(487871,1633632)) {
 258         clog << "zeta(-1) + zeta(-2) + ... + zeta(-15) erroneously returned "
 259              << e << " instead of 487871/1633632" << endl;
 260         ++result;
 261     }
 262
 263     return result;
 264 }
 265
 266 unsigned inifcns_consist(void)
 267 {
 268     unsigned result = 0;
 269
 270     cout << "checking consistency of symbolic functions..." << flush;
 271     clog << "---------consistency of symbolic functions:" << endl;
 272
 273     result += inifcns_consist_sin();
 274     result += inifcns_consist_cos();
 275     result += inifcns_consist_trans();
 276     result += inifcns_consist_gamma();
 277     result += inifcns_consist_psi();
 278     result += inifcns_consist_zeta();
 279
 280     if (!result) {
 281         cout << " passed ";
 282         clog << "(no output)" << endl;
 283     } else {
 284         cout << " failed ";
 285     }
 286
 287     return result;
 288 }