-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.FIR.PolyInterp
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Polynomial interpolators.  Taken from:
--
-- Olli Niemitalo (ollinie\@freenet.hut.fi), "Polynomial Interpolators for
-- High-Quality Resampling of Oversampled Audio" Search for "deip.pdf" with
-- Google and you will find it.
--
-----------------------------------------------------------------------------

-- TODO: limit the export list

-- TODO: figure out better way to create the coeficeints where you don't
-- have to explicitly state the number of interpolation points.

module DSP.Filter.FIR.PolyInterp where

import Data.Array

import Polynomial.Basic

-- | 'mkcoef' takes the continuous impluse response function (one of the
-- functions below, @f@) and number of points in the interpolation, @p@, time
-- shifts it by @x@, samples it, and creates an array with the interpolation
-- coeficients that can be used as a FIR filter.

mkcoef :: (Num a, Ix b, Integral b) => (a -> a) -- ^ f
       -> b -- ^ p
       -> a -- ^ x
       -> Array b a -- ^ h[n]

mkcoef f p x = listArray (0,p-1) $ map f [ x - fromIntegral i | i <- [p1..p2] ]
    where p1 = -(p `div` 2 - 1)
	  p2 = p `div` 2

---------------------------------------------------------------------------------

-- The impulse responses, centered around zero

-- The following functions are named like

-- blah_ApBo or optimal_ApBoCx

-- A = number of points in the interpolation
-- B = the polynomial order
-- C = the oversampling rate that the function is designed for

---------------------------------------------------------------------------------

-- B-Splines

bspline_1p0o :: (Ord a, Fractional a) => a -> a
bspline_1p0o x | 0 <= x && x < 1 = polyeval [ 1 ] x
               | otherwise       = 0

bspline_2p1o :: (Ord a, Fractional a) => a -> a
bspline_2p1o x | 0 <= x && x < 1 = polyeval [ 1, -1 ] x
               | 1 <= x          = 0
               | otherwise       = bspline_2p1o (-x)

bspline_4p3o :: (Ord a, Fractional a) => a -> a
bspline_4p3o x | 0 <= x && x < 1 = polyeval [ 2/3,  0, -1,  1/2 ] x
               | 1 <= x && x < 2 = polyeval [ 4/3, -2,  1, -1/6 ] x
               | 2 <= x          = 0
               | otherwise       = bspline_4p3o (-x)

bspline_6p5o :: (Ord a, Fractional a) => a -> a
bspline_6p5o x | 0 <= x && x < 1 = polyeval [ 11/20,     0, -1/2,    0,  1/4,  -1/12 ] x
               | 1 <= x && x < 2 = polyeval [ 17/40,   5/8, -7/4,  5/4, -3/8,   1/24 ] x
               | 2 <= x && x < 3 = polyeval [ 81/40, -27/8,  9/4, -3/4,  1/8, -1/120 ] x
               | 3 <= x          = 0
               | otherwise       = bspline_6p5o (-x)

---------------------------------------------------------------------------------

-- Lagrange polynomials

lagrange_4p3o :: (Ord a, Fractional a) => a -> a
lagrange_4p3o x | 0 <= x && x < 1 = polyeval [ 1,  -1/2, -1,  1/2 ] x
                | 1 <= x && x < 2 = polyeval [ 1, -11/6,  1, -1/6 ] x
                | 2 <= x          = 0
		| otherwise       = lagrange_4p3o (-x)

lagrange_6p5o :: (Ord a, Fractional a) => a -> a
lagrange_6p5o x | 0 <= x && x < 1 = polyeval [ 1,    -1/3, -5/4,   5/12,  1/4,  -1/12 ] x
                | 1 <= x && x < 2 = polyeval [ 1,  -13/12, -5/8,  25/24, -3/8,   1/24 ] x
                | 2 <= x && x < 3 = polyeval [ 1, -137/60, 15/8, -17/24,  1/8, -1/120 ] x
                | 3 <= x          = 0
		| otherwise       = lagrange_6p5o (-x)

---------------------------------------------------------------------------------

-- Hermite (1st-order-osculating) polynomials

hermite_4p3o :: (Ord a, Fractional a) => a -> a
hermite_4p3o x | 0 <= x && x < 1 = polyeval [ 1,  0, -5/2,  3/2 ] x
               | 1 <= x && x < 2 = polyeval [ 2, -4,  5/2, -1/2 ] x
               | 2 <= x          = 0
	       | otherwise       = hermite_4p3o (-x)

