#include "Math/SpecFunc.h"
#include <iostream>
#include <cmath>

using std::cout; 
using std::endl;

// input for wolframalpha.com:  Table[{N[x/10], N[CosIntegral[x/10], 18]}, {x, -100, 100}],
// delete imaginary part for x<0, remove entry for x=0
double civalues[200][2]=
{{-10.,-0.045456433004455373},
{-9.9,-0.036763956296836376},
{-9.8,-0.027519181109809388},
{-9.7,-0.017804097705837473},
{-9.6,-0.007707036058533522},
{-9.5,0.002678058835650657},
{-9.4,0.013252418669884902},
{-9.3,0.02391330446927572},
{-9.2,0.03455491341975988},
{-9.1,0.045069332542612144},
{-9.,0.055347531333133607},
{-8.9,0.06528038500431541},
{-8.8,0.074759719566173204},
{-8.7,0.083679369633444155},
{-8.6,0.091936239592257399},
{-8.5,0.099431358573421916},
{-8.4,0.106070919578643913},
{-8.3,0.11176729308811125},
{-8.2,0.116440005544566539},
{-8.1,0.120016673260596569},
{-8.,0.122433882532009557},
{-7.9,0.123638007059717843},
{-7.8,0.123585954183608557},
{-7.7,0.122245831911846316},
{-7.6,0.119597529284565836},
{-7.5,0.11563320323793427},
{-7.4,0.110357665828378243},
{-7.3,0.103788666432027628},
{-7.2,0.095957064345180793},
{-7.1,0.08690688807134784},
{-7.,0.076695278482184518},
{-6.9,0.065392313975951504},
{-6.8,0.053080716720199198},
{-6.7,0.039855440047043456},
{-6.6,0.025823138061263264},
{-6.5,0.011101519514930109},
{-6.4,-0.004181411011335063},
{-6.3,-0.019888220609842217},
{-6.2,-0.035873019273454991},
{-6.1,-0.051982528980021969},
{-6.,-0.068057243893247126},
{-5.9,-0.083932674118556494},
{-5.8,-0.099440664689378585},
{-5.7,-0.114410780761679062},
{-5.6,-0.128671749369807815},
{-5.5,-0.142052947551519255},
{-5.4,-0.154385926190724433},
{-5.3,-0.165505958558927253},
{-5.2,-0.175253602265659485},
{-5.1,-0.183476263159298894},
{-5.,-0.190029749656643879},
{-4.9,-0.194779806026237225},
{-4.8,-0.197603613309935236},
{-4.7,-0.198391246842472852},
{-4.6,-0.197047079722356195},
{-4.5,-0.193491122101738757},
{-4.4,-0.187660286800440685},
{-4.3,-0.179509572512633253},
{-4.2,-0.169013156767156739},
{-4.1,-0.15616539182812111},
{-4.,-0.140981697886930412},
{-3.9,-0.123499349207815143},
{-3.8,-0.103778150356897706},
{-3.7,-0.08190100128429845},
{-3.6,-0.057974351859800879},
{-3.5,-0.032128548512481116},
{-3.4,-0.004518077930741954},
{-3.3,0.024678284607958114},
{-3.2,0.055257411719942492},
{-3.1,0.086991831195536989},
{-3.,0.119629786008000328},
{-2.9,0.152895324159588311},
{-2.8,0.186488389643175768},
{-2.7,0.220084878632961618},
{-2.6,0.253336616062584192},
{-2.5,0.285871196365383495},
{-2.4,0.317291617436697984},
{-2.3,0.347175617540316224},
{-2.2,0.375074599049832154},
{-2.1,0.40051198784439639},
{-2.,0.422980828774864996},
{-1.9,0.441940349681598846},
{-1.8,0.456811129418336893},
{-1.7,0.46696836417695463},
{-1.6,0.471732516931877803},
{-1.5,0.470356317195399887},
{-1.4,0.462006585094677276},
{-1.3,0.445738567528534523},
{-1.2,0.420459182894240503},
{-1.1,0.384873377424650815},
{-1.,0.337403922900968135},
{-0.9,0.27606783046777286},
{-0.8,0.198278615952467177},
{-0.7,0.100514707008897833},
{-0.6,-0.022270706959279763},
{-0.5,-0.177784078806612901},
{-0.4,-0.378809346425244332},
{-0.3,-0.649172932971161745},
{-0.2,-1.042205595672781975},
{-0.1,-1.727868386657296639},
{0.1,-1.727868386657296639},
{0.2,-1.042205595672781975},
{0.3,-0.649172932971161745},
{0.4,-0.3788093464252443321},
{0.5,-0.1777840788066129013},
{0.6,-0.02227070695927976253},
{0.7,0.1005147070088978327},
{0.8,0.198278615952467177},
{0.9,0.2760678304677728602},
{1.,0.3374039229009681347},
{1.1,0.3848733774246508155},
{1.2,0.4204591828942405027},
{1.3,0.445738567528534523},
{1.4,0.4620065850946772763},
{1.5,0.4703563171953998867},
{1.6,0.4717325169318778034},
{1.7,0.4669683641769546303},
{1.8,0.4568111294183368931},
{1.9,0.4419403496815988459},
{2.,0.4229808287748649957},
{2.1,0.4005119878443963905},
{2.2,0.375074599049832154},
{2.3,0.3471756175403162244},
{2.4,0.3172916174366979837},
{2.5,0.2858711963653834954},
{2.6,0.2533366160625841922},
{2.7,0.2200848786329616183},
{2.8,0.1864883896431757677},
{2.9,0.1528953241595883109},
{3.,0.1196297860080003276},
{3.1,0.08699183119553698867},
{3.2,0.05525741171994249172},
{3.3,0.02467828460795811365},
{3.4,-0.004518077930741953583},
{3.5,-0.03212854851248111562},
{3.6,-0.05797435185980087899},
{3.7,-0.08190100128429845044},
{3.8,-0.103778150356897706},
{3.9,-0.1234993492078151427},
{4.,-0.1409816978869304116},
{4.1,-0.1561653918281211096},
{4.2,-0.1690131567671567387},
{4.3,-0.1795095725126332535},
{4.4,-0.1876602868004406848},
{4.5,-0.1934911221017387574},
{4.6,-0.1970470797223561954},
{4.7,-0.1983912468424728519},
{4.8,-0.1976036133099352364},
{4.9,-0.1947798060262372251},
{5.,-0.1900297496566438786},
{5.1,-0.1834762631592988937},
{5.2,-0.1752536022656594845},
{5.3,-0.1655059585589272527},
{5.4,-0.1543859261907244327},
{5.5,-0.1420529475515192553},
{5.6,-0.1286717493698078152},
{5.7,-0.1144107807616790621},
{5.8,-0.09944066468937858453},
{5.9,-0.08393267411855649376},
{6.,-0.0680572438932471262},
{6.1,-0.05198252898002196938},
{6.2,-0.03587301927345499141},
{6.3,-0.01988822060984221721},
{6.4,-0.004181411011335062901},
{6.5,0.01110151951493010868},
{6.6,0.0258231380612632636},
{6.7,0.03985544004704345647},
{6.8,0.05308071672019919828},
{6.9,0.06539231397595150445},
{7.,0.07669527848218451838},
{7.1,0.08690688807134783966},
{7.2,0.09595706434518079339},
{7.3,0.1037886664320276284},
{7.4,0.1103576658283782431},
{7.5,0.1156332032379342704},
{7.6,0.1195975292845658358},
{7.7,0.1222458319118463165},
{7.8,0.1235859541836085566},
{7.9,0.1236380070597178432},
{8.,0.1224338825320095573},
{8.1,0.1200166732605965689},
{8.2,0.1164400055445665387},
{8.3,0.1117672930881112503},
{8.4,0.1060709195786439133},
{8.5,0.09943135857342191604},
{8.6,0.09193623959225739866},
{8.7,0.08367936963344415541},
{8.8,0.07475971956617320393},
{8.9,0.06528038500431541014},
{9.,0.05534753133313360709},
{9.1,0.04506933254261214352},
{9.2,0.03455491341975987981},
{9.3,0.02391330446927571965},
{9.4,0.01325241866988490234},
{9.5,0.002678058835650657419},
{9.6,-0.007707036058533521531},
{9.7,-0.01780409770583747337},
{9.8,-0.02751918110980938821},
{9.9,-0.03676395629683637606},
{10.,-0.04545643300445537263}};

