GSL 0.4.23

//
// A rust binding for the GSL library by Guillaume Gomez (guillaume1.gomez@gmail.com)
//

#![allow(non_upper_case_globals)]

use ffi;
use enums;
use std::intrinsics::{fabsf64, logf64, powf64, sinf64, cosf64};
use num::Float;
use c_vec::CSlice;

static XI : [f64; 33] = [
  -1f64, -0.99518472667219688624f64, -0.98078528040323044912f64,
  -0.95694033573220886493f64, -0.92387953251128675612f64,
  -0.88192126434835502970f64, -0.83146961230254523708f64,
  -0.77301045336273696082f64, -0.70710678118654752440f64,
  -0.63439328416364549822f64, -0.55557023301960222475f64,
  -0.47139673682599764857f64, -0.38268343236508977173f64,
  -0.29028467725446236764f64, -0.19509032201612826785f64,
  -0.098017140329560601995f64, 0f64, 0.098017140329560601995f64,
  0.19509032201612826785f64, 0.29028467725446236764f64, 0.38268343236508977173f64,
  0.47139673682599764857f64, 0.55557023301960222475f64, 0.63439328416364549822f64,
  0.70710678118654752440f64, 0.77301045336273696082f64, 0.83146961230254523708f64,
  0.88192126434835502970f64, 0.92387953251128675612f64, 0.95694033573220886493f64,
  0.98078528040323044912f64, 0.99518472667219688624f64, 1f64
];

static bee : [f64; 68] = [
  0.00000000000000e+00f64, 2.28868854108532e-01f64, 0.00000000000000e+00f64,
  -8.15740215243451e-01f64, 0.00000000000000e+00f64, 5.31212715259731e-01f64,
  0.00000000000000e+00f64, 1.38538036812454e-02f64, 0.00000000000000e+00f64,
  3.74405228908818e-02f64, 0.00000000000000e+00f64, 2.12224115039342e-01f64,
  0.00000000000000e+00f64, -8.16362644507898e-01f64, 0.00000000000000e+00f64,
  5.35648426691481e-01f64, 0.00000000000000e+00f64, 1.52417902753662e-03f64,
  0.00000000000000e+00f64, 2.63058840550873e-03f64, 0.00000000000000e+00f64,
  4.15292106318904e-03f64, 0.00000000000000e+00f64, 6.97106011119775e-03f64,
  0.00000000000000e+00f64, 1.35535708431058e-02f64, 0.00000000000000e+00f64,
  3.52132898424856e-02f64, 0.00000000000000e+00f64, 2.06946714741884e-01f64,
  0.00000000000000e+00f64, -8.15674251283876e-01f64, 0.00000000000000e+00f64,
  5.38841175520580e-01f64, 0.00000000000000e+00f64, 1.84909689577590e-04f64,
  0.00000000000000e+00f64, 2.90936325007499e-04f64, 0.00000000000000e+00f64,
  3.84877750950089e-04f64, 0.00000000000000e+00f64, 4.86436656735046e-04f64,
  0.00000000000000e+00f64, 6.08688640346879e-04f64, 0.00000000000000e+00f64,
  7.66732830740331e-04f64, 0.00000000000000e+00f64, 9.82753336104205e-04f64,
  0.00000000000000e+00f64, 1.29359957505615e-03f64, 0.00000000000000e+00f64,
  1.76616363801885e-03f64, 0.00000000000000e+00f64, 2.53323433039089e-03f64,
  0.00000000000000e+00f64, 3.88872172121956e-03f64, 0.00000000000000e+00f64,
  6.58635106468291e-03f64, 0.00000000000000e+00f64, 1.30326736343254e-02f64,
  0.00000000000000e+00f64, 3.44353850696714e-02f64, 0.00000000000000e+00f64,
  2.05025409531915e-01f64, 0.00000000000000e+00f64, -8.14985893995401e-01f64,
  0.00000000000000e+00f64, 5.40679930965238e-01f64
];

static Lalpha : [f64; 33] = [
  5.77350269189626e-01f64, 5.16397779494322e-01f64, 5.07092552837110e-01f64,
  5.03952630678970e-01f64, 5.02518907629606e-01f64, 5.01745206004255e-01f64,
  5.01280411827603e-01f64, 5.00979432868120e-01f64, 5.00773395667191e-01f64,
  5.00626174321759e-01f64, 5.00517330712619e-01f64, 5.00434593736979e-01f64,
  5.00370233297676e-01f64, 5.00319182924304e-01f64, 5.00278009473803e-01f64,
  5.00244319584578e-01f64, 5.00216403386025e-01f64, 5.00193012939056e-01f64,
  5.00173220168024e-01f64, 5.00156323280355e-01f64, 5.00141783641018e-01f64,
  5.00129182278347e-01f64, 5.00118189340972e-01f64, 5.00108542278496e-01f64,
  5.00100030010004e-01f64, 5.00092481273333e-01f64, 5.00085755939229e-01f64,
  5.00079738458365e-01f64, 5.00074332862969e-01f64, 5.00069458915387e-01f64,
  5.00065049112355e-01f64, 5.00061046334395e-01f64, 5.00057401986298e-01f64
];

static Lgamma : [f64; 33] = [
  0f64, 0f64, 5.16397779494322e-01f64, 5.07092552837110e-01f64, 5.03952630678970e-01f64,
  5.02518907629606e-01f64, 5.01745206004255e-01f64, 5.01280411827603e-01f64,
  5.00979432868120e-01f64, 5.00773395667191e-01f64, 5.00626174321759e-01f64,
  5.00517330712619e-01f64, 5.00434593736979e-01f64, 5.00370233297676e-01f64,
  5.00319182924304e-01f64, 5.00278009473803e-01f64, 5.00244319584578e-01f64,
  5.00216403386025e-01f64, 5.00193012939056e-01f64, 5.00173220168024e-01f64,
  5.00156323280355e-01f64, 5.00141783641018e-01f64, 5.00129182278347e-01f64,
  5.00118189340972e-01f64, 5.00108542278496e-01f64, 5.00100030010003e-01f64,
  5.00092481273333e-01f64, 5.00085755939229e-01f64, 5.00079738458365e-01f64,
  5.00074332862969e-01f64, 5.00069458915387e-01f64, 5.00065049112355e-01f64,
  5.00061046334395e-01f64
];

static V1inv : [f64; 25] = [
 0.47140452079103168293e-1f64, 0.37712361663282534635f64, 0.56568542494923801952f64,
 0.37712361663282534635f64,0.47140452079103168293e-1f64,
  -0.81649658092772603273e-1f64, -0.46188021535170061160f64, 0f64,
 0.46188021535170061160f64,0.81649658092772603273e-1f64,0.15058465048420853962f64,
 0.12046772038736683169f64, -0.54210474174315074262f64,0.12046772038736683169f64,
 0.15058465048420853962f64, -0.21380899352993950775f64,0.30237157840738178177f64, -0f64,
  -0.30237157840738178177f64,0.21380899352993950775f64,0.10774960475223581324f64,
  -0.21549920950447162648f64,0.21549920950447162648f64, -0.21549920950447162648f64,
 0.10774960475223581324f64
];

static V2inv : [f64; 81] = [
 0.11223917161691230546e-1f64,0.10339219839658349826f64,0.19754094204576565761f64,
 0.25577315077753587922f64,0.27835314560994251755f64,0.25577315077753587922f64,
 0.19754094204576565761f64,0.10339219839658349826f64,0.11223917161691230546e-1f64,
  -0.19440394783993476970e-1f64, -0.16544884625069155470f64, -0.24193725566041460608f64,
  -0.16953338808305493604f64, 0f64,0.16953338808305493604f64,0.24193725566041460608f64,
 0.16544884625069155470f64,0.19440394783993476970e-1f64,0.26466393115406349388e-1f64,
 0.17766815796285469394f64,0.11316664642449611462f64, -0.16306601003711325980f64,
  -0.30847037493128779631f64, -0.16306601003711325980f64,0.11316664642449611462f64,
 0.17766815796285469394f64,0.26466393115406349388e-1f64,
  -0.32395302049990834508e-1f64, -0.15521142532414866547f64,
 0.88573492664788602740e-1f64,0.29570405784974857322f64, 0f64,
  -0.29570405784974857322f64, -0.88573492664788602740e-1f64,0.15521142532414866547f64,
 0.32395302049990834508e-1f64,0.41442155673936851246e-1f64,
 0.98186757907405608245e-1f64, -0.23056908429499411784f64,
  -0.68047008326360625520e-1f64,0.31797435808002456774f64,
  -0.68047008326360625520e-1f64, -0.23056908429499411784f64,
 0.98186757907405608245e-1f64,0.41442155673936851246e-1f64,
  -0.49981120317798783134e-1f64, -0.24861810572835756217e-1f64,
 0.23561326072010832539f64, -0.24472785656448415351f64, 0f64,0.24472785656448415351f64,
  -0.23561326072010832539f64,0.24861810572835756217e-1f64,
 0.49981120317798783134e-1f64,0.79691635865674781228e-1f64,
  -0.95725617891693941833e-1f64, -0.57957553356854386344e-1f64,
 0.21164072460540271452f64, -0.27529837844505833514f64,0.21164072460540271452f64,
  -0.57957553356854386344e-1f64, -0.95725617891693941833e-1f64,
 0.79691635865674781228e-1f64,
  -0.10894869830716590913f64,0.20131094491947531782f64, -0.15407672674888869038f64,
 0.83385723639789791384e-1f64, 0f64, -0.83385723639789791384e-1f64,
 0.15407672674888869038f64, -0.20131094491947531782f64,0.10894869830716590913f64,
 0.54581057089643838221e-1f64, -0.10916211417928767644f64,0.10916211417928767644f64,
  -0.10916211417928767644f64,0.10916211417928767644f64, -0.10916211417928767644f64,
 0.10916211417928767644f64, -0.10916211417928767644f64,0.54581057089643838221e-1f64
];

static V3inv : [f64; 289] = [
 0.27729677693590098996e-2f64,0.26423663180333065153e-1f64,
 0.53374068493933898312e-1f64,0.77007854739523195947e-1f64,
 0.98257061072911596869e-1f64,0.11538049741786835604f64,0.12832134344120884559f64,
 0.13612785914022865001f64,0.13888293186236181317f64,0.13612785914022865001f64,
 0.12832134344120884559f64,0.11538049741786835604f64,0.98257061072911596869e-1f64,
 0.77007854739523195947e-1f64,0.53374068493933898312e-1f64,
 0.26423663180333065153e-1f64,0.27729677693590098996e-2f64,
  -0.48029210642807413690e-2f64, -0.44887724635478800254e-1f64,
  -0.85409520147301089416e-1f64, -0.11090267822061423050f64, -0.12033983162705862441f64,
  -0.11102786862182788886f64, -0.85054870109799336515e-1f64,
  -0.45998467987742225160e-1f64, 0f64,0.45998467987742225160e-1f64,
 0.85054870109799336515e-1f64,0.11102786862182788886f64,0.12033983162705862441f64,
 0.11090267822061423050f64,0.85409520147301089416e-1f64,0.44887724635478800254e-1f64,
 0.48029210642807413690e-2f64,0.62758546879582030087e-2f64,
 0.55561297093529155869e-1f64,
 0.93281491021051539742e-1f64,0.92320151237493695139e-1f64,
 0.55077987469605684531e-1f64,
  -0.96998141716497488255e-2f64, -0.80285961895427405567e-1f64,
  -0.13496839655913850224f64,
  -0.15512521776684524331f64, -0.13496839655913850224f64, -0.80285961895427405567e-1f64,
  -0.96998141716497488255e-2f64,0.55077987469605684531e-1f64,
 0.92320151237493695139e-1f64,0.93281491021051539742e-1f64,
 0.55561297093529155869e-1f64,0.62758546879582030087e-2f64,
  -0.74850969394858555939e-2f64, -0.61751608943839234096e-1f64,
  -0.82974150437304275958e-1f64, -0.38437763431942633378e-1f64,
 0.45745502025779701366e-1f64,0.12369235652734542162f64,0.14720439712852868239f64,
 0.98768034347019704401e-1f64, 0f64,
  -0.98768034347019704401e-1f64, -0.14720439712852868239f64, -0.12369235652734542162f64,
  -0.45745502025779701366e-1f64,0.38437763431942633378e-1f64,
 0.82974150437304275958e-1f64,0.61751608943839234096e-1f64,
 0.74850969394858555939e-2f64,0.86710099994384056338e-2f64,
 0.64006230103659573344e-1f64,0.58517426396091675690e-1f64,
  -0.29743410528985802680e-1f64,
  -0.11934127779157114754f64, -0.12686773515361299409f64, -0.30729137153877447035e-1f64,
 0.97307836256600731568e-1f64,0.15635811574451401023f64,0.97307836256600731568e-1f64,
  -0.30729137153877447035e-1f64, -0.12686773515361299409f64, -0.11934127779157114754f64,
  -0.29743410528985802680e-1f64,0.58517426396091675690e-1f64,
 0.64006230103659573344e-1f64,0.86710099994384056338e-2f64,
  -0.97486395666294840165e-2f64, -0.62995604908060224672e-1f64,
  -0.24373234450275529219e-1f64,0.87760984413626872730e-1f64,
 0.12205204576993351394f64,
 0.16216004196864002088e-1f64, -0.12422320942156845775f64, -0.13682714580929614678f64,
  0f64,0.13682714580929614678f64,0.12422320942156845775f64,
  -0.16216004196864002088e-1f64, -0.12205204576993351394f64,
  -0.87760984413626872730e-1f64,0.24373234450275529219e-1f64,
 0.62995604908060224672e-1f64,0.97486395666294840165e-2f64,
 0.10956271233750488468e-1f64,0.58613204255294358939e-1f64,
  -0.13306063940736618859e-1f64, -0.11606666444978454399f64,
  -0.52059598001115805639e-1f64,0.10868540217796151849f64,0.12594452879014618005f64,
  -0.44678658254872910434e-1f64, -0.15617684362128533405f64,
  -0.44678658254872910434e-1f64,0.12594452879014618005f64,0.10868540217796151849f64,
  -0.52059598001115805639e-1f64, -0.11606666444978454399f64,
  -0.13306063940736618859e-1f64,0.58613204255294358939e-1f64,
 0.10956271233750488468e-1f64, -0.12098893000863087230e-1f64,
  -0.51626244709126208453e-1f64,0.48919433304746979330e-1f64,
 0.10467644465949427090f64,
  -0.48729879523084673782e-1f64, -0.13668732103524749234f64,0.28190838706814496438e-1f64,
 0.15434223333238741600f64, 0f64, -0.15434223333238741600f64,
  -0.28190838706814496438e-1f64,0.13668732103524749234f64,
 0.48729879523084673782e-1f64, -0.10467644465949427090f64,
  -0.48919433304746979330e-1f64,0.51626244709126208453e-1f64,
 0.12098893000863087230e-1f64,0.13542668300437944822e-1f64,
 0.41712033418258689308e-1f64,
  -0.76190463272803434388e-1f64, -0.58303943170068132010e-1f64,0.12158068748245606853f64,
 0.42121099930651007882e-1f64, -0.14684425840766337756f64,
  -0.16108203535058647043e-1f64,0.15698075850757976092f64,
  -0.16108203535058647043e-1f64, -0.14684425840766337756f64,
 0.42121099930651007882e-1f64,0.12158068748245606853f64,
  -0.58303943170068132010e-1f64, -0.76190463272803434388e-1f64,
 0.41712033418258689308e-1f64,0.13542668300437944822e-1f64,
  -0.14939634995117694417e-1f64, -0.30047246373341564039e-1f64,
 0.91624635082546425678e-1f64, -0.79133374319110026377e-2f64,
  -0.12292558212072233355f64,0.90013382617762643524e-1f64,
 0.84013717196539593395e-1f64, -0.14813033309980695856f64, 0f64,
 0.14813033309980695856f64, -0.84013717196539593395e-1f64,
  -0.90013382617762643524e-1f64,
 0.12292558212072233355f64,0.79133374319110026377e-2f64, -0.91624635082546425678e-1f64,
 0.30047246373341564039e-1f64,0.14939634995117694417e-1f64,
 0.16986031342807474208e-1f64,
 0.15760203882617033601e-1f64, -0.91494054040950941996e-1f64,
 0.70082459207876130806e-1f64,
 0.53390713710144539104e-1f64, -0.14340746778352039430f64,0.84048122493418898508e-1f64,
 0.72456667788091316868e-1f64, -0.15564535320096811360f64,
 0.72456667788091316868e-1f64,0.84048122493418898508e-1f64,
  -0.14340746778352039430f64,0.53390713710144539104e-1f64,
 0.70082459207876130806e-1f64, -0.91494054040950941996e-1f64,
 0.15760203882617033601e-1f64,
 0.16986031342807474208e-1f64, -0.18994065631858742028e-1f64,
  -0.82901821370405592927e-3f64,0.77239669773015192888e-1f64,
  -0.10850735431039424680f64,0.47524484622086496464e-1f64,
 0.69148184871588737021e-1f64, -0.14829314646228194928f64,0.11992057742398672066f64,
  0f64, -0.11992057742398672066f64,0.14829314646228194928f64,
  -0.69148184871588737021e-1f64, -0.47524484622086496464e-1f64,
 0.10850735431039424680f64, -0.77239669773015192888e-1f64,
 0.82901821370405592927e-3f64,0.18994065631858742028e-1f64,
 0.22761703826371535132e-1f64, -0.17728848711449643358e-1f64,
  -0.47496371572480503788e-1f64,0.10659958402328690063f64, -0.11696013966166296514f64,
 0.63073750910894244526e-1f64,0.32928881123602721303e-1f64,
  -0.12280950532497593683f64,0.15926189077282729505f64, -0.12280950532497593683f64,
 0.32928881123602721303e-1f64,0.63073750910894244526e-1f64,
  -0.11696013966166296514f64,0.10659958402328690063f64, -0.47496371572480503788e-1f64,
  -0.17728848711449643358e-1f64,0.22761703826371535132e-1f64,
  -0.26493215276042203434e-1f64,0.35579780856128386192e-1f64,
 0.10447309718398935122e-1f64, -0.68616154085314996709e-1f64,
 0.11775363082763954214f64, -0.13918901977011837274f64,0.12312819418827395690f64,
  -0.72053565748259077905e-1f64, 0f64,0.72053565748259077905e-1f64,
  -0.12312819418827395690f64,0.13918901977011837274f64, -0.11775363082763954214f64,
 0.68616154085314996709e-1f64, -0.10447309718398935122e-1f64,
  -0.35579780856128386192e-1f64,
 0.26493215276042203434e-1f64,0.40742523354399706918e-1f64,
  -0.73124912999529117195e-1f64,0.49317266444153837821e-1f64,
  -0.13686605413876015320e-1f64, -0.28342624942191100464e-1f64,
 0.70371855298258216249e-1f64, -0.10600251632853603875f64,0.12981016288391131812f64,
  -0.13817029659318161476f64,0.12981016288391131812f64, -0.10600251632853603875f64,
 0.70371855298258216249e-1f64, -0.28342624942191100464e-1f64,
  -0.13686605413876015320e-1f64,
 0.49317266444153837821e-1f64, -0.73124912999529117195e-1f64,
 0.40742523354399706918e-1f64, -0.54944368958699908688e-1f64,
 0.10777725663147408190f64, -0.10152395581538265428f64,0.91369146312596428468e-1f64,
  -0.77703071757424700773e-1f64,0.61050911730999815031e-1f64,
  -0.42052599404498348871e-1f64,0.21438229266251454773e-1f64, 0f64,
  -0.21438229266251454773e-1f64,0.42052599404498348871e-1f64,
  -0.61050911730999815031e-1f64,0.77703071757424700773e-1f64,
  -0.91369146312596428468e-1f64,
 0.10152395581538265428f64, -0.10777725663147408190f64,0.54944368958699908688e-1f64,
 0.27485608464748840573e-1f64, -0.54971216929497681146e-1f64,
 0.54971216929497681146e-1f64,
  -0.54971216929497681146e-1f64,0.54971216929497681146e-1f64,
  -0.54971216929497681146e-1f64,0.54971216929497681146e-1f64,
  -0.54971216929497681146e-1f64,0.54971216929497681146e-1f64,
  -0.54971216929497681146e-1f64,0.54971216929497681146e-1f64,
  -0.54971216929497681146e-1f64,0.54971216929497681146e-1f64,
  -0.54971216929497681146e-1f64,0.54971216929497681146e-1f64,
  -0.54971216929497681146e-1f64,0.27485608464748840573e-1f64
];

