
# -- num/qag.lua.in
# -- 
# -- Copyright (C) 2009-2011 Francesco Abbate
# -- 
# -- This program is free software; you can redistribute it and/or modify
# -- it under the terms of the GNU General Public License as published by
# -- the Free Software Foundation; either version 3 of the License, or (at
# -- your option) any later version.
# -- 
# -- This program is distributed in the hope that it will be useful, but
# -- WITHOUT ANY WARRANTY; without even the implied warranty of
# -- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# -- General Public License for more details.
# -- 
# -- You should have received a copy of the GNU General Public License
# -- along with this program; if not, write to the Free Software
# -- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
# --

# -- Adapted from the GSL Library, version 1.14

# -- template parameters are 'limit' and 'order'

# half_order = (order - 1) / 2

# GSL_DBL_EPSILON = 2.2204460492503131e-16
# GSL_DBL_MIN     = 2.2250738585072014e-308
# GSL_DBL_MAX     = 1.7976931348623157e+308

local abs, min, max, pow = math.abs, math.min, math.max, math.pow

local ffi = require "ffi"

$(include 'gauss-kronrod-x-wgs')

local ws_size = 0
local ws_nrmax
local ws_i
local ws_maximum_level = 0

local ws_alist = ffi.new('double[$(limit)]')
local ws_blist = ffi.new('double[$(limit)]')
local ws_rlist = ffi.new('double[$(limit)]')
local ws_elist = ffi.new('double[$(limit)]')

local ws_order = ffi.new('unsigned int[$(limit)]')
local ws_level = ffi.new('unsigned int[$(limit)]')

local fv1 = ffi.new('unsigned int[$(half_order+1)]')
local fv2 = ffi.new('unsigned int[$(half_order+1)]')

local function rescale_error (err, result_abs, result_asc)
   err = abs(err)

   if result_asc ~= 0 and err ~= 0 then
      local scale = pow((200 * err / result_asc), 1.5)
      scale = min(scale, 1)
      err = result_asc * scale
   end

   if result_abs > 2.2250738585072e-308 / (50 * 2.2204460492503e-16) then
      local min_err = 50 * 2.2204460492503e-16 * result_abs
      err = max(err, min_err)
   end
   
   return err
end

local function q_integrator (f, a, b)
   local center = 0.5 * (a + b)
   local half_length = 0.5 * (b - a)
   local abs_half_length = abs(half_length)
   local f_center = f(center)

   local result_gauss = 0
   local result_kronrod = f_center * wgk[$(half_order)]

   local result_abs = abs(result_kronrod)

#  if half_order % 2 == 0 then
   local result_gauss = 0
#  else
   local result_gauss = f_center * wg[$((half_order-1)/2)]
#  end

   for j = 0, $(half_order/2 - 1) do
      local jtw = j * 2 + 1        -- j=1,2,3 jtw=2,4,6
      local abscissa = half_length * xgk[jtw]
      local fval1 = f(center - abscissa)
      local fval2 = f(center + abscissa)
      local fsum = fval1 + fval2
      fv1[jtw] = fval1
      fv2[jtw] = fval2
      result_gauss = result_gauss + wg[j] * fsum
      result_kronrod = result_kronrod + wgk[jtw] * fsum
      result_abs = result_abs + wgk[jtw] * (abs(fval1) + abs(fval2))
   end

   for j = 0, $((half_order+1)/2 - 1) do
      local jtwm1 = j * 2
      local abscissa = half_length * xgk[jtwm1]
      local fval1 = f(center - abscissa)
      local fval2 = f(center + abscissa)
      fv1[jtwm1] = fval1
      fv2[jtwm1] = fval2
      result_kronrod = result_kronrod + wgk[jtwm1] * (fval1 + fval2)
      result_abs = result_abs + wgk[jtwm1] * (abs(fval1) + abs(fval2))
   end

   local mean = result_kronrod * 0.5

   local result_asc = wgk[$(half_order)] * abs(f_center - mean)

   for j = 0, $(half_order-1) do
      result_asc = result_asc + wgk[j] * (abs(fv1[j] - mean) + abs(fv2[j] - mean))
   end

   -- scale by the width of the integration region

   local err = (result_kronrod - result_gauss) * half_length

   result_kronrod = result_kronrod * half_length
   result_abs = result_abs * abs_half_length
   result_asc = result_asc * abs_half_length

   local abserr = rescale_error (err, result_abs, result_asc)

   return result_kronrod, abserr, result_abs, result_asc