hermite_6p3o :: (Ord a, Fractional a) => a -> a
hermite_6p3o x | 0 <= x && x < 1 = polyeval [ 1,        0, -7/3,   4/3 ] x
               | 1 <= x && x < 2 = polyeval [ 5/2, -59/12,    3, -7/12 ] x
               | 2 <= x && x < 3 = polyeval [ -3/2,   7/4, -2/3,  1/12 ] x
               | 3 <= x          = 0
               | otherwise       = hermite_6p3o (-x)

hermite_6p5o :: (Ord a, Fractional a) => a -> a
hermite_6p5o x | 0 <= x && x < 1 = polyeval [ 1,     0, -25/12,   5/12, 13/12, -5/12 ] x
               | 1 <= x && x < 2 = polyeval [ 1,  5/12,  -35/8,   35/8, -13/8,  5/24 ] x
               | 2 <= x && x < 3 = polyeval [ 3, -29/4, 155/24, -65/24, 13/24, -1/24 ] x
               | 3 <= x          = 0
               | otherwise       = hermite_6p5o (-x)

---------------------------------------------------------------------------------

-- 2nd-order-osculating polynomials

sndosc_4p5o :: (Ord a, Fractional a) => a -> a
sndosc_4p5o x | 0 <= x && x < 1 = polyeval [  1, 0,   -1, -9/2,  15/2, -3 ] x
              | 1 <= x && x < 2 = polyeval [ -4, 18, -29, 43/2, -15/2,  1 ] x
              | 2 <= x          = 0
	      | otherwise       = sndosc_4p5o (-x)

sndosc_6p5o :: (Ord a, Fractional a) => a -> a
sndosc_6p5o x | 0 <= x && x < 1 = polyeval [  1,      0,   -5/4,  -35/12,  21/4, -25/12 ] x
              | 1 <= x && x < 2 = polyeval [ -4,   75/4, -245/8,  545/24, -63/8,  25/24 ] x
              | 2 <= x && x < 3 = polyeval [ 18, -153/4,  255/8, -313/24,  21/8,  -5/24 ] x
              | 3 <= x          = 0
              | otherwise       = sndosc_6p5o (-x)

---------------------------------------------------------------------------------

-- Misc

watte_4p2o :: (Ord a, Fractional a) => a -> a
watte_4p2o x | 0 <= x && x < 1 = polyeval [ 1, -1/2, -1/2 ] x
             | 1 <= x && x < 2 = polyeval [ 1, -3/2,  1/2 ] x
             | 2 <= x          = 0
	     | otherwise       = watte_4p2o (-x)

parabolic2x_4p2o :: (Ord a, Fractional a) => a -> a
parabolic2x_4p2o x | 0 <= x && x < 1 = polyeval [ 1/2, 0, -1/4 ] x
                   | 1 <= x && x < 2 = polyeval [ 1,  -1,  1/4 ] x
                   | 2 <= x          = 0
		   | otherwise       = parabolic2x_4p2o (-x)

---------------------------------------------------------------------------------

-- Optimal designs

optimal_2p3o2x :: (Ord a, Fractional a) => a -> a
optimal_2p3o2x x | 0 <= x && x < 1 = polyeval [ 0.80607906469176971, 0.17594740788514596,
						-2.35977550974341630, 1.57015627178718420 ] x
                 | 1 <= x          = 0
		 | otherwise       = optimal_2p3o2x (-x)

optimal_2p3o4x :: (Ord a, Fractional a) => a -> a
optimal_2p3o4x x | 0 <= x && x < 1 = polyeval [ 0.88207975731800936, -0.10012219395448523,
						-1.99054787320203810, 1.32598918957298410 ] x
                 | 1 <= x          = 0
		 | otherwise       = optimal_2p3o4x (-x)

optimal_2p3o8x :: (Ord a, Fractional a) => a -> a
optimal_2p3o8x x | 0 <= x && x < 1 = polyeval [ 0.94001491168487883, -0.51213628865925998,
					        -1.10319974084152170, 0.73514591836770027 ] x
                 | 1 <= x          = 0
		 | otherwise       = optimal_2p3o8x (-x)