static V4inv : [f64; 1089] = [
 0.69120897476690862600e-3f64,0.66419939766331555194e-2f64,
 0.13600665164323186111e-1f64,0.20122785860913684493e-1f64,
 0.26583214101668429944e-1f64,0.32712713318999268739e-1f64,
 0.38576221976287138036e-1f64,0.44033030938268925133e-1f64,
 0.49092709529622799673e-1f64,0.53657949874312515646e-1f64,
 0.57724533144734311859e-1f64,0.61219564530655179096e-1f64,
 0.64138907503837875026e-1f64,0.66427905189318792009e-1f64,
 0.68088956652280022887e-1f64,0.69083051391555695878e-1f64,
 0.69422738116739271449e-1f64,0.69083051391555695878e-1f64,
 0.68088956652280022887e-1f64,0.66427905189318792009e-1f64,
 0.64138907503837875026e-1f64,0.61219564530655179096e-1f64,
 0.57724533144734311859e-1f64,0.53657949874312515646e-1f64,
 0.49092709529622799673e-1f64,0.44033030938268925133e-1f64,
 0.38576221976287138036e-1f64,0.32712713318999268739e-1f64,
 0.26583214101668429944e-1f64,0.20122785860913684493e-1f64,
 0.13600665164323186111e-1f64,0.66419939766331555194e-2f64,
 0.69120897476690862600e-3f64, -0.11972090629438798134e-2f64,
  -0.11448874821643225573e-1f64, -0.23104401104002905904e-1f64,
  -0.33352899418646530133e-1f64, -0.42538626424075425908e-1f64,
  -0.49969730733911825941e-1f64, -0.55555454015360728353e-1f64,
  -0.58955533624852604918e-1f64, -0.60126044219122513907e-1f64,
  -0.58959430451175833624e-1f64, -0.55546925396227130606e-1f64,
  -0.49984739749347973762e-1f64, -0.42513009141170294365e-1f64,
  -0.33399140950669746346e-1f64, -0.23007690803851790829e-1f64,
  -0.11728275717520066169e-1f64, 0f64,0.11728275717520066169e-1f64,
 0.23007690803851790829e-1f64,0.33399140950669746346e-1f64,
 0.42513009141170294365e-1f64,0.49984739749347973762e-1f64,
 0.55546925396227130606e-1f64,0.58959430451175833624e-1f64,
 0.60126044219122513907e-1f64,0.58955533624852604918e-1f64,
 0.55555454015360728353e-1f64,0.49969730733911825941e-1f64,
 0.42538626424075425908e-1f64,0.33352899418646530133e-1f64,
 0.23104401104002905904e-1f64,0.11448874821643225573e-1f64,
 0.11972090629438798134e-2f64,0.15501585012936019146e-2f64,
 0.14628781502199620482e-1f64,0.28684915921474815271e-1f64,
 0.39299396074628048026e-1f64,0.46393418975496284204e-1f64,
 0.48756902531094699526e-1f64,0.46331333488337494692e-1f64,
 0.39012645376980228775e-1f64,0.27452795421085791153e-1f64,
 0.12430953621169863781e-1f64, -0.47682978056024928800e-2f64,
  -0.22825828045428973853e-1f64,
  -0.40195512090720278312e-1f64, -0.55503004262826221955e-1f64,
  -0.67424537752827046308e-1f64, -0.75020199300113606452e-1f64,
  -0.77607844312483656131e-1f64, -0.75020199300113606452e-1f64,
  -0.67424537752827046308e-1f64, -0.55503004262826221955e-1f64,
  -0.40195512090720278312e-1f64, -0.22825828045428973853e-1f64,
  -0.47682978056024928800e-2f64,0.12430953621169863781e-1f64,
 0.27452795421085791153e-1f64,0.39012645376980228775e-1f64,
 0.46331333488337494692e-1f64,0.48756902531094699526e-1f64,
 0.46393418975496284204e-1f64,0.39299396074628048026e-1f64,
 0.28684915921474815271e-1f64,0.14628781502199620482e-1f64,
 0.15501585012936019146e-2f64, -0.18377757558949194214e-2f64,
  -0.17050470050949761565e-1f64, -0.31952119564923250836e-1f64,
  -0.40197423449026348155e-1f64,
  -0.41205649520281371624e-1f64, -0.33909965817492272248e-1f64,
  -0.19393664422115332144e-1f64,0.56661049630886784692e-3f64,
 0.22948272173686561721e-1f64,0.44489719570904738207e-1f64,
 0.61790363672287920596e-1f64,0.72121014727028013894e-1f64,
 0.73627151185287858579e-1f64,0.65784665375961398923e-1f64,
 0.49369676372333667559e-1f64,0.26444326317059715065e-1f64, 0f64,
  -0.26444326317059715065e-1f64, -0.49369676372333667559e-1f64,
  -0.65784665375961398923e-1f64, -0.73627151185287858579e-1f64,
  -0.72121014727028013894e-1f64, -0.61790363672287920596e-1f64,
  -0.44489719570904738207e-1f64, -0.22948272173686561721e-1f64,
  -0.56661049630886784692e-3f64,0.19393664422115332144e-1f64,
 0.33909965817492272248e-1f64,0.41205649520281371624e-1f64,
 0.40197423449026348155e-1f64,0.31952119564923250836e-1f64,
 0.17050470050949761565e-1f64,0.18377757558949194214e-2f64,
 0.20942714740729767769e-2f64,0.18935902405146518232e-1f64,
 0.33335840852491735126e-1f64,0.36770680999102286065e-1f64,
 0.28873194534132768509e-1f64,0.10267303017729535513e-1f64,
  -0.14607738306201572890e-1f64, -0.40139568545572305818e-1f64,
  -0.59808326733858291561e-1f64, -0.68528358823372627506e-1f64,
  -0.63306535387619244879e-1f64, -0.44508601817574921056e-1f64,
  -0.15449116105605395357e-1f64,0.17941083795006546367e-1f64,
 0.48747356011657242123e-1f64,0.70329553984201665523e-1f64,
 0.78106117292526169663e-1f64,0.70329553984201665523e-1f64,
 0.48747356011657242123e-1f64,0.17941083795006546367e-1f64,
  -0.15449116105605395357e-1f64, -0.44508601817574921056e-1f64,
  -0.63306535387619244879e-1f64, -0.68528358823372627506e-1f64,
  -0.59808326733858291561e-1f64,
  -0.40139568545572305818e-1f64, -0.14607738306201572890e-1f64,
 0.10267303017729535513e-1f64,0.28873194534132768509e-1f64,
 0.36770680999102286065e-1f64,0.33335840852491735126e-1f64,
 0.18935902405146518232e-1f64,0.20942714740729767769e-2f64,
  -0.23245285491878278419e-2f64, -0.20401404737639389919e-1f64,
  -0.33019548231022514097e-1f64, -0.29709828426463720091e-1f64,
  -0.11760070922697422156e-1f64,0.15987584743850393793e-1f64,
 0.43619012891472813485e-1f64,0.61177322409671487721e-1f64,
 0.61144030218486655594e-1f64,
 0.41895377620089086167e-1f64,0.80232011820644308033e-2f64,
  -0.30574701186675900915e-1f64,
  -0.62072243008844865848e-1f64, -0.76336186183574765586e-1f64,
  -0.68435466095345537115e-1f64, -0.40237669208466966207e-1f64, 0f64,
 0.40237669208466966207e-1f64,0.68435466095345537115e-1f64,
 0.76336186183574765586e-1f64,0.62072243008844865848e-1f64,
 0.30574701186675900915e-1f64, -0.80232011820644308033e-2f64,
  -0.41895377620089086167e-1f64, -0.61144030218486655594e-1f64,
  -0.61177322409671487721e-1f64, -0.43619012891472813485e-1f64,
  -0.15987584743850393793e-1f64,0.11760070922697422156e-1f64,
 0.29709828426463720091e-1f64,0.33019548231022514097e-1f64,
 0.20401404737639389919e-1f64,0.23245285491878278419e-2f64,
 0.25451717261579269307e-2f64,0.21480418595666878775e-1f64,
 0.31177212469293007998e-1f64,0.19816333607013379373e-1f64,
  -0.72439496274458793681e-2f64, -0.38404203906598342397e-1f64,
  -0.57633632255322221046e-1f64, -0.54070547403585392952e-1f64,
  -0.26249823354368866005e-1f64,0.15643058212336881516e-1f64,
 0.54539832735118677194e-1f64,0.73283028002473989724e-1f64,
 0.62835303524135936213e-1f64,0.26175977027801048141e-1f64,
  -0.22193636309998606610e-1f64, -0.62597049956093311234e-1f64,
  -0.78206986173170212505e-1f64, -0.62597049956093311234e-1f64,
  -0.22193636309998606610e-1f64,0.26175977027801048141e-1f64,
 0.62835303524135936213e-1f64,
 0.73283028002473989724e-1f64,0.54539832735118677194e-1f64,
 0.15643058212336881516e-1f64,
  -0.26249823354368866005e-1f64, -0.54070547403585392952e-1f64,
  -0.57633632255322221046e-1f64, -0.38404203906598342397e-1f64,
  -0.72439496274458793681e-2f64,0.19816333607013379373e-1f64,
 0.31177212469293007998e-1f64,0.21480418595666878775e-1f64,
 0.25451717261579269307e-2f64, -0.27506573922483820005e-2f64,
  -0.22224442095099251870e-1f64, -0.27949927254215773020e-1f64,
  -0.80918481053370034987e-2f64,0.25121859354449306916e-1f64,
 0.51563535009373061074e-1f64,0.51936965107145960512e-1f64,
 0.22146626648171527753e-1f64,
  -0.24172689882103382748e-1f64, -0.61731229104853568296e-1f64,
  -0.68477262429344201201e-1f64, -0.38311232728303704742e-1f64,
 0.14160578713659552679e-1f64,0.61248813427564184033e-1f64,
 0.77136328841293031805e-1f64,0.52514801765183697988e-1f64, 0f64,
  -0.52514801765183697988e-1f64, -0.77136328841293031805e-1f64,
  -0.61248813427564184033e-1f64, -0.14160578713659552679e-1f64,
 0.38311232728303704742e-1f64,
 0.68477262429344201201e-1f64,0.61731229104853568296e-1f64,
 0.24172689882103382748e-1f64,
  -0.22146626648171527753e-1f64, -0.51936965107145960512e-1f64,
  -0.51563535009373061074e-1f64, -0.25121859354449306916e-1f64,
 0.80918481053370034987e-2f64,0.27949927254215773020e-1f64,
 0.22224442095099251870e-1f64,0.27506573922483820005e-2f64,
 0.29562461131654311467e-2f64,0.22630271480554450613e-1f64,
 0.23547399831373800971e-1f64, -0.43964593440902476642e-2f64,
  -0.39055315767504970597e-1f64, -0.52369643937940066804e-1f64,
  -0.28506131614971613422e-1f64,0.19906048093338832322e-1f64,
 0.60408880866392420279e-1f64,0.62493397473656883090e-1f64,
 0.21391278377641297859e-1f64, -0.37302864786623254746e-1f64,
  -0.73665127933539496872e-1f64, -0.61706142476854010202e-1f64,
  -0.78065168882546327888e-2f64,0.52335307373945544428e-1f64,
 0.78278746279419264777e-1f64,0.52335307373945544428e-1f64,
  -0.78065168882546327888e-2f64, -0.61706142476854010202e-1f64,
  -0.73665127933539496872e-1f64, -0.37302864786623254746e-1f64,
 0.21391278377641297859e-1f64,0.62493397473656883090e-1f64,
 0.60408880866392420279e-1f64,0.19906048093338832322e-1f64,
  -0.28506131614971613422e-1f64, -0.52369643937940066804e-1f64,
  -0.39055315767504970597e-1f64, -0.43964593440902476642e-2f64,
 0.23547399831373800971e-1f64,0.22630271480554450613e-1f64,
 0.29562461131654311467e-2f64, -0.31515718415504761303e-2f64,
  -0.22739451096655080673e-1f64, -0.18157123602272119779e-1f64,
 0.16496480897167303621e-1f64,0.46921166788569301124e-1f64,
 0.40644395739978416354e-1f64, -0.46275803430732216900e-2f64,
  -0.52883375891308909486e-1f64, -0.61116483226324111734e-1f64,
  -0.17411698764545629853e-1f64,0.44773430013166822765e-1f64,
 0.73441577962383869198e-1f64,0.42127368371995472815e-1f64,
  -0.25504645957196772465e-1f64, -0.74126818045972742488e-1f64,
  -0.62780077864719287317e-1f64, 0f64,0.62780077864719287317e-1f64,
 0.74126818045972742488e-1f64,0.25504645957196772465e-1f64,
  -0.42127368371995472815e-1f64, -0.73441577962383869198e-1f64,
  -0.44773430013166822765e-1f64,0.17411698764545629853e-1f64,
 0.61116483226324111734e-1f64,0.52883375891308909486e-1f64,
 0.46275803430732216900e-2f64, -0.40644395739978416354e-1f64,
  -0.46921166788569301124e-1f64, -0.16496480897167303621e-1f64,
 0.18157123602272119779e-1f64,0.22739451096655080673e-1f64,
 0.31515718415504761303e-2f64,0.33536559294882188208e-2f64,
 0.22535348942792006185e-1f64,
 0.12048629300953560767e-1f64, -0.27166076791299493403e-1f64,
  -0.47492745604230978367e-1f64, -0.19246623430993153174e-1f64,
 0.36231297307556299322e-1f64,0.61713617181636122004e-1f64,
 0.25928029734266134490e-1f64, -0.40478700752883602818e-1f64,
  -0.71053889866326412049e-1f64, -0.31870824482961751482e-1f64,
 0.41515251100219081281e-1f64,0.76481960760098381651e-1f64,
 0.36726509155999912440e-1f64, -0.40090067032627055969e-1f64,
  -0.78270742903374539397e-1f64, -0.40090067032627055969e-1f64,
 0.36726509155999912440e-1f64,0.76481960760098381651e-1f64,
 0.41515251100219081281e-1f64, -0.31870824482961751482e-1f64,
  -0.71053889866326412049e-1f64, -0.40478700752883602818e-1f64,
 0.25928029734266134490e-1f64,0.61713617181636122004e-1f64,
 0.36231297307556299322e-1f64, -0.19246623430993153174e-1f64,
  -0.47492745604230978367e-1f64, -0.27166076791299493403e-1f64,
 0.12048629300953560767e-1f64,0.22535348942792006185e-1f64,
 0.33536559294882188208e-2f64,
  -0.35481220456925318865e-2f64, -0.22062913693073191150e-1f64,
  -0.54487362861834144999e-2f64,0.35438821865804087489e-1f64,
 0.40733077820527411302e-1f64, -0.67403098138950720914e-2f64,
  -0.55559584405239171054e-1f64, -0.42417050790865158745e-1f64,
 0.24499901971884704925e-1f64,0.68721232891705409302e-1f64,
 0.34086082787461126592e-1f64, -0.43441000373118474002e-1f64,
  -0.73878085292669148950e-1f64, -0.18846995664706657127e-1f64,
 0.59827776178286834498e-1f64,0.70644634584085901794e-1f64, 0f64,
  -0.70644634584085901794e-1f64, -0.59827776178286834498e-1f64,
 0.18846995664706657127e-1f64,0.73878085292669148950e-1f64,
 0.43441000373118474002e-1f64, -0.34086082787461126592e-1f64,
  -0.68721232891705409302e-1f64, -0.24499901971884704925e-1f64,
 0.42417050790865158745e-1f64,0.55559584405239171054e-1f64,
 0.67403098138950720914e-2f64, -0.40733077820527411302e-1f64,
  -0.35438821865804087489e-1f64,0.54487362861834144999e-2f64,
 0.22062913693073191150e-1f64,0.35481220456925318865e-2f64,
 0.37554176816665075631e-2f64,0.21297045781589919482e-1f64,
  -0.13327293083183431816e-2f64,
  -0.40635299172764596484e-1f64, -0.27659860508374175359e-1f64,
 0.31089232744083445986e-1f64,0.56113781541334176109e-1f64,
 0.37577840643257763400e-2f64, -0.60511227350664590865e-1f64,
  -0.46670556446129053853e-1f64,0.33263195878575888247e-1f64,
 0.72757324720645228775e-1f64,0.15011712351692283635e-1f64,
  -0.65601212994924119078e-1f64, -0.60016855838843789772e-1f64,
 0.26220858553188665966e-1f64,0.78322776605833552980e-1f64,
 0.26220858553188665966e-1f64, -0.60016855838843789772e-1f64,
  -0.65601212994924119078e-1f64,
 0.15011712351692283635e-1f64,0.72757324720645228775e-1f64,
 0.33263195878575888247e-1f64,
  -0.46670556446129053853e-1f64, -0.60511227350664590865e-1f64,
 0.37577840643257763400e-2f64,0.56113781541334176109e-1f64,
 0.31089232744083445986e-1f64, -0.27659860508374175359e-1f64,
  -0.40635299172764596484e-1f64, -0.13327293083183431816e-2f64,
 0.21297045781589919482e-1f64,0.37554176816665075631e-2f64,
  -0.39566995305720591229e-2f64, -0.20291873414438919995e-1f64,
 0.80617453830770930551e-2f64,0.42270189157016547906e-1f64,
 0.10332624526759093004e-1f64, -0.48054759547616142024e-1f64,
  -0.37678032941171643972e-1f64,
 0.36617192625732482394e-1f64,0.61009425973424865714e-1f64,
  -0.95589113168026591466e-2f64,
  -0.71023202645076922361e-1f64, -0.25097788086808784456e-1f64,
 0.62406621963267050244e-1f64,0.56907293171100693511e-1f64,
  -0.36435383083882206257e-1f64, -0.75790105119208756348e-1f64, 0f64,
 0.75790105119208756348e-1f64,0.36435383083882206257e-1f64,
  -0.56907293171100693511e-1f64, -0.62406621963267050244e-1f64,
 0.25097788086808784456e-1f64,0.71023202645076922361e-1f64,
 0.95589113168026591466e-2f64,
  -0.61009425973424865714e-1f64, -0.36617192625732482394e-1f64,
 0.37678032941171643972e-1f64,0.48054759547616142024e-1f64,
  -0.10332624526759093004e-1f64, -0.42270189157016547906e-1f64,
  -0.80617453830770930551e-2f64,0.20291873414438919995e-1f64,
 0.39566995305720591229e-2f64,0.41776092289182138591e-2f64,
 0.19013221163904414395e-1f64, -0.14420609729849899876e-1f64,
  -0.40259160586844441220e-1f64,0.86327811113710831649e-2f64,
 0.53564430703021034399e-1f64,0.65469185402150431933e-2f64,
  -0.60383116311280629856e-1f64,
  -0.25657793784058876939e-1f64,0.58745680576829226900e-1f64,
 0.45649937869034420296e-1f64,
  -0.49167932056844167772e-1f64, -0.62696614328552187977e-1f64,
 0.32540234556426699997e-1f64,0.74280410383464269758e-1f64,
  -0.11425672633410999870e-1f64, -0.78280649404686404903e-1f64,
  -0.11425672633410999870e-1f64,0.74280410383464269758e-1f64,
 0.32540234556426699997e-1f64, -0.62696614328552187977e-1f64,
  -0.49167932056844167772e-1f64,0.45649937869034420296e-1f64,
 0.58745680576829226900e-1f64, -0.25657793784058876939e-1f64,
  -0.60383116311280629856e-1f64,0.65469185402150431933e-2f64,
 0.53564430703021034399e-1f64,
 0.86327811113710831649e-2f64, -0.40259160586844441220e-1f64,
  -0.14420609729849899876e-1f64,0.19013221163904414395e-1f64,
 0.41776092289182138591e-2f64, -0.43935502082478059199e-2f64,
  -0.17528761237509401631e-1f64,0.20208915249153872535e-1f64,
 0.34734743119040669109e-1f64, -0.26275910172353637955e-1f64,
  -0.46368003346018878786e-1f64,
 0.26800056330709381025e-1f64,0.56681476464606609921e-1f64,
  -0.24749011438127255898e-1f64,
  -0.64934612189056658992e-1f64,0.20333742247679279535e-1f64,
 0.71429299070059318651e-1f64,
  -0.14452513210428671266e-1f64, -0.75793341281736586582e-1f64,
 0.74717094137184935270e-2f64,0.78034921554757317374e-1f64, 0f64,
  -0.78034921554757317374e-1f64, -0.74717094137184935270e-2f64,
 0.75793341281736586582e-1f64,0.14452513210428671266e-1f64,
  -0.71429299070059318651e-1f64, -0.20333742247679279535e-1f64,
 0.64934612189056658992e-1f64,0.24749011438127255898e-1f64,
  -0.56681476464606609921e-1f64,
  -0.26800056330709381025e-1f64,0.46368003346018878786e-1f64,
 0.26275910172353637955e-1f64,
  -0.34734743119040669109e-1f64, -0.20208915249153872535e-1f64,
 0.17528761237509401631e-1f64,0.43935502082478059199e-2f64,
 0.46379089482818671473e-2f64,0.15791188144791287229e-1f64,
  -0.25134290048737455284e-1f64, -0.26249795071946841205e-1f64,
 0.39960457575789924651e-1f64,0.28111892450146525404e-1f64,
  -0.51026476400767918226e-1f64,
  -0.27266747278681831364e-1f64,0.60708796647861610865e-1f64,
 0.23532306960642115854e-1f64,
  -0.68169639871532441111e-1f64, -0.18204924701958312032e-1f64,
 0.73822890510656128485e-1f64,0.11373392486424717019e-1f64,
  -0.77133324017644609416e-1f64, -0.39295877480342619961e-2f64,
 0.78351902829418987960e-1f64, -0.39295877480342619961e-2f64,
  -0.77133324017644609416e-1f64,0.11373392486424717019e-1f64,
 0.73822890510656128485e-1f64, -0.18204924701958312032e-1f64,
  -0.68169639871532441111e-1f64,0.23532306960642115854e-1f64,
 0.60708796647861610865e-1f64, -0.27266747278681831364e-1f64,
  -0.51026476400767918226e-1f64,0.28111892450146525404e-1f64,
 0.39960457575789924651e-1f64, -0.26249795071946841205e-1f64,
  -0.25134290048737455284e-1f64,0.15791188144791287229e-1f64,
 0.46379089482818671473e-2f64, -0.48780095920069827068e-2f64,
  -0.13886961667516983541e-1f64,0.29071311049368895844e-1f64,
 0.15480559452075811600e-1f64, -0.47527977686242313065e-1f64,
  -0.31929089844361042178e-2f64,0.58015667638415922967e-1f64,
  -0.14547915466597622925e-1f64, -0.61067668299848923244e-1f64,
 0.35093678009090186851e-1f64,0.55378399159800654657e-1f64,
  -0.54277226474891610385e-1f64, -0.42023830782434076509e-1f64,
 0.69197384645944912066e-1f64,0.22610783557709586445e-1f64,
  -0.77269275900637030185e-1f64, 0f64,0.77269275900637030185e-1f64,
  -0.22610783557709586445e-1f64,
  -0.69197384645944912066e-1f64,0.42023830782434076509e-1f64,
 0.54277226474891610385e-1f64,
  -0.55378399159800654657e-1f64, -0.35093678009090186851e-1f64,
 0.61067668299848923244e-1f64,0.14547915466597622925e-1f64,
  -0.58015667638415922967e-1f64,0.31929089844361042178e-2f64,
 0.47527977686242313065e-1f64, -0.15480559452075811600e-1f64,
  -0.29071311049368895844e-1f64,0.13886961667516983541e-1f64,
 0.48780095920069827068e-2f64,0.51591759101720291381e-2f64,
 0.11747497650231330965e-1f64, -0.31777863364694653331e-1f64,
  -0.34555825499804605557e-2f64,0.47914131921157015198e-1f64,
  -0.22573685920142225247e-1f64, -0.45320344390022666738e-1f64,
 0.49660630547172186418e-1f64,0.25707858143963615736e-1f64,
  -0.68132707341917233933e-1f64,0.67534860185243140399e-2f64,
 0.69268150370037450063e-1f64, -0.41585011920451477177e-1f64,
  -0.51622397460510041271e-1f64,0.68408139576363036148e-1f64,
 0.18981259024768933323e-1f64, -0.78265472429342305554e-1f64,
 0.18981259024768933323e-1f64,0.68408139576363036148e-1f64,
  -0.51622397460510041271e-1f64,
  -0.41585011920451477177e-1f64,0.69268150370037450063e-1f64,
 0.67534860185243140399e-2f64,
  -0.68132707341917233933e-1f64,0.25707858143963615736e-1f64,
 0.49660630547172186418e-1f64,
  -0.45320344390022666738e-1f64, -0.22573685920142225247e-1f64,
 0.47914131921157015198e-1f64, -0.34555825499804605557e-2f64,
  -0.31777863364694653331e-1f64,0.11747497650231330965e-1f64,
 0.51591759101720291381e-2f64, -0.54365757412741340377e-2f64,
  -0.94862516619529080191e-2f64,0.33240472093448190877e-1f64,
  -0.88698898099681552229e-2f64,
  -0.40973252097216337576e-1f64,0.42995673349795657065e-1f64,
 0.17320914507876958783e-1f64,
  -0.62201292691914856803e-1f64,0.24726274174637346693e-1f64,
 0.51320859246515407288e-1f64,
  -0.62882063373810501763e-1f64, -0.11003569131725622672e-1f64,
 0.73842261324108943465e-1f64, -0.39240120294802923208e-1f64,
  -0.49293966443941122807e-1f64,0.73552644778818223475e-1f64, 0f64,
  -0.73552644778818223475e-1f64,0.49293966443941122807e-1f64,
 0.39240120294802923208e-1f64, -0.73842261324108943465e-1f64,
 0.11003569131725622672e-1f64,0.62882063373810501763e-1f64,
  -0.51320859246515407288e-1f64,
  -0.24726274174637346693e-1f64,0.62201292691914856803e-1f64,
  -0.17320914507876958783e-1f64, -0.42995673349795657065e-1f64,
 0.40973252097216337576e-1f64,0.88698898099681552229e-2f64,
  -0.33240472093448190877e-1f64,0.94862516619529080191e-2f64,
 0.54365757412741340377e-2f64,0.57750194549356126240e-2f64,
 0.69981166020044116791e-2f64, -0.33274982140403110792e-1f64,
 0.20297071020698356116e-1f64,0.27898517839646066582e-1f64,
  -0.53368678853282030262e-1f64,0.16656482990394548343e-1f64,
 0.46342901447260614255e-1f64,
  -0.60536796508149003365e-1f64,0.29109107483842596340e-2f64,
 0.63224486124385124504e-1f64,
  -0.59028872851312033411e-1f64, -0.14783105962696191734e-1f64,
 0.74269399241069253865e-1f64, -0.49053677339382384625e-1f64,
  -0.33525466624811186739e-1f64,0.78397349622515386647e-1f64,
  -0.33525466624811186739e-1f64, -0.49053677339382384625e-1f64,
 0.74269399241069253865e-1f64, -0.14783105962696191734e-1f64,
  -0.59028872851312033411e-1f64,
 0.63224486124385124504e-1f64,0.29109107483842596340e-2f64,
  -0.60536796508149003365e-1f64,
 0.46342901447260614255e-1f64,0.16656482990394548343e-1f64,
  -0.53368678853282030262e-1f64,
 0.27898517839646066582e-1f64,0.20297071020698356116e-1f64,
  -0.33274982140403110792e-1f64,
 0.69981166020044116791e-2f64,0.57750194549356126240e-2f64,
  -0.61100308370519200637e-2f64, -0.44383614355738148616e-2f64,
 0.32011283412619094811e-1f64, -0.29965011866372897633e-1f64,
  -0.10560682331349193348e-1f64,0.51110336443392506342e-1f64,
  -0.45012284729681775492e-1f64, -0.94236825555873320102e-2f64,
 0.60860695783141264746e-1f64,
  -0.55014628647083368926e-1f64, -0.73474782382499482121e-2f64,
 0.66640148475243034781e-1f64, -0.62533116045749887988e-1f64,
  -0.38650525912400102585e-2f64,0.68429769005837003777e-1f64,
  -0.66984505412544901945e-1f64, 0f64,0.66984505412544901945e-1f64,
  -0.68429769005837003777e-1f64,0.38650525912400102585e-2f64,
 0.62533116045749887988e-1f64, -0.66640148475243034781e-1f64,
 0.73474782382499482121e-2f64,
 0.55014628647083368926e-1f64, -0.60860695783141264746e-1f64,
 0.94236825555873320102e-2f64,
 0.45012284729681775492e-1f64, -0.51110336443392506342e-1f64,
 0.10560682331349193348e-1f64,
 0.29965011866372897633e-1f64, -0.32011283412619094811e-1f64,
 0.44383614355738148616e-2f64,
 0.61100308370519200637e-2f64,0.65409373892036191538e-2f64,
 0.16350101107071157065e-2f64, -0.29301957285983144319e-1f64,
 0.36838667173388832579e-1f64, -0.81922703976491586393e-2f64,
  -0.36955670021050133434e-1f64,0.58374851095540469865e-1f64,
  -0.31977016246946181856e-1f64, -0.25311073698658094646e-1f64,
 0.66674413950106952577e-1f64,
  -0.54865713324521039571e-1f64, -0.39797027891537985440e-2f64,
 0.62830285264808449064e-1f64, -0.72226313251296100676e-1f64,
 0.22560232697133353980e-1f64,0.46455784709904033738e-1f64,
  -0.78200930751070349956e-1f64,0.46455784709904033738e-1f64,
 0.22560232697133353980e-1f64, -0.72226313251296100676e-1f64,
 0.62830285264808449064e-1f64, -0.39797027891537985440e-2f64,
  -0.54865713324521039571e-1f64,0.66674413950106952577e-1f64,
  -0.25311073698658094646e-1f64, -0.31977016246946181856e-1f64,
 0.58374851095540469865e-1f64, -0.36955670021050133434e-1f64,
  -0.81922703976491586393e-2f64,0.36838667173388832579e-1f64,
  -0.29301957285983144319e-1f64,0.16350101107071157065e-2f64,
 0.65409373892036191538e-2f64, -0.69686180931868703196e-2f64,
 0.11849538727632789870e-2f64,0.25452286414610537766e-1f64,
  -0.40522480651713943230e-1f64,0.25694679053362813183e-1f64,
 0.14057118113748390637e-1f64, -0.52037614725803488893e-1f64,
 0.58849342223684035589e-1f64,
  -0.25075229077361409271e-1f64, -0.29559771094034181083e-1f64,
 0.68296746944165720199e-1f64, -0.62890462146423984955e-1f64,
 0.14457636466274596445e-1f64,0.45787612031322361496e-1f64,
  -0.77231759014655809742e-1f64,0.57881203613910543657e-1f64, 0f64,
  -0.57881203613910543657e-1f64,0.77231759014655809742e-1f64,
  -0.45787612031322361496e-1f64, -0.14457636466274596445e-1f64,
 0.62890462146423984955e-1f64,
  -0.68296746944165720199e-1f64,0.29559771094034181083e-1f64,
 0.25075229077361409271e-1f64,
  -0.58849342223684035589e-1f64,0.52037614725803488893e-1f64,
  -0.14057118113748390637e-1f64, -0.25694679053362813183e-1f64,
 0.40522480651713943230e-1f64, -0.25452286414610537766e-1f64,
  -0.11849538727632789870e-2f64,0.69686180931868703196e-2f64,
 0.75611653617520254845e-2f64, -0.43290610418608409141e-2f64,
  -0.20277062025115566914e-1f64,
 0.40362947027704828926e-1f64, -0.38938808024132120254e-1f64,
 0.11831186195916702262e-1f64,
 0.28476667401744525357e-1f64, -0.59320969056617684621e-1f64,
 0.61101629747436200186e-1f64,
  -0.29514834848355389223e-1f64, -0.20668001885001084821e-1f64,
 0.62923592802445122793e-1f64, -0.73558456263588833115e-1f64,
 0.45314556330160999776e-1f64,0.79031645918426015574e-2f64,
  -0.58136953576334689357e-1f64,0.78538474524006405758e-1f64,
  -0.58136953576334689357e-1f64,0.79031645918426015574e-2f64,
 0.45314556330160999776e-1f64, -0.73558456263588833115e-1f64,
 0.62923592802445122793e-1f64, -0.20668001885001084821e-1f64,
  -0.29514834848355389223e-1f64,0.61101629747436200186e-1f64,
  -0.59320969056617684621e-1f64,0.28476667401744525357e-1f64,
 0.11831186195916702262e-1f64, -0.38938808024132120254e-1f64,
 0.40362947027704828926e-1f64, -0.20277062025115566914e-1f64,
  -0.43290610418608409141e-2f64,0.75611653617520254845e-2f64,
  -0.81505692478987769484e-2f64,0.74297333588288568430e-2f64,
 0.14314212513540223314e-1f64, -0.36711242251332751607e-1f64,
 0.46240027755503814626e-1f64, -0.34921532671769023773e-1f64,
 0.46930051972353714773e-2f64,
 0.32842770336385381562e-1f64, -0.61317813706529588466e-1f64,
 0.67000809902468893103e-1f64,
  -0.45337449655535622885e-1f64,0.35794459576271920867e-2f64,
 0.41830061526027213385e-1f64,
  -0.72091371931944711708e-1f64,0.74150028530317793195e-1f64,
  -0.46487632538609942002e-1f64, 0f64,0.46487632538609942002e-1f64,
  -0.74150028530317793195e-1f64,0.72091371931944711708e-1f64,
  -0.41830061526027213385e-1f64, -0.35794459576271920867e-2f64,
 0.45337449655535622885e-1f64, -0.67000809902468893103e-1f64,
 0.61317813706529588466e-1f64, -0.32842770336385381562e-1f64,
  -0.46930051972353714773e-2f64,0.34921532671769023773e-1f64,
  -0.46240027755503814626e-1f64,0.36711242251332751607e-1f64,
  -0.14314212513540223314e-1f64, -0.74297333588288568430e-2f64,
 0.81505692478987769484e-2f64,0.90693182942442189743e-2f64,
  -0.11121000903959576737e-1f64, -0.71308296141317458546e-2f64,
 0.29219439765986671645e-1f64, -0.45820286629778129593e-1f64,
 0.49088381175879124421e-1f64, -0.35614888785023038938e-1f64,
 0.78906970900092777895e-2f64,
 0.26262843038404929480e-1f64, -0.56143674270125757857e-1f64,
 0.71700220472378350694e-1f64,
  -0.66963544500697307945e-1f64,0.42215091779892228883e-1f64,
  -0.41338867413966866997e-2f64, -0.36164891772995367321e-1f64,
 0.66584367783847858225e-1f64, -0.77874712365070098328e-1f64,
 0.66584367783847858225e-1f64, -0.36164891772995367321e-1f64,
  -0.41338867413966866997e-2f64,0.42215091779892228883e-1f64,
  -0.66963544500697307945e-1f64,
 0.71700220472378350694e-1f64, -0.56143674270125757857e-1f64,
 0.26262843038404929480e-1f64,
 0.78906970900092777895e-2f64, -0.35614888785023038938e-1f64,
 0.49088381175879124421e-1f64,
  -0.45820286629778129593e-1f64,0.29219439765986671645e-1f64,
  -0.71308296141317458546e-2f64, -0.11121000903959576737e-1f64,
 0.90693182942442189743e-2f64, -0.99848472706332791043e-2f64,
 0.14701271465939718856e-1f64, -0.32917820356048383366e-3f64,
  -0.19201195309873585230e-1f64,0.38409681836626963278e-1f64,
  -0.51647324405878909521e-1f64,0.54522171113149311354e-1f64,
  -0.45040302741689006270e-1f64,0.24183738595685990149e-1f64,
 0.42204134165479735097e-2f64, -0.34317295181348742251e-1f64,
 0.59542472465494579941e-1f64, -0.74135115907618101263e-1f64,
 0.74491937840566532596e-1f64, -0.60042604725161994304e-1f64,
 0.33437677409000083169e-1f64, 0f64,
  -0.33437677409000083169e-1f64,0.60042604725161994304e-1f64,
  -0.74491937840566532596e-1f64,0.74135115907618101263e-1f64,
  -0.59542472465494579941e-1f64,0.34317295181348742251e-1f64,
  -0.42204134165479735097e-2f64, -0.24183738595685990149e-1f64,
 0.45040302741689006270e-1f64, -0.54522171113149311354e-1f64,
 0.51647324405878909521e-1f64, -0.38409681836626963278e-1f64,
 0.19201195309873585230e-1f64,0.32917820356048383366e-3f64,
  -0.14701271465939718856e-1f64,0.99848472706332791043e-2f64,
 0.11775579274769383373e-1f64, -0.19892153937316935880e-1f64,
 0.95335114477449041055e-2f64,0.57661528440359081617e-2f64,
  -0.23382690532380910781e-1f64,0.40237257037170725321e-1f64,
  -0.53280289903551636474e-1f64,0.59974361806023689068e-1f64,
  -0.58701684061992853224e-1f64,0.49033407111597129616e-1f64,
  -0.31818835267847249219e-1f64,0.90800541261162098886e-2f64,
 0.16272906819312603838e-1f64, -0.40863896581186229487e-1f64,
 0.61346046297517367703e-1f64,
  -0.74896047554167268919e-1f64,0.79632642148310325817e-1f64,
  -0.74896047554167268919e-1f64,0.61346046297517367703e-1f64,
  -0.40863896581186229487e-1f64,0.16272906819312603838e-1f64,
 0.90800541261162098886e-2f64, -0.31818835267847249219e-1f64,
 0.49033407111597129616e-1f64, -0.58701684061992853224e-1f64,
 0.59974361806023689068e-1f64, -0.53280289903551636474e-1f64,
 0.40237257037170725321e-1f64, -0.23382690532380910781e-1f64,
 0.57661528440359081617e-2f64,0.95335114477449041055e-2f64,
  -0.19892153937316935880e-1f64,
 0.11775579274769383373e-1f64, -0.13562702617218467450e-1f64,
 0.24885419969649845849e-1f64, -0.18368693901908875583e-1f64,
 0.81673147806084084638e-2f64,0.47890591326129587131e-2f64,
  -0.19313752945227974024e-1f64,0.34065953398362954708e-1f64,
  -0.47667045133463415672e-1f64,0.58820377816690514309e-1f64,
  -0.66424139824618415970e-1f64,
 0.69667606260856092515e-1f64, -0.68102459384364543253e-1f64,
 0.61683024923302547971e-1f64,
  -0.50771943476441639136e-1f64,0.36110771847327189215e-1f64,
  -0.18758028464284563358e-1f64, 0f64,0.18758028464284563358e-1f64,
  -0.36110771847327189215e-1f64,0.50771943476441639136e-1f64,
  -0.61683024923302547971e-1f64,0.68102459384364543253e-1f64,
  -0.69667606260856092515e-1f64,0.66424139824618415970e-1f64,
  -0.58820377816690514309e-1f64,0.47667045133463415672e-1f64,
  -0.34065953398362954708e-1f64,0.19313752945227974024e-1f64,
  -0.47890591326129587131e-2f64, -0.81673147806084084638e-2f64,
 0.18368693901908875583e-1f64, -0.24885419969649845849e-1f64,
 0.13562702617218467450e-1f64,0.20576545037980523979e-1f64,
  -0.40093155172981004337e-1f64,0.36954083167944054826e-1f64,
  -0.31856506837591907746e-1f64,0.24996323181546255126e-1f64,
  -0.16637165210473614136e-1f64,0.71002706773325085237e-2f64,
 0.32478629093205201133e-2f64,
  -0.14009562579050569518e-1f64,0.24771262248780618922e-1f64,
  -0.35119395835433647559e-1f64,0.44656290368574753171e-1f64,
  -0.53015448339647394161e-1f64,0.59875631995693046782e-1f64,
  -0.64973208326045193862e-1f64,0.68112280331082143373e-1f64,
  -0.69172215234062186994e-1f64,0.68112280331082143373e-1f64,
  -0.64973208326045193862e-1f64,0.59875631995693046782e-1f64,
  -0.53015448339647394161e-1f64,0.44656290368574753171e-1f64,
  -0.35119395835433647559e-1f64,0.24771262248780618922e-1f64,
  -0.14009562579050569518e-1f64,0.32478629093205201133e-2f64,
 0.71002706773325085237e-2f64, -0.16637165210473614136e-1f64,
 0.24996323181546255126e-1f64, -0.31856506837591907746e-1f64,
 0.36954083167944054826e-1f64, -0.40093155172981004337e-1f64,
 0.20576545037980523979e-1f64, -0.27584914609096156163e-1f64,
 0.54904171411058497973e-1f64, -0.54109756419563083153e-1f64,
 0.52794234894345577483e-1f64, -0.50970276026831042415e-1f64,
 0.48655445537990983379e-1f64,
  -0.45872036510847994332e-1f64,0.42646854695899611372e-1f64,
  -0.39010960357087507670e-1f64,0.34999369144476467749e-1f64,
  -0.30650714874402762189e-1f64,0.26006877464703437057e-1f64,
  -0.21112579608213651273e-1f64,0.16014956068786763273e-1f64,
  -0.10763099747751940252e-1f64,0.54075888924374485533e-2f64, 0f64,
  -0.54075888924374485533e-2f64,0.10763099747751940252e-1f64,
  -0.16014956068786763273e-1f64,
 0.21112579608213651273e-1f64, -0.26006877464703437057e-1f64,
 0.30650714874402762189e-1f64,
  -0.34999369144476467749e-1f64,0.39010960357087507670e-1f64,
  -0.42646854695899611372e-1f64,0.45872036510847994332e-1f64,
  -0.48655445537990983379e-1f64,0.50970276026831042415e-1f64,
  -0.52794234894345577483e-1f64,0.54109756419563083153e-1f64,
  -0.54904171411058497973e-1f64,0.27584914609096156163e-1f64,
 0.13794141262469565740e-1f64, -0.27588282524939131481e-1f64,
 0.27588282524939131481e-1f64, -0.27588282524939131481e-1f64,
 0.27588282524939131481e-1f64, -0.27588282524939131481e-1f64,
 0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.27588282524939131481e-1f64,
  -0.27588282524939131481e-1f64,0.13794141262469565740e-1f64
];

