/* Numerical testers - all confirmed with Matlab 6r12.0 */

f():=[];  //[//Echo(CurrentLine()+1);];

BuiltinPrecisionSet(10); //Echo("BIP set to ",10);

NumericEqual(N(Sqrt(2),6), 1.41421,6);f();
NumericEqual(N(N(1+Pi,20)-Pi,20),1,20);f();		// "N" should have "HoldArgument" in some way, so inner "N" is evaluated with outer precision 20

/* Got the first digits of Pi from the following page:
   http://www.cecm.sfu.ca/projects/ISC/dataB/isc/C/pi10000.txt
  Just checking that Yacas agrees.
  First, however, we need to set BuiltinPrecision way higher than 10 !!
*/
BuiltinPrecisionSet(90);//Echo("BIP set to ",90);
NumericEqual(N(Pi,70),3.141592653589793238462643383279502884197169399375105820974944592307816,70);f();
BuiltinPrecisionSet(10);//Echo("BIP set to ",10);

NumericEqual( N(Sec(2),10), -2.402997962, 9);f();
NumericEqual( N(Csc(2),9),  1.09975017,9);f();
NumericEqual( N(Cot(2),9), -0.457657554, 9);f();
NumericEqual( N(Sinh(2),10), 3.626860408,10); f();// matter of discussion whether rounding should be to nearest

NumericEqual( N(ArcSin(2), 9), Complex(1.570796327,1.316957897),9);f();
NumericEqual( N(ArcCos(2),9), Complex(0,-1.316957897),9);f();
NumericEqual( N(ArcTan(2*I), 12), N(Complex(1.5707963267950,0.54930614433405),12),11); f();// calculating to precision+1 because RoundToPrecision rounds... cluttering the last digit with round-off
NumericEqual( N(ArcSinh(2), 10), 1.443635475,9);f();
NumericEqual( N(ArcCosh(2), 10), 1.316957897,9);f();
NumericEqual( N(ArcCosh(-2), 8), Complex(-1.3169579,3.14159265),8);f();
NumericEqual( N(ArcTanh(2), 9), Complex(0.549306144,1.570796327),9);f();

/* Numerical tests - all confirmed with Maple */
BuiltinPrecisionSet(55);//Echo("BIP set to ",55);
// NOTE (hso,100311) BuiltinPrecisionSet MUST specify a value higher
//      than any of the precisions to be used in tests below.  
//      Otherwise, false errors are reported.  I have changed
//      the value from 50 to 55, to satisfy this requirement.

NumericEqual(
RoundToPrecision(N(Pi), 50)
, 3.14159265358979323846264338327950288419716939937511
, 50);f();

NumericEqual(
RoundToPrecision(N(Sin(2.0)), 49)
, 0.9092974268256816953960198659117448427022549714479
,48);f();

NumericEqual(
RoundToPrecision(N(Sin(2.0)), 50)
, 0.90929742682568169539601986591174484270225497144789
,49);f();

NumericEqual(
RoundToPrecision(N(Sin(2.0)), 51)
, 0.90929742682568169539601986591174484270225497144789
,50);f();

NumericEqual(
RoundToPrecision(N(Cos(20.0)), 49)
, 0.4080820618133919860622678609276449570992995103163
, 48); f();

NumericEqual(
RoundToPrecision(N(Tan(20.0)), 49)
, 2.2371609442247422652871732477303491783724839749188
, 48); f(); // to here. hso

NumericEqual(
RoundToPrecision(N(Exp(10.32),54), 54)
, 30333.2575962246035600343483350109621778376486335450125
,48); f();  // This one rounds off the wrong direction (125 rounded to 12 iso 13). But alas, change was needed because new interpretation means the required precision was actually higher (not number of decimals after point, but total number of digits were meant).

NumericEqual(
RoundToPrecision(N(Ln(10.32/4.07)), 50)
, 0.93044076059891305468974486564632598071134270468002
, 50); f();  //  fixed 20101229

NumericEqual(
RoundToPrecision(N(1.3^10.32), 48)
, 14.99323664825717956473936947123246987802978985306
, 48); f();

NumericEqual(
RoundToPrecision(N(Sqrt(5.3),51), 51)
, 2.302172886644267644194841586420201850185830282633675
,51); f(); // increased to 51 digits so round-off is obviously downwards (previous rounding was defendably wrong)

// this failed in gmp due to broken SqrtN()
NumericEqual(
RoundToPrecision(N(Sqrt(25.3)), 50)
, 5.0299105359837166789719353820984438468186649281130
,50);f();

// this failed due to broken RoundToPrecision()
NumericEqual(
RoundToPrecision(PowerN(13, -23), 50)
, 0.23949855470974885180294666343025235387321690490245e-25
, 50);f();