// input for wolframalpha.com: Table[{N[x/10], N[SinIntegral[x/10], 18]}, {x, -100, 100}]
double sivalues[201][2]=
{{-10.,-1.658347594218874049},
{-9.9,-1.663384056595864698},
{-9.8,-1.667569616851386723},
{-9.7,-1.67084456972736341},
{-9.6,-1.673156980105444773},
{-9.5,-1.674463342281433091},
{-9.4,-1.674729172532594998},
{-9.3,-1.673929528316134798},
{-9.2,-1.672049447994015827},
{-9.1,-1.669084305598515127},
{-9.,-1.665040075829602495},
{-8.9,-1.659933505204107197},
{-8.8,-1.653792186051813822},
{-8.7,-1.646654530868780776},
{-8.6,-1.638569645386548639},
{-8.5,-1.629597099590385592},
{-8.4,-1.619806596812884436},
{-8.3,-1.609277541933421013},
{-8.2,-1.59809851062137326},
{-8.1,-1.586366622463643107},
{-8.,-1.574186821706942052},
{-7.9,-1.561671070214550172},
{-7.8,-1.548937458077995749},
{-7.7,-1.53610923812865961},
{-7.6,-1.523313791355258246},
{-7.5,-1.510681530943385878},
{-7.4,-1.498344753306055538},
{-7.3,-1.486436445063168036},
{-7.2,-1.475089055447246115},
{-7.1,-1.46443324405734104},
{-7.,-1.454596614248093591},
{-6.9,-1.445702442722501157},
{-6.8,-1.437868416091684784},
{-6.7,-1.43120538527026366},
{-6.6,-1.425816148589978626},
{-6.5,-1.421794274435881687},
{-6.4,-1.419222974038433272},
{-6.3,-1.418174034791727003},
{-6.2,-1.418706824114094167},
{-6.1,-1.420867373424620212},
{-6.,-1.424687551280506536},
{-5.9,-1.430184334109366196},
{-5.8,-1.437359182281820432},
{-5.7,-1.446197528508234603},
{-5.6,-1.456668384714831391},
{-5.5,-1.468724072665098669},
{-5.4,-1.482300082649290016},
{-5.3,-1.497315063575331138},
{-5.2,-1.513670946766480634},
{-5.1,-1.531253204712921712},
{-5.,-1.549931244944674137},
{-4.9,-1.569558938100652056},
{-4.8,-1.589975278172365569},
{-4.7,-1.611005171809781381},
{-4.6,-1.63246035250034989},
{-4.5,-1.654140414379243984},
{-4.4,-1.675833959408374161},
{-4.3,-1.697319850682468461},
{-4.2,-1.718368563690868599},
{-4.1,-1.738743626491768926},
{-4.,-1.758203138949053058},
{-3.9,-1.776501360447805439},
{-3.8,-1.793390354849570171},
{-3.7,-1.80862168087845377},
{-3.6,-1.821948115649503541},
{-3.5,-1.833125398665997048},
{-3.4,-1.841913983326143036},
{-3.3,-1.848080782795211419},
{-3.2,-1.851400897018440279},
{-3.1,-1.851659307674519836},
{-3.,-1.848652527999468256},
{-2.9,-1.842190194645858628},
{-2.8,-1.832096589081322327},
{-2.7,-1.818212076470866054},
{-2.6,-1.80039445052677016},
{-2.5,-1.778520173443826642},
{-2.4,-1.752485500761767386},
{-2.3,-1.722207481805503392},
{-2.2,-1.68762482724109852},
{-2.1,-1.64869863624441878},
{-2.,-1.605412976802694849},
{-1.9,-1.55777531374881852},
{-1.8,-1.505816780255578572},
{-1.7,-1.449592289683321133},
{-1.6,-1.389180485870438428},
{-1.5,-1.32468353117211968},
{-1.4,-1.256226732779217943},
{-1.3,-1.183958009076062919},
{-1.2,-1.10804719901371859},
{-1.1,-1.028685218673733522},
{-1.,-0.9460830703671830149},
{-0.9,-0.8604707107452929328},
{-0.8,-0.7720957854819965603},
{-0.7,-0.6812222391166113109},
{-0.6,-0.5881288096080800669},
{-0.5,-0.4931074180430666892},
{-0.4,-0.396461464751372883},
{-0.3,-0.2985040438070431614},
{-0.2,-0.1995560885262338214},
{-0.1,-0.09994446110827695016},
{0,0},
{0.1,0.09994446110827695016},
{0.2,0.1995560885262338214},
{0.3,0.2985040438070431614},
{0.4,0.396461464751372883},
{0.5,0.4931074180430666892},
{0.6,0.5881288096080800669},
{0.7,0.6812222391166113109},
{0.8,0.7720957854819965603},
{0.9,0.8604707107452929328},
{1.,0.9460830703671830149},
{1.1,1.028685218673733522},
{1.2,1.10804719901371859},
{1.3,1.183958009076062919},
{1.4,1.256226732779217943},
{1.5,1.32468353117211968},
{1.6,1.389180485870438428},
{1.7,1.449592289683321133},
{1.8,1.505816780255578572},
{1.9,1.55777531374881852},
{2.,1.605412976802694849},
{2.1,1.64869863624441878},
{2.2,1.68762482724109852},
{2.3,1.722207481805503392},
{2.4,1.752485500761767386},
{2.5,1.778520173443826642},
{2.6,1.80039445052677016},
{2.7,1.818212076470866054},
{2.8,1.832096589081322327},
{2.9,1.842190194645858628},
{3.,1.848652527999468256},
{3.1,1.851659307674519836},
{3.2,1.851400897018440279},
{3.3,1.848080782795211419},
{3.4,1.841913983326143036},
{3.5,1.833125398665997048},
{3.6,1.821948115649503541},
{3.7,1.80862168087845377},
{3.8,1.793390354849570171},
{3.9,1.776501360447805439},
{4.,1.758203138949053058},
{4.1,1.738743626491768926},
{4.2,1.718368563690868599},
{4.3,1.697319850682468461},
{4.4,1.675833959408374161},
{4.5,1.654140414379243984},
{4.6,1.63246035250034989},
{4.7,1.611005171809781381},
{4.8,1.589975278172365569},
{4.9,1.569558938100652056},
{5.,1.549931244944674137},
{5.1,1.531253204712921712},
{5.2,1.513670946766480634},
{5.3,1.497315063575331138},
{5.4,1.482300082649290016},
{5.5,1.468724072665098669},
{5.6,1.456668384714831391},
{5.7,1.446197528508234603},
{5.8,1.437359182281820432},
{5.9,1.430184334109366196},
{6.,1.424687551280506536},
{6.1,1.420867373424620212},
{6.2,1.418706824114094167},
{6.3,1.418174034791727003},
{6.4,1.419222974038433272},
{6.5,1.421794274435881687},
{6.6,1.425816148589978626},
{6.7,1.43120538527026366},
{6.8,1.437868416091684784},
{6.9,1.445702442722501157},
{7.,1.454596614248093591},
{7.1,1.46443324405734104},
{7.2,1.475089055447246115},
{7.3,1.486436445063168036},
{7.4,1.498344753306055538},
{7.5,1.510681530943385878},
{7.6,1.523313791355258246},
{7.7,1.53610923812865961},
{7.8,1.548937458077995749},
{7.9,1.561671070214550172},
{8.,1.574186821706942052},
{8.1,1.586366622463643107},
{8.2,1.59809851062137326},
{8.3,1.609277541933421013},
{8.4,1.619806596812884436},
{8.5,1.629597099590385592},
{8.6,1.638569645386548639},
{8.7,1.646654530868780776},
{8.8,1.653792186051813822},
{8.9,1.659933505204107197},
{9.,1.665040075829602495},
{9.1,1.669084305598515127},
{9.2,1.672049447994015827},
{9.3,1.673929528316134798},
{9.4,1.674729172532594998},
{9.5,1.674463342281433091},
{9.6,1.673156980105444773},
{9.7,1.67084456972736341},
{9.8,1.667569616851386723},
{9.9,1.663384056595864698},
{10.,1.658347594218874049}};