static Tleft : [f64; 1089] = [
  1f64, -0.86602540378443864678f64, 0f64,0.33071891388307382381f64, 0f64,
  -0.20728904939721249057f64, 0f64,0.15128841196122722208f64, 0f64,
  -0.11918864298744029244f64, 0f64,0.98352013661686631224e-1f64, 0f64,
  -0.83727065404940845733e-1f64, 0f64,0.72893399403505841203e-1f64, 0f64,
  -0.64544632643375022436e-1f64, 0f64,0.57913170372415565639e-1f64, 0f64,
  -0.52518242575729562263e-1f64, 0f64,0.48043311993977520457e-1f64, 0f64,
  -0.44271433659733990243e-1f64, 0f64,0.41048928022856771981e-1f64, 0f64,
  -0.38263878662008271459e-1f64, 0f64,0.35832844026365304501e-1f64, 0f64, 0f64,
 0.50000000000000000000f64, -0.96824583655185422130f64,0.57282196186948000082f64,
 0.21650635094610966169f64, -0.35903516540862679125f64, -0.97578093724974971969e-1f64,
 0.26203921611325660506f64,0.55792409597991015609e-1f64, -0.20644078533943456204f64,
  -0.36172381205961199479e-1f64,0.17035068468874958194f64,
 0.25371838001497225980e-1f64, -0.14501953125000000000f64,
  -0.18786835250972344757e-1f64,0.12625507130328301066f64,
 0.14473795929590520582e-1f64, -0.11179458309419422675f64,
  -0.11494434254897626155e-1f64,0.10030855351241635862f64,
 0.93498556820544479096e-2f64, -0.90964264465390582629e-1f64,
  -0.77546391824364392762e-2f64,0.83213457337452292745e-1f64,
 0.65358085945588638605e-2f64, -0.76680372422574234569e-1f64,
  -0.55835321940047427169e-2f64,0.71098828931825789428e-1f64,
 0.48253327982967591019e-2f64, -0.66274981937248958553e-1f64,
  -0.42118078245337801387e-2f64,0.62064306433355646267e-1f64,
 0.37083386598903548973e-2f64, 0f64, 0f64,0.25000000000000000000f64,
  -0.73950997288745200531f64,0.83852549156242113615f64, -0.23175620272173946716f64,
  -0.37791833195149451496f64,0.25710129174850522325f64,0.21608307321780204633f64,
  -0.22844049245646009157f64, -0.14009503000335388415f64,0.19897685605518413847f64,
 0.98264706042471226893e-1f64, -0.17445445004279014046f64,
  -0.72761100054958328401e-1f64,0.15463589893742108388f64,
 0.56056770591708784481e-1f64, -0.13855313872640495158f64,
  -0.44517752443294564781e-1f64,0.12534277657695128850f64,
 0.36211835346039665762e-1f64, -0.11434398255136139683f64,
  -0.30033588409423828125e-1f64,0.10506705408753910481f64,
 0.25313077840725783008e-1f64, -0.97149327637744872155e-1f64,
  -0.21624927200393328444e-1f64,0.90319582367202122625e-1f64,
 0.18688433567711780666e-1f64, -0.84372291635345108584e-1f64,
  -0.16312261561845420752e-1f64,0.79149526894804751586e-1f64,
 0.14362333871852474757e-1f64, 0f64, 0f64, 0f64,0.12500000000000000000f64,
  -0.49607837082461073572f64,0.82265291131801144317f64, -0.59621200088559103072f64,
  -0.80054302859059362371e-1f64,0.42612156697795759420f64,
  -0.90098145270865592887e-1f64, -0.29769623255090078484f64,0.13630307904779758221f64,
 0.21638835185708931831f64, -0.14600247270306082052f64, -0.16348801804014290453f64,
 0.14340708728599057249f64,0.12755243353979286190f64, -0.13661523715071346961f64,
  -0.10215585947881057394f64,0.12864248070157166547f64,0.83592528025348693602e-1f64,
  -0.12066728689302565222f64, -0.69633728678718053052e-1f64,0.11314245177331919532f64,
 0.58882939251410088028e-1f64, -0.10621835858758221487f64,
  -0.50432266865187597572e-1f64,0.99916834723527771581e-1f64,
 0.43672094283057258509e-1f64, -0.94206380251950852413e-1f64,
  -0.38181356812697746418e-1f64,0.89035739656537771225e-1f64,
 0.33661934598216332678e-1f64, 0f64, 0f64, 0f64, 0f64,0.62500000000000000000e-1f64,
  -0.31093357409581873586f64,0.67604086414949799246f64, -0.75644205980613611039f64,
 0.28990586430124175741f64,0.30648508196770360914f64, -0.35801372616842500052f64,
  -0.91326869828709014708e-1f64,0.31127929687500000000f64,
  -0.90915752838698393094e-2f64, -0.25637381283965534330f64,
 0.57601077850322797594e-1f64,0.21019685709225757945f64,
  -0.81244992138514014256e-1f64, -0.17375078516720988858f64,
 0.92289437277967051125e-1f64,0.14527351914265391374f64,
  -0.96675340792832019889e-1f64, -0.12289485697108543415f64,
 0.97448175340011084006e-1f64,0.10511755943298339844f64,
  -0.96242247086378239657e-1f64, -0.90822942272780513537e-1f64,
 0.93966350452322132384e-1f64,0.79189411876493712558e-1f64,
  -0.91139307067989309325e-1f64, -0.69613039934383197265e-1f64,
 0.88062491671135767870e-1f64,0.61646331729340817494e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64,
 0.31250000000000000000e-1f64, -0.18684782411095934408f64,0.50176689760410660236f64,
  -0.74784031498626095398f64,0.56472001151566251186f64,0.14842464993721351203e-1f64,
  -0.41162920273003120936f64,0.20243071230196532282f64,0.23772054897172750436f64,
  -0.24963810923972235950f64, -0.12116179938394678936f64,0.24330535483519110663f64,
 0.47903849781124471359e-1f64, -0.22133299683101224293f64,
  -0.20542915138527200983e-2f64,0.19653465717678146728f64,
  -0.26818172626509178444e-1f64, -0.17319122357631210944f64,
 0.45065391411065545445e-1f64,0.15253391395444065941f64,
  -0.56543897711725408302e-1f64, -0.13469154928743585367f64,
 0.63632471400208840155e-1f64,0.11941684923913523817f64,
  -0.67828850207933293098e-1f64, -0.10636309084510652670f64,
 0.70095786922999181504e-1f64,0.95187373095150709082e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64,0.15625000000000000000e-1f64, -0.10909562534194485289f64,
 0.34842348626527747318f64, -0.64461114561628111443f64,0.69382480527334683659f64,
  -0.29551102358528827763f64, -0.25527584713978439819,0.38878771718544715394f64,
  -0.82956185835347407489e-2f64, -0.31183177761966943912f64,0.12831420840372374767f64,
 0.22067618205599434368f64, -0.17569196937129496961f64, -0.14598057000132284135f64,
 0.18864406621763419484f64,0.89921002550386645767e-1f64, -0.18571835020187122114f64,
  -0.48967672227195481777e-1f64,0.17584685670380332798f64,
 0.19267984545067426324e-1f64, -0.16335437520503462738f64,
 0.22598055455032407594e-2f64,0.15032800884170631129f64,
  -0.17883358353754640871e-1f64, -0.13774837869432209951f64,
 0.29227555960587143675e-1f64,0.12604194747513151053f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64,0.78125000000000000000e-2f64, -0.62377810244809812496e-1f64,
 0.23080781467370883845f64, -0.50841310636012325368f64,0.69834547012574056043f64,
  -0.52572723156526459672f64,0.11464215704954976471e-1f64,0.38698869011491210342f64,
  -0.26125646622255207507, -0.16951698812361607510f64,0.29773875898928782269f64,
 0.20130501202570367491e-1f64, -0.26332493149159310198f64,
 0.67734613690401207009e-1f64,0.21207315477103762715f64, -0.11541543390889415193f64,
  -0.16249634759782417533f64,0.13885887405041735068f64,0.11996491328010275427f64,
  -0.14810432001630926895f64, -0.85177658352556243411e-1f64,0.14918860659904380587f64,
 0.57317789510444151564e-1f64, -0.14569827645586660151f64,
  -0.35213090145965327390e-1f64,0.13975998126844578198f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64,0.39062500000000000000e-2f64, -0.35101954600803571207e-1f64,
 0.14761284084133737720f64, -0.37655033076080192966f64,0.62410290231517322776,
  -0.64335622317683389875f64,0.28188168266139524244f64,0.22488495672137010675f64,
  -0.39393811089283576186f64,0.75184777995770096714e-1f64,0.28472023119398293003,
  -0.20410910833705899572f64, -0.15590046962908511750f64,0.23814567544617953125f64,
 0.54442805556829031204e-1f64, -0.22855930338589720954f64,
 0.16303223615756629897e-1f64,0.20172722433875559213f64,
  -0.62723406421217419404e-1f64, -0.17012230831020922010f64,
 0.91754642766136561612e-1f64,0.13927644821381121197, -0.10886600968068418181f64,
  -0.11139075654373395292f64,0.11797455976331702879f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64,0.19531250000000000000e-2f64, -0.19506820659607596598e-1f64,
 0.91865676095362231937e-1f64, -0.26604607809696493849f64,0.51425874205091288223,
  -0.66047561132505329292f64,0.48660109511591303851f64, -0.17575661168678285615e-1f64,
  -0.36594333408055703366f64,0.29088854695378694533f64,0.11318677346656537927f64,
  -0.31110645235730182168f64,0.60733219161008787341e-1f64,0.24333848233620420826,
  -0.15254312332655419708f64, -0.15995968483455388613f64,0.19010344455215289289f64,
 0.86040636766440260000e-1f64, -0.19652589954665259945f64,
  -0.27633388517205837713e-1f64,0.18660848552712880387f64,
  -0.15942583868416775867e-1f64, -0.16902042462382064786f64,
 0.47278526495327740646e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
 0.97656250000000000000e-3f64, -0.10731084460857378207e-1f64,
 0.55939644713816406331e-1f64, -0.18118487371914493668f64,0.39914857299829864263f64,
  -0.60812322949933902435f64,0.60011887183061967583f64, -0.26002695805835928795f64,
  -0.20883922404786010096f64,0.38988130966114638081f64, -0.11797833550782589082f64,
  -0.25231824756239520077f64,0.24817859972953934712f64,0.90516417677868996417e-1f64,
  -0.26079073291293066798f64,0.30259468817169480161e-1f64,0.22178195264114178432f64,
  -0.10569877864302048175f64, -0.16679648389266977455f64,0.14637718550245050850f64,
 0.11219272032739559870f64, -0.16359363640525750353f64, -0.64358194509092101393e-1f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.48828125000000000000e-3f64,
  -0.58542865274813470967e-2f64,0.33461741635290096452e-1f64,
  -0.11979993155896201271f64,0.29580223766987206958f64, -0.51874761979436016742f64,
 0.62861483498014306968f64, -0.44868895761051453296f64,0.12567502628371529386e-1f64,
 0.35040366183235474275f64, -0.30466868455569500886f64, -0.70903913601490112666e-1f64,
 0.30822791893032512740f64, -0.11969443264190207736f64, -0.20764760317621313946f64,
 0.20629838355452128532f64,0.95269702915334718507e-1f64, -0.22432624768705133300f64,
  -0.33103381593477797101e-2f64,0.20570036048155716333f64,
  -0.62208282720094518964e-1f64, -0.17095309330441436348f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.24414062500000000000e-3f64,
  -0.31714797501871532475e-2f64,0.19721062526127334100e-1f64,
  -0.77311181185536498246e-1f64,0.21124871792841566575f64, -0.41777980401893650886f64,
 0.59401977834943551650f64, -0.56132417807488349048f64,0.23433675061367565951f64,
 0.20222775295220942126, -0.38280372496506190127f64,0.14443804214023095767f64,
 0.22268950939178466797, -0.27211314150777981984f64, -0.34184876506180717313e-1f64,
 0.26006498895669734842f64, -0.97650425186005090107e-1f64, -0.19024527660129101293f64,
 0.16789164198044635671f64,0.10875811641651905252f64, -0.19276785058805921298f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.12207031250000000000e-3f64,
  -0.17078941137247586143e-2f64,0.11477733754843910060e-1f64,
  -0.48887017020924625462e-1f64,0.14634927241421789683f64, -0.32156282683019547854f64,
 0.52165811920227223937, -0.60001958466396926460f64,0.41208501541480733755f64,
  -0.11366945503190350975e-2f64, -0.33968093962672089159f64,0.30955190935923386766f64,
 0.40657421856578262210e-1f64, -0.29873400409871531764,0.16094481791768257440f64,
 0.16876122436206497694f64, -0.23650217045022161255f64, -0.33070260090574765012e-1f64,
 0.22985258456375907796f64, -0.68645651043827097771e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.61035156250000000000e-4f64,
  -0.91501857608428649078e-3f64,0.66085179496951987952e-2f64,
  -0.30383171695850355404e-1f64,0.98840838845366876117e-1f64,
  -0.23855447246420318989f64,0.43322017468145613917f64, -0.58049033744876107191,
 0.52533893203742699346f64, -0.20681056202371946180, -0.20180000924562504384f64,
 0.37503922291962681797, -0.15988102869837429062f64, -0.19823558102762374094f64,
 0.28393023878803799622f64, -0.11188133439357510403e-1f64, -0.24730368377168229255f64,
 0.14731529061377942839f64,0.14878558042884266021f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.30517578125000000000e-4f64,
  -0.48804277318479845551e-3f64,0.37696080990601968396e-2f64,
  -0.18603912108994738255e-1f64,0.65325006755649582964e-1f64,
  -0.17162960707938819795f64,0.34411527956476971322f64, -0.52289350347082497959f64,
 0.57319653625674910592f64, -0.37662253421045430413f64, -0.14099055105384663902e-1f64,
 0.33265570610216904208f64, -0.30921265572647566661f64, -0.19911390594166455281e-1f64,
 0.28738590811031797718f64, -0.18912130469738472647f64, -0.13235936203215819193f64,
 0.25076406142356675279f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64,0.15258789062500000000e-4f64, -0.25928719280954633249e-3f64,
 0.21327398937568540428e-2f64, -0.11244626133630732010e-1f64,
 0.42375605740664331966e-1f64, -0.12031130345907846211,0.26352562258934426830f64,
  -0.44590628258512682078f64,0.56682835613700749379f64, -0.49116715128261660395f64,
 0.17845943097110339078f64,0.20541650677432497477f64, -0.36739803642257458221f64,
 0.16776034069210108273f64,0.17920950989905112908f64, -0.28867732805385066532f64,
 0.46473465543376206337e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64,0.76293945312500000000e-5f64, -0.13727610943181290891e-3f64,
 0.11979683091449349286e-2f64, -0.67195313034570709806e-2f64,
 0.27044920779931968175e-1f64, -0.82472196498517457862e-1f64,
 0.19570475044896150093f64, -0.36391620788543817693f64,0.52241392782736588032f64,
  -0.54727504974907879912f64,0.34211551468813581183f64,0.31580472732719957762e-1f64,
  -0.32830006549176759667f64,0.30563797665254420769f64,0.64905014620683140120e-2f64,
  -0.27642986248995073032f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.38146972656250000000e-5f64,
  -0.72454147007837596854e-4f64,0.66859847582761390285e-3f64,
  -0.39751311980366118437e-2f64,0.17015198650201528366e-1f64,
  -0.55443621868993855715e-1f64,0.14157060481641692131f64, -0.28641242619559616836f64,
 0.45610665490966615415f64, -0.55262786406029265394f64,0.45818352706035500108f64,
  -0.14984403004611673047f64, -0.21163807462970713245f64,0.36007252928843413718f64,
  -0.17030961385712954159f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.19073486328125000000e-5f64,
  -0.38135049864067468562e-4f64,0.37101393638555730015e-3f64,
  -0.23305339886279723213e-2f64,0.10569913448297127219e-1f64,
  -0.36640175162216897547e-1f64,0.10010476414320235508f64, -0.21860074212675559892f64,
 0.38124757096345313719, -0.52020999209879669177f64,0.52172632730659212045f64,
  -0.30841620620308814614f64, -0.50322546186721500184e-1f64,0.32577618885114899053f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64,0.95367431640625000000e-6f64, -0.20021483206955925244e-4f64,
 0.20481807322420625431e-3f64, -0.13553476938058909882e-2f64,
 0.64919676350791905019e-2f64, -0.23848725425069251903e-1f64,
 0.69384632678886421292e-1f64, -0.16249711393618776934,0.30736618106830314788f64,
  -0.46399909601971539157f64,0.53765031034002467225f64, -0.42598991476520183929f64,
 0.12130445348350215652f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.47683715820312500000e-6f64,
  -0.10487707828484902486e-4f64,0.11254146162337528943e-3f64,
  -0.78248929534271987118e-3f64,0.39468337145306794566e-2f64,
  -0.15313546659475671763e-1f64,0.47249070825218564146e-1f64,
  -0.11804374107101480543f64,0.24031796927792491122f64, -0.39629215049166341285f64,
 0.51629108968402548545f64, -0.49622372075429782915f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
 0.23841857910156250000e-6f64, -0.54823314130625337326e-5f64,
 0.61575377321535518154e-4f64, -0.44877834366497538134e-3f64,
 0.23774612048621955857e-2f64, -0.97136347645161687796e-2f64,
 0.31671599547606636717e-1f64, -0.84028665767000747480e-1f64,
 0.18298487576742964949f64, -0.32647878537696945218f64,0.46970971486488895077f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64,0.11920928955078125000e-6f64, -0.28604020001177375838e-5f64,
 0.33559227978295551013e-4f64, -0.25583821662860610560e-3f64,
 0.14201552747787302339e-2f64, -0.60938046986874414969e-2f64,
 0.20930869247951926793e-1f64, -0.58745021125678072911e-1f64,
 0.13613725780285953720f64, -0.26083988356030237586f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
 0.59604644775390625000e-7f64, -0.14898180663526043291e-5f64,
 0.18224991282807693921e-4f64, -0.14504433444608833821e-3f64,
 0.84184722720281809548e-3f64, -0.37846965430000478789e-2f64,
 0.13656355548211376864e-1f64, -0.40409541997718853934e-1f64,
 0.99226988101858325902e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
 0.29802322387695312500e-7f64, -0.77471708843445529468e-6f64,
 0.98649879372606876995e-5f64, -0.81814934772838523887e-4f64,
 0.49554483992403011328e-3f64, -0.23290922072351413938e-2f64,
 0.88068134250844034186e-2f64, -0.27393666952485719070e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64,0.14901161193847656250e-7f64, -0.40226235946098233685e-6f64,
 0.53236418690561306700e-5f64, -0.45933829691164002269e-4f64,
 0.28982005232838857913e-3f64, -0.14212974043211018374e-2f64,
 0.56192363087488842264e-2f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
 0.74505805969238281250e-8f64, -0.20858299254133430408e-6f64,
 0.28648457300134381744e-5f64, -0.25677535898258910850e-4f64,
 0.16849420429491355445e-3f64, -0.86062824010315834002e-3f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64,0.37252902984619140625e-8f64, -0.10801736017613096861e-6f64,
 0.15376606719887104015e-5f64, -0.14296523739727437959e-4f64,
 0.97419023656050887203e-4f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
 0.18626451492309570312e-8f64, -0.55871592916438890146e-7f64,
 0.82331193828137454068e-6f64, -0.79302250528382787666e-5f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.93132257461547851562e-9f64,
  -0.28867244235852488244e-7f64,0.43982811713864556957e-6f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.46566128730773925781e-9f64,
  -0.14899342093408253335e-7f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
  0f64, 0f64,0.23283064365386962891e-9f64
];