optimal_2p3o16x :: (Ord a, Fractional a) => a -> a
optimal_2p3o16x x | 0 <= x && x < 1 = polyeval [ 0.96964782067188493, -0.74617479745643256,
						 -0.57923093055631791, 0.38606621963374965 ] x
                  | 1 <= x          = 0
		  | otherwise       = optimal_2p3o16x (-x)

optimal_2p3o32x :: (Ord a, Fractional a) => a -> a
optimal_2p3o32x x | 0 <= x && x < 1 = polyeval [ 0.98472017575676363, -0.87053863725307623,
					         -0.29667081825572522, 0.19775766248673177 ] x
                  | 1 <= x          = 0
	          | otherwise       = optimal_2p3o32x (-x)

optimal_4p2o2x :: (Ord a, Fractional a) => a -> a
optimal_4p2o2x x | 0 <= x && x < 1 = polyeval [ 0.50061662213752656, -0.04782068534965925,
					        -0.21343978756177684 ] x
                 | 1 <= x && x < 2 = polyeval [ 0.92770135528027386, -0.88689658749623701,
					        0.21303593243799016  ] x
                 | 2 <= x          = 0
	         | otherwise       = optimal_4p2o2x (-x)

optimal_4p2o4x :: (Ord a, Fractional a) => a -> a
optimal_4p2o4x x | 0 <= x && x < 1 = polyeval [ 0.33820365736567115, 0.2114449807519728,
					        -0.22865399531858188  ] x
                 | 1 <= x && x < 2 = polyeval [ 1.12014639874555470, -1.01414466618792900,
					        0.22858390767180370  ] x
                 | 2 <= x          = 0
	         | otherwise       = optimal_4p2o4x (-x)

optimal_4p2o8x :: (Ord a, Fractional a) => a -> a
optimal_4p2o8x x | 0 <= x && x < 1 = polyeval [ 0.09224718574204172, 0.59257579283164508,
					        -0.24005206207889518  ] x
                 | 1 <= x && x < 2 = polyeval [ 1.38828036063664320, -1.17126532964206100,
					        0.24004281672637814  ] x
                 | 2 <= x          = 0
	         | otherwise       = optimal_4p2o8x (-x)

optimal_4p2o16x :: (Ord a, Fractional a) => a -> a
optimal_4p2o16x x | 0 <= x && x < 1 = polyeval [ -0.41849525763976203, 1.36361593203840510,
					         -0.24506117865474364  ] x
                  | 1 <= x && x < 2 = polyeval [ 1.90873339502208310, -1.44144384373471430,
					         0.24506002360805534  ] x
                  | 2 <= x          = 0
	          | otherwise       = optimal_4p2o16x (-x)

optimal_4p2o32x :: (Ord a, Fractional a) => a -> a
optimal_4p2o32x x | 0 <= x && x < 1 = polyeval [ -1.42170796824052890, 2.87083485132510450,
					         -0.24755243839713828 ] x
                  | 1 <= x && x < 2 = polyeval [ 2.91684291662070860, -1.95043794419108290,
					        0.24755229501840223 ] x
                  | 2 <= x          = 0
	          | otherwise       = optimal_4p2o32x (-x)

optimal_4p3o2x :: (Ord a, Fractional a) => a -> a
optimal_4p3o2x x | 0 <= x && x < 1 = polyeval [ 0.59244492420272321, 0.03573669883299365,
					        -0.78664888597764893, 0.36030925263849456 ] x
                 | 1 <= x && x < 2 = polyeval [ 1.20220428331406090, -1.60101160971478710,
					        0.70401463131621556, -0.10174985775982505 ] x
                 | 2 <= x          = 0
	         | otherwise       = optimal_4p3o2x (-x)

optimal_4p3o4x :: (Ord a, Fractional a) => a -> a
optimal_4p3o4x x | 0 <= x && x < 1 = polyeval [ 0.60304009430474115, 0.05694012453786401,
					        -0.89223007211175309, 0.42912649274763925 ] x
                 | 1 <= x && x < 2 = polyeval [ 1.31228823423882930, -1.85072890189700660,
					        0.87687351895686727, -0.13963062613760227 ] x
                 | 2 <= x          = 0
	         | otherwise       = optimal_4p3o4x (-x)

optimal_4p3o8x :: (Ord a, Fractional a) => a -> a
optimal_4p3o8x x | 0 <= x && x < 1 = polyeval [ 0.60658368706046584, 0.07280793921972525,
					        -0.95149675410360302, 0.46789242171187317 ] x
                 | 1 <= x && x < 2 = polyeval [ 1.35919815911169020, -1.95618744839533010,
					        0.94949311590826524, -0.15551896027602030 ] x
                 | 2 <= x          = 0
	         | otherwise       = optimal_4p3o8x (-x)