int testSiCi() {
  int fail = 0;
  double xmaxrel = -1, maxrel = 0, xmaxabs=-1, maxabs=0;
  for (int i = 0; i < 201; ++i) {
    double x = sivalues[i][0];
    double SiMathematica = sivalues[i][1];
    double SiRoot = ROOT::Math::sinint(x);
    if (std::abs(SiRoot-SiMathematica) > maxabs) {
      maxabs = std::fabs(SiRoot-SiMathematica);
      xmaxabs = x;
    }
    if (std::abs(SiRoot/SiMathematica-1) > maxrel) {
      maxrel = std::abs(SiRoot/SiMathematica-1);
      xmaxrel = x;
    }
    
//     cout << "x=" << x << ": Si(x) = " << SiMathematica
//          << ", sinint(x) = " << SiRoot
//          << ", rel. diff = " << SiRoot/SiMathematica-1
//          << endl;
  }
  cout << "Si: Maximum relative deviation: " << maxrel
       << " for x=" << xmaxrel 
       << (maxrel > 1E-15 ? " -> FAIL" : " -> pass") << endl;
  cout << "Si: Maximum absolute deviation: " << maxabs
       << " for x=" << xmaxabs 
       << (maxabs > 1E-15 ? " -> FAIL" : " -> pass") << endl;
  if (maxrel > 1E-15) fail += 1;
  if (maxabs > 1E-15) fail += 2;
  
  xmaxrel = -1; maxrel = 0; xmaxabs=-1; maxabs=0;
  for (int i = 0; i < 200; ++i) {
    double x = civalues[i][0];
    double CiMathematica = civalues[i][1];
    double CiRoot = ROOT::Math::cosint(x);
    if (std::fabs(CiRoot-CiMathematica) > maxabs) {
      maxabs = std::fabs(CiRoot-CiMathematica);
      xmaxabs = x;
    }
    if (std::fabs(CiRoot/CiMathematica-1) > maxrel) {
      maxrel = std::fabs(CiRoot/CiMathematica-1);
      xmaxrel = x;
    }
//     cout << "x=" << x << ": Ci(x) = " << CiMathematica
//          << ", cosint(x) = " << CiRoot
//          << ", rel. diff = " << CiRoot/CiMathematica-1
//          << endl;
  }
  cout << "Ci: Maximum relative deviation: " << maxrel
       << " for x=" << xmaxrel
       << (maxrel > 1E-13 ? " -> FAIL" : " -> pass") <<endl;
  cout << "Ci: Maximum absolute deviation: " << maxabs
       << " for x=" << xmaxabs 
       << (maxabs > 1E-15 ? " -> FAIL" : " -> pass")<< endl;
   
  if (maxrel > 1E-13) fail += 4;
  if (maxabs > 1E-15) fail += 8;
  return fail;
}

int testSpecFuncSiCi() { 
   int status = testSiCi(); 
   std::cerr << "Test Special Function Sin integral and Cos integral :\t";
   if (status)  std::cout << "FAILED !!!" << std::endl;
   else 
      std::cerr << "OK " << std::endl;
   return status; 
}

int main() { 
   return testSpecFuncSiCi();
}