static Tright : [f64; 1089] = [
  1f64,0.86602540378443864678f64, 0f64, -0.33071891388307382381f64, 0f64,
   0.20728904939721249057f64, 0f64, -0.15128841196122722208f64, 0f64,
   0.11918864298744029244f64, 0f64, -0.98352013661686631224e-1f64, 0f64,
   0.83727065404940845733e-1f64, 0f64, -0.72893399403505841203e-1f64, 0f64,
   0.64544632643375022436e-1f64, 0f64, -0.57913170372415565639e-1f64, 0f64,
   0.52518242575729562263e-1f64, 0f64, -0.48043311993977520457e-1f64, 0f64,
   0.44271433659733990243e-1f64, 0f64, -0.41048928022856771981e-1f64, 0f64,
   0.38263878662008271459e-1f64, 0f64, -0.35832844026365304501e-1f64, 0f64, 0f64,
   0.50000000000000000000f64,0.96824583655185422130f64,0.57282196186948000082f64,
    -0.21650635094610966169f64, -0.35903516540862679125f64,0.97578093724974971969e-1f64,
   0.26203921611325660506f64, -0.55792409597991015609e-1f64, -0.20644078533943456204f64,
   0.36172381205961199479e-1f64,0.17035068468874958194f64,
    -0.25371838001497225980e-1f64, -0.14501953125000000000f64,
   0.18786835250972344757e-1f64,0.12625507130328301066f64,
    -0.14473795929590520582e-1f64, -0.11179458309419422675f64,
   0.11494434254897626155e-1f64,0.10030855351241635862f64,
    -0.93498556820544479096e-2f64, -0.90964264465390582629e-1f64,
   0.77546391824364392762e-2f64,0.83213457337452292745e-1f64,
    -0.65358085945588638605e-2f64, -0.76680372422574234569e-1f64,
   0.55835321940047427169e-2f64,0.71098828931825789428e-1f64,
    -0.48253327982967591019e-2f64, -0.66274981937248958553e-1f64,
   0.42118078245337801387e-2f64,0.62064306433355646267e-1f64,
    -0.37083386598903548973e-2f64, 0f64, 0f64,0.25000000000000000000f64,
   0.73950997288745200531f64,0.83852549156242113615f64,0.23175620272173946716f64,
    -0.37791833195149451496f64, -0.25710129174850522325f64,0.21608307321780204633f64,
   0.22844049245646009157f64, -0.14009503000335388415f64, -0.19897685605518413847f64,
   0.98264706042471226893e-1f64,0.17445445004279014046f64,
    -0.72761100054958328401e-1f64, -0.15463589893742108388f64,
   0.56056770591708784481e-1f64,0.13855313872640495158f64,
    -0.44517752443294564781e-1f64, -0.12534277657695128850f64,
   0.36211835346039665762e-1f64,0.11434398255136139683f64,
    -0.30033588409423828125e-1f64, -0.10506705408753910481f64,
   0.25313077840725783008e-1f64,0.97149327637744872155e-1f64,
    -0.21624927200393328444e-1f64, -0.90319582367202122625e-1f64,
   0.18688433567711780666e-1f64,0.84372291635345108584e-1f64,
    -0.16312261561845420752e-1f64, -0.79149526894804751586e-1f64,
   0.14362333871852474757e-1f64, 0f64, 0f64, 0f64,0.12500000000000000000f64,
   0.49607837082461073572f64,0.82265291131801144317f64,0.59621200088559103072f64,
    -0.80054302859059362371e-1f64, -0.42612156697795759420f64,
    -0.90098145270865592887e-1f64,0.29769623255090078484f64,0.13630307904779758221f64,
    -0.21638835185708931831f64, -0.14600247270306082052f64,0.16348801804014290453f64,
   0.14340708728599057249f64, -0.12755243353979286190f64, -0.13661523715071346961f64,
   0.10215585947881057394f64,0.12864248070157166547f64, -0.83592528025348693602e-1f64,
    -0.12066728689302565222f64,0.69633728678718053052e-1f64,0.11314245177331919532f64,
    -0.58882939251410088028e-1f64, -0.10621835858758221487f64,
   0.50432266865187597572e-1f64,0.99916834723527771581e-1f64,
    -0.43672094283057258509e-1f64, -0.94206380251950852413e-1f64,
   0.38181356812697746418e-1f64,0.89035739656537771225e-1f64,
    -0.33661934598216332678e-1f64, 0f64, 0f64, 0f64, 0f64,0.62500000000000000000e-1f64,
   0.31093357409581873586f64,0.67604086414949799246f64,0.75644205980613611039f64,
   0.28990586430124175741f64, -0.30648508196770360914f64, -0.35801372616842500052f64,
   0.91326869828709014708e-1f64,0.31127929687500000000f64,0.90915752838698393094e-2f64,
    -0.25637381283965534330f64, -0.57601077850322797594e-1f64,0.21019685709225757945f64,
   0.81244992138514014256e-1f64, -0.17375078516720988858f64,
    -0.92289437277967051125e-1f64,0.14527351914265391374f64,
   0.96675340792832019889e-1f64, -0.12289485697108543415f64,
    -0.97448175340011084006e-1f64,0.10511755943298339844f64,
   0.96242247086378239657e-1f64, -0.90822942272780513537e-1f64,
    -0.93966350452322132384e-1f64,0.79189411876493712558e-1f64,
   0.91139307067989309325e-1f64, -0.69613039934383197265e-1f64,
    -0.88062491671135767870e-1f64,0.61646331729340817494e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64,
   0.31250000000000000000e-1f64,0.18684782411095934408f64,0.50176689760410660236f64,
   0.74784031498626095398f64,0.56472001151566251186f64, -0.14842464993721351203e-1f64,
    -0.41162920273003120936f64, -0.20243071230196532282f64,0.23772054897172750436f64,
   0.24963810923972235950f64, -0.12116179938394678936f64, -0.24330535483519110663f64,
   0.47903849781124471359e-1f64,0.22133299683101224293f64,
    -0.20542915138527200983e-2f64, -0.19653465717678146728f64,
    -0.26818172626509178444e-1f64,0.17319122357631210944f64,
   0.45065391411065545445e-1f64, -0.15253391395444065941f64,
    -0.56543897711725408302e-1f64,0.13469154928743585367f64,
   0.63632471400208840155e-1f64, -0.11941684923913523817f64,
    -0.67828850207933293098e-1f64,0.10636309084510652670f64,
   0.70095786922999181504e-1f64, -0.95187373095150709082e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64,0.15625000000000000000e-1f64,0.10909562534194485289f64,
   0.34842348626527747318f64,0.64461114561628111443f64,0.69382480527334683659f64,
   0.29551102358528827763f64, -0.25527584713978439819, -0.38878771718544715394f64,
    -0.82956185835347407489e-2f64,0.31183177761966943912f64,0.12831420840372374767f64,
    -0.22067618205599434368f64, -0.17569196937129496961f64,0.14598057000132284135f64,
   0.18864406621763419484f64, -0.89921002550386645767e-1f64, -0.18571835020187122114f64,
   0.48967672227195481777e-1f64,0.17584685670380332798f64,
    -0.19267984545067426324e-1f64, -0.16335437520503462738f64,
    -0.22598055455032407594e-2f64,0.15032800884170631129f64,
   0.17883358353754640871e-1f64, -0.13774837869432209951f64,
    -0.29227555960587143675e-1f64,0.12604194747513151053f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64,0.78125000000000000000e-2f64,0.62377810244809812496e-1f64,
   0.23080781467370883845f64,0.50841310636012325368f64,0.69834547012574056043f64,
   0.52572723156526459672f64,0.11464215704954976471e-1f64, -0.38698869011491210342f64,
    -0.26125646622255207507,0.16951698812361607510f64,0.29773875898928782269f64,
    -0.20130501202570367491e-1f64, -0.26332493149159310198f64,
    -0.67734613690401207009e-1f64,0.21207315477103762715f64,0.11541543390889415193f64,
    -0.16249634759782417533f64, -0.13885887405041735068f64,0.11996491328010275427f64,
   0.14810432001630926895f64, -0.85177658352556243411e-1f64, -0.14918860659904380587f64,
   0.57317789510444151564e-1f64,0.14569827645586660151f64,
    -0.35213090145965327390e-1f64, -0.13975998126844578198f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64,0.39062500000000000000e-2f64,0.35101954600803571207e-1f64,
   0.14761284084133737720f64,0.37655033076080192966f64,0.62410290231517322776,
   0.64335622317683389875f64,0.28188168266139524244f64, -0.22488495672137010675f64,
    -0.39393811089283576186f64, -0.75184777995770096714e-1f64,0.28472023119398293003,
   0.20410910833705899572f64, -0.15590046962908511750f64, -0.23814567544617953125f64,
   0.54442805556829031204e-1f64,0.22855930338589720954f64,0.16303223615756629897e-1f64,
    -0.20172722433875559213f64, -0.62723406421217419404e-1f64,0.17012230831020922010f64,
   0.91754642766136561612e-1f64, -0.13927644821381121197, -0.10886600968068418181f64,
   0.11139075654373395292f64,0.11797455976331702879f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64,0.19531250000000000000e-2f64,0.19506820659607596598e-1f64,
   0.91865676095362231937e-1f64,0.26604607809696493849f64,0.51425874205091288223,
   0.66047561132505329292f64,0.48660109511591303851f64,0.17575661168678285615e-1f64,
    -0.36594333408055703366f64, -0.29088854695378694533f64,0.11318677346656537927f64,
   0.31110645235730182168f64,0.60733219161008787341e-1f64, -0.24333848233620420826,
    -0.15254312332655419708f64,0.15995968483455388613f64,0.19010344455215289289f64,
    -0.86040636766440260000e-1f64, -0.19652589954665259945f64,
   0.27633388517205837713e-1f64,0.18660848552712880387f64,0.15942583868416775867e-1f64,
    -0.16902042462382064786f64, -0.47278526495327740646e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64,0.97656250000000000000e-3f64,0.10731084460857378207e-1f64,
   0.55939644713816406331e-1f64,0.18118487371914493668f64,0.39914857299829864263f64,
   0.60812322949933902435f64,0.60011887183061967583f64,0.26002695805835928795f64,
    -0.20883922404786010096f64, -0.38988130966114638081f64, -0.11797833550782589082f64,
   0.25231824756239520077f64,0.24817859972953934712f64, -0.90516417677868996417e-1f64,
    -0.26079073291293066798f64, -0.30259468817169480161e-1f64,0.22178195264114178432f64,
   0.10569877864302048175f64, -0.16679648389266977455f64, -0.14637718550245050850f64,
   0.11219272032739559870f64,0.16359363640525750353f64, -0.64358194509092101393e-1f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.48828125000000000000e-3f64,
   0.58542865274813470967e-2f64,0.33461741635290096452e-1f64,0.11979993155896201271f64,
   0.29580223766987206958f64,0.51874761979436016742f64,0.62861483498014306968f64,
   0.44868895761051453296f64,0.12567502628371529386e-1f64, -0.35040366183235474275f64,
    -0.30466868455569500886f64,0.70903913601490112666e-1f64,0.30822791893032512740f64,
   0.11969443264190207736f64, -0.20764760317621313946f64, -0.20629838355452128532f64,
   0.95269702915334718507e-1f64,0.22432624768705133300f64,
    -0.33103381593477797101e-2f64, -0.20570036048155716333f64,
    -0.62208282720094518964e-1f64,0.17095309330441436348f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.24414062500000000000e-3f64,
   0.31714797501871532475e-2f64,0.19721062526127334100e-1f64,
   0.77311181185536498246e-1f64,0.21124871792841566575f64,0.41777980401893650886f64,
   0.59401977834943551650f64,0.56132417807488349048f64,0.23433675061367565951f64,
    -0.20222775295220942126, -0.38280372496506190127f64, -0.14443804214023095767f64,
   0.22268950939178466797,0.27211314150777981984f64, -0.34184876506180717313e-1f64,
    -0.26006498895669734842f64, -0.97650425186005090107e-1f64,0.19024527660129101293f64,
   0.16789164198044635671f64, -0.10875811641651905252f64, -0.19276785058805921298f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.12207031250000000000e-3f64,
   0.17078941137247586143e-2f64,0.11477733754843910060e-1f64,
   0.48887017020924625462e-1f64,0.14634927241421789683f64,0.32156282683019547854f64,
   0.52165811920227223937,0.60001958466396926460f64,0.41208501541480733755f64,
   0.11366945503190350975e-2f64, -0.33968093962672089159f64, -0.30955190935923386766f64,
   0.40657421856578262210e-1f64,0.29873400409871531764,0.16094481791768257440f64,
    -0.16876122436206497694f64, -0.23650217045022161255f64,0.33070260090574765012e-1f64,
   0.22985258456375907796f64,0.68645651043827097771e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.61035156250000000000e-4f64,
   0.91501857608428649078e-3f64,0.66085179496951987952e-2f64,
   0.30383171695850355404e-1f64,0.98840838845366876117e-1f64,0.23855447246420318989f64,
   0.43322017468145613917f64,0.58049033744876107191,0.52533893203742699346f64,
   0.20681056202371946180, -0.20180000924562504384f64, -0.37503922291962681797,
    -0.15988102869837429062f64,0.19823558102762374094f64,0.28393023878803799622f64,
   0.11188133439357510403e-1f64, -0.24730368377168229255f64, -0.14731529061377942839f64,
   0.14878558042884266021f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64,0.30517578125000000000e-4f64,0.48804277318479845551e-3f64,
   0.37696080990601968396e-2f64,0.18603912108994738255e-1f64,
   0.65325006755649582964e-1f64,0.17162960707938819795f64,0.34411527956476971322f64,
   0.52289350347082497959f64,0.57319653625674910592f64,0.37662253421045430413f64,
    -0.14099055105384663902e-1f64, -0.33265570610216904208f64, -0.30921265572647566661f64,
   0.19911390594166455281e-1f64,0.28738590811031797718f64,0.18912130469738472647f64,
    -0.13235936203215819193f64, -0.25076406142356675279f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.15258789062500000000e-4f64,
   0.25928719280954633249e-3f64,0.21327398937568540428e-2f64,
   0.11244626133630732010e-1f64,0.42375605740664331966e-1f64,0.12031130345907846211,
   0.26352562258934426830f64,0.44590628258512682078f64,0.56682835613700749379f64,
   0.49116715128261660395f64,0.17845943097110339078f64, -0.20541650677432497477f64,
    -0.36739803642257458221f64, -0.16776034069210108273f64,0.17920950989905112908f64,
   0.28867732805385066532f64,0.46473465543376206337e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.76293945312500000000e-5f64,
   0.13727610943181290891e-3f64,0.11979683091449349286e-2f64,
   0.67195313034570709806e-2f64,0.27044920779931968175e-1f64,
   0.82472196498517457862e-1f64,0.19570475044896150093f64,0.36391620788543817693f64,
   0.52241392782736588032f64,0.54727504974907879912f64,0.34211551468813581183f64,
    -0.31580472732719957762e-1f64, -0.32830006549176759667f64, -0.30563797665254420769f64,
   0.64905014620683140120e-2f64,0.27642986248995073032f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.38146972656250000000e-5f64,
   0.72454147007837596854e-4f64,0.66859847582761390285e-3f64,
   0.39751311980366118437e-2f64,0.17015198650201528366e-1f64,
   0.55443621868993855715e-1f64,0.14157060481641692131f64,0.28641242619559616836f64,
   0.45610665490966615415f64,0.55262786406029265394f64,0.45818352706035500108f64,
   0.14984403004611673047f64, -0.21163807462970713245f64, -0.36007252928843413718f64,
    -0.17030961385712954159f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.19073486328125000000e-5f64,
   0.38135049864067468562e-4f64,0.37101393638555730015e-3f64,
   0.23305339886279723213e-2f64,0.10569913448297127219e-1f64,
   0.36640175162216897547e-1f64,0.10010476414320235508f64,0.21860074212675559892f64,
   0.38124757096345313719,0.52020999209879669177f64,0.52172632730659212045f64,
   0.30841620620308814614f64, -0.50322546186721500184e-1f64, -0.32577618885114899053f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64,0.95367431640625000000e-6f64,0.20021483206955925244e-4f64,
   0.20481807322420625431e-3f64,0.13553476938058909882e-2f64,
   0.64919676350791905019e-2f64,0.23848725425069251903e-1f64,
   0.69384632678886421292e-1f64,0.16249711393618776934,0.30736618106830314788f64,
   0.46399909601971539157f64,0.53765031034002467225f64,0.42598991476520183929f64,
   0.12130445348350215652f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.47683715820312500000e-6f64,
   0.10487707828484902486e-4f64,0.11254146162337528943e-3f64,
   0.78248929534271987118e-3f64,0.39468337145306794566e-2f64,
   0.15313546659475671763e-1f64,0.47249070825218564146e-1f64,0.11804374107101480543f64,
   0.24031796927792491122f64,0.39629215049166341285f64,0.51629108968402548545f64,
   0.49622372075429782915f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.23841857910156250000e-6f64,
   0.54823314130625337326e-5f64,0.61575377321535518154e-4f64,
   0.44877834366497538134e-3f64,0.23774612048621955857e-2f64,
   0.97136347645161687796e-2f64,0.31671599547606636717e-1f64,
   0.84028665767000747480e-1f64,0.18298487576742964949f64,0.32647878537696945218f64,
   0.46970971486488895077f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.11920928955078125000e-6f64,
   0.28604020001177375838e-5f64,0.33559227978295551013e-4f64,
   0.25583821662860610560e-3f64,0.14201552747787302339e-2f64,
   0.60938046986874414969e-2f64,0.20930869247951926793e-1f64,
   0.58745021125678072911e-1f64,0.13613725780285953720f64,0.26083988356030237586f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.59604644775390625000e-7f64,
   0.14898180663526043291e-5f64,0.18224991282807693921e-4f64,
   0.14504433444608833821e-3f64,0.84184722720281809548e-3f64,
   0.37846965430000478789e-2f64,0.13656355548211376864e-1f64,
   0.40409541997718853934e-1f64,0.99226988101858325902e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64,0.29802322387695312500e-7f64,0.77471708843445529468e-6f64,
   0.98649879372606876995e-5f64,0.81814934772838523887e-4f64,
   0.49554483992403011328e-3f64,0.23290922072351413938e-2f64,
   0.88068134250844034186e-2f64,0.27393666952485719070e-1f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64,0.14901161193847656250e-7f64,0.40226235946098233685e-6f64,
   0.53236418690561306700e-5f64,0.45933829691164002269e-4f64,
   0.28982005232838857913e-3f64,0.14212974043211018374e-2f64,
   0.56192363087488842264e-2f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
   0.74505805969238281250e-8f64,0.20858299254133430408e-6f64,
   0.28648457300134381744e-5f64,0.25677535898258910850e-4f64,
   0.16849420429491355445e-3f64,0.86062824010315834002e-3f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64,0.37252902984619140625e-8f64,0.10801736017613096861e-6f64,
   0.15376606719887104015e-5f64,0.14296523739727437959e-4f64,
   0.97419023656050887203e-4f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
   0.18626451492309570312e-8f64,0.55871592916438890146e-7f64,
   0.82331193828137454068e-6f64,0.79302250528382787666e-5f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.93132257461547851562e-9f64,
   0.28867244235852488244e-7f64,0.43982811713864556957e-6f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,0.46566128730773925781e-9f64,
   0.14899342093408253335e-7f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64, 0f64,
    0f64, 0f64,0.23283064365386962891e-9f64
];

/// The QAG algorithm is a simple adaptive integration procedure. The integration region is divided into subintervals, and on each iteration
/// the subinterval with the largest estimated error is bisected. This reduces the overall error rapidly, as the subintervals become concentrated
/// around local difficulties in the integrand. These subintervals are managed by a gsl_integration_workspace struct, which handles the memory
/// for the subinterval ranges, results and error estimates.
pub struct IntegrationWorkspace {
    w: *mut ffi::gsl_integration_workspace
}

struct InternParam<'r, T:'r> {
    func: ::function<T>,
    param: &'r mut T,
    p2: f64
}

impl IntegrationWorkspace {
    /// This function allocates a workspace sufficient to hold n double precision intervals, their integration results and error estimates. One
    /// workspace may be used multiple times as all necessary reinitialization is performed automatically by the integration routines.
    pub fn new(n: u64) -> Option<IntegrationWorkspace> {
        let tmp = unsafe { ffi::gsl_integration_workspace_alloc(n) };

        if tmp.is_null() {
            None
        } else {
            Some(IntegrationWorkspace {
                w: tmp
            })
        }
    }

    /// This function applies an integration rule adaptively until an estimate of the integral of f over (a,b) is achieved within the desired
    /// absolute and relative error limits, epsabs and epsrel. The function returns the final approximation, result, and an estimate of the
    /// absolute error, abserr. The integration rule is determined by the value of key, which should be chosen from the following symbolic names,
    /// 
    /// GSL_INTEG_GAUSS15  (key = 1)
    /// 
    /// GSL_INTEG_GAUSS21  (key = 2)
    /// 
    /// GSL_INTEG_GAUSS31  (key = 3)
    /// 
    /// GSL_INTEG_GAUSS41  (key = 4)
    /// 
    /// GSL_INTEG_GAUSS51  (key = 5)
    /// 
    /// GSL_INTEG_GAUSS61  (key = 6)
    /// 
    /// corresponding to the 15f64, 21f64, 31f64, 41f64, 51 and 61 point Gauss-Kronrod rules. The higher-order rules give better accuracy for smooth functions,
    /// while lower-order rules save time when the function contains local difficulties, such as discontinuities.
    /// 
    /// On each iteration the adaptive integration strategy bisects the interval with the largest error estimate. The subintervals and their
    /// results are stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the
    /// allocated size of the workspace.
    pub fn qag<T>(&self, f: ::function<T>, arg: &mut T, a: f64, b: f64, epsabs: f64, epsrel: f64, limit: u64, key: enums::gauss_konrod_rule::GaussKonrodRule,
        result: &mut f64, abserr: &mut f64) -> ::Value {
        match key {
            ::GaussKonrodRule::Gauss15 => {
                intern_qag(f, arg, a, b, epsabs, epsrel, limit, self, result, abserr,
                    ::integration::qk15)
            }
            ::GaussKonrodRule::Gauss21 => {
                intern_qag(f, arg, a, b, epsabs, epsrel, limit, self, result, abserr,
                    ::integration::qk21)
            }
            ::GaussKonrodRule::Gauss31 => {
                intern_qag(f, arg, a, b, epsabs, epsrel, limit, self, result, abserr,
                    ::integration::qk31)
            }
            ::GaussKonrodRule::Gauss41 => {
                intern_qag(f, arg, a, b, epsabs, epsrel, limit, self, result, abserr,
                    ::integration::qk41)
            }
            ::GaussKonrodRule::Gauss51 => {
                intern_qag(f, arg, a, b, epsabs, epsrel, limit, self, result, abserr,
                    ::integration::qk51)
            }
            ::GaussKonrodRule::Gauss61 => {
                intern_qag(f, arg, a, b, epsabs, epsrel, limit, self, result, abserr,
                    ::integration::qk61)
            }
            /*_ => {
                let file = file!();
                "value of key does specify a known integration rule".with_c_str(|c_str|{
                    file.with_c_str(|c_file|{
                        unsafe { ffi::gsl_error(c_str, c_file, line!() as i32, ::Value::Inval as i32) }
                    });
                });
                // this line is not used but just for compilation...
                ::integration::qk15
            }*/
        }
    }

    /// This function applies the Gauss-Kronrod 21-point integration rule adaptively until an estimate of the integral of f over (a,b) is achieved
    /// within the desired absolute and relative error limits, epsabs and epsrel. The results are extrapolated using the epsilon-algorithm, which
    /// accelerates the convergence of the integral in the presence of discontinuities and integrable singularities. The function returns the
    /// final approximation from the extrapolation, result, and an estimate of the absolute error, abserr. The subintervals and their results are
    /// stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated size
    /// of the workspace.
    pub fn qags<T>(&self, f: ::function<T>, arg: &mut T, a: f64, b: f64, epsabs: f64, epsrel: f64, limit: u64, result: &mut f64,
        abserr: &mut f64) -> ::Value {
        unsafe { intern_qags(f, arg, a, b, epsabs, epsrel, limit, self, result, abserr, ::integration::qk21) }
    }

    /// This function applies the adaptive integration algorithm QAGS taking account of the user-supplied locations of singular points. The array
    /// pts of length npts should contain the endpoints of the integration ranges defined by the integration region and locations of the singularities.
    /// For example, to integrate over the region (a,b) with break-points at x_1, x_2, x_3 (where a < x_1 < x_2 < x_3 < b) the following pts array
    /// should be used
    /// 
    /// pts[0] = a
    /// pts[1] = x_1
    /// pts[2] = x_2
    /// pts[3] = x_3
    /// pts[4] = b
    /// with npts = 5.
    /// 
    /// If you know the locations of the singular points in the integration region then this routine will be faster than QAGS.
    pub fn qagp<T>(&self, f: ::function<T>, arg: &mut T, pts: &mut [f64], epsabs: f64, epsrel: f64, limit: u64, result: &mut f64,
        abserr: &mut f64) -> ::Value {
        unsafe { intern_qagp(f, arg, pts, epsabs, epsrel, limit, self, result, abserr, ::integration::qk21) }
    }

    /// This function computes the integral of the function f over the infinite interval (-\infty,+\infty). The integral is mapped onto the
    /// semi-open interval (0,1] using the transformation x = (1-t)/t,
    /// 
    /// \int_{-\infty}^{+\infty} dx f(x) = 
    ///      \int_0^1 dt (f((1-t)/t) + f((-1+t)/t))/t^2.
    /// 
    /// It is then integrated using the QAGS algorithm. The normal 21-point Gauss-Kronrod rule of QAGS is replaced by a 15-point rule, because
    /// the transformation can generate an integrable singularity at the origin. In this case a lower-order rule is more efficient.
    pub fn qagi<T>(&self, f: ::function<T>, arg: &mut T, epsabs: f64, epsrel: f64, limit: u64, result: &mut f64, abserr: &mut f64) -> ::Value {
        let mut s = InternParam{func: f, param: arg, p2: 0f64};

        unsafe { intern_qags(i_transform, &mut s, 0f64, 1f64, epsabs, epsrel, limit, self, result, abserr, ::integration::qk15) }
    }

    /// This function computes the integral of the function f over the semi-infinite interval (a,+\infty). The integral is mapped onto the
    /// semi-open interval (0,1] using the transformation x = a + (1-t)/t,
    /// 
    /// \int_{a}^{+\infty} dx f(x) = 
    ///      \int_0^1 dt f(a + (1-t)/t)/t^2
    /// 
    /// and then integrated using the QAGS algorithm.
    pub fn qagiu<T>(&self, f: ::function<T>, arg: &mut T, a: f64, epsabs: f64, epsrel: f64, limit: u64, result: &mut f64, abserr: &mut f64) -> ::Value {
        let mut s = InternParam{func: f, param: arg, p2: a};

        unsafe { intern_qags(iu_transform, &mut s, 0f64, 1f64, epsabs, epsrel, limit, self, result, abserr, ::integration::qk15) }
    }

    /// This function computes the integral of the function f over the semi-infinite interval (-\infty,b). The integral is mapped onto the semi-open
    /// interval (0,1] using the transformation x = b - (1-t)/t,
    /// 
    /// \int_{-\infty}^{b} dx f(x) = 
    ///      \int_0^1 dt f(b - (1-t)/t)/t^2
    /// 
    /// and then integrated using the QAGS algorithm.
    pub fn qagil<T>(&self, f: ::function<T>, arg: &mut T, b: f64, epsabs: f64, epsrel: f64, limit: u64, result: &mut f64, abserr: &mut f64) -> ::Value {
        let mut s = InternParam{func: f, param: arg, p2: b};

        unsafe { intern_qags(il_transform, &mut s, 0f64, 1f64, epsabs, epsrel, limit, self, result, abserr, ::integration::qk15) }
    }

