root/trunk/mathextra/testconsistency.d

/**
 Test consistency of all of the statistical functions which have inverses.
 These tests are really tough -- most will fail on machines which don't have
 80-bit reals.
*/
module mathextra.testconsistency;

import std.stdio;
import mathextra.realtest;

import mathextra.mathstat;
import mathextra.mathspecial;


version (testetcgamma) {
import mathextra.etcgamma;
alias mathextra.etcgamma.tgamma gamma;
alias mathextra.etcgamma.lgamma loggamma;
} else {
import std.math; // for lgamma, tgamma
alias std.math.tgamma gamma;
alias std.math.lgamma loggamma;
}

import mathextra.complex;
// do this just so it will compile.
 real cosh(real x) { return std.math.cosh(x); }
 real sinh(real x) { return std.math.sinh(x); }
 real sin(real x) { return std.math.sin(x); }
 real cos(real x) { return std.math.cos(x); }
 real exp(real x) { return std.math.exp(x); }
 real log(real x) { return std.math.log(x); }
 real acos(real x) { return std.math.acos(x); }
 real asin(real x) { return std.math.asin(x); }
 real atan(real x) { return std.math.atan(x); }

// -------------------------------
// Tests of standard mathematical functions
// -------------------------------

unittest { // lgamma and tgamma

    writefln("CONSISTENCY TESTS OF STANDARD FUNCTIONS");
    
    static real logabsgamma(real x)
    {
     // For poles, gamma(x) returns nan, but loggamma() returns infinity.
      if (x<0 && x==floor(x)) return real.infinity;
      return log(fabs(gamma(x)));
    }

  // Check that loggamma(x) = log(gamma(x)) provided x is not -1, -2, -3, ...

// Test some high-precision values (50 decimal digits)
    const real SQRT_PI = 1.77245385090551602729816748334114518279754945612238L;
    assert(feqrel(gamma(0.5L), SQRT_PI) == real.mant_dig); 

//  assert(consistencyTwoFuncs(&loggamma, &logabsgamma, -1000, 1700) > real.mant_dig-7);
  assert(
    consistencyTwoFuncs(
        (real x) { return exp(loggamma(x)); },
        (real x) { if (x<0 && x==floor(x)) return real.infinity;
                   return fabs(gamma(x)); },
       -2000, real.infinity)
       > real.mant_dig-16);
}

// x*exp(x). Used only for testing
real productExp(real x)
{
 return x*exp(x);
}

unittest { // productLog, productLog_1
 assert(isNegZero(productLog(-0.0)));
 assert(isPosZero(productLog(0.0)));
 assert( consistencyRealInverse(&productLog, &productExp, 
        -1/E, double.max) >= double.mant_dig);
 assert( consistencyRealInverse(&productExp, &productLog,-1, 11000)
        >= real.mant_dig-3);
/+        
 assert( consistencyRealInverse(&productLog_1, &productExp, 
        -1/E, -real.epsilon) >= double.mant_dig);
 assert( consistencyRealInverse(&productExp, &productLog_1, 
    log(real.epsilon), -1) >= double.mant_dig);
+/    
}

unittest { // cosh, acosh
  assert( consistencyRealUnity( (real x){ return cosh(acosh(x)); }, 1, double.max)
             >= double.mant_dig);
  assert( consistencyRealUnity( (real x){ return acosh(cosh(x)); }, 1, log(real.max)*(1-real.epsilon))
             >= real.mant_dig-1);
}

unittest { // sinh, asinh
     assert( consistencyRealUnity((real x){ return sinh(asinh(x)); }, -double.max, double.max) >= double.mant_dig);
     assert( consistencyRealUnity((real x){ return sinh(asinh(x)); },-1, 1) >= real.mant_dig-2);
     assert( consistencyRealUnity((real x){ return asinh(sinh(x));},
             -log(real.max)*(1+real.epsilon), log(real.max)*(1-real.epsilon)
             )>= double.mant_dig);
}

unittest { // tanh, atanh
 assert( consistencyRealInverse(&atanh, &tanh, -1, 1)
         >= real.mant_dig-2);
 assert( consistencyRealInverse(&tanh, &atanh, -1,1)
         >= real.mant_dig-2);        
}


unittest { // tan, atan
 assert( consistencyRealInverse(&atan, &tan, -PI_2, PI_2)
         >= real.mant_dig-1);
 assert( consistencyRealInverse(&tan, &atan, -PI_2*(1-real.epsilon), PI_2*(1-real.epsilon))
         >= real.mant_dig-2);
}
/+
// Sin and cos require absolute accuracy tests.

unittest { // cos, acos
writefln("cos");
//writefln(-PI_2, "  ", cos(-PI_2), "  ", acos(cos(PI_2*2)));
writefln(-PI_2, "  %a", cos(0x1p-63), "  ", cos(acos(PI)));
//writefln(-PI_2, "  ", acos(-2), "  ", acos(cos(-1)));
 assert( consistencyRealInverse(&acos, &cos, 0x1p-64, 1)
         >= real.mant_dig-1);
         
/*         
 assert( consistencyRealInverse(&cos, &acos, -PI_2*(1-real.epsilon), PI_2*(1-real.epsilon))
         >= real.mant_dig-2);
*/         
}
+/