end

local function qpsrt()
   local last = ws_size - 1
   local i_nrmax = ws_nrmax
   local i_maxerr = ws_order[i_nrmax] 
  
   -- Check whether the list contains more than two error estimates

   if last < 2 then
      ws_order[0] = 0 
      ws_order[1] = 1 
      ws_i = i_maxerr 
      return 
   end

   local errmax = ws_elist[i_maxerr] 

   -- This part of the routine is only executed if, due to a difficult
   -- integrand, subdivision increased the error estimate. In the normal
   -- case the insert procedure should start after the nrmax-th largest
   -- error estimate.

   while i_nrmax > 0 and errmax > ws_elist[ws_order[i_nrmax - 1]] do
      ws_order[i_nrmax] = ws_order[i_nrmax - 1] 
      i_nrmax = i_nrmax - 1
   end

   -- Compute the number of elements in the list to be maintained in
   -- descending ws_order. This number depends on the number of
   -- subdivisions still allowed.
  
   local top = last < $(limit/2 + 2) and last or $(limit) - last + 1
  
   -- Insert errmax by traversing the list top-down, starting
   -- comparison from the element ws_elist(ws_order(i_nrmax+1)).
  
   local i = i_nrmax + 1 
  
   -- The ws_order of the tests in the following line is important to
   -- prevent a segmentation fault

   while i < top and errmax < ws_elist[ws_order[i]] do
      ws_order[i-1] = ws_order[i] 
      i = i+1 
   end
  
   ws_order[i-1] = i_maxerr 
  
   -- Insert errmin by traversing the list bottom-up
  
   local errmin = ws_elist[last] 
  
   local k = top - 1 
  
   while k > i - 2 and errmin >= ws_elist[ws_order[k]] do
      ws_order[k+1] = ws_order[k] 
      k = k - 1 
   end
  
   ws_order[k+1] = last 

   -- Set i_max and e_max
   i_maxerr = ws_order[i_nrmax] 
  
   ws_i = i_maxerr 
   ws_nrmax = i_nrmax 
end



local function update(a1, b1, area1, error1, a2, b2, area2, error2)
   local i_max = ws_i
   local i_new = ws_size

   local new_level = ws_level[i_max] + 1

   -- append the newly-created intervals to the list
  
   if error2 > error1 then
      ws_alist[i_max] = a2         -- blist[maxerr] is already == b2
      ws_rlist[i_max] = area2
      ws_elist[i_max] = error2
      ws_level[i_max] = new_level
      
      ws_alist[i_new] = a1
      ws_blist[i_new] = b1
      ws_rlist[i_new] = area1
      ws_elist[i_new] = error1
      ws_level[i_new] = new_level
   else
      ws_blist[i_max] = b1        -- alist[maxerr] is already == a1
      ws_rlist[i_max] = area1
      ws_elist[i_max] = error1
      ws_level[i_max] = new_level
      
      ws_alist[i_new] = a2
      ws_blist[i_new] = b2
      ws_rlist[i_new] = area2
      ws_elist[i_new] = error2
      ws_level[i_new] = new_level
   end
  
   ws_size = ws_size + 1

   if new_level > ws_maximum_level then
      ws_maximum_level = new_level
   end

   qpsrt ()
end

local function initialise(a, b)
  ws_size = 0
  ws_nrmax = 0
  ws_i = 0
  ws_alist[0] = a
  ws_blist[0] = b
  ws_rlist[0] = 0
  ws_elist[0] = 0
  ws_order[0] = 0
  ws_level[0] = 0

  ws_maximum_level = 0
end