    /// This function computes the Cauchy principal value of the integral of f over (a,b), with a singularity at c,
    /// 
    /// I = \int_a^b dx f(x) / (x - c)
    /// 
    /// The adaptive bisection algorithm of QAG is used, with modifications to ensure that subdivisions do not occur at the singular point x = c.
    /// When a subinterval contains the point x = c or is close to it then a special 25-point modified Clenshaw-Curtis rule is used to control
    /// the singularity. Further away from the singularity the algorithm uses an ordinary 15-point Gauss-Kronrod integration rule.
    #[allow(unused_assignments)]
    pub fn qawc<T>(&self, f: ::function<T>, arg: &mut T, a: f64, b: f64, c: f64, epsabs: f64, epsrel: f64, limit: u64,
        result: &mut f64, abserr: &mut f64) -> ::Value {
        let mut result0 = 0f64;
        let mut abserr0 = 0f64;
        let mut roundoff_type1 = 0i32;
        let mut roundoff_type2 = 0i32;
        let mut error_type = 0i32;
        let mut err_reliable = 0i32;
        let mut sign = 1f64;
        let mut lower = 0f64;
        let mut higher = 0f64;

        /* Initialize results */
        *result = 0f64;
        *abserr = 0f64;

        unsafe {
            if limit > (*self.w).limit {
                rgsl_error!("iteration limit exceeds available workspace", ::Value::Inval);
            }

            if b < a {
                lower = b;
                higher = a;
                sign = -1f64;
            } else {
                lower = a;
                higher = b;
            }

            self.initialise(lower, higher);

            if epsabs <= 0f64 && (epsrel < 50f64 * ::DBL_EPSILON || epsrel < 0.5e-28f64) {
                rgsl_error!("tolerance cannot be acheived with given epsabs and epsrel", ::Value::BadTol);
            }

            if c == a || c == b {
                rgsl_error!("cannot integrate with singularity on endpoint", ::Value::Inval);
            }      

            /* perform the first integration */
            qc25c(f, arg, lower, higher, c, &mut result0, &mut abserr0, &mut err_reliable);

            self.set_initial_result(result0, abserr0);

            /* Test on accuracy, use 0.01 relative error as an extra safety margin on the first iteration (ignored for subsequent iterations) */
            let mut tolerance = epsabs.max(epsrel * fabsf64(result0));

            if abserr0 < tolerance && abserr0 < 0.01f64 * fabsf64(result0) {
                *result = sign * result0;
                *abserr = abserr0;

                return ::Value::Success;
            } else if limit == 1 {
                *result = sign * result0;
                *abserr = abserr0;

                rgsl_error!("a maximum of one iteration was insufficient", ::Value::MaxIter);
            }

            let mut area = result0;
            let mut errsum = abserr0;

            let mut iteration = 1;

            loop {
                let mut a_i = 0f64;
                let mut b_i = 0f64;
                let mut r_i = 0f64;
                let mut e_i = 0f64;
                let mut area1 = 0f64;
                let mut area2 = 0f64;
                let mut error1 = 0f64;
                let mut error2 = 0f64;
                let mut err_reliable1 = 0i32;
                let mut err_reliable2 = 0i32;

                /* Bisect the subinterval with the largest error estimate */
                self.retrieve(&mut a_i, &mut b_i, &mut r_i, &mut e_i);

                let a1 = a_i; 
                let mut b1 = 0.5 * (a_i + b_i);
                let mut a2 = b1;
                let b2 = b_i;

                if c > a1 && c <= b1 {
                    b1 = 0.5f64 * (c + b2);
                    a2 = b1;
                } else if c > b1 && c < b2 {
                    b1 = 0.5f64 * (a1 + c);
                    a2 = b1;
                }

                qc25c(f, arg, a1, b1, c, &mut area1, &mut error1, &mut err_reliable1);
                qc25c(f, arg, a2, b2, c, &mut area2, &mut error2, &mut err_reliable2);

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

                errsum += error12 - e_i;
                area += area12 - r_i;

                if err_reliable1 != 0 && err_reliable2 != 0 {
                    let delta = r_i - area12;

                    if fabsf64(delta) <= 1.0e-5f64 * fabsf64(area12) && error12 >= 0.99f64 * e_i {
                        roundoff_type1 += 1;
                    }
                    if iteration >= 10 && error12 > e_i {
                        roundoff_type2 += 1;
                    }
                }

                tolerance = epsabs.max(epsrel * fabsf64(area));

                if errsum > tolerance {
                    if roundoff_type1 >= 6 || roundoff_type2 >= 20 {
                        /* round off error */
                        error_type = 2;
                    }

                    /* set error flag in the case of bad integrand behaviour at a point of the integration range */
                    if ::util::subinterval_too_small(a1, a2, b2) {
                        error_type = 3;
                    }
                }

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

                self.retrieve(&mut a_i, &mut b_i, &mut r_i, &mut e_i);

                iteration += 1;

                if iteration < limit && error_type == 0 && errsum > tolerance {
                } else {
                    break;
                }
            }

            *result = sign * self.sum_results();
            *abserr = errsum;

            if errsum <= tolerance {
              ::Value::Success
            } else if error_type == 2 {
                rgsl_error!("roundoff error prevents tolerance from being achieved", ::Value::Round);
                ::Value::Round
            } else if error_type == 3 {
                rgsl_error!("bad integrand behavior found in the integration interval", ::Value::Sing);
                ::Value::Sing
            } else if iteration == limit {
                rgsl_error!("maximum number of subdivisions reached", ::Value::MaxIter);
                ::Value::MaxIter
            } else {
                rgsl_error!("could not integrate function", ::Value::Failed);
                ::Value::Failed
            }
        }
    }

    pub fn sort_results(&self) {
        unsafe {
            let nint = (*self.w).size as usize;
            let mut t_elist = CSlice::new((*self.w).elist, nint);
            let mut t_order = CSlice::new((*self.w).order, nint);
            let elist = t_elist.as_mut();
            let order = t_order.as_mut();

            for i in 0usize..nint {
                let i1 = order[i] as usize;
                let mut e1 = elist[i1];
                let mut i_max = i1;

                for j in (i + 1)..nint {
                    let i2 = order[j] as usize;
                    let e2 = elist[i2];

                    if e2 >= e1 {
                        i_max = i2;
                        e1 = e2;
                    }
                }

                if i_max != i1 {
                    order[i] = order[i_max];
                    order[i_max] = i1 as u64;
                }
            }
            (*self.w).i = order[0];
        }
    }

    pub fn qpsrt(&self) {
        let w = self.w;

        unsafe {
            let mut order = CSlice::new((*w).order, (*w).size as usize + 1usize);
            let last = (*w).size - 1;
            let limit = (*w).limit;

            let mut i_nrmax = (*w).nrmax;
            let mut i_maxerr = order.as_ref()[i_nrmax as usize];

            // Check whether the list contains more than two error estimates
            if last < 2 {
                order.as_mut()[0usize] = 0;
                order.as_mut()[1usize] = 1;
                (*w).i = i_maxerr;
                return ;
            }

            // Compute the number of elements in the list to be maintained in descending order. This number depends on the number of
            // subdivisions still allowed.
            let top =  if last < (limit / 2 + 2) {
                last
            } else {
                limit - last + 1
            };

            let elist = CSlice::new((*w).elist, top as usize + 1usize);
            let errmax = elist.as_ref()[i_maxerr as usize];

            // 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 && errmax > elist.as_ref()[order.as_ref()[i_nrmax as usize - 1usize] as usize] {
                order.as_mut()[i_nrmax as usize] = order.as_ref()[i_nrmax as usize - 1usize];
                i_nrmax -= 1;
            }

            // Insert errmax by traversing the list top-down, starting comparison from the element elist(order(i_nrmax+1)).
            let mut i = i_nrmax + 1;

            // The order of the tests in the following line is important to prevent a segmentation fault
            while i < top && errmax < elist.as_ref()[order.as_ref()[i as usize] as usize] {
                order.as_mut()[i as usize - 1usize] = order.as_ref()[i as usize];
                i += 1;
            }

            order.as_mut()[i as usize - 1usize] = i_maxerr;

            // Insert errmin by traversing the list bottom-up
            let errmin = elist.as_ref()[last as usize];

            let mut k = top - 1;
            while k > i - 2 && errmin >= elist.as_ref()[order.as_ref()[k as usize] as usize] {
                order.as_mut()[k as usize + 1usize] = order.as_ref()[k as usize];
                k -= 1;
            }
      
            order.as_mut()[k as usize + 1usize] = last;

            // Set i_max and e_max
            i_maxerr = order.as_ref()[i_nrmax as usize] ;

            (*w).i = i_maxerr ;
            (*w).nrmax = i_nrmax ;
        }
    }

    pub fn sum_results(&self) -> f64 {
        unsafe {
            let f_w = self.w;
            let mut result_sum = 0f64;
            let v_rlist = CSlice::new((*f_w).rlist, (*f_w).size as usize);
            let rlist = v_rlist.as_ref();

            for k in 0usize..((*f_w).size as usize) {
                result_sum += rlist[k];
            }

            result_sum
        }
    }

    pub fn retrieve(&self, a: &mut f64, b: &mut f64, r: &mut f64, e: &mut f64) {
        unsafe {
            let w = self.w;
            let alist = CSlice::new((*w).alist, (*w).i as usize + 1usize);
            let blist = CSlice::new((*w).blist, (*w).i as usize + 1usize);
            let rlist = CSlice::new((*w).rlist, (*w).i as usize + 1usize);
            let elist = CSlice::new((*w).elist, (*w).i as usize + 1usize);

            let i = (*w).i as usize;

            *a = alist.as_ref()[i];
            *b = blist.as_ref()[i];
            *r = rlist.as_ref()[i];
            *e = elist.as_ref()[i];
        }
    }

    pub fn update(&self, a1: f64, b1: f64, area1: f64, error1: f64, a2: f64, b2: f64, area2: f64, error2: f64) {
        let w = self.w;

        unsafe {
            let i_max = (*w).i as usize;
            let i_new = (*w).size as usize;
            let tmp = if i_max > i_new {
                i_max + 1usize
            } else {
                i_new + 1usize
            };
            let mut alist = CSlice::new((*w).alist, tmp);
            let mut blist = CSlice::new((*w).blist, tmp);
            let mut rlist = CSlice::new((*w).rlist, tmp);
            let mut elist = CSlice::new((*w).elist, tmp);
            let mut level = CSlice::new((*w).level, tmp);

            let new_level = level.as_ref()[i_max] + 1;

            /* append the newly-created intervals to the list */

            if error2 > error1 {
                // blist[maxerr] is already == b2
                alist.as_mut()[i_max] = a2;
                rlist.as_mut()[i_max] = area2;
                elist.as_mut()[i_max] = error2;
                level.as_mut()[i_max] = new_level;

                alist.as_mut()[i_new] = a1;
                blist.as_mut()[i_new] = b1;
                rlist.as_mut()[i_new] = area1;
                elist.as_mut()[i_new] = error1;
                level.as_mut()[i_new] = new_level;
            } else {
                // alist[maxerr] is already == a1
                blist.as_mut()[i_max] = b1;
                rlist.as_mut()[i_max] = area1;
                elist.as_mut()[i_max] = error1;
                level.as_mut()[i_max] = new_level;

                alist.as_mut()[i_new] = a2;
                blist.as_mut()[i_new] = b2;
                rlist.as_mut()[i_new] = area2;
                elist.as_mut()[i_new] = error2;
                level.as_mut()[i_new] = new_level;
            }

            (*w).size += 1;

            if new_level > (*w).maximum_level {
                (*w).maximum_level = new_level;
            }

            self.qpsrt();
        }
    }

    pub fn set_initial_result(&self, result: f64, error: f64) {
        unsafe {
            let mut rlist = CSlice::new((*self.w).rlist, 1);
            let mut elist = CSlice::new((*self.w).elist, 1);

            (*self.w).size = 1;
            rlist.as_mut()[0] = result;
            elist.as_mut()[0] = error;
        }
    }

    pub fn initialise(&self, a: f64, b: f64) {
        let w = self.w;

        unsafe {
            let mut alist = CSlice::new((*w).alist, 1);
            let mut blist = CSlice::new((*w).blist, 1);
            let mut rlist = CSlice::new((*w).rlist, 1);
            let mut elist = CSlice::new((*w).elist, 1);
            let mut order = CSlice::new((*w).order, 1);
            let mut level = CSlice::new((*w).level, 1);

            (*w).size = 0;
            (*w).nrmax = 0;
            (*w).i = 0;
            alist.as_mut()[0] = a;
            blist.as_mut()[0] = b;
            rlist.as_mut()[0] = 0f64;
            elist.as_mut()[0] = 0f64;
            order.as_mut()[0] = 0u64;
            level.as_mut()[0] = 0u64;

            (*w).maximum_level = 0;
        }
    }

    fn append_interval(&self, a1: f64, b1: f64, area1: f64, error1: f64) {
        unsafe {
            let w = self.w;
            let i_new = (*w).size as usize;
            let mut alist = CSlice::new((*w).alist, i_new + 1usize);
            let mut blist = CSlice::new((*w).blist, i_new + 1usize);
            let mut rlist = CSlice::new((*w).rlist, i_new + 1usize);
            let mut elist = CSlice::new((*w).elist, i_new + 1usize);
            let mut order = CSlice::new((*w).order, i_new + 1usize);
            let mut level = CSlice::new((*w).level, i_new + 1usize);

            alist.as_mut()[i_new] = a1;
            blist.as_mut()[i_new] = b1;
            rlist.as_mut()[i_new] = area1;
            elist.as_mut()[i_new] = error1;
            order.as_mut()[i_new] = i_new as u64;
            level.as_mut()[i_new] = 0;

            (*w).size += 1;
        }
    }

    pub fn limit(&self) -> u64 {
        unsafe { (*self.w).limit }
    }

    pub fn size(&self) -> u64 {
        unsafe { (*self.w).size }
    }
}

impl Drop for IntegrationWorkspace {
    fn drop(&mut self) {
        unsafe { ffi::gsl_integration_workspace_free(self.w) };
        self.w = ::std::ptr::null_mut();
    }
}

impl ffi::FFI<ffi::gsl_integration_workspace> for IntegrationWorkspace {
    fn wrap(w: *mut ffi::gsl_integration_workspace) -> IntegrationWorkspace {
        IntegrationWorkspace {
            w: w
        }
    }

    fn unwrap(w: &IntegrationWorkspace) -> *mut ffi::gsl_integration_workspace {
        w.w
    }
}

/// The QAWS algorithm is designed for integrands with algebraic-logarithmic singularities at the end-points of an integration region. In order
/// to work efficiently the algorithm requires a precomputed table of Chebyshev moments.
pub struct IntegrationQawsTable {
    w: *mut ffi::gsl_integration_qaws_table
}

impl IntegrationQawsTable {
    /// This function allocates space for a gsl_integration_qaws_table struct describing a singular weight function W(x) with the parameters
    /// (\alpha, \beta, \mu, \nu),
    /// 
    /// W(x) = (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x)
    /// 
    /// where \alpha > -1f64, \beta > -1f64, and \mu = 0, 1, \nu = 0, 1. The weight function can take four different forms depending on the values
    /// of \mu and \nu,
    /// 
    /// W(x) = (x-a)^alpha (b-x)^beta                   (mu = 0, nu = 0)
    /// W(x) = (x-a)^alpha (b-x)^beta log(x-a)          (mu = 1, nu = 0)
    /// W(x) = (x-a)^alpha (b-x)^beta log(b-x)          (mu = 0, nu = 1)
    /// W(x) = (x-a)^alpha (b-x)^beta log(x-a) log(b-x) (mu = 1, nu = 1)
    /// 
    /// The singular points (a,b) do not have to be specified until the integral is computed, where they are the endpoints of the integration
    /// range.
    /// 
    /// The function returns a pointer to the newly allocated table gsl_integration_qaws_table if no errors were detected, and 0 in the case
    /// of error.
    pub fn new(alpha: f64, beta: f64, mu: i32, nu: i32) -> Option<IntegrationQawsTable> {
        let tmp = unsafe { ffi::gsl_integration_qaws_table_alloc(alpha, beta, mu, nu) };

        if tmp.is_null() {
            None
        } else {
            Some(IntegrationQawsTable {
                w: tmp
            })
        }
    }

    /// This function modifies the parameters (\alpha, \beta, \mu, \nu)
    pub fn set(&self, alpha: f64, beta: f64, mu: i32, nu: i32) -> ::Value {
        unsafe { ffi::gsl_integration_qaws_table_set(self.w, alpha, beta, mu, nu) }
    }

    /// This function computes the integral of the function f(x) over the interval (a,b) with the singular weight function (x-a)^\alpha
    /// (b-x)^\beta \log^\mu (x-a) \log^\nu (b-x). The parameters of the weight function (\alpha, \beta, \mu, \nu) are taken from the
    /// table self. The integral is,
    /// 
    /// I = \int_a^b dx f(x) (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x).
    /// 
    /// The adaptive bisection algorithm of QAG is used. When a subinterval contains one of the endpoints then a special 25-point modified
    /// Clenshaw-Curtis rule is used to control the singularities. For subintervals which do not include the endpoints an ordinary 15-point
    /// Gauss-Kronrod integration rule is used.
    #[allow(unused_assignments)]
    pub fn qaws<T>(&self, f: ::function<T>, arg: &mut T, a: f64, b: f64, epsabs: f64, epsrel: f64, limit: u64, workspace: &IntegrationWorkspace,
        result: &mut f64, abserr: &mut f64) -> ::Value {
        let mut result0 = 0f64;
        let mut abserr0 = 0f64;
        let mut roundoff_type1 = 0i32;
        let mut roundoff_type2 = 0i32;
        let mut error_type = 0i32;

        /* Initialize results */
        workspace.initialise(a, b);

        *result = 0f64;
        *abserr = 0f64;

        unsafe {
            if limit > (*workspace.w).limit {
                rgsl_error!("iteration limit exceeds available workspace", ::Value::Inval);
            }

            if b <= a {
                rgsl_error!("limits must form an ascending sequence, a < b", ::Value::Inval);
            }

            if epsabs <= 0f64 && (epsrel < 50f64 * ::DBL_EPSILON || epsrel < 0.5e-28f64) {
                rgsl_error!("tolerance cannot be acheived with given epsabs and epsrel", ::Value::BadTol);
            }

            /* perform the first integration */
            {
                let mut area1 = 0f64;
                let mut area2 = 0f64;
                let mut error1 = 0f64;
                let mut error2 = 0f64;
                let mut err_reliable1 = false;
                let mut err_reliable2 = false;
                let a1 = a;
                let b1 = 0.5f64 * (a + b);
                let a2 = b1;
                let b2 = b;

                qc25s(f, arg, a, b, a1, b1, self.w, &mut area1, &mut error1, &mut err_reliable1);
                qc25s(f, arg, a, b, a2, b2, self.w, &mut area2, &mut error2, &mut err_reliable2);

                if error1 > error2 {
                    workspace.append_interval(a1, b1, area1, error1);
                    workspace.append_interval(a2, b2, area2, error2);
                } else {
                    workspace.append_interval(a2, b2, area2, error2);
                    workspace.append_interval(a1, b1, area1, error1);
                }

                result0 = area1 + area2;
                abserr0 = error1 + error2;
            }

            /* Test on accuracy */
            let mut tolerance = epsabs.max(epsrel * fabsf64(result0));

            // Test on accuracy, use 0.01 relative error as an extra safety margin on the first iteration (ignored for subsequent iterations)
            if abserr0 < tolerance && abserr0 < 0.01 * fabsf64(result0) {
                *result = result0;
                *abserr = abserr0;

                return ::Value::Success;
            } else if limit == 1 {
                *result = result0;
                *abserr = abserr0;

                rgsl_error!("a maximum of one iteration was insufficient", ::Value::MaxIter);
            }

            let mut area = result0;
            let mut errsum = abserr0;

            let mut iteration = 2;

            loop {
                let mut a_i = 0f64;
                let mut b_i = 0f64;
                let mut r_i = 0f64;
                let mut e_i = 0f64;
                let mut area1 = 0f64;
                let mut area2 = 0f64;
                let mut error1 = 0f64;
                let mut error2 = 0f64;
                let mut err_reliable1 = false;
                let mut err_reliable2 = false;

                /* Bisect the subinterval with the largest error estimate */
                workspace.retrieve(&mut a_i, &mut b_i, &mut r_i, &mut e_i);

                let a1 = a_i; 
                let b1 = 0.5f64 * (a_i + b_i);
                let a2 = b1;
                let b2 = b_i;

                qc25s(f, arg, a, b, a1, b1, self.w, &mut area1, &mut error1, &mut err_reliable1);
                qc25s(f, arg, a, b, a2, b2, self.w, &mut area2, &mut error2, &mut err_reliable2);

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

                errsum += error12 - e_i;
                area += area12 - r_i;

                if err_reliable1 && err_reliable2 {
                    let delta = r_i - area12;

                    if fabsf64(delta) <= 1.0e-5f64 * fabsf64(area12) && error12 >= 0.99 * e_i {
                        roundoff_type1 += 1;
                    }
                    if iteration >= 10 && error12 > e_i {
                        roundoff_type2 += 1;
                    }
                }

                tolerance = epsabs.max(epsrel * fabsf64(area));

                if errsum > tolerance {
                    if roundoff_type1 >= 6 || roundoff_type2 >= 20 {
                        /* round off error */
                        error_type = 2;
                    }

                    /* set error flag in the case of bad integrand behaviour at a point of the integration range */
                    if ::util::subinterval_too_small(a1, a2, b2) {
                        error_type = 3;
                    }
                }

                workspace.update(a1, b1, area1, error1, a2, b2, area2, error2);
                workspace.retrieve(&mut a_i, &mut b_i, &mut r_i, &mut e_i);

                iteration += 1;

                if iteration < limit && error_type == 0 && errsum > tolerance {
                } else {
                    break;
                }
            }

            *result = workspace.sum_results();
            *abserr = errsum;

            if errsum <= tolerance {
                ::Value::Success
            } else if error_type == 2 {
                rgsl_error!("roundoff error prevents tolerance from being achieved", ::Value::Round);
                ::Value::Round
            } else if error_type == 3 {
                rgsl_error!("bad integrand behavior found in the integration interval", ::Value::Sing);
                ::Value::Sing
            } else if iteration == limit {
                rgsl_error!("maximum number of subdivisions reached", ::Value::MaxIter);
                ::Value::MaxIter
            } else {
                rgsl_error!("could not integrate function", ::Value::Failed);
                ::Value::Failed
            }
        }
    }
}

impl Drop for IntegrationQawsTable {
    fn drop(&mut self) {
        unsafe { ffi::gsl_integration_qaws_table_free(self.w) };
        self.w = ::std::ptr::null_mut();
    }
}

impl ffi::FFI<ffi::gsl_integration_qaws_table> for IntegrationQawsTable {
    fn wrap(w: *mut ffi::gsl_integration_qaws_table) -> IntegrationQawsTable {
        IntegrationQawsTable {
            w: w
        }
    }

    fn unwrap(w: &IntegrationQawsTable) -> *mut ffi::gsl_integration_qaws_table {
        w.w
    }
}

/// The QAWO algorithm is designed for integrands with an oscillatory factor, \sin(\omega x) or \cos(\omega x). In order to work efficiently
/// the algorithm requires a table of Chebyshev moments which must be pre-computed with calls to the functions below.
pub struct IntegrationQawoTable {
    w: *mut ffi::gsl_integration_qawo_table
}

impl IntegrationQawoTable {
    /// This function allocates space for a gsl_integration_qawo_table struct and its associated workspace describing a sine or cosine weight
    /// function W(x) with the parameters (\omega, L),
    /// 
    /// W(x) = sin(omega x)
    /// W(x) = cos(omega x)
    /// 
    /// The parameter L must be the length of the interval over which the function will be integrated L = b - a. The choice of sine or cosine
    /// is made with the parameter sine which should be chosen from one of the two following symbolic values:
    /// 
    /// ::Cosine
    /// ::IntegrationQawo::Sine
    /// 
    /// The gsl_integration_qawo_table is a table of the trigonometric coefficients required in the integration process. The parameter n determines
    /// the number of levels of coefficients that are computed. Each level corresponds to one bisection of the interval L, so that n levels are
    /// sufficient for subintervals down to the length L/2^n. The integration routine gsl_integration_qawo returns the error ::Table if the
    /// number of levels is insufficient for the requested accuracy.
    pub fn new(omega: f64, l: f64, sine: ::IntegrationQawo, n: u64) -> Option<IntegrationQawoTable> {
        let tmp = unsafe { ffi::gsl_integration_qawo_table_alloc(omega, l, sine, n) };

        if tmp.is_null() {
            None
        } else {
            Some(IntegrationQawoTable {
                w: tmp
            })
        }
    }

    /// This function changes the parameters omega, L and sine of the existing self workspace.
    pub fn set(&self, omega: f64, l: f64, sine: ::IntegrationQawo) -> ::Value {
        unsafe { ffi::gsl_integration_qawo_table_set(self.w, omega, l, sine) }
    }

    /// This function allows the length parameter l of the self workspace to be changed.
    pub fn set_length(&self, l: f64) -> ::Value {
        unsafe { ffi::gsl_integration_qawo_table_set_length(self.w, l) }
    }