optimal_4p3o16x :: (Ord a, Fractional a) => a -> a
optimal_4p3o16x x | 0 <= x && x < 1 = polyeval [ 0.60844825096346644, 0.07980169577604959,
					         -0.97894238166068270, 0.48601256046234864 ] x
                  | 1 <= x && x < 2 = polyeval [ 1.37724137476464990, -1.99807048591354810,
					         0.97870442828560433, -0.16195131297091253 ] x
                  | 2 <= x          = 0
	          | otherwise       = optimal_4p3o16x (-x)

optimal_4p3o32x :: (Ord a, Fractional a) => a -> a
optimal_4p3o32x x | 0 <= x && x < 1 = polyeval [ 0.60908264223655417, 0.08298544053689563,
					         -0.99052586766084594, 0.49369595780454456 ] x
                  | 1 <= x && x < 2 = polyeval [ 1.38455689452848450, -2.01496368680360890,
					         0.99049753216621961, -0.16455902278580614 ] x
                  | 2 <= x          = 0
	          | otherwise       = optimal_4p3o32x (-x)

optimal_4p4o2x :: (Ord a, Fractional a) => a -> a
optimal_4p4o2x x | 0 <= x && x < 1 = polyeval [ 0.58448510036125145, 0.04442540676862300,
					        -0.7586487041827807, 0.29412762852131868,
					        0.04252164479749607 ] x
                 | 1 <= x && x < 2 = polyeval [ 1.06598379704160570, -1.16581445347275190,
					        0.21256821036268256, 0.13781898240764315,
					        -0.04289144034653719 ] x
                 | 2 <= x          = 0
	         | otherwise       = optimal_4p4o2x (-x)

optimal_4p4o4x :: (Ord a, Fractional a) => a -> a
optimal_4p4o4x x | 0 <= x && x < 1 = polyeval [ 0.61340295990566229, 0.06128937679587994,
					        -0.94057832565094635, 0.44922093286355397,
					        0.00986988334359864 ] x
                 | 1 <= x && x < 2 = polyeval [ 1.30835018075821670, -1.82814511658458520,
					        0.81943257721092366, -0.09642760567543440,
					        -0.00989340017126506 ] x
                 | 2 <= x          = 0
	         | otherwise       = optimal_4p4o4x (-x)

optimal_4p4o8x :: (Ord a, Fractional a) => a -> a
optimal_4p4o8x x | 0 <= x && x < 1 = polyeval [ 0.62095991632974834, 0.06389302461261143,
					       -0.98489647972932193, 0.48698871865064902,
					        0.00255074537015887 ] x
                 | 1 <= x && x < 2 = polyeval [ 1.35943398999940390, -1.97277963497287720,
					        0.95410568622888214, -0.14868053358928229,
					       -0.00255226912537286 ] x
                 | 2 <= x          = 0
	         | otherwise       = optimal_4p4o8x (-x)

optimal_4p4o16x :: (Ord a, Fractional a) => a -> a
optimal_4p4o16x x | 0 <= x && x < 1 = polyeval [ 0.62293049365660191, 0.06443376638262904,
					        -0.99620011474430481, 0.49672182806667398,
					         0.00064264050033187 ] x
                  | 1 <= x && x < 2 = polyeval [ 1.37216269878963180, -2.00931632449031920,
					         0.98847675044522398, -0.16214364417487748,
					        -0.00064273459469381 ] x
                  | 2 <= x          = 0
	          | otherwise       = optimal_4p4o16x (-x)

optimal_4p4o32x :: (Ord a, Fractional a) => a -> a
optimal_4p4o32x x | 0 <= x && x < 1 = polyeval [ 0.62342449465938121, 0.06456923251842608,
					        -0.99904509583176049, 0.49917660509564427,
					         0.00016095224137360 ] x
                  | 1 <= x && x < 2 = polyeval [ 1.37534629142898650, -2.01847637982642340,
					         0.99711292321092770, -0.16553360612350931,
					        -0.00016095810460478 ] x
                  | 2 <= x          = 0
	          | otherwise       = optimal_4p4o32x (-x)

