/** Tests of integration from NumericalIntegration.frink */

use NumericalIntegration.frink

f = {|x| x exp[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| sin[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]

// Tough integral because it oscillates to infinite frequency at 0
// Should be
// 1.530545039439096744804606825931883291397090653672480521828229153575456594997211465903528121943039515
f = {|x| sin[1/x]}
test[f, 1/100, 3, 1.53054503943909674480460682593]

// EllipticF with additional (constant) parameter m.
x = 0.9
m = 1
exp = 1.05392311374108557671836881250653549463789148650801692985338393194531897
f = {|t,m| 1 / sqrt[1 - m sin[t]^2] }
testData[f,0,x,m,exp]

// Alternate formulation of EllipticF with additional (constant) parameter m.
f = {|t,m| 1 / (sqrt[1-t^2] sqrt[1 - m t^2])}
testData[f,0,sin[x],m,exp]

// Gamma function at 1.1
f = {|x, g| x^(g-1) exp[-x]}
g = 1.1
testData[f,0,"inf", g, 0.9513507698668731836292487177]

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

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

/*
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] :=
{
   v = NIntegrate[f,a,b]
   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] :=
{
   v = NIntegrateData[f,a,b,data]
   if target == 0
      err = abs[target] - v
   else
      err = (abs[target] - v) / target

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