    /// This function uses an adaptive algorithm to compute the integral of f over (a,b) with the weight function \sin(\omega x) or \cos(\omega x)
    /// defined by the table wf,
    /// 
    /// I = \int_a^b dx f(x) sin(omega x)
    /// I = \int_a^b dx f(x) cos(omega x)
    /// 
    /// The results are extrapolated using the epsilon-algorithm to accelerate the convergence of the integral. The function returns the final
    /// approximation from the extrapolation, result, and an estimate of the absolute error, abserr. The subintervals and their results are
    /// stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated
    /// size of the workspace.
    /// 
    /// Those subintervals with “large” widths d where d\omega > 4 are computed using a 25-point Clenshaw-Curtis integration rule, which handles
    /// the oscillatory behavior. Subintervals with a “small” widths where d\omega < 4 are computed using a 15-point Gauss-Kronrod integration.
    pub fn qawo<T>(&self, f: ::function<T>, arg: &mut T, a: f64, epsabs: f64, epsrel: f64, limit: u64, workspace: &IntegrationWorkspace,
        result: &mut f64, abserr: &mut f64) -> ::Value {
        let mut result0 = 0f64;
        let mut abserr0 = 0f64;
        let mut resabs0 = 0f64;
        let mut resasc0 = 0f64;

        let mut ertest = 0f64;
        let mut error_over_large_intervals = 0f64;
        let mut reseps = 0f64;
        let mut abseps = 0f64;
        let mut correc = 0f64;
        let mut ktmin = 0u64;
        let mut roundoff_type1 = 0i32;
        let mut roundoff_type2 = 0i32;
        let mut roundoff_type3 = 0i32;
        let mut error_type = 0i32;
        let mut error_type2 = 0i32;

        let mut extrapolate = false;
        let mut extall = false;
        let mut disallow_extrapolation = false;

        unsafe {
            let mut table : ffi::extrapolation_table = ::std::mem::zeroed();

            let b = a + (*self.w).L ;
            let abs_omega = fabsf64((*self.w).omega) ;

            /* Initialize results */
            workspace.initialise(a, b);

            *result = 0f64;
            *abserr = 0f64;

            if limit > (*workspace.w).limit {
                rgsl_error!("iteration limit exceeds available workspace", ::Value::Inval);
            }

            /* Test on accuracy */
            if epsabs <= 0f64 && (epsrel < 50f64 * ::DBL_EPSILON || epsrel < 0.5e-28f64) {
                rgsl_error!("tolerance cannot be acheived with given epsabs and epsrel", ::Value::BadTol);
            }

            /* Perform the first integration */
            qc25f(f, arg, a, b, self.w, 0, &mut result0, &mut abserr0, &mut resabs0, &mut resasc0);

            workspace.set_initial_result(result0, abserr0);

            let mut tolerance = epsabs.max(epsrel * fabsf64(result0));

            if abserr0 <= 100f64 * ::DBL_EPSILON * resabs0 && abserr0 > tolerance {
                *result = result0;
                *abserr = abserr0;

                rgsl_error!("cannot reach tolerance because of roundoff error on first attempt", ::Value::Round);
            } else if (abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0f64 {
                *result = result0;
                *abserr = abserr0;

                return ::Value::Success;
            } else if limit == 1 {
                *result = result0;
                *abserr = abserr0;

                rgsl_error!("a maximum of one iteration was insufficient", ::Value::MaxIter);
            }

            /* Initialization */
            initialise_table(&mut table);

            if 0.5f64 * abs_omega * fabsf64(b - a) <= 2f64 {
                append_table(&mut table, result0);
                extall = true;
            }

            let mut area = result0;
            let mut errsum = abserr0;

            let mut res_ext = result0;
            let mut err_ext = ::DBL_MAX;

            let positive_integrand = test_positivity(result0, resabs0);

            let mut iteration = 1;

            loop {
                let mut a_i = 0f64;
                let mut b_i = 0f64;
                let mut r_i = 0f64;
                let mut e_i = 0f64;
                let mut area1 = 0f64;
                let mut area2 = 0f64;
                let mut error1 = 0f64;
                let mut error2 = 0f64;
                let mut label70 = false;
                let mut resasc1 = 0f64;
                let mut resasc2 = 0f64;
                let mut resabs1 = 0f64;
                let mut resabs2 = 0f64;

                /* Bisect the subinterval with the largest error estimate */
                workspace.retrieve(&mut a_i, &mut b_i, &mut r_i, &mut e_i);

                let level = CSlice::new((*workspace.w).level, (*workspace.w).i as usize + 1);

                let current_level = level.as_ref()[(*workspace.w).i as usize] + 1;

                if current_level >= (*self.w).n {
                    /* exceeded limit of table */
                    error_type = -1;
                    break;
                }

                let a1 = a_i;
                let b1 = 0.5f64 * (a_i + b_i);
                let a2 = b1;
                let b2 = b_i;

                iteration += 1;

                qc25f(f, arg, a1, b1, self.w, current_level, &mut area1, &mut error1, &mut resabs1, &mut resasc1);
                qc25f(f, arg, a2, b2, self.w, current_level, &mut area2, &mut error2, &mut resabs2, &mut resasc2);

                let area12 = area1 + area2;
                let error12 = error1 + error2;
                let last_e_i = e_i;

                /* Improve previous approximations to the integral and test for accuracy.

                We write these expressions in the same way as the original QUADPACK code so that the rounding errors are the same, which
                makes testing easier. */
                errsum = errsum + error12 - e_i;
                area = area + area12 - r_i;

                tolerance = epsabs.max(epsrel * fabsf64(area));

                if resasc1 != error1 && resasc2 != error2 {
                    let delta = r_i - area12;

                    if fabsf64(delta) <= 1.0e-5f64 * fabsf64(area12) && error12 >= 0.99f64 * e_i {
                        if !extrapolate {
                            roundoff_type1 += 1;
                        } else {
                            roundoff_type2 += 1;
                        }
                    }
                    if iteration > 10 && error12 > e_i {
                        roundoff_type3 += 1;
                    }
                }

                /* Test for roundoff and eventually set error flag */
                if roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20 {
                    /* round off error */
                    error_type = 2;
                }

                if roundoff_type2 >= 5 {
                    error_type2 = 1;
                }

                /* set error flag in the case of bad integrand behaviour at a point of the integration range */
                if ::util::subinterval_too_small(a1, a2, b2) {
                    error_type = 4;
                }

                /* append the newly-created intervals to the list */
                workspace.update(a1, b1, area1, error1, a2, b2, area2, error2);

                if errsum <= tolerance {
                    *result = workspace.sum_results();
                    *abserr = errsum;
                    return return_error(error_type);
                }

                if error_type != 0 {
                    break;
                }

                if iteration >= limit - 1 {
                    error_type = 1;
                    break;
                }

                /* set up variables on first iteration */
                if iteration == 2 && extall {
                    error_over_large_intervals = errsum;
                    ertest = tolerance;
                    append_table(&mut table, area);
                    continue;
                }

                if disallow_extrapolation {
                    continue;
                }

                if extall {
                    error_over_large_intervals += -last_e_i;

                    if current_level < (*workspace.w).maximum_level {
                        error_over_large_intervals += error12;
                    }

                    if extrapolate {
                        label70 = true;
                    }
                }

                if label70 == false {
                    if large_interval(workspace.w) {
                        continue;
                    }

                    if extall {
                        extrapolate = true;
                        (*workspace.w).nrmax = 1;
                    } else {
                        /* test whether the interval to be bisected next is the smallest interval. */
                        let i = (*workspace.w).i as usize;
                        let blist = CSlice::new((*workspace.w).blist, i + 1usize);
                        let alist = CSlice::new((*workspace.w).alist, i + 1usize);
                        let width = blist.as_ref()[i] - alist.as_ref()[i];

                        if 0.25f64 * fabsf64(width) * abs_omega > 2f64 {
                            continue;
                        }

                        extall = true;
                        error_over_large_intervals = errsum;
                        ertest = tolerance;
                        continue;
                    }
                }

                if error_type2 == 0 && error_over_large_intervals > ertest {
                    if increase_nrmax(workspace.w) {
                        continue;
                    }
                }

                /* Perform extrapolation */
                append_table(&mut table, area);

                if table.n < 3 {
                    reset_nrmax(workspace.w);
                    extrapolate = false;
                    error_over_large_intervals = errsum;
                    continue;
                }

                intern_qelg(&mut table, &mut reseps, &mut abseps);

                ktmin += 1;

                if ktmin > 5 && err_ext < 0.001f64 * errsum {
                    error_type = 5;
                }

                if abseps < err_ext {
                    ktmin = 0;
                    err_ext = abseps;
                    res_ext = reseps;
                    correc = error_over_large_intervals;
                    ertest = epsabs.max(epsrel * fabsf64(reseps));
                    if err_ext <= ertest {
                        break;
                    }
                }

                /* Prepare bisection of the smallest interval. */
                if table.n == 1 {
                    disallow_extrapolation = true;
                }

                if error_type == 5 {
                    break;
                }

                /* work on interval with largest error */
                reset_nrmax(workspace.w);
                extrapolate = false;
                error_over_large_intervals = errsum;

                if iteration >= limit {
                    break;
                }
            }

            *result = res_ext;
            *abserr = err_ext;

            if err_ext == ::DBL_MAX {
                *result = workspace.sum_results();
                *abserr = errsum;
                return return_error(error_type);
            }

            if error_type != 0 || error_type2 != 0 {
                if error_type2 != 0 {
                    err_ext += correc;
                }

                if error_type == 0 {
                    error_type = 3;
                }

                if *result != 0f64 && area != 0f64 {
                    if err_ext / fabsf64(res_ext) > errsum / fabsf64(area) {
                        *result = workspace.sum_results();
                        *abserr = errsum;
                        return return_error(error_type);
                    }
                } else if err_ext > errsum {
                    *result = workspace.sum_results();
                    *abserr = errsum;
                    return return_error(error_type);
                } else if area == 0f64 {
                    return return_error(error_type);
                }
            }

            /*  Test on divergence. */
            {
                let max_area = fabsf64(res_ext).max(fabsf64(area));

                if !positive_integrand && max_area < 0.01f64 * resabs0 {
                    return return_error(error_type);
                }
            }

            {
                let ratio = res_ext / area;

                if ratio < 0.01f64 || ratio > 100f64 || errsum > fabsf64(area) {
                    error_type = 6;
                }
            }
            return_error(error_type)
        }
    }
}

impl Drop for IntegrationQawoTable {
    fn drop(&mut self) {
        unsafe { ffi::gsl_integration_qawo_table_free(self.w) };
        self.w = ::std::ptr::null_mut();
    }
}

impl ffi::FFI<ffi::gsl_integration_qawo_table> for IntegrationQawoTable {
    fn wrap(w: *mut ffi::gsl_integration_qawo_table) -> IntegrationQawoTable {
        IntegrationQawoTable {
            w: w
        }
    }

    fn unwrap(w: &IntegrationQawoTable) -> *mut ffi::gsl_integration_qawo_table {
        w.w
    }
}

/// CQUAD is a new doubly-adaptive general-purpose quadrature routine which can handle most types of singularities, non-numerical function
/// values such as Inf or NaN, as well as some divergent integrals. It generally requires more function evaluations than the integration
/// routines in QUADPACK, yet fails less often for difficult integrands.
/// 
/// The underlying algorithm uses a doubly-adaptive scheme in which Clenshaw-Curtis quadrature rules of increasing degree are used to compute
/// the integral in each interval. The L_2-norm of the difference between the underlying interpolatory polynomials of two successive rules
/// is used as an error estimate. The interval is subdivided if the difference between two successive rules is too large or a rule of maximum
/// degree has been reached.
pub struct CquadWorkspace {
    w: *mut ffi::gsl_integration_cquad_workspace
}

impl CquadWorkspace {
    /// This function allocates a workspace sufficient to hold the data for n intervals. The number n is not the maximum number of intervals
    /// that will be evaluated. If the workspace is full, intervals with smaller error estimates will be discarded. A minimum of 3 intervals
    /// is required and for most functions, a workspace of size 100 is sufficient.
    pub fn new(n: u64) -> Option<CquadWorkspace> {
        let tmp = unsafe { ffi::gsl_integration_cquad_workspace_alloc(n) };

        if tmp.is_null() {
            None
        } else {
            Some(CquadWorkspace {
                w: tmp
            })
        }
    }

    /// This function computes the integral of f over (a,b) within the desired absolute and relative error limits, epsabs and epsrel using
    /// the CQUAD algorithm. The function returns the final approximation, result, an estimate of the absolute error, abserr, and the number
    /// of function evaluations required, nevals.
    /// 
    /// The CQUAD algorithm divides the integration region into subintervals, and in each iteration, the subinterval with the largest estimated
    /// error is processed. The algorithm uses Clenshaw-Curits quadrature rules of degree 4, 8, 16 and 32 over 5, 9, 17 and 33 nodes respectively.
    /// Each interval is initialized with the lowest-degree rule. When an interval is processed, the next-higher degree rule is evaluated and
    /// an error estimate is computed based on the L_2-norm of the difference between the underlying interpolating polynomials of both rules.
    /// If the highest-degree rule has already been used, or the interpolatory polynomials differ significantly, the interval is bisected.
    /// 
    /// The subintervals and their results are stored in the memory provided by workspace. If the error estimate or the number of function
    /// evaluations is not needed, the pointers abserr and nevals can be set to NULL (not in rgsl).
    #[allow(unused_assignments)]
    pub fn cquad<T>(&self, f: ::function<T>, arg: &mut T, a: f64, b: f64, epsabs: f64, epsrel: f64, result: &mut f64, abserr: &mut f64,
        nevals: &mut u64) -> ::Value {
        /* Some constants that we will need. */
        static n : [i32; 4] = [4i32, 8i32, 16i32, 32i32];
        static skip : [i32; 4] = [8i32, 4i32, 2i32, 1i32];
        static idx : [i32; 4] = [0i32, 5i32, 14i32, 31i32];
        let w = ::std::f64::consts::SQRT_2 / 2f64;
        let ndiv_max = 20i32;

        /* Actual variables (as opposed to constants above). */
        let mut temp = 0f64;
        let mut neval = 0i32;
        let mut nans : [i32; 32] = [0i32; 32];

        /* Check the input arguments. */
        /*if (f == NULL)
            GSL_ERROR ("function pointer shouldn't be NULL", GSL_EINVAL);
        if (result == NULL)
            GSL_ERROR ("result pointer shouldn't be NULL", GSL_EINVAL);
        if (ws == NULL)
            GSL_ERROR ("workspace pointer shouldn't be NULL", GSL_EINVAL);*/

        unsafe {
            /* Check for unreasonable accuracy demands */
            if epsabs < 0f64 || epsrel < 0f64 {
                rgsl_error!("tolerances may not be negative", ::Value::BadTol);
            }
            if epsabs <= 0f64 && epsrel < ::DBL_EPSILON {
                rgsl_error!("unreasonable accuracy requirement", ::Value::BadTol);
            }

            /* Create the first interval. */
            let mut iv = (*self.w).ivals;
            let mut m = (a + b) / 2f64;
            let mut h = (b - a) / 2f64;
            let mut nnans = 0usize;

            for i in 0usize..(n[3] as usize + 1usize) {
                (*iv).fx[i] = f(m + XI[i] * h, arg);
                neval += 1;
                if !(*iv).fx[i].is_finite() {
                    nans[nnans] = i as i32;
                    nnans += 1;
                    (*iv).fx[i] = 0.0;
                }
            }
            Vinvfx(&(*iv).fx, &mut ((*iv).c), 0);
            let mut tmp_c = CSlice::new((::std::mem::transmute::<&mut f64, *mut f64>(&mut (*iv).c[0])).offset(3), (*iv).c.len() - 3);
            Vinvfx(&(*iv).fx, tmp_c.as_mut(), 3);
            let mut tmp_c = CSlice::new((::std::mem::transmute::<&mut f64, *mut f64>(&mut (*iv).c[0])).offset(2), (*iv).c.len() - 2);
            Vinvfx(&(*iv).fx, tmp_c.as_mut(), 2);

            for i in 0usize..nnans {
                (*iv).fx[nans[i] as usize] = ::NAN;
            }

            (*iv).a = a;
            (*iv).b = b;
            (*iv).depth = 3;
            (*iv).rdepth = 1;
            (*iv).ndiv = 0;
            (*iv).igral = 2f64 * h * (*iv).c[idx[3] as usize] * w;
            let mut nc = 0f64;

            for i in (n[2] as usize)..(n[3] as usize + 1usize) {
                temp = (*iv).c[idx[3] as usize + i];
                nc += temp * temp;
            }

            let mut ncdiff = nc;
            for i in 0usize..(n[2] as usize + 1usize) {
                temp = (*iv).c[idx[2] as usize + i] - (*iv).c[idx[3] as usize + i];
                ncdiff += temp * temp;
                nc += (*iv).c[idx[3] as usize + i] * (*iv).c[idx[3] as usize + i];
            }

            ncdiff = ncdiff.sqrt();
            nc = nc.sqrt();
            (*iv).err = ncdiff * 2f64 * h;
            if ncdiff / nc > 0.1f64 && (*iv).err < 2f64 * h * nc {
                (*iv).err = 2f64 * h * nc;
            }

            let mut t_heap = Vec::from_raw_buf((*self.w).heap, (*self.w).size as usize + 1usize);
            let heap = t_heap.as_mut();
            let mut t_ivals = Vec::from_raw_buf((*self.w).ivals, (*self.w).size as usize + 1usize);
            let ivals = t_ivals.as_mut();

            /* Initialize the heaps. */
            for i in 0usize..((*self.w).size as usize) {
                heap[i] = i as u64;
            }

            /* Initialize some global values. */
            let mut igral = (*iv).igral;
            let mut err = (*iv).err;
            let mut nivals = 132;
            let mut igral_final = 0f64;
            let mut err_final = 0f64;

            /* Main loop. */
            while nivals > 0 && err > 0f64 && !(err <= fabsf64(igral) * epsrel || err <= epsabs)
                && !(err_final > fabsf64(igral) * epsrel && err - err_final < fabsf64(igral) * epsrel)
                && !(err_final > epsabs && err - err_final < epsabs) {

                /* Put our finger on the interval with the largest error. */
                iv = (*self.w).ivals.offset(heap[0] as isize);
                m = ((*iv).a + (*iv).b) / 2f64;
                h = ((*iv).b - (*iv).a) / 2f64;

                /* printf
                   ("cquad: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
                   ws->heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
                */
                /* Should we try to increase the degree? */
                let split = if (*iv).depth < 3 {
                    /* Keep tabs on some variables. */
                    (*iv).depth += 1;
                    let d = (*iv).depth;

                    /* Get the new (missing) function values */
                    let mut it = skip[d as usize] as usize;
                    while it < 33 {
                        (*iv).fx[it] = f(m + XI[it] * h, arg);
                        neval += 1;
                        it += 2 * skip[d as usize] as usize;
                    }
                    nnans = 0;
                    it = 0;
                    while it < 33 {
                        if !(*iv).fx[it].is_finite() {
                            nans[nnans] = it as i32;
                            nnans += 1;
                            (*iv).fx[it] = 0f64;
                        }
                        it += skip[d as usize] as usize;
                    }

                    /* Compute the new coefficients. */
                    let size = idx[d as usize] as usize;
                    let mut tmp_c = CSlice::new((::std::mem::transmute::<&mut f64, *mut f64>(&mut (*iv).c[0])).offset(size as isize), (*iv).c.len() - size);
                    Vinvfx(&(*iv).fx, tmp_c.as_mut(), d);

                    /* Downdate any NaNs. */
                    if nnans > 0 {
                        downdate(tmp_c.as_mut(), n[d as usize], d, &mut nans, nnans as i32);
                        for i in 0usize..nnans {
                            (*iv).fx[nans[i] as usize] = ::NAN;
                        }
                    }

                    /* Compute the error estimate. */
                    nc = 0f64;
                    for i in (n[d as usize - 1usize] as usize)..(n[d as usize] as usize + 1usize) {
                        temp = (*iv).c[idx[d as usize] as usize + i];
                        nc += temp * temp;
                    }

                    ncdiff = nc;
                    for i in 0usize..(n[d as usize - 1usize] as usize + 1usize) {
                        temp = (*iv).c[idx[d as usize - 1usize] as usize + i] - (*iv).c[idx[d as usize] as usize + i];
                        ncdiff += temp * temp;
                        nc += (*iv).c[idx[d as usize] as usize + i] * (*iv).c[idx[d as usize] as usize + i];
                    }

                    ncdiff = ncdiff.sqrt();
                    nc = nc.sqrt();
                    (*iv).err = ncdiff * 2f64 * h;

                    /* Compute the local integral. */
                    (*iv).igral = 2f64 * h * w * (*iv).c[idx[d as usize] as usize];

                    /* Split the interval prematurely? */
                    (nc > 0f64 && ncdiff / nc > 0.1f64)
                } else {
                    /* Maximum degree reached, just split. */
                    true
                };

                /* Should we drop this interval? */
                if (m + h * XI[0]) >= (m + h * XI[1]) || (m + h * XI[31]) >= (m + h * XI[32])
                    || (*iv).err < fabsf64((*iv).igral) * ::DBL_EPSILON * 10f64 {
                    /* printf
                       ("cquad: dumping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
                       ws->heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
                       iv->depth);
                    */
                    /* Keep this interval's contribution */
                    err_final += (*iv).err;
                    igral_final += (*iv).igral;

                    /* Swap with the last element on the heap */
                    let mut t = heap[nivals - 1];
                    heap[nivals - 1] = heap[0];
                    heap[0] = t;
                    nivals -= 1;

                    /* Fix up the heap */
                    let mut it = 0;
                    while 2 * it + 1 < nivals {
                        /* Get the kids */
                        let mut j = 2 * it + 1;

                        /* If the j+1st entry exists and is larger than the jth, use it instead. */
                        if j + 1 < nivals && ivals[heap[j + 1] as usize].err >= ivals[heap[j] as usize].err {
                            j += 1;
                        }

                        /* Do we need to move the ith entry up? */
                        if ivals[heap[j] as usize].err <= ivals[heap[it] as usize].err {
                            break;
                        } else {
                            t = heap[j];
                            heap[j] = heap[it];
                            heap[it] = t;
                            it = j;
                        }
                    }
                } else if split {
                    /* Do we need to split this interval? */
                    /* Some values we will need often... */
                    let d = (*iv).depth;

                    /* Generate the interval on the left */
                    let ivl = (*self.w).ivals.offset(heap[nivals] as isize);
                    nivals += 1;
                    (*ivl).a = (*iv).a;
                    (*ivl).b = m;
                    (*ivl).depth = 0;
                    (*ivl).rdepth = (*iv).rdepth + 1;
                    (*ivl).fx[0] = (*iv).fx[0];
                    (*ivl).fx[32] = (*iv).fx[16];

                    let mut it = skip[0] as usize;
                    while it < 32usize {
                        (*ivl).fx[it] = f(((*ivl).a + (*ivl).b) / 2f64 + XI[it] * h / 2f64, arg);
                        neval += 1;
                        it += skip[0] as usize;
                    }

                    nnans = 0;
                    it = 0usize;
                    while it < 33 {
                        if !(*ivl).fx[it].is_finite() {
                            nans[nnans] = it as i32;
                            nnans += 1;
                            (*ivl).fx[it] = 0f64;
                        }
                        it += skip[0] as usize;
                    }
                    Vinvfx(&(*ivl).fx, &mut (*ivl).c, 0);
                    if nnans > 0 {
                        downdate(&mut (*ivl).c, n[0], 0, &mut nans, nnans as i32);
                        for i in 0usize..nnans {
                            (*ivl).fx[nans[i] as usize] = ::NAN;
                        }
                    }
                    for i in 0usize..(n[d as usize] as usize + 1usize) {
                        (*ivl).c[idx[d as usize] as usize + i] = 0f64;
                        for j in i..(n[d as usize] as usize + 1usize) {
                            (*ivl).c[idx[d as usize] as usize + i] += Tleft[i * 33 + j] * (*iv).c[idx[d as usize] as usize + j];
                        }
                    }
                    ncdiff = 0f64;
                    for i in 0usize..(n[0] as usize + 1usize) {
                        temp = (*ivl).c[i] - (*ivl).c[idx[d as usize] as usize + i];
                        ncdiff += temp * temp;
                    }
                    for i in (n[0] as usize + 1usize)..(n[d as usize] as usize + 1usize) {
                        temp = (*ivl).c[idx[d as usize] as usize + i];
                        ncdiff += temp * temp;
                    }
                    ncdiff = ncdiff.sqrt();
                    (*ivl).err = ncdiff * h;

                    /* Check for divergence. */
                    (*ivl).ndiv = (*iv).ndiv + (fabsf64((*iv).c[0]) > 0f64 && (*ivl).c[0] / (*iv).c[0] > 2f64) as i32;
                    if (*ivl).ndiv > ndiv_max && 2 * (*ivl).ndiv > (*ivl).rdepth {
                        /* need copysign(INFINITY, igral) */
                        *result = if igral >= 0f64 {
                            ::POSINF
                        } else {
                            ::NEGINF
                        };
                        *nevals = neval as u64;
                        return ::Value::Diverge;
                    }

                    /* Compute the local integral. */
                    (*ivl).igral = h * w * (*ivl).c[0];


                    /* Generate the interval on the right */
                    let ivr = (*self.w).ivals.offset(heap[nivals] as isize);
                    nivals += 1;
                    (*ivr).a = m;
                    (*ivr).b = (*iv).b;
                    (*ivr).depth = 0;
                    (*ivr).rdepth = (*iv).rdepth + 1;
                    (*ivr).fx[0] = (*iv).fx[16];
                    (*ivr).fx[32] = (*iv).fx[32];

                    it = skip[0] as usize;
                    while it < 32usize {
                        (*ivr).fx[it] = f(((*ivr).a + (*ivr).b) / 2f64 + XI[it] * h / 2f64, arg);
                        neval += 1;
                        it += skip[0] as usize;
                    }
                    nnans = 0;
                    it = 0usize;
                    while it < 33usize {
                        if !(*ivr).fx[it].is_finite() {
                          nans[nnans] = it as i32;
                          nnans += 1;
                          (*ivr).fx[it] = 0f64;
                        }
                        it += skip[0] as usize;
                    }
                    Vinvfx(&(*ivr).fx, &mut (*ivr).c, 0);
                    if nnans > 0 {
                        downdate(&mut (*ivr).c, n[0], 0, &mut nans, nnans as i32);
                        for i in 0usize..nnans {
                            (*ivr).fx[nans[i] as usize] = ::NAN;
                        }
                    }

                    for i in 0usize..(n[d as usize] as usize + 1usize) {
                        (*ivr).c[idx[d as usize] as usize + i] = 0f64;
                        for j in i..(n[d as usize] as usize + 1usize) {
                            (*ivr).c[idx[d as usize] as usize + i] += Tright[i * 33 + j] * (*iv).c[idx[d as usize] as usize + j];
                        }
                    }
                    ncdiff = 0f64;
                    for i in 0usize..(n[0] as usize + 1usize) {
                        temp = (*ivr).c[i] - (*ivr).c[idx[d as usize] as usize + i];
                        ncdiff += temp * temp;
                    }
                    for i in (n[0] as usize + 1usize)..(n[d as usize] as usize + 1usize) {
                        temp = (*ivr).c[idx[d as usize] as usize + i];
                        ncdiff += temp * temp;
                    }
                    ncdiff = ncdiff.sqrt();
                    (*ivr).err = ncdiff * h;

                    /* Check for divergence. */
                    (*ivr).ndiv = (*iv).ndiv + (fabsf64((*iv).c[0]) > 0f64 && (*ivr).c[0] / (*iv).c[0] > 2f64) as i32;
                
                    if (*ivr).ndiv > ndiv_max && 2i32 * (*ivr).ndiv > (*ivr).rdepth {
                        /* need copysign(INFINITY, igral) */
                        *result = if igral >= 0f64 {
                            ::POSINF
                        } else {
                            ::NEGINF
                        };
                        //if (nevals != NULL)
                        *nevals = neval as u64;
                        return ::Value::Diverge;
                    }

                    /* Compute the local integral. */
                    (*ivr).igral = h * w * (*ivr).c[0];

                    /* Fix-up the heap: we now have one interval on top that we don't need any more and two new, unsorted ones at the bottom. */

                    /* Flip the last interval to the top of the heap and sift down. */
                    let mut t = heap[nivals - 1];
                    heap[nivals - 1] = heap[0];
                    heap[0] = t;
                    nivals -= 1;

                    /* Shift this interval back down the heap. */
                    it = 0;
                    while 2 * it + 1 < nivals - 1 {
                        let mut j = 2 * it + 1;
                        if j + 1 < nivals - 1 && ivals[heap[j + 1] as usize].err >= ivals[heap[j] as usize].err {
                            j += 1;
                        }
                        if ivals[heap[j] as usize].err <= ivals[heap[it] as usize].err {
                            break;
                        } else {
                            t = heap[j];
                            heap[j] = heap[it];
                            heap[it] = t;
                            it = j;
                        }
                    }

                    /* Now grab the last interval and sift it up the heap. */
                    it = nivals - 1;
                    while it > 0 {
                        let j = (it - 1) / 2;
                        if ivals[heap[j] as usize].err < ivals[heap[it] as usize].err {
                            t = heap[j];
                            heap[j] = heap[it];
                            heap[it] = t;
                            it = j;
                        } else {
                            break;
                        }
                    }
                } else {
                    /* Otherwise, just fix-up the heap. */
                    let mut it = 0usize;
                    while 2 * it + 1 < nivals {
                        let mut j = 2 * it + 1;
                        if j + 1 < nivals && ivals[heap[j + 1] as usize].err >= ivals[heap[j] as usize].err {
                            j += 1;
                        }
                        if ivals[heap[j] as usize].err <= ivals[heap[it] as usize].err {
                            break;
                        } else {
                            let t = heap[j];
                            heap[j] = heap[it];
                            heap[it] = t;
                            it = j;
                        }
                    }
                }

                /* If the heap is about to overflow, remove the last two intervals. */
                while nivals > (*self.w).size as usize - 2usize {
                    iv = (*self.w).ivals.offset(heap[nivals - 1usize] as isize);

                    /* printf
                       ("cquad: dumping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
                       (*self.w).heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
                       iv->depth);
                    */
                    err_final += (*iv).err;
                    igral_final += (*iv).igral;
                    nivals -= 1;
                }

                /* Collect the value of the integral and error. */
                igral = igral_final;
                err = err_final;
                for i in 0usize..nivals {
                    let tmp = heap[i] as usize;
                    igral += ivals[tmp].igral;
                    err += ivals[heap[i] as usize].err;
                }
            }

            /* Dump the contents of the heap. */
            /*  for (i = 0; i < nivals; i++)
            {
                iv = &(ws->ivals[ws->heap[i]]);
                printf
                    ("cquad: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i\n",
                    i, ws->heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth,
                    iv->rdepth);
            }
            */
            /* Clean up and present the results. */
            *result = igral;
            //if (abserr != NULL)
            *abserr = err;
            //if (nevals != NULL)
            *nevals = neval as u64;

            /* All is well that ends well. */
            ::Value::Success
        }
    }
}

impl Drop for CquadWorkspace {
    fn drop(&mut self) {
        unsafe { ffi::gsl_integration_cquad_workspace_free(self.w) };
        self.w = ::std::ptr::null_mut();
    }
}

impl ffi::FFI<ffi::gsl_integration_cquad_workspace> for CquadWorkspace {
    fn wrap(w: *mut ffi::gsl_integration_cquad_workspace) -> CquadWorkspace {
        CquadWorkspace {
            w: w
        }
    }