NumericEqual(
RoundToPrecision([Local(x);x:=Newton(x*Exp(x)-4,x,1,10^(-49)); N(x*Exp(x));], 49)
, 4
,49); f();
Verify(Newton(x^2+1,x,1,0.1,-3,3), Fail); f();
NumericEqual(Newton(x^2-1,x,1,0.1,-3,3), 1,BuiltinPrecisionGet()); f();

NumericEqual(
RoundToPrecision(N(ArcSin(0.32)), 49)
, 0.3257294872946301593103199105324500784354180998122808003
,49); f();

NumericEqual(
RoundToPrecision(N(Sin(N(ArcSin(0.1234567)))), 49)
, 0.1234567
,49); f();

/* ArcSin(x) for x close to 1 */

NumericEqual(
RoundToPrecision(N( (1-Sin(N(ArcSin(1-10^(-25)))))*10^25), 25)
, 1
, 25); f();

NumericEqual(
N(ArcSin(N(Sin(1.234567),50)),50)
, N(1.234567,50)
, 49); f(); // calculating to precision+1 because RoundToPrecision rounds... cluttering the last digit with round-off

NumericEqual(
RoundToPrecision(N(ArcCos(0.32)), 49)
, 1.2450668395002664599210017811073013636631665998753
, 49); f();

NumericEqual(
RoundToPrecision(N(ArcTan(0.32)), 49)
, 0.3097029445424561999173808103924156700884366304804
, 49); f();

NumericEqual(
RoundToPrecision(N(Cos(N(ArcCos(0.1234567)))), 49)
, 0.1234567
, 49);f();

NumericEqual(
RoundToPrecision(N(ArcCos(N(Cos(1.234567)))), 49)
, 1.234567
, 49);f();

NumericEqual(
RoundToPrecision(N(Tan(N(ArcTan(20)))), 46)	// large roundoff error on Tan() calculation due to subtraction from Pi/2 -- unavoidable loss of precision
, 20
, 46);f();
//KnownFailure(
NumericEqual(
RoundToPrecision(N(Tan(N(ArcTan(500000)))), 38)
, 500000
//)
, 38);f();

BuiltinPrecisionSet(60); //Echo("BIP set to ",60); // obviously, 50 is not enough for the following
//KnownFailure(
NumericEqual(
RoundToPrecision(N((Pi/2-ArcTan(N(Tan(N(Pi/2)-10^(-24)))))*10^24 ), 25)
, 1
//)
, 25);f();

/// special functions
BuiltinPrecisionSet(50); //Echo("BIP set to ",50); // needs a pretty high value for Gamma

TestMathPiper(
Gamma(10.5)
, (654729075*Sqrt(Pi))/1024
);f();

TestMathPiper(
Gamma(9/2)
, (105*Sqrt(Pi))/16
);f();

TestMathPiper(
Gamma(-10.5)
, (-2048*Sqrt(Pi))/13749310575
);f();

TestMathPiper(
Gamma(-7/2)
, (16*Sqrt(Pi))/105
);f();