optimal_6p4o2x :: (Ord a, Fractional a) => a -> a
optimal_6p4o2x x | 0 <= x && x < 1 = polyeval [ 0.42640922432669054, -0.0052558029434142,
					       -0.20486985491012843, 0.00255494211547300,
					        0.03134095684084392 ] x
                 | 1 <= x && x < 2 = polyeval [ 0.30902529029941583, 0.37868437559565432,
					       -0.70564644117967990, 0.31182026815653541,
					       -0.04385804833432710 ] x
                 | 2 <= x && x < 3 = polyeval [ 1.51897639740576910, -1.83761742915820410,
					        0.83217835730406542, -0.16695522597587154,
					        0.01249475765486819 ] x
                 | 3 <= x          = 0
	         | otherwise       = optimal_6p4o2x (-x)

optimal_6p4o4x :: (Ord a, Fractional a) => a -> a
optimal_6p4o4x x | 0 <= x && x < 1 = polyeval [ 0.20167941634921072, -0.06119274485321008,
					        0.56468711069379207, -0.42059475673758634,
					        0.02881527997393852 ] x
                 | 1 <= x && x < 2 = polyeval [ -0.64579641436229407, 2.33580825807694700,
					        -1.85350543411307390, 0.51926458031522660,
					        -0.04250898918476453 ] x
                 | 2 <= x && x < 3 = polyeval [ 2.76228852293285200, -3.09936092833253300,
					        1.27147464005834010, -0.22283280665600644,
					        0.01369173779618459 ] x
                 | 3 <= x          = 0
	         | otherwise       = optimal_6p4o4x (-x)

optimal_6p4o8x :: (Ord a, Fractional a) => a -> a
optimal_6p4o8x x | 0 <= x && x < 1 = polyeval [ -0.17436452172055789, -0.15190225510786248,
					         1.87551558979819120, -1.15976496200057480,
					         0.03401038103941584 ] x
                 | 1 <= x && x < 2 = polyeval [ -2.26955357035241170, 5.73320660746477540,
					        -3.92391712129699590, 0.93463067895166918,
					        -0.05090907029392906 ] x
                 | 2 <= x && x < 3 = polyeval [ 4.84834508915762540, -5.25661448354449060,
					        2.04584149450148180, -0.32814290420019698,
					        0.01689861603514873 ] x
                 | 3 <= x          = 0
	         | otherwise       = optimal_6p4o8x (-x)

optimal_6p4o16x :: (Ord a, Fractional a) => a -> a
optimal_6p4o16x x | 0 <= x && x < 1 = polyeval [ -0.94730014688427577, -0.33649680079382827,
					          4.53807483241466340, -2.64598691215356660,
					          0.03755086455339280 ] x
                  | 1 <= x && x < 2 = polyeval [ -5.55035312316726960, 12.52871168241192600,
					         -7.98288364772738750, 1.70665858343069510,
					         -0.05631219122315393 ] x
                  | 2 <= x && x < 3 = polyeval [ 8.94785524286246310, -9.37021675593126700,
					         3.44447036756440590, -0.49470749109917245,
					         0.01876132424143207 ] x
                  | 3 <= x          = 0
	          | otherwise       = optimal_6p4o16x (-x)

optimal_6p4o32x :: (Ord a, Fractional a) => a -> a
optimal_6p4o32x x | 0 <= x && x < 1 = polyeval [ -2.44391738331193720, -0.69468212315980082,
					          9.67889243081689440, -5.50592307590218160,
					          0.03957507923965987 ] x
                  | 1 <= x && x < 2 = polyeval [ -11.87524595267807600, 25.58633277328986500,
					         -15.73068663442630400, 3.15288929279855570,
					         -0.05936083498715066 ] x
                  | 2 <= x && x < 3 = polyeval [ 16.79403235763479100, -17.17264148794549100,
					         6.05175140696421730, -0.79053754554850286,
					         0.01978575568000696 ] x
                  | 3 <= x          = 0
	          | otherwise       = optimal_6p4o32x (-x)