    fn unwrap(w: &CquadWorkspace) -> *mut ffi::gsl_integration_cquad_workspace {
        w.w
    }
}

/// The fixed-order Gauss-Legendre integration routines are provided for fast integration of smooth functions with known polynomial order.
/// The n-point Gauss-Legendre rule is exact for polynomials of order 2*n-1 or less. For example, these rules are useful when integrating
/// basis functions to form mass matrices for the Galerkin method. Unlike other numerical integration routines within the library, these
/// routines do not accept absolute or relative error bounds.
pub struct GLFixedTable {
    w: *mut ffi::gsl_integration_glfixed_table
}

impl GLFixedTable {
    /// This function determines the Gauss-Legendre abscissae and weights necessary for an n-point fixed order integration scheme. If possible,
    /// high precision precomputed coefficients are used. If precomputed weights are not available, lower precision coefficients are computed
    /// on the fly.
    pub fn new(n: u64) -> Option<GLFixedTable> {
        let tmp = unsafe { ffi::gsl_integration_glfixed_table_alloc(n) };

        if tmp.is_null() {
            None
        } else {
            Some(GLFixedTable {
                w: tmp
            })
        }
    }

    /// For i in [0, …, t->n - 1], this function obtains the i-th Gauss-Legendre point xi and weight wi on the interval [a,b]. The points
    /// and weights are ordered by increasing point value. A function f may be integrated on [a,b] by summing wi * f(xi) over i.
    pub fn point(&self, a: f64, b: f64, i: u64, xi: &mut f64, wi: &mut f64) -> ::Value {
        unsafe { ffi::gsl_integration_glfixed_point(a, b, i, xi, wi, self.w) }
    }

    /// This function applies the Gauss-Legendre integration rule contained in table self and returns the result.
    pub fn glfixed<T>(&self, f: ::function<T>, arg: &mut T, a: f64, b: f64) -> f64 {
        unsafe {
            let n = (*self.w).n;
            let mut s = 0f64;

            let m = (n + 1) >> 1;
            let t_w = CSlice::new((*self.w).w, m as usize + 1usize);
            let w = t_w.as_ref();
            let t_x = CSlice::new((*self.w).x, m as usize + 1usize);
            let x = t_x.as_ref();
            let A = 0.5 * (b - a);
            let B = 0.5 * (b + a);

            /* n - odd */
            if n % 2 != 0 {
                s = w[0] * f(B, arg);

                for i in 1usize..(m as usize) {
                    let Ax = A * x[i];
                    s += w[i] * (f(B + Ax, arg) + f(B - Ax, arg));
                }
            } else {
                /* n - even */
                for i in 0usize..(m as usize) {
                    let Ax = A * x[i];
                    s += w[i] * (f(B + Ax, arg) + f(B - Ax, arg));
                }
            }
            A * s
        }
    }
}

impl Drop for GLFixedTable {
    fn drop(&mut self) {
        unsafe { ffi::gsl_integration_glfixed_table_free(self.w) };
        self.w = ::std::ptr::null_mut();
    }
}

impl ffi::FFI<ffi::gsl_integration_glfixed_table> for GLFixedTable {
    fn wrap(w: *mut ffi::gsl_integration_glfixed_table) -> GLFixedTable {
        GLFixedTable {
            w: w
        }
    }