NumericEqual(RoundToPlace(N( Internal'GammaNum(10.5) ), 13), 1133278.3889487855673, 13);f();

NumericEqual(RoundToPlace(N( Internal'GammaNum(-11.5) ), 20), 0.00000002295758104824, 20);f();

NumericEqual(RoundToPlace(N( Internal'GammaNum(-12.5) ), 20), -0.00000000183660648386, 20);f();

// Check for one example that N(Gamma(x)) returns the same as Internal'GammaNum
NumericEqual(RoundToPlace(N( Gamma(10.5) ), 13), 1133278.3889487855673, 13);f();



NumericEqual(	// lost 2 digits b/c of imprecise arithmetic
RoundToPlace(N( Zeta(-11.5) ), 18)
, 0.020396978715942792
,18);f();

TestMathPiper(
Zeta(40)
, (261082718496449122051*Pi^40)/20080431172289638826798401128390556640625
);f();

TestMathPiper(
Zeta(-11)
, 691/32760
);f();

TestMathPiper(
Zeta(-12)
, 0
);f();

NumericEqual(
RoundToPrecision(N(Zeta(40)), 19)
, 1.0000000000009094948
,19);f();

NumericEqual(
RoundToPrecision(N(Zeta(1.5)), 19)
, 2.6123753486854883433
,19);f();

// test correctness of Zeta(3)
NumericEqual(
RoundToPrecision(Internal'ZetaNum(3)-N(Zeta(3)), 20)
, 0
,20);f();

TestMathPiper(
Bernoulli(40)
, -261082718496449122051/13530
);f();

Verify(
ContFracList(355/113)
, {3,7,16}
);f();

Verify(
ContFracList(-24, 4)
, {-24}
);f();

Verify(
ContFracList(-355/113)
, {-4,1,6,16}
);f();

//BuiltinPrecisionSet(7);//Echo("BIP set to ",7);

Verify(
NearRational(N(Pi),3)
, 201/64
);f();

/*
  For the NearRational test, perhaps better would be a real test that
  checks that the result is correct up to the required number of digits
  accuracy.
*/
BuiltinPrecisionSet(10);//Echo("BIP set to ",10);
Verify(
NearRational(N(Pi))
, 355/113, 
);f();

// Lambert's W function
BuiltinPrecisionSet(10);

NumericEqual(
N(LambertW(-0.24),BuiltinPrecisionGet())
, -0.3357611648
, BuiltinPrecisionGet());f();

NumericEqual(
N(LambertW(10),BuiltinPrecisionGet())
, 1.7455280027
, BuiltinPrecisionGet());f();

// Bessel Functions
// These results are from GNU bc, matlab seems to suck.
BuiltinPrecisionSet(50);//Echo("BIP set to ",50);
NumericEqual( N(BesselJ(0,.5)), RoundToPrecision(.93846980724081290422840467359971262556892679709682,50),50 );f();
NumericEqual( N(BesselJ(0,.9)), RoundToPrecision(.80752379812254477730240904228745534863542363027564,50),50 );f();
NumericEqual( N(BesselJ(0,.99999)), RoundToPrecision(.76520208704756659155313775543958045290339472808482,50),50 );f();
NumericEqual( N(BesselJ(10,.75)), RoundToPrecision(.000000000014962171311759681469871248362168283485781647202136,50),50 );f();
NumericEqual( N(BesselJ(5,1)), RoundToPrecision(.00024975773021123443137506554098804519815836777698007,50),50 );f();
NumericEqual( N(BesselJ(4,2)), RoundToPrecision(.033995719807568434145759211288531044714832968346313,50),50 );f();
NumericEqual( N(BesselJ(10,3)), RoundToPrecision( .000012928351645715883777534530802580170743420832844164,50),50 );f();

NumericEqual( N(BesselJ(11,11)), RoundToPrecision( .20101400990926940339478738551009382430831534125484,50),50 );f();
NumericEqual( N(BesselJ(-11,11)), RoundToPrecision( -.20101400990926940339478738551009382430831534125484,50),50 );f();
NumericEqual( RoundToPrecision(N(BesselJ(1,10)),50), RoundToPrecision( .043472746168861436669748768025859288306272867118594, 50),50 );f();
NumericEqual( N(BesselJ(10,10)), RoundToPrecision( .20748610663335885769727872351875342803274461128682, 50 ),50 );f();
NumericEqual( RoundToPrecision(N(BesselJ(1,3.6)),50), RoundToPrecision( .095465547177876403845706744226060986019432754908851, 50 ),50) ;f();

BuiltinPrecisionSet(20);//Echo("BIP set to ",20);
Verify( RoundToPrecision(N(Erf(Sqrt(0.8)),20),19), 
RoundToPrecision(.79409678926793169113034892342, 19)
);f();

Verify( RoundToPrecision(N(Erf(50*I+20)/10^910,22),19), 
RoundToPrecision(1.09317119002909585408+I*0.00475463306931818955275, 19)
);f();

// testing GammaConstNum against Maple
Testing("Gamma constant");
BuiltinPrecisionSet(40);//Echo("BIP set to ",40);
NumericEqual(Internal'gamma()+0, 0.5772156649015328606065120900824024310422,BuiltinPrecisionGet());f();
BuiltinPrecisionSet(20);//Echo("BIP set to ",20);

Verify(gamma,ToAtom("gamma"));f();

NumericEqual(RoundToPrecision(Internal'gamma()+0,19), 0.5772156649015328606,19);f();

NumericEqual(RoundToPrecision(N(1/2+gamma+Pi), 19), 4.2188083184913260991,19);f();

// From GSL 1.0
//NumericEqual( N(PolyLog(2,-0.001),20), -0.00099975011104865108, 20 );
// PolyLog I didn't write PolyLog, but it seems to not always calculate correctly up to the last digit.
Verify( RoundToPrecision(N(PolyLog(2,-0.001)+0.00099975011104865108,20),20),0);f();

// Round-off errors
N([
  Local(a,b);
  a:= 77617;
  b:= 33096;
  // this expression gives a wrong answer on any hardware floating-point platform
  NumericEqual( 333.75*b^6 + a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8 +a/(2*b), -0.827396,6);
],40);f();