optimal_6p5o2x :: (Ord a, Fractional a) => a -> a
optimal_6p5o2x x | 0 <= x && x < 1 = polyeval [ 0.48217702203158502, -0.00127577239632662,
					       -0.3267507171395277, -0.02014846731685776,
					        0.14640674192652170, -0.04317950185225609 ] x
                 | 1 <= x && x < 2 = polyeval [ 0.35095903476754237, 0.53534756396439365,
					       -1.22477236472789920, 0.74995484587342742,
					       -0.19234043023690772, 0.01802814255926417 ] x
                 | 2 <= x && x < 3 = polyeval [ 1.62814578813495040, -2.26168360510917840,
					        1.22220278720010690, -0.31577407091450355,
					        0.03768876199398620, -0.00152170021558204 ] x
                 | 3 <= x          = 0
	         | otherwise       = optimal_6p5o2x (-x)

optimal_6p5o4x :: (Ord a, Fractional a) => a -> a
optimal_6p5o4x x | 0 <= x && x < 1 = polyeval [ 0.50164509338655083, -0.00256790184606694,
					       -0.36229943140977111, -0.04512026308730401,
					        0.20620318519804220, -0.06607747864416924 ] x
                 | 1 <= x && x < 2 = polyeval [ 0.30718330223223800, 0.78336433172501685,
					       -1.66940481896969310, 1.08365113099941970,
					       -0.30560854964737405, 0.03255079211953620 ] x
                 | 2 <= x && x < 3 = polyeval [ 2.05191571792256240, -3.19403437421534920,
					        1.99766476840488070, -0.62765808573554227,
					        0.09909173357642603, -0.00628989632244913 ] x
		 | 3 <= x          = 0
	         | otherwise       = optimal_6p5o4x (-x)

optimal_6p5o8x :: (Ord a, Fractional a) => a -> a
optimal_6p5o8x x | 0 <= x && x < 1 = polyeval [ 0.50513183702821474, -0.00368143670114908,
					       -0.36434084624989699, -0.06070462616102962,
					        0.22942797169644802, -0.07517133281176167 ] x
                 | 1 <= x && x < 2 = polyeval [ 0.28281884957695946, 0.88385964850687193,
					       -1.82581238657617080, 1.19588167464050650,
					       -0.34363487882262922, 0.03751837438141215 ] x
                 | 2 <= x && x < 3 = polyeval [ 2.15756386503245070, -3.42137079071284810,
					        2.18592382088982260, -0.70370361187427199,
					        0.11419603882898799, -0.00747588873055296 ] x
                 | 3 <= x          = 0
	         | otherwise       = optimal_6p5o8x (-x)

optimal_6p5o16x :: (Ord a, Fractional a) => a -> a
optimal_6p5o16x x | 0 <= x && x < 1 = polyeval [ 0.50819303579369868, -0.00387117789818541,
					        -0.36990908725555449, -0.06616250180411522,
					         0.24139298776307896, -0.07990500783668089 ] x
                  | 1 <= x && x < 2 = polyeval [ 0.27758734130911511, 0.91870010875159547,
					        -1.89281840112089440, 1.24834464824612510,
					        -0.36203450650610985, 0.03994519162531633   ] x
                  | 2 <= x && x < 3 = polyeval [ 2.19284545406407450, -3.50786533926449100,
					         2.26228244623301580, -0.73559668875725392,
					         0.12064126711558003, -0.00798609327859495   ] x
                  | 3 <= x          = 0
	          | otherwise       = optimal_6p5o16x (-x)

optimal_6p5o32x :: (Ord a, Fractional a) => a -> a
optimal_6p5o32x x | 0 <= x && x < 1 = polyeval [ 0.52558916128536759, 0.00010896283126635,
					        -0.42682321682847008, -0.04095676092513167,
					         0.25041444762720882, -0.08349799235675044 ] x
                  | 1 <= x && x < 2 = polyeval [ 0.33937904183610190, 0.80946953063234006,
					        -1.86228986389877100, 1.27215033630638800,
					        -0.37562266426589430, 0.04174912841630993 ] x
                  | 2 <= x && x < 3 = polyeval [ 2.13606003964474490, -3.48774662195185850,
					         2.28912105276248390, -0.75510203509083995,
					         0.12520821766375972, -0.00834987866042734  ] x
                  | 3 <= x          = 0
	          | otherwise       = optimal_6p5o32x (-x)

---------------------------------------------------------------------------------