    fn unwrap(w: &GLFixedTable) -> *mut ffi::gsl_integration_glfixed_table {
        w.w
    }
}

fn i_transform<T>(t: f64, params: &mut InternParam<T>) -> f64 {
    let f = params.func;
    let x = (1f64 - t) / t;
    let y = f(x, params.param) + f(-x, params.param);
  
    (y / t) / t
}

fn iu_transform<T>(t: f64, p: &mut InternParam<T>) -> f64 {
    let a = p.p2;
    let f = p.func;
    let x = a + (1f64 - t) / t;
    let y = f(x, p.param);

    (y / t) / t
}

fn il_transform<T>(t: f64, p: &mut InternParam<T>) -> f64 {
    let b = p.p2;
    let f = p.func;
    let x = b - (1f64 - t) / t;
    let y = f(x, p.param);

    (y / t) / t
}

fn intern_qag<T>(f: ::function<T>, arg: &mut T, a: f64, b: f64, epsabs: f64, epsrel: f64, limit: u64, f_w: &IntegrationWorkspace,
    result: &mut f64, abserr: &mut f64, q: ::integration_function<T>) -> ::Value {
    let w = f_w.w;
    let mut roundoff_type1 = 0i32;
    let mut roundoff_type2 = 0i32;
    let mut error_type = 0i32;

    *result = 0f64;
    *abserr = 0f64;

    // Initialize results
    f_w.initialise(a, b);

    if unsafe { limit > (*w).limit } {
        rgsl_error!("iteration limit exceeds available workspace", ::Value::Inval);
    }
    if epsabs <= 0f64 && (epsrel < 50f64 * ::DBL_EPSILON || epsrel < 0.5e-28f64) {
        rgsl_error!("tolerance cannot be achieved with given epsabs and epsrel", ::Value::BadTol);
    }

    // perform the first integration
    let mut result0 = 0f64;
    let mut abserr0 = 0f64;
    let mut resabs0 = 0f64;
    let mut resasc0 = 0f64;
    q(f, arg, a, b, &mut result0, &mut abserr0, &mut resabs0, &mut resasc0);

    f_w.set_initial_result(result0, abserr0);

    // Test on accuracy
    let mut tolerance = unsafe { epsabs.max(epsrel * fabsf64(result0)) };

    // need IEEE rounding here to match original quadpack behavior
    let round_off = 50f64 * ::DBL_EPSILON * resabs0;

    if abserr0 <= round_off && abserr0 > tolerance {
        *result = result0;
        *abserr = abserr0;

        rgsl_error!("cannot reach tolerance because of roundoff error on first attempt", ::Value::Round);
    } else if (abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0f64 {
        *result = result0;
        *abserr = abserr0;

        return ::Value::Success;
    } else if limit == 1 {
        *result = result0;
        *abserr = abserr0;

        rgsl_error!("a maximum of one iteration was insufficient", ::Value::MaxIter);
    }

    let mut area = result0;
    let mut errsum = abserr0;

    let iteration = 1u64;

    loop {
        let mut area1 = 0f64;
        let mut area2 = 0f64;
        let mut error1 = 0f64;
        let mut error2 = 0f64;
        let mut a_i = 0f64;
        let mut b_i = 0f64;
        let mut r_i = 0f64;
        let mut e_i = 0f64;
        let mut resabs1 = 0f64;
        let mut resasc1 = 0f64;
        let mut resabs2 = 0f64;
        let mut resasc2 = 0f64;

        // Bisect the subinterval with the largest error estimate
        f_w.retrieve(&mut a_i, &mut b_i, &mut r_i, &mut e_i);

        let a1 = a_i; 
        let b1 = 0.5 * (a_i + b_i);
        let a2 = b1;
        let b2 = b_i;

        q(f, arg, a1, b1, &mut area1, &mut error1, &mut resabs1, &mut resasc1);
        q(f, arg, a2, b2, &mut area2, &mut error2, &mut resabs2, &mut resasc2);

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

        errsum += error12 - e_i;
        area += area12 - r_i;

        if resasc1 != error1 && resasc2 != error2 {
            let delta = r_i - area12;

            if unsafe { fabsf64(delta) <= 1.0e-5f64 * fabsf64(area12) && error12 >= 0.99f64 * e_i } {
                roundoff_type1 += 1;
            }
            if iteration >= 10 && error12 > e_i {
                roundoff_type2 += 1;
            }
        }

        tolerance = unsafe { epsabs.max(epsrel * fabsf64(area)) };

        if errsum > tolerance {
            if roundoff_type1 >= 6 || roundoff_type2 >= 20 {
                error_type = 2;   /* round off error */
            }

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

            if ::util::subinterval_too_small(a1, a2, b2) {
                error_type = 3;
            }
        }

        if iteration < limit && error_type == 0i32 && errsum > tolerance {
        } else {
            break;
        }
    }
    *result = f_w.sum_results();
    *abserr = errsum;

    if errsum <= tolerance {
        ::Value::Success
    } else if error_type == 2 {
        rgsl_error!("roundoff error prevents tolerance from being achieved", ::Value::Round);
        ::Value::Round
    } else if error_type == 3 {
        rgsl_error!("bad integrand behavior found in the integration interval", ::Value::Sing);
        ::Value::Sing
    } else if iteration == limit {
        rgsl_error!("maximum number of subdivisions reached", ::Value::MaxIter);
        ::Value::MaxIter
    } else {
        rgsl_error!("could not integrate function", ::Value::Failed);
        ::Value::Failed
    }
}

#[doc(hidden)]
pub unsafe fn initialise_table(table: *mut ffi::extrapolation_table) {
    (*table).n = 0;
    (*table).nres = 0;
}

#[doc(hidden)]
pub unsafe fn append_table(table: &mut ffi::extrapolation_table, y: f64) {
    let n = (*table).n as usize;

    (*table).rlist2[n] = y;
    (*table).n += 1;
}

#[doc(hidden)]
pub fn append_interval(w: &IntegrationWorkspace, a: f64, b: f64, area: f64, error: f64) {
    w.append_interval(a, b, area, error)
}

#[doc(hidden)]
pub unsafe fn intern_qelg(table: &mut ffi::extrapolation_table, result: &mut f64, abserr: &mut f64) {
    let epstab = &mut (*table).rlist2;//Vec::from_raw_buf((*table).rlist2 as *mut f64, (*table).n as usize + 3);
    let res3la = &mut (*table).res3la;//Vec::from_raw_buf((*table).res3la as *mut f64, 3u);
    let n = (*table).n as usize - 1usize;

    let current = (*epstab)[n];

    let mut absolute = ::DBL_MAX;
    let mut relative = 5f64 * ::DBL_EPSILON * fabsf64(current);

    let newelm = n / 2usize;
    let n_orig = n;
    let mut n_final = n;

    let nres_orig = (*table).nres;

    *result = current;
    *abserr = ::DBL_MAX;

    if n < 2 {
        *result = current;
        *abserr = absolute.max(relative);
        return;
    }

    epstab[n + 2] = epstab[n];
    epstab[n] = ::DBL_MAX;

    for i in 0..newelm {
        let mut res = epstab[n - 2 * i + 2];
        let e0 = epstab[n - 2 * i - 2];
        let e1 = epstab[n - 2 * i - 1];
        let e2 = res;

        let e1abs = fabsf64(e1);
        let delta2 = e2 - e1;
        let err2 = fabsf64(delta2);
        let tol2 = fabsf64(e2).max(e1abs) * ::DBL_EPSILON;
        let delta3 = e1 - e0;
        let err3 = fabsf64(delta3);
        let tol3 = e1abs.max(fabsf64(e0)) * ::DBL_EPSILON;

        if err2 <= tol2 && err3 <= tol3 {
            /* If e0, e1 and e2 are equal to within machine accuracy, convergence is assumed.  */
            *result = res;
            absolute = err2 + err3;
            relative = 5f64 * ::DBL_EPSILON * fabsf64(res);
            *abserr = absolute.max(relative);
            return;
        }

        let e3 = epstab[n - 2 * i];
        epstab[n - 2 * i] = e1;
        let delta1 = e1 - e3;
        let err1 = fabsf64(delta1);
        let tol1 = e1abs.max(fabsf64(e3)) * ::DBL_EPSILON;

        /* If two elements are very close to each other, omit a part of the table by adjusting the value of n */
        if err1 <= tol1 || err2 <= tol2 || err3 <= tol3 {
            n_final = 2 * i;
            break;
        }

        let ss = (1f64 / delta1 + 1f64 / delta2) - 1f64 / delta3;

        /* Test to detect irregular behaviour in the table, and eventually omit a part of the table by adjusting the value of n. */

        if fabsf64(ss * e1) <= 0.0001f64 {
            n_final = 2 * i;
            break;
        }

        /* Compute a new element and eventually adjust the value of result. */

        res = e1 + 1f64 / ss;
        epstab[n - 2 * i] = res;

        {
            let error = err2 + fabsf64(res - e2) + err3;

            if error <= *abserr {
                *abserr = error;
                *result = res;
            }
        }
    }

    /* Shift the table */
    {
        let limexp = 49u64;

        if n_final == limexp as usize {
            n_final = 2usize * (limexp as usize / 2usize);
        }
    }

    if n_orig % 2 == 1 {
        for i in 0..(newelm + 1) {
          epstab[1 + i * 2] = epstab[i * 2 + 3];
        }
    } else {
        for i in 0..(newelm + 1) {
            epstab[i * 2] = epstab[i * 2 + 2];
        }
    }

    if n_orig != n_final {
        for i in 0..(n_final + 1) {
            epstab[i] = epstab[n_orig - n_final + i];
        }
    }

    (*table).n = n_final as u64 + 1;

    if nres_orig < 3 {
        res3la.as_mut()[nres_orig as usize] = *result;
        *abserr = ::DBL_MAX;
    } else {
        /* Compute error estimate */
        *abserr = fabsf64(*result - res3la[2]) + fabsf64(*result - res3la[1]) + fabsf64(*result - res3la[0]);

        res3la[0] = res3la[1];
        res3la[1] = res3la[2];
        res3la[2] = *result;
    }

    /* In QUADPACK the variable table->nres is incremented at the top of qelg, so it increases on every call. This leads to the array
       res3la being accessed when its elements are still undefined, so I have moved the update to this point so that its value more
       useful. */

    (*table).nres = nres_orig + 1;  

    *abserr = (*abserr).max(5f64 * ::DBL_EPSILON * fabsf64(*result));
}

unsafe fn test_positivity(result: f64, resabs: f64) -> bool {
    (fabsf64(result) >= (1f64 - 50f64 * ::DBL_EPSILON) * resabs)
}

#[allow(unused_variables)]
unsafe fn increase_nrmax(workspace: *mut ffi::gsl_integration_workspace) -> bool {
    let id = (*workspace).nrmax;

    let t_order = CSlice::new((*workspace).order, (*workspace).nrmax as usize + 1usize);
    let order = t_order.as_ref();
    let t_level = CSlice::new((*workspace).level, order[(*workspace).nrmax as usize] as usize + 1usize);
    let level = t_level.as_ref();

    let limit = (*workspace).limit;
    let last = (*workspace).size - 1;

    let jupbnd = if last > (1 + limit / 2) {
        limit + 1 - last
    } else {
        last
    };
  
    for k in id..(jupbnd + 1) {
        let i_max = order[(*workspace).nrmax as usize];
      
        (*workspace).i = i_max ;
        if level[i_max as usize] < (*workspace).maximum_level {
            return true;
        }
        (*workspace).nrmax += 1;
    }
    false
}

unsafe fn large_interval(workspace: *mut ffi::gsl_integration_workspace) -> bool {
    let i = (*workspace).i ;
    let level = CSlice::new((*workspace).level, i as usize + 1usize);
  
    if level.as_ref()[i as usize] < (*workspace).maximum_level {
        true
    } else {
        false
    }
}

unsafe fn reset_nrmax(workspace: *mut ffi::gsl_integration_workspace) {
    (*workspace).nrmax = 0;
    (*workspace).i = *(*workspace).order;
}

#[doc(hidden)]
pub unsafe fn return_error(t_error_type: i32) -> ::Value {
    let error_type = if t_error_type > 2 {
        t_error_type - 1
    } else {
        t_error_type
    };

    match error_type {
        0 => ::Value::Success,
        1 => {
            rgsl_error!("number of iterations was insufficient", ::Value::MaxIter);
            ::Value::MaxIter
        }
        2 => {
            rgsl_error!("cannot reach tolerance because of roundoff error", ::Value::Round);
            ::Value::Round
        }
        3 => {
            rgsl_error!("bad integrand behavior found in the integration interval", ::Value::Sing);
            ::Value::Sing
        }
        4 => {
            rgsl_error!("roundoff error detected in the extrapolation table", ::Value::Round);
            ::Value::Round
        }
        5 => {
            rgsl_error!("integral is divergent, or slowly convergent", ::Value::Round);
            ::Value::Round
        }
        _ => {
            rgsl_error!("could not integrate function", ::Value::Failed);
            ::Value::Failed
        }
    }
}

unsafe fn intern_qags<T>(f: ::function<T>, arg: &mut T, a: f64, b: f64, epsabs: f64, epsrel: f64, limit: u64, f_w: &IntegrationWorkspace,
    result: &mut f64, abserr: &mut f64, q: ::integration_function<T>) -> ::Value {
    let w = f_w.w;
    let mut ertest = 0f64;
    let mut error_over_large_intervals = 0f64;
    let mut reseps = 0f64;
    let mut abseps = 0f64;
    let mut correc = 0f64;
    let mut ktmin = 0u64;
    let mut roundoff_type1 = 0i32;
    let mut roundoff_type2 = 0i32;
    let mut roundoff_type3 = 0i32;
    let mut error_type = 0i32;
    let mut error_type2 = 0i32;
    let mut result0 = 0f64;
    let mut abserr0 = 0f64;
    let mut resabs0 = 0f64;
    let mut resasc0 = 0f64;

    let mut extrapolate = 0i32;
    let mut disallow_extrapolation = 0i32;

    let mut table : ffi::extrapolation_table = ::std::mem::zeroed();

    /* Initialize results */
    f_w.initialise(a, b);

    *result = 0f64;
    *abserr = 0f64;

    if limit > (*w).limit {
        rgsl_error!("iteration limit exceeds available workspace", ::Value::Inval);
    }

    /* Test on accuracy */
    if epsabs <= 0f64 && (epsrel < 50f64 * ::DBL_EPSILON || epsrel < 0.5e-28f64) {
        rgsl_error!("tolerance cannot be acheived with given epsabs and epsrel", ::Value::BadTol);
    }

    /* Perform the first integration */
    q(f, arg, a, b, &mut result0, &mut abserr0, &mut resabs0, &mut resasc0);

    f_w.set_initial_result(result0, abserr0);

    let mut tolerance = epsabs.max(epsrel * fabsf64(result0));

    if abserr0 <= 100f64 * ::DBL_EPSILON * resabs0 && abserr0 > tolerance {
        *result = result0;
        *abserr = abserr0;

        rgsl_error!("cannot reach tolerance because of roundoff error on first attempt", ::Value::Round);
    } else if (abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0f64 {
        *result = result0;
        *abserr = abserr0;

        return ::Value::Success;
    } else if limit == 1 {
        *result = result0;
        *abserr = abserr0;

        rgsl_error!("a maximum of one iteration was insufficient", ::Value::MaxIter);
    }

    /* Initialization */
    initialise_table(&mut table);
    append_table(&mut table, result0);

    let mut area = result0;
    let mut errsum = abserr0;

    let mut res_ext = result0;
    let mut err_ext = ::DBL_MAX;

    let positive_integrand = test_positivity(result0, resabs0);

    let mut iteration = 1u64;

    loop {
        let mut a_i = 0f64;
        let mut b_i = 0f64;
        let mut r_i = 0f64;
        let mut e_i = 0f64;
        let mut area1 = 0f64;
        let mut area2 = 0f64;
        let mut error1 = 0f64;
        let mut error2 = 0f64;
        let mut resasc1 = 0f64;
        let mut resasc2 = 0f64;
        let mut resabs1 = 0f64;
        let mut resabs2 = 0f64;

        /* Bisect the subinterval with the largest error estimate */
        f_w.retrieve(&mut a_i, &mut b_i, &mut r_i, &mut e_i);

        let t_level = CSlice::new((*w).level, (*w).i as usize + 1);
        let current_level = t_level.as_ref()[(*w).i as usize] + 1;

        let a1 = a_i;
        let b1 = 0.5f64 * (a_i + b_i);
        let a2 = b1;
        let b2 = b_i;

        iteration += 1;

        q(f, arg, a1, b1, &mut area1, &mut error1, &mut resabs1, &mut resasc1);
        q(f, arg, a2, b2, &mut area2, &mut error2, &mut resabs2, &mut resasc2);

        let area12 = area1 + area2;
        let error12 = error1 + error2;
        let last_e_i = e_i;

        /* Improve previous approximations to the integral and test for accuracy.

        We write these expressions in the same way as the original QUADPACK code so that the rounding errors are the same, which
        makes testing easier. */

        errsum = errsum + error12 - e_i;
        area = area + area12 - r_i;

        tolerance = epsabs.max(epsrel * fabsf64(area));

        if resasc1 != error1 && resasc2 != error2 {
            let delta = r_i - area12;

            if fabsf64(delta) <= 1.0e-5f64 * fabsf64(area12) && error12 >= 0.99f64 * e_i {
                if extrapolate == 0 {
                  roundoff_type1 += 1;
                } else {
                  roundoff_type2 += 1;
                }
            }
            if iteration > 10 && error12 > e_i {
                roundoff_type3 += 1;
            }
        }

        /* Test for roundoff and eventually set error flag */
        if roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20 {
            /* round off error */
            error_type = 2;
        }

        if roundoff_type2 >= 5 {
            error_type2 = 1;
        }

        /* set error flag in the case of bad integrand behaviour at a point of the integration range */
        if ::util::subinterval_too_small(a1, a2, b2) {
            error_type = 4;
        }

        /* append the newly-created intervals to the list */
        f_w.update(a1, b1, area1, error1, a2, b2, area2, error2);

        if errsum <= tolerance {
            *result = f_w.sum_results();
            *abserr = errsum;
            return return_error(error_type);
        }

        if error_type != 0 {
            break;
        }

        if iteration >= limit - 1 {
            error_type = 1;
            break;
        }

        /* set up variables on first iteration */
        if iteration == 2 {
            error_over_large_intervals = errsum;
            ertest = tolerance;
            append_table(&mut table, area);
            continue;
        }

        if disallow_extrapolation != 0 {
            continue;
        }

        error_over_large_intervals += -last_e_i;

        if current_level < (*w).maximum_level {
            error_over_large_intervals += error12;
        }

        if extrapolate == 0 {
            /* test whether the interval to be bisected next is the
             smallest interval. */

            if large_interval(w) {
                continue;
            }

            extrapolate = 1;
            (*w).nrmax = 1;
        }

        if error_type2 == 0 && error_over_large_intervals > ertest {
            if increase_nrmax(w) {
                continue;
            }
        }

        /* Perform extrapolation */
        append_table(&mut table, area);

        intern_qelg(&mut table, &mut reseps, &mut abseps);

        ktmin += 1;

        if ktmin > 5 && err_ext < 0.001f64 * errsum {
            error_type = 5;
        }

        if abseps < err_ext {
            ktmin = 0;
            err_ext = abseps;
            res_ext = reseps;
            correc = error_over_large_intervals;
            ertest = epsabs.max(epsrel * fabsf64(reseps));
            if err_ext <= ertest {
                break;
            }
        }

        /* Prepare bisection of the smallest interval. */
        if table.n == 1 {
            disallow_extrapolation = 1;
        }

        if error_type == 5 {
            break;
        }

        /* work on interval with largest error */
        reset_nrmax(w);
        extrapolate = 0;
        error_over_large_intervals = errsum;
        if iteration >= limit {
            break;
        }
    }

    *result = res_ext;
    *abserr = err_ext;

    if err_ext == ::DBL_MAX {
        *result = f_w.sum_results();
        *abserr = errsum;
        return return_error(error_type);
    }

    if error_type != 0 || error_type2 != 0 {
        if error_type2 != 0 {
            err_ext += correc;
        }

        if error_type == 0 {
            error_type = 3;
        }

        if res_ext != 0f64 && area != 0f64 {
            if err_ext / fabsf64(res_ext) > errsum / fabsf64(area) {
                *result = f_w.sum_results();
                *abserr = errsum;
                return return_error(error_type);
            }
        } else if err_ext > errsum {
            *result = f_w.sum_results();
            *abserr = errsum;
            return return_error(error_type);
        } else if area == 0f64 {
            return return_error(error_type);
        }
    }

    /*  Test on divergence. */
    let max_area = fabsf64(res_ext).max(fabsf64(area));

    if !positive_integrand && max_area < 0.01f64 * resabs0 {
        return return_error(error_type);
    }

    let ratio = res_ext / area;
    if ratio < 0.01f64 || ratio > 100f64 || errsum > fabsf64(area) {
        error_type = 6;
    }
    return_error(error_type)
}

unsafe fn intern_qagp<T>(f: ::function<T>, arg: &mut T, pts: &mut [f64], epsabs: f64, epsrel: f64, limit: u64, f_w: &IntegrationWorkspace,
    result: &mut f64, abserr: &mut f64, q: ::integration_function<T>) -> ::Value {
    let w = f_w.w;
    let mut reseps = 0f64;
    let mut abseps = 0f64;
    let mut correc = 0f64;
    let mut ktmin = 0u64;
    let mut roundoff_type1 = 0i32;
    let mut roundoff_type2 = 0i32;
    let mut roundoff_type3 = 0i32;
    let mut error_type = 0i32;
    let mut error_type2 = 0i32;

    let mut extrapolate = 0i32;
    let mut disallow_extrapolation = 0i32;

    let mut table : ffi::extrapolation_table = ::std::mem::zeroed();

    /* number of intervals */
    let nint = pts.len() as u64 - 1u64;

    /* temporarily alias ndin to level */
    let mut t_ndin = Vec::from_raw_buf((*w).level, pts.len());
    let ndin = t_ndin.as_mut();

    /* Initialize results */
    *result = 0f64;
    *abserr = 0f64;

    /* Test on validity of parameters */
    if limit > (*w).limit {
        rgsl_error!("iteration limit exceeds available workspace", ::Value::Inval);
    }

    if pts.len() as u64 > (*w).limit {
        rgsl_error!("pts length exceeds size of workspace", ::Value::Inval);
    }

    if epsabs <= 0f64 && (epsrel < 50f64 * ::DBL_EPSILON || epsrel < 0.5e-28f64) {
        rgsl_error!("tolerance cannot be acheived with given epsabs and epsrel", ::Value::BadTol);
    }

    /* Check that the integration range and break points are an ascending sequence */
    for i in 0usize..(nint as usize) {
        if pts[i + 1] < pts[i] {
            rgsl_error!("points are not in an ascending sequence", ::Value::Inval);
        }
    }

    /* Perform the first integration */
    let mut result0 = 0f64;
    let mut abserr0 = 0f64;
    let mut resabs0 = 0f64;

    f_w.initialise(0f64, 0f64);

    for i in 0usize..(nint as usize) {
        let mut area1 = 0f64;
        let mut error1 = 0f64;
        let mut resabs1 = 0f64;
        let mut resasc1 = 0f64;
        let a1 = pts[i];
        let b1 = pts[i + 1];

        q(f, arg, a1, b1, &mut area1, &mut error1, &mut resabs1, &mut resasc1);

        result0 = result0 + area1;
        abserr0 = abserr0 + error1;
        resabs0 = resabs0 + resabs1;

        f_w.append_interval(a1, b1, area1, error1);

        if error1 == resasc1 && error1 != 0f64 {
            ndin[i] = 1;
        } else {
            ndin[i] = 0;
        }
    }

    /* Compute the initial error estimate */
    let mut errsum = 0f64;
    let mut t_elist = CSlice::new((*w).elist, nint as usize);
    let elist = t_elist.as_mut();
    let mut t_level = CSlice::new((*w).level, nint as usize);
    let level = t_level.as_mut();

    for i in 0usize..(nint as usize) {
        if ndin[i] != 0 {
            elist[i] = abserr0;
        }
        errsum = errsum + elist[i];
    }

    for i in 0usize..(nint as usize) {
        level[i] = 0u64;
    }

    /* Sort results into order of decreasing error via the indirection array order[] */
    f_w.sort_results();

    /* Test on accuracy */
    let mut tolerance = epsabs.max(epsrel * fabsf64(result0));

    if abserr0 <= 100f64 * ::DBL_EPSILON * resabs0 && abserr0 > tolerance {
        *result = result0;
        *abserr = abserr0;

        rgsl_error!("cannot reach tolerance because of roundoff error on first attempt", ::Value::Round);
    } else if abserr0 <= tolerance {
        *result = result0;
        *abserr = abserr0;

        return ::Value::Success;
    } else if limit == 1 {
        *result = result0;
        *abserr = abserr0;

        rgsl_error!("a maximum of one iteration was insufficient", ::Value::MaxIter);
    }

    /* Initialization */
    initialise_table(&mut table);
    append_table(&mut table, result0);

    let mut area = result0;

    let mut res_ext = result0;
    let mut err_ext = ::DBL_MAX;

    let mut error_over_large_intervals = errsum;
    let mut ertest = tolerance;

    let positive_integrand = test_positivity(result0, resabs0);

    let mut iteration = nint - 1; 

    loop {
        let mut a_i = 0f64;
        let mut b_i = 0f64;
        let mut r_i = 0f64;
        let mut e_i = 0f64;
        let mut area1 = 0f64;
        let mut area2 = 0f64;
        let mut error1 = 0f64;
        let mut error2 = 0f64;
        let mut resasc1 = 0f64;
        let mut resasc2 = 0f64;
        let mut resabs1 = 0f64;
        let mut resabs2 = 0f64;

        /* Bisect the subinterval with the largest error estimate */
        f_w.retrieve(&mut a_i, &mut b_i, &mut r_i, &mut e_i);

        let current_level = level[(*w).i as usize] + 1u64;

        let a1 = a_i;
        let b1 = 0.5f64 * (a_i + b_i);
        let a2 = b1;
        let b2 = b_i;

        iteration += 1;

        q(f, arg, a1, b1, &mut area1, &mut error1, &mut resabs1, &mut resasc1);
        q(f, arg, a2, b2, &mut area2, &mut error2, &mut resabs2, &mut resasc2);

        let area12 = area1 + area2;
        let error12 = error1 + error2;
        let last_e_i = e_i;

        /* Improve previous approximations to the integral and test for accuracy.

         We write these expressions in the same way as the original QUADPACK code so that the rounding errors are the same, which
         makes testing easier. */

        errsum = errsum + error12 - e_i;
        area = area + area12 - r_i;

        tolerance = epsabs.max(epsrel * fabsf64(area));

        if resasc1 != error1 && resasc2 != error2 {
            let delta = r_i - area12;

            if fabsf64(delta) <= 1.0e-5f64 * fabsf64(area12) && error12 >= 0.99f64 * e_i {
                if extrapolate == 0 {
                    roundoff_type1 += 1;
                } else {
                  roundoff_type2 += 1;
                }
            }

            if nint + 1 > 10 && error12 > e_i {
                roundoff_type3 += 1;
            }
        }

        /* Test for roundoff and eventually set error flag */
        if roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20 {
            /* round off error */
            error_type = 2;
        }

        if roundoff_type2 >= 5 {
            error_type2 = 1;
        }

        /* set error flag in the case of bad integrand behaviour at a point of the integration range */
        if ::util::subinterval_too_small(a1, a2, b2) {
            error_type = 4;
        }

        /* append the newly-created intervals to the list */
        f_w.update(a1, b1, area1, error1, a2, b2, area2, error2);

        if errsum <= tolerance {
            *result = f_w.sum_results();
            *abserr = errsum;
            return return_error(error_type);
        }

        if error_type != 0 {
            break;
        }

        if iteration >= limit - 1 {
            error_type = 1;
            break;
        }

        if disallow_extrapolation != 0 {
            continue;
        }

        error_over_large_intervals += -last_e_i;

        if current_level < (*w).maximum_level {
            error_over_large_intervals += error12;
        }

        if extrapolate == 0 {
            /* test whether the interval to be bisected next is the smallest interval. */
            if large_interval(w) {
                continue;
            }

            extrapolate = 1;
            (*w).nrmax = 1;
        }

        /* The smallest interval has the largest error.  Before bisecting decrease the sum of the errors over the larger
         intervals (error_over_large_intervals) and perform extrapolation. */
        if error_type2 == 0 && error_over_large_intervals > ertest {
            if increase_nrmax(w) {
                continue;
            }
        }

        /* Perform extrapolation */
        append_table (&mut table, area);

        if table.n > 2 {
            intern_qelg(&mut table, &mut reseps, &mut abseps);

            ktmin += 1;

            if ktmin > 5 && err_ext < 0.001f64 * errsum {
                error_type = 5;
            }

            if abseps < err_ext {
                ktmin = 0;
                err_ext = abseps;
                res_ext = reseps;
                correc = error_over_large_intervals;
                ertest = epsabs.max(epsrel * fabsf64(reseps));
                if err_ext <= ertest {
                    break;
                }
            }

            /* Prepare bisection of the smallest interval. */
            if table.n == 1 {
                disallow_extrapolation = 1;
            }

            if error_type == 5 {
                break;
            }
        }

        reset_nrmax(w);
        extrapolate = 0;
        error_over_large_intervals = errsum;
        if iteration >= limit {
            break;
        }
    }

    *result = res_ext;
    *abserr = err_ext;

    if err_ext == ::DBL_MAX {
        *result = f_w.sum_results();
        *abserr = errsum;
        return return_error(error_type);
    }

    if error_type != 0 || error_type2 != 0 {
        if error_type2 != 0 {
            err_ext += correc;
        }

        if error_type == 0 {
            error_type = 3;
        }

        if *result != 0f64 && area != 0f64 {
            if err_ext / fabsf64(res_ext) > errsum / fabsf64(area) {
                *result = f_w.sum_results();
                *abserr = errsum;
                return return_error(error_type);
            }
        } else if err_ext > errsum {
            *result = f_w.sum_results();
            *abserr = errsum;
            return return_error(error_type);
        } else if area == 0f64 {
            return return_error(error_type);
        }
    }

    /*  Test on divergence. */
    {
        let max_area = fabsf64(res_ext).max(fabsf64(area));

        if !positive_integrand && max_area < 0.01f64 * resabs0 {
            return return_error(error_type);
        }
    }

    {
        let ratio = res_ext / area;

        if ratio < 0.01f64 || ratio > 100f64 || errsum > fabsf64(area) {
            error_type = 6;
        }
    }

    return_error(error_type)
}

unsafe fn qc25c<T>(f: ::function<T>, arg: &mut T, a: f64, b: f64, c: f64, result: &mut f64, abserr: &mut f64, err_reliable: &mut i32) {
    let cc = (2f64 * c - b - a) / (b - a);

    if fabsf64(cc) > 1.1f64 {
        let mut resabs = 0f64;
        let mut resasc = 0f64;

        let mut tmp_arg = InternParam{func: f, param: arg, p2: c};

        ::integration::qk15(fn_cauchy, &mut tmp_arg, a, b, result, abserr, &mut resabs, &mut resasc);

        if *abserr == resasc {
            *err_reliable = 0;
        } else  {
            *err_reliable = 1;
        }
    } else {
        let mut cheb12 : [f64; 13] = [0f64; 13];
        let mut cheb24 : [f64; 25] = [0f64; 25];
        let mut moment : [f64; 25] = [0f64; 25];
        let mut res12 = 0f64;
        let mut res24 = 0f64;

        gsl_integration_qcheb(f, arg, a, b, &mut cheb12, &mut cheb24);
        compute_moments(cc, &mut moment);

        for i in 0..13 {
            res12 += cheb12[i] * moment[i];
        }

        for i in 0..25 {
            res24 += cheb24[i] * moment[i];
        }

        *result = res24;
        *abserr = fabsf64(res24 - res12) ;
        *err_reliable = 0;
    }
}

fn fn_cauchy<T>(x: f64, p: &mut InternParam<T>) -> f64 {
    let c = p.p2;
    let f = p.func;

    f(x, p.param) / (x - c)
}

unsafe fn compute_moments(cc: f64, moment: &mut [f64]) {
    let mut a0 = logf64(fabsf64((1f64 - cc) / (1f64 + cc)));
    let mut a1 = 2f64 + a0 * cc;

    moment[0] = a0;
    moment[1] = a1;

    for k in 2usize..25usize {
        let a2 = if (k % 2) == 0 {
            2f64 * cc * a1 - a0
        } else {
            let km1 = k as f64 - 1f64;

            2f64 * cc * a1 - a0 - 4.0 / (km1 * km1 - 1f64)
        };

        moment[k] = a2;

        a0 = a1;
        a1 = a2;
    }
}

// This function computes the 12-th order and 24-th order Chebyshev approximations to f(x) on [a,b]
fn gsl_integration_qcheb<T>(f: ::function<T>, arg: &mut T, a: f64, b: f64, cheb12: &mut [f64], cheb24: &mut [f64]) {
    let mut fval : [f64; 25] = [0f64; 25];
    let mut v : [f64; 12] = [0f64; 12];

    /* These are the values of cos(pi*k/24) for k=1..11 needed for the Chebyshev expansion of f(x) */
    let x : [f64; 11] = [
        0.9914448613738104f64,     
        0.9659258262890683f64,
        0.9238795325112868f64,     
        0.8660254037844386f64,
        0.7933533402912352f64,     
        0.7071067811865475f64,
        0.6087614290087206f64,     
        0.5000000000000000f64,
        0.3826834323650898f64,     
        0.2588190451025208f64,
        0.1305261922200516f64];
  
    let center = 0.5f64 * (b + a);
    let half_length =  0.5f64 * (b - a);
  
    fval[0] = 0.5f64 * f(b, arg);
    fval[12] = f(center, arg);
    fval[24] = 0.5f64 * f(a, arg);

    for i in 1usize..12usize {
        let j = 24usize - i;
        let u = half_length * x[i - 1usize];

        fval[i] = f(center + u, arg);
        fval[j] = f(center - u, arg);
    }

    for i in 1usize..12usize {
        let j = 24usize - i;

        v[i] = fval[i] - fval[j];
        fval[i] = fval[i] + fval[j];
    }

    {
        let alam1 = v[0] - v[8];
        let alam2 = x[5] * (v[2] - v[6] - v[10]);

        cheb12[3] = alam1 + alam2;
        cheb12[9] = alam1 - alam2;
    }

    {
        let alam1 = v[1] - v[7] - v[9];
        let alam2 = v[3] - v[5] - v[11];
        {
            let alam = x[2] * alam1 + x[8] * alam2;

            cheb24[3] = cheb12[3] + alam;
            cheb24[21] = cheb12[3] - alam;
        }

        {
            let alam = x[8] * alam1 - x[2] * alam2;
            cheb24[9] = cheb12[9] + alam;
            cheb24[15] = cheb12[9] - alam;
        }
    }

    {
        let part1 = x[3] * v[4];
        let part2 = x[7] * v[8];
        let part3 = x[5] * v[6];
        
        {
            let alam1 = v[0] + part1 + part2;
            let alam2 = x[1] * v[2] + part3 + x[9] * v[10];
          
            cheb12[1] = alam1 + alam2;
            cheb12[11] = alam1 - alam2;
        }
        
        {
            let alam1 = v[0] - part1 + part2;
            let alam2 = x[9] * v[2] - part3 + x[1] * v[10];

            cheb12[5] = alam1 + alam2;
            cheb12[7] = alam1 - alam2;
        }
    }

    {
        let alam = x[0] * v[1] + x[2] * v[3] + x[4] * v[5] + x[6] * v[7] + x[8] * v[9] + x[10] * v[11];

        cheb24[1] = cheb12[1] + alam;
        cheb24[23] = cheb12[1] - alam;
    }

    {
        let alam = x[10] * v[1] - x[8] * v[3] + x[6] * v[5] - x[4] * v[7] + x[2] * v[9] - x[0] * v[11];

        cheb24[11] = cheb12[11] + alam;
        cheb24[13] = cheb12[11] - alam;
    }

    {
        let alam = x[4] * v[1] - x[8] * v[3] - x[0] * v[5] - x[10] * v[7] + x[2] * v[9] + x[6] * v[11];

        cheb24[5] = cheb12[5] + alam;
        cheb24[19] = cheb12[5] - alam;
    }

    {
        let alam = x[6] * v[1] - x[2] * v[3] - x[10] * v[5] + x[0] * v[7] - x[8] * v[9] - x[4] * v[11];

        cheb24[7] = cheb12[7] + alam;
        cheb24[17] = cheb12[7] - alam;
    }

    for i in 0usize..6usize {
        let j = 12usize - i;

        v[i] = fval[i] - fval[j];
        fval[i] = fval[i] + fval[j];
    }

    {
        let alam1 = v[0] + x[7] * v[4];
        let alam2 = x[3] * v[2];

        cheb12[2] = alam1 + alam2;
        cheb12[10] = alam1 - alam2;
    }

    cheb12[6] = v[0] - v[4];

    {
        let alam = x[1] * v[1] + x[5] * v[3] + x[9] * v[5];

        cheb24[2] = cheb12[2] + alam;
        cheb24[22] = cheb12[2] - alam;
    }

    {
        let alam = x[5] * (v[1] - v[3] - v[5]);

        cheb24[6] = cheb12[6] + alam;
        cheb24[18] = cheb12[6] - alam;
    }

    {
        let alam = x[9] * v[1] - x[5] * v[3] + x[1] * v[5];

        cheb24[10] = cheb12[10] + alam;
        cheb24[14] = cheb12[10] - alam;
    }

    for i in 0usize..3usize {
        let j = 6 - i;

        v[i] = fval[i] - fval[j];
        fval[i] = fval[i] + fval[j];
    }

    cheb12[4] = v[0] + x[7] * v[2];
    cheb12[8] = fval[0] - x[7] * fval[2];

    {
        let alam = x[3] * v[1];

        cheb24[4] = cheb12[4] + alam;
        cheb24[20] = cheb12[4] - alam;
    }

    {
        let alam = x[7] * fval[1] - fval[3];

        cheb24[8] = cheb12[8] + alam;
        cheb24[16] = cheb12[8] - alam;
    }

    cheb12[0] = fval[0] + fval[2];

    {
        let alam = fval[1] + fval[3];

        cheb24[0] = cheb12[0] + alam;
        cheb24[24] = cheb12[0] - alam;
    }

    cheb12[12] = v[0] - v[2];
    cheb24[12] = cheb12[12];

    let mut tmp = 1f64 / 6f64;

    for i in 0usize..12usize {
        cheb12[i] *= tmp;
    }

    tmp = 1f64 / 12f64;

    cheb12[0] *= tmp;
    cheb12[12] *= tmp;

    for i in 1usize..24usize {
        cheb24[i] *= tmp;
    }

    tmp = 1f64 / 24f64;
    cheb24[0] *= tmp;
    cheb24[24] *= tmp;
}

struct fn_qaws_params<'r, T:'r> {
    function: ::function<T>,
    a: f64,
    b: f64,
    table: *mut ffi::gsl_integration_qaws_table,
    arg: &'r mut T
}

unsafe fn qc25s<T>(f: ::function<T>, arg: &mut T, a: f64, b: f64, a1: f64, b1: f64, t: *mut ffi::gsl_integration_qaws_table, result: &mut f64,
    abserr: &mut f64, err_reliable: &mut bool)
{
    let mut fn_params = fn_qaws_params{function: f, a: a, b: b, table: t, arg: arg};

    if a1 == a && ((*t).alpha != 0f64 || (*t).mu != 0) {
        let mut cheb12 : [f64; 13] = [0f64; 13];
        let mut cheb24 : [f64; 25] = [0f64; 25];

        let factor = powf64(0.5f64 * (b1 - a1), (*t).alpha + 1f64);

        gsl_integration_qcheb(fn_qaws_R, &mut fn_params, a1, b1, &mut cheb12, &mut cheb24);

        if (*t).mu == 0 {
            let mut res12 = 0f64;
            let mut res24 = 0f64;
            let u = factor;

            qc25s_compute_result(&(*t).ri, &cheb12, &cheb24, &mut res12, &mut res24);

            *result = u * res24;
            *abserr = fabsf64(u * (res24 - res12));
        } else {
            let mut res12a = 0f64;
            let mut res24a = 0f64;
            let mut res12b = 0f64;
            let mut res24b = 0f64;

            let u = factor * logf64(b1 - a1);
            let v = factor;

            qc25s_compute_result(&(*t).ri, &cheb12, &cheb24, &mut res12a, &mut res24a);
            qc25s_compute_result(&(*t).rg, &cheb12, &cheb24, &mut res12b, &mut res24b);

            *result = u * res24a + v * res24b;
            *abserr = fabsf64(u * (res24a - res12a)) + fabsf64(v * (res24b - res12b));
        }

        *err_reliable = false;
    } else if b1 == b && ((*t).beta != 0.0 || (*t).nu != 0) {
        let mut cheb12 : [f64; 13] = [0f64; 13];
        let mut cheb24 : [f64; 25] = [0f64; 25];
        let factor = powf64(0.5f64 * (b1 - a1), (*t).beta + 1f64);

        gsl_integration_qcheb(fn_qaws_L, &mut fn_params, a1, b1, &mut cheb12, &mut cheb24);

        if (*t).nu == 0 {
            let mut res12 = 0f64;
            let mut res24 = 0f64;
            let u = factor;

            qc25s_compute_result(&(*t).rj, &cheb12, &cheb24, &mut res12, &mut res24);

            *result = u * res24;
            *abserr = fabsf64(u * (res24 - res12));
        } else {
            let mut res12a = 0f64;
            let mut res24a = 0f64;
            let mut res12b = 0f64;
            let mut res24b = 0f64;

            let u = factor * logf64(b1 - a1);
            let v = factor;

            qc25s_compute_result(&(*t).rj, &cheb12, &cheb24, &mut res12a, &mut res24a);
            qc25s_compute_result(&(*t).rh, &cheb12, &cheb24, &mut res12b, &mut res24b);

            *result = u * res24a + v * res24b;
            *abserr = fabsf64(u * (res24a - res12a)) + fabsf64(v * (res24b - res12b));
        }
        *err_reliable = false;
    }
    else {
        let mut resabs = 0f64;
        let mut resasc = 0f64;

        ::integration::qk15(fn_qaws, &mut fn_params, a1, b1, result, abserr, &mut resabs, &mut resasc);

        if *abserr == resasc {
            *err_reliable = false;
        } else {
            *err_reliable = true;
        }
    }
}

fn fn_qaws<T>(x: f64, p: &mut fn_qaws_params<T>) -> f64 {
    let f = p.function;
    let t = p.table;

    let mut factor = 1f64;
  
    unsafe {
        if (*t).alpha != 0f64 {
            factor *= powf64(x - p.a, (*t).alpha);
        }

        if (*t).beta != 0f64 {
            factor *= powf64(p.b - x, (*t).beta);
        }

        if (*t).mu == 1 {
            factor *= logf64(x - p.a);
        }

        if (*t).nu == 1 {
            factor *= logf64(p.b - x);
        }

        factor * f(x, p.arg)
    }
}

fn fn_qaws_R<T>(x: f64, p: &mut fn_qaws_params<T>) -> f64 {
    let f = p.function;
    let t = p.table;

    let mut factor = 1f64;
  
    unsafe {
        if (*t).beta != 0f64 {
            factor *= powf64(p.b - x, (*t).beta);
        }

        if (*t).nu == 1 {
            factor *= logf64(p.b - x);
        }

        factor * f(x, p.arg)
    }
}

fn fn_qaws_L<T>(x: f64, p: &mut fn_qaws_params<T>) -> f64 {
    let f = p.function;
    let t = p.table;

    let mut factor = 1f64;
  
    unsafe {
        if (*t).alpha != 0f64 {
            factor *= powf64(x - p.a, (*t).alpha);
        }

        if (*t).mu == 1 {
            factor *= logf64(x - p.a);
        }

        factor * f(x, p.arg)
    }
}

unsafe fn qc25s_compute_result(r: &[f64], cheb12: &[f64], cheb24: &[f64], result12: &mut f64, result24: &mut f64) {
    let mut res12 = 0f64;
    let mut res24 = 0f64;

    for i in 0usize..13usize {
        res12 += r[i] * cheb12[i];
    }

    for i in 0usize..25usize {
        res24 += r[i] * cheb24[i];
    }

    *result12 = res12;
    *result24 = res24;
}

struct fn_fourier_params<'r, T:'r> {
    omega: f64,
    function: ::function<T>,
    arg: &'r mut T
}

unsafe fn qc25f<T>(f: ::function<T>, arg: &mut T, a: f64, b: f64, wf: *mut ffi::gsl_integration_qawo_table, level: u64, result: &mut f64,
    abserr: &mut f64, resabs: &mut f64, resasc: &mut f64) {
    let center = 0.5f64 * (a + b);
    let half_length = 0.5f64 * (b - a);
    let omega = (*wf).omega ;
  
    let par = omega * half_length;

    if fabsf64(par) < 2f64 {
        let mut fn_params = fn_fourier_params{omega: omega, function: f, arg: arg};

        ::integration::qk15(if (*wf).sine == ::IntegrationQawo::Sine {
                fn_sin
            } else {
                fn_cos
            }, &mut fn_params, a, b, result, abserr, resabs, resasc);
    } else {
        let mut cheb12 : [f64; 13] = [0f64; 13];
        let mut cheb24 : [f64; 25] = [0f64; 25];

        gsl_integration_qcheb(f, arg, a, b, &mut cheb12, &mut cheb24);

        if level >= (*wf).n {
            /* table overflow should not happen, check before calling */
            rgsl_error!("table overflow in internal function", ::Value::Sanity);
        }

        /* obtain moments from the table */
        let t_moment = CSlice::new((*wf).chebmo.offset(25isize * level as isize), 25);
        let moment = t_moment.as_ref();

        let mut res12_cos = cheb12[12] * moment[12];
        let mut res12_sin = 0f64;

        for i in 0usize..6usize {
            let k = 10usize - 2usize * i;

            res12_cos += cheb12[k] * moment[k];
            res12_sin += cheb12[k + 1] * moment[k + 1];
        }

        let mut res24_cos = cheb24[24] * moment[24];
        let mut res24_sin = 0f64;

        let mut result_abs = fabsf64(cheb24[24]) ;

        for i in 0usize..12usize {
            let k = 22usize - 2usize * i;

            res24_cos += cheb24[k] * moment[k];
            res24_sin += cheb24[k + 1] * moment[k + 1];
            result_abs += fabsf64(cheb24[k]) + fabsf64(cheb24[k + 1]);
        }

        let est_cos = fabsf64(res24_cos - res12_cos);
        let est_sin = fabsf64(res24_sin - res12_sin);

        let c = half_length * cosf64(center * omega);
        let s = half_length * sinf64(center * omega);

        if (*wf).sine == ::IntegrationQawo::Sine {
            *result = c * res24_sin + s * res24_cos;
            *abserr = fabsf64(c * est_sin) + fabsf64(s * est_cos);
        } else {
            *result = c * res24_cos - s * res24_sin;
            *abserr = fabsf64(c * est_cos) + fabsf64(s * est_sin);
        }

        *resabs = result_abs * half_length;
        *resasc = ::DBL_MAX;
    }
}

fn fn_sin<T>(x: f64, p: &mut fn_fourier_params<T>) -> f64 {
    let f = p.function;
    let w = p.omega;
    let wx = w * x;
    let sinwx = unsafe { sinf64(wx) };

    f(x, p.arg) * sinwx
}

fn fn_cos<T>(x: f64, p: &mut fn_fourier_params<T>) -> f64 {
    let f = p.function;
    let w = p.omega;
    let wx = w * x;
    let coswx = unsafe { cosf64(wx) };

    f(x, p.arg) * coswx
}

/* Compute the product of the fx with one of the inverse Vandermonde-like matrices. */
fn Vinvfx(fx: &[f64], c: &mut [f64], d: i32) {
    match d {
        0i32 => {
            for i in 0usize..5usize {
                c[i] = 0f64;
                for j in 0usize..5usize {
                    c[i] += V1inv[i * 5 + j] * fx[j * 8];
                }
            }
        }
        1i32 => {
            for i in 0usize..9usize {
                c[i] = 0f64;
                for j in 0usize..9usize {
                    c[i] += V2inv[i * 9 + j] * fx[j * 4];
                }
            }
        }
        2i32 => {
            for i in 0usize..17usize {
                c[i] = 0f64;
                for j in 0usize..17usize {
                    c[i] += V3inv[i * 17 + j] * fx[j * 2];
                }
            }
        }
        3i32 => {
            for i in 0usize..33usize {
                c[i] = 0f64;
                for j in 0usize..33usize {
                    c[i] += V4inv[i * 33 + j] * fx[j];
                }
            }
        }
        _ => {}
    }
}

/* Downdate the interpolation given by the n coefficients c by removing the nodes with indices in nans. */
fn downdate(c: &mut [f64], t_n: i32, d: i32, nans: &mut [i32], nnans: i32) {
    static bidx : [i32; 4] = [0i32, 6i32, 16i32, 34i32];
    let mut b_new : [f64; 34] = [0f64; 34];
    let mut n = t_n;

    for i in 0..(n + 2) {
        b_new[i as usize] = bee[bidx[d as usize] as usize + i as usize];
    }
    for i in 0usize..(nnans as usize) {
        let tmp = nans[i] as usize;
        b_new[n as usize + 1usize] = b_new[n as usize + 1usize] / Lalpha[n as usize];
        b_new[n as usize] = (b_new[n as usize] + XI[tmp] * b_new[n as usize + 1usize]) / Lalpha[n as usize - 1usize];
        for j in (n - 1)..0 {
            b_new[j as usize] = (b_new[j as usize] + XI[nans[i] as usize] * b_new[j as usize + 1usize] - Lgamma[j as usize + 1usize]
                * b_new[j as usize + 2usize]) / Lalpha[j as usize - 1usize];
        }
        for j in 0..(n + 1) {
            b_new[j as usize] = b_new[j as usize + 1usize];
        }

        let alpha = c[n as usize] / b_new[n as usize];
        for j in 0..n {
            c[j as usize] -= alpha * b_new[j as usize];
        }
        c[n as usize] = 0f64;
        n -= 1;
    }

}