/** Tests of integration from NumericalIntegration.frink

    at about 80 digits, Frink:The Next Generation is about 24 times faster, so
    you should be using that.
*/

use NumericalIntegrationArbitrary.frink

// We need to make this a global definition but this lets us use the builtin
// 2-arguments sqrt[x, digits] with the default working precision by just
// calling arbitrarySqrt[x]
arbitrarySqrt[x, digits=getPrecision[]] := sqrt[x, digits]

// You shouldn't trust this above 90 because that's about how many digits of
// the "right" answer there are for some of these integrals.
digits = 80
setPrecision[digits]

pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
f = {|x| x arbitraryExp[x^2]}
test[f,-1,1,0]

f = {|x| 4 / (1+x^2) }     // Integrates to pi.  (4 arctan[1])
test[f, 0, 1., pi]

f = {|x| arbitrarySin[x]}           // Integrates to 0
test[f, 0, 2 pi, 0]

f = {|x| 1}
test[f, 0, 1, 1]

f = {|x| x}
test[f, -1, 1, 0]

f = {|x| x^2}
test[f, -1, 1, 2/3]

f = {|x| x^3}
test[f, -1, 1, 0]

// This won't be exact for Simpson's rule since degree > 3
// but adaptive Simpson will make it exact for degree up to 5.
f = {|x| x^4}
test[f, -1, 1, 2/5]

// This won't be exact for Simpson's rule since degree > 3
// but adaptive Simpson will make it exact for degree up to 5.
f = {|x| x^5}
test[f, 0, 1, 1/6]

// This won't be exact for Simpson's rule since degree > 3
// nor adaptive Simpson which is exact up to degree 5
f = {|x| x^6}
test[f, 0, 1, 1/7]

f = {|x| x^25}      // Integrates to 1/26
test[f, 0, 1, 1/26]

// EllipticF with additional (constant) parameter m.  This converges better.
x = 0.9
m = 1
exp = 1.0539231137410855767183688125065354946378914865080169298533839319453189726158212071060197769755030468304078734878239296454851973733534166863938
f = {|t,m| 1 / arbitrarySqrt[1 - m arbitrarySin[t]^2] }
testData[f,0,x,m,exp]

// Alternate formulation of EllipticF with additional (constant) parameter m.
// Use the above instead because it converges better.
f = {|t,m| 1 / (arbitrarySqrt[1-t^2] arbitrarySqrt[1 - m t^2])}
testData[f,0,arbitrarySin[x],m,exp]

// Gamma function at 1.1
f = {|x, g| arbitraryPow[x, (g-1)] arbitraryExp[-x]}
g = 1.1
testData[f,0,"inf", g, 0.9513507698668731836292487177265402192550578626088377343050000770434265383322821011537163794266447209797395258921102153439811044984082990793952]

// Gamma function at 100
g = 100
testData[f,0,"inf", g, (g-1)!]

// Gamma function at 1000.  This does not converge very well and gives only
// about 4 correct digits.
g = 1000
testData[f,0,"inf", g, (g-1)!]

// Tough integral because it oscillates to infinite frequency at 0
// WARNING:  This is really slow if you have rational bounds!
f = {|x| arbitrarySin[1/x]}
test[f, 0.01, 3, 1.530545039439096744804606825931883291397090653672480521828229153575456594997211465903528121943039515]

// Euler-Mascheroni constant, a few different ways
em = 0.577215664901532860606512090082402431042159335939923598805767234884867726777664670936947063291746749514631447249807082480960504014486542836224174

// below gives 78/80 digits, best
f = {|x| -arbitraryLn[-arbitraryLn[x]] } 
test[f, 0, 1, em]

// below gives 28/80 digits, slow convergence
// f = {|x| -x arbitraryExp[x-arbitraryExp[x]] } 
// test[f, "-inf", "inf", em]

// Below gives 57/80 digits, medium convergence
// f = {|x| -arbitraryLn[x] / arbitraryExp[x] } 
// test[f, 0, "inf", em]

/*
f = {|x| 1 m}    // Has units of measure in function
test[f, 0, 1, 1 m]

f = {|x| 1 }    // Has units of measure in integration bounds but not function
test[f, 0 m, 1 m, 1 m]

f = {|x| 1 m}    // Has units of measure in both integration bounds and function
test[f, 0 m, 1 m, 1 m^2]
*/

// Test integration against the expected target
test[f, a, b, target, eps=10^-digits] :=
{
   v = NIntegrate[f,a,b,eps]   // TODO:  Pass in digits?
   if target == 0
      err = abs[target] - v
   else
      err = (abs[target] - v) / target
   
   println[inputForm[f] + ": $v, err=$err"]
   println[]
}

// Test integration with additional parameters against the expected target
testData[f, a, b, data, target, eps=10^-digits] :=
{
   v = NIntegrateData[f,a,b,data,eps]  // TODO:  Pass in digits?
   if target == 0
      err = abs[target] - v
   else
      err = (abs[target] - v) / target

   println[inputForm[f] + ": $v, err=$err"]
   println[]
}