{-------------------

Test routines

y = [ sin $ 0.345 + 0.1234 * fromIntegral i | i <- [0..10] ]

h1 = mkcoef bspline_4p3o     4 0.2
h2 = mkcoef hermite_4p3o     4 0.2
h3 = mkcoef lagrange_4p3o    4 0.2
h4 = mkcoef hermite_4p3o     4 0.2
h5 = mkcoef sndosc_4p5o      4 0.2
h6 = mkcoef watte_4p2o       4 0.2
h7 = mkcoef parabolic2x_4p2o 4 0.2

h8  = mkcoef bspline_6p5o  6 0.2
h9  = mkcoef lagrange_6p5o 6 0.2
h10 = mkcoef hermite_6p3o  6 0.2
h11 = mkcoef hermite_6p5o  6 0.2
h12 = mkcoef sndosc_6p5o   6 0.2

h2p3o2x  = mkcoef optimal_2p3o2x  2 0.2
h2p3o4x  = mkcoef optimal_2p3o4x  2 0.2
h2p3o8x  = mkcoef optimal_2p3o8x  2 0.2
h2p3o16x = mkcoef optimal_2p3o16x 2 0.2
h2p3o32x = mkcoef optimal_2p3o32x 2 0.2

h4p2o2x  = mkcoef optimal_4p2o2x  4 0.2
h4p2o4x  = mkcoef optimal_4p2o4x  4 0.2
h4p2o8x  = mkcoef optimal_4p2o8x  4 0.2
h4p2o16x = mkcoef optimal_4p2o16x 4 0.2
h4p2o32x = mkcoef optimal_4p2o32x 4 0.2

h4p3o2x  = mkcoef optimal_4p3o2x  4 0.2
h4p3o4x  = mkcoef optimal_4p3o4x  4 0.2
h4p3o8x  = mkcoef optimal_4p3o8x  4 0.2
h4p3o16x = mkcoef optimal_4p3o16x 4 0.2
h4p3o32x = mkcoef optimal_4p3o32x 4 0.2

h4p4o2x  = mkcoef optimal_4p4o2x  4 0.2
h4p4o4x  = mkcoef optimal_4p4o4x  4 0.2
h4p4o8x  = mkcoef optimal_4p4o8x  4 0.2
h4p4o16x = mkcoef optimal_4p4o16x 4 0.2
h4p4o32x = mkcoef optimal_4p4o32x 4 0.2

h6p4o2x  = mkcoef optimal_6p4o2x  4 0.2
h6p4o4x  = mkcoef optimal_6p4o4x  4 0.2
h6p4o8x  = mkcoef optimal_6p4o8x  4 0.2
h6p4o16x = mkcoef optimal_6p4o16x 4 0.2
h6p4o32x = mkcoef optimal_6p4o32x 4 0.2

h6p5o2x  = mkcoef optimal_6p5o2x  4 0.2
h6p5o4x  = mkcoef optimal_6p5o4x  4 0.2
h6p5o8x  = mkcoef optimal_6p5o8x  4 0.2
h6p5o16x = mkcoef optimal_6p5o16x 4 0.2
h6p5o32x = mkcoef optimal_6p5o32x 4 0.2

interpolate y h = sum $ zipWith (*) y (elems h)

x1  = sin $ 0.345 + 0.1234 * 1.2
x1' = map (interpolate y) [ h1, h2, h3, h4, h5, h6, h7 ]

x2  = sin $ 0.345 + 0.1234 * 2.2
x2' = map (interpolate y) [ h8, h9, h10, h11, h12 ]

The values of all these lists should be one, or nearly one.  They
aren't for the 6p4o optimal designs, but I'm not sure why.  Olli's
paper states that these are a little screwy, though.

h_test = map (sum . elems) [ h1, h2, h3, h4, h5, h6, h7, h8, h9, h10, h11, h12 ]
h2p3o_test = map (sum . elems) [ h2p3o2x, h2p3o4x, h2p3o8x, h2p3o16x, h2p3o32x ]
h4p2o_test = map (sum . elems) [ h4p2o2x, h4p2o4x, h4p2o8x, h4p2o16x, h4p2o32x ]
h4p3o_test = map (sum . elems) [ h4p4o2x, h4p4o4x, h4p4o8x, h4p4o16x, h4p4o32x ]
h4p4o_test = map (sum . elems) [ h4p4o2x, h4p4o4x, h4p4o8x, h4p4o16x, h4p4o32x ]
h6p4o_test = map (sum . elems) [ h6p4o2x, h6p4o4x, h6p4o8x, h6p4o16x, h6p4o32x ]
h6p5o_test = map (sum . elems) [ h6p5o2x, h6p5o4x, h6p5o8x, h6p5o16x, h6p5o32x ]

-------------------}