begin
public pli, inv, diffh, diffvh, integral

use math: nan, floor

function pli(x0,d,a)
   n = len(a)
   undefined = nan if n==0 else nan*a[0]
   return fn|x|
      k = int(floor((x-x0)/d))
      if k<0 or k+1>=n
         return undefined
      else
         return (a[k+1]-a[k])/d*(x-x0-k*d)+a[k]
      end
   end
end

function inv(f,x,a,b)
   a1 = a; b1 = b
   s = sgn(f(b)-f(a))
   if s==0 then s=1 end

   for k in 60
      m = (a+b)/2
      if s*(f(m)-x)<0
         a = m
      else
         b = m
      end
   end
   if abs(m-a1)<1E-8
      return nan
   elif abs(m-b1)<1E-8
      return nan
   else
      return m
   end
end

function diffh(argm={})
   {h=0.001,order=false,fast=false} = argm
   if order
      use math.na.difftab: tab_diff
      return tab_diff(h)
   elif fast
      return fn diff|f,x|
         return (f(x+h)-f(x-h))/(2*h)
      end
   else
      return fn diff|f,x|
         return (f(x-2*h)-8*f(x-h)+8*f(x+h)-f(x+2*h))/(12*h)
      end
   end
end

function diffvh(argm={})
   {h=0.001,order=false,fast=false} = argm
   if order
      use math.na.difftab: tab_diffv
      return tab_diffv(h)
   elif fast
      return fn diffv|v,f,x|
         hv = h*v
         return (f(x+hv)-f(x-hv))/(2*h)
      end
   else
      return fn diffv|v,f,x|
         hv = h*v
         return (f(x-2*hv)-8*f(x-hv)+8*f(x+hv)-f(x+2*hv))/(12*h)
      end
   end
end


GL32 = [
   [-0.9972638618494816, 0.007018610009470141],
   [-0.9856115115452684, 0.01627439473090563],
   [-0.9647622555875064, 0.02539206530926208],
   [-0.9349060759377397, 0.03427386291302148],
   [-0.8963211557660522, 0.04283589802222668],
   [-0.8493676137325699, 0.05099805926237622],
   [-0.7944837959679425, 0.05868409347853551],
   [-0.7321821187402897, 0.06582222277636184],
   [-0.6630442669302152, 0.07234579410884850],
   [-0.5877157572407623, 0.07819389578707028],
   [-0.5068999089322295, 0.08331192422694673],
   [-0.4213512761306354, 0.08765209300440380],
   [-0.3318686022821277, 0.09117387869576390],
   [-0.2392873622521371, 0.09384439908080457],
   [-0.1444719615827965, 0.09563872007927489],
   [-0.04830766568773831,0.09654008851472778],
   [ 0.04830766568773831,0.09654008851472778],
   [ 0.1444719615827965, 0.09563872007927489],
   [ 0.2392873622521371, 0.09384439908080457],
   [ 0.3318686022821277, 0.09117387869576390],
   [ 0.4213512761306354, 0.08765209300440380],
   [ 0.5068999089322295, 0.08331192422694673],
   [ 0.5877157572407623, 0.07819389578707028],
   [ 0.6630442669302152, 0.07234579410884850],
   [ 0.7321821187402897, 0.06582222277636184],
   [ 0.7944837959679425, 0.05868409347853551],
   [ 0.8493676137325699, 0.05099805926237622],
   [ 0.8963211557660522, 0.04283589802222668],
   [ 0.9349060759377397, 0.03427386291302148],
   [ 0.9647622555875064, 0.02539206530926208],
   [ 0.9856115115452684, 0.01627439473090563],
   [ 0.9972638618494816, 0.007018610009470141]
]

function gauss(gl)
   return fn integral|a,b,f,n=1|
      h = (b-a)/n
      p = 0.5*h
      s = 0
      for j in 0..n-1
         q = p+a+j*h
         sj = 0
         for t in gl
            sj += t[1]*f(p*t[0]+q)
         end
         s += p*sj
      end
      return s
   end
end

integral = gauss(GL32)

end