local function subinterval_too_small(a1, a2, b2)
   local tmp = (1 + 100 * $(GSL_DBL_EPSILON)) * (abs(a2) + 1000 * $(GSL_DBL_MIN))
   return abs(a1) <= tmp and abs(b2) <= tmp
end

local function sum_results()
   local result_sum = 0
   local n = ws_size
   for k = 0, n-1 do
      result_sum = result_sum + ws_rlist[k]
   end
   return result_sum
end

local function qag (f, a, b, epsabs, epsrel, limit)
   local error_type = 0

   limit = limit or $(limit)
   initialise (a, b)

   if limit > $(limit) then
      error "iteration limit exceeds available workspace" 
   end

   if epsabs <= 0 and (epsrel < 50 * $(GSL_DBL_EPSILON) or epsrel < 0.5e-28) then
      error "tolerance cannot be acheived with given epsabs and epsrel"
   end

   -- perform the first integration

   local result0, abserr0, resabs0, resasc0 = q_integrator (f, a, b)

   ws_size = 1
   ws_rlist[0] = result0
   ws_elist[0] = abserr0

   -- Test on accuracy

   local tolerance = max(epsabs, epsrel * abs(result0))

   -- need IEEE rounding here to match original quadpack behavior

   round_off = 50 * $(GSL_DBL_EPSILON) * resabs0

   if abserr0 <= round_off and abserr0 > tolerance then
      error "cannot reach tolerance because of roundoff error on first attempt"
   elseif (abserr0 <= tolerance and abserr0 ~= resasc0) or abserr0 == 0 then
      return result0, abserr0
   elseif limit == 1 then
      error "a maximum of one iteration was insufficient"
   end

   local area, errsum = result0, abserr0
   local iteration = 1

   repeat
      local roundoff_type1, roundoff_type2 = 0, 0
      -- Bisect the subinterval with the largest error estimate

      local a_i, b_i = ws_alist[ws_i], ws_blist[ws_i]
      local r_i, e_i = ws_rlist[ws_i], ws_elist[ws_i]

      local m_i = 0.5 * (a_i + b_i)
      local a1, b1, a2, b2 = a_i, m_i, m_i, b_i;

      local area1, error1, resabs1, resasc1 = q_integrator (f, a1, b1)
      local area2, error2, resabs2, resasc2 = q_integrator (f, a2, b2)

      local area12 = area1 + area2
      local error12 = error1 + error2;

      errsum = errsum + (error12 - e_i)
      area = area + (area12 - r_i)

      if resasc1 ~= error1 and resasc2 ~= error2 then
	 local delta = r_i - area12

	 if abs(delta) <= 1.0e-5 * abs(area12) and error12 >= 0.99 * e_i then
	    roundoff_type1 = roundoff_type1 + 1
	 end
	 if iteration >= 10 and error12 > e_i then
	    roundoff_type2 = roundoff_type2 + 1
	 end
      end

      tolerance = max(epsabs, epsrel * abs(area))

      if errsum > tolerance then
	 if roundoff_type1 >= 6 or roundoff_type2 >= 20 then
	    error_type = 2   -- round off error
	 end

	 -- set error flag in the case of bad integrand behaviour at
         -- a point of the integration range

	 if subinterval_too_small(a1, a2, b2) then
	    error_type = 3
	 end
      end

      update (a1, b1, area1, error1, a2, b2, area2, error2)

      -- the following steps seems to be not useful
      -- a_i = ws_alist[ws_i]
      -- b_i = ws_blist[ws_i]
      -- r_i = ws_rlist[ws_i]
      -- e_i = ws_elist[ws_i]

      iteration = iteration + 1
   until iteration >= limit or error_type ~= 0 or errsum <= tolerance

   if errsum <= tolerance then
      return sum_results(), errsum
   elseif error_type == 2 then
      error "roundoff error prevents tolerance from being achieved"
   elseif error_type == 3 then
      error "bad integrand behavior found in the integration interval"
   elseif iteration == limit then
      error "maximum number of subdivisions reached"
   else
      error "could not integrate function" -- je peux pas
   end
end

return qag
