utlMath.h

Go to the documentation of this file.

#ifndef __utlMath_h
#define __utlMath_h

#include <tr1/cmath>
#include "utlCore.h"
#include "utlConstants.h"

namespace utl
{

typedef double (*Func1To1)(double); 

const static double BesselJPrimeZerosTable[] = 
  {
   // 10 in order 0
   0.0,3.8317059702075125,7.0155866698156055,10.173468135062722,13.323691936314223,16.47063005087763,19.615858510468243,22.760084380592772,25.903672087618382,29.046828534916855,           
   // 10 in order 1
   1.841183781340659,5.3314427735250325,8.536316366346284,11.706004902592063,14.863588633909034,18.015527862681804,21.16436985918879,24.311326857210776,27.457050571059245,30.601922972669094,
   // 10 in order 2
   3.0542369282271404,6.706133194158457,9.969467823087596,13.170370856016122,16.347522318321783,19.512912782488204,22.671581772477424,25.826037141785264,28.977672772993678,32.127327020443474,
   // 10 in order 3
   4.201188941210528,8.015236598375951,11.345924310742964,14.585848286167028,17.78874786606647,20.9724769365377,24.144897432909264,27.310057930204348,30.470268806290424,33.62694918279668,
   // 10 in order 4
   5.317553126083994,9.282396285241617,12.68190844263889,15.964107037731551,19.196028800048904,22.401032267689004,25.589759681386735,28.767836217666503,31.938539340972785,35.10391667734677,
   // 10 in order 5
   6.41561637570024,10.519860873772291,13.9871886301403,17.312842487884627,20.57551452138689,23.803581476593862,27.01030789777772,30.20284907898166,33.38544390101012,36.56077768688036,
   // 10 in order 6
   7.5012661446841475,11.734935953042752,15.268181461097873,18.637443009666203,21.931715017802237,25.183925599499627,28.409776362510083,31.617875716105036,34.81339298429743,37.9996408977153,
   // 10 in order 7
   8.57783648971416,12.93238623709143,16.529365884366943,19.941853366527344,23.26805292645757,26.545032061823576,29.790748583196613,33.015178641375144,36.22438054878716,39.42227457893926,
   // 10 in order 8
   9.647421651997242,14.11551890789478,17.774012366915255,21.229062622853125,24.58719748631768,27.889269427955092,31.155326556188324,34.396628554272176,37.620078044197086,40.83017868182204,
   // 10 in order 9
   10.711433970699948,15.286737667333524,19.004593537946054,22.501398726777285,25.891277276839137,29.21856349993608,32.50524735237553,35.7637929288088,39.00190281151422,42.22463843075328,
   // 10 in order 10
   11.770876674955586,16.447852748486536,20.223031412681703,23.760715860327448,27.182021527190532,30.534504754007074,33.84196577513571,37.118000423665606,40.37106890533389,43.60676490137951,
  };

const static double BesselJPrimeZerosOrder1[60] = 
  {
 1.841183781340659, 5.331442773525033, 8.536316366346325, 11.706004902592063, 14.863588633909032, 18.015527862681804, 21.164369859188792, 
 24.311326857210776, 27.457050571059245, 30.601922972669094, 33.746182898667385, 36.88998740923681, 40.03344405335068, 43.17662896544882, 
 46.319597561173914, 49.46239113970275, 52.60504111155669, 55.74757179225101, 58.8900022991857, 62.03234787066199, 65.17462080254445, 
 68.31683112595181, 71.458987105851, 74.60109561345641, 77.74316240819677, 80.88519235387844, 84.02718958629353, 87.16915764454028, 
 90.31109957490342, 93.45301801376003, 96.59491525429114, 99.73679330057391, 102.87865391175445, 106.0204986383608, 109.16232885234086, 
 112.30414577205505, 115.44595048318557, 118.58774395631993, 121.7295270618102, 124.87130058238789, 128.01306522392028, 131.15482162462013, 
 134.29657036296211, 137.43831196451296, 140.58004690784534, 143.72177562967596, 146.86349852934376, 150.0052159727252, 153.14692829566764, 
 156.28863580700812, 159.43033879123553, 162.57203751084367, 165.71373220841684, 168.85542310848228, 171.99711041915972, 175.13879433363329, 
 178.28047503146777, 181.42215267978807, 184.56382743433812, 187.70549944043358
  };


const static double GammaHalfIntegerTable[30] =
  {
1.7724538509055160272981674833411452, 
0.88622692545275801364908374167057259, 
1.3293403881791370204736256125058589, 
3.3233509704478425511840640312646472, 
11.631728396567448929144224109426265, 
52.342777784553520181149008492418194, 
287.88527781504436099631954670830007, 
1871.2543057977883464760770536039504, 
14034.407293483412598570577902029628, 
119292.46199460900708784991216725184, 
1.1332783889487855673345741655888925*1e6, 
1.1899423083962248457013028738683371*1e7, 
1.3684336546556585725564983049485877*1e8, 
1.7105420683195732156956228811857346*1e9, 
2.3092317922314238411890908896007417*1e10, 
3.3483860987355645697241817899210754*1e11, 
5.1899984530401250830724817743776669*1e12, 
8.5634974475162063870695949277231504*1e13, 
1.4986120533153361177371791123515513*1e15, 
2.7724322986333718178137813578503700*1e16, 
5.4062429823350750447368736478082214*1e17, 
1.1082798113786903841710590978006854*1e19, 
2.3828015944641843259677770602714736*1e20, 
5.3613035875444147334274983856108156*1e21, 
1.2599063430729374623554621206185417*1e23, 
3.0867705405286967827708821955154271*1e24, 
7.8712648783481767960657495985643390*1e25, 
2.0858851927622668509574236436195498*1e27, 
5.7361842800962338401329150199537621*1e28, 
1.6348125198274266444378807806868222*1e30
  };

  
const static unsigned long FactorialTable[21]=
  {
  1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 
  39916800, 479001600, 6227020800, 87178291200, 1307674368000, 20922789888000, 
  355687428096000, 6402373705728000, 121645100408832000,2432902008176640000
  // , 
  // 51090942171709440000, 1124000727777607680000, 
  // 25852016738884976640000, 620448401733239439360000, 
  // 15511210043330985984000000, 403291461126605635584000000, 
  // 10888869450418352160768000000, 304888344611713860501504000000, 
  // 8841761993739701954543616000000, 265252859812191058636308480000000
  };

inline double
ExpNegtiveLUT(const double dist, const double distMax=30.0, const int precision=1000 )
{
  if (dist >= distMax-1) 
    return 0.0;

  static int LUT_LENGTH = (int) rintf( distMax * (double) precision);
  static std::vector<double> EXP_LUT(LUT_LENGTH,-1.0);
  
  static bool is_firstTime = true;
  if (is_firstTime)
    {
    for (int i=0; i< LUT_LENGTH; i++)   
      EXP_LUT[i] = std::exp( - (double) i / (double)precision);
    is_firstTime = false;
    }

  double distPrecision = dist* (double)precision;
  int x = (int) std::floor(distPrecision);
  return EXP_LUT[x] + (EXP_LUT[x+1]-EXP_LUT[x])*(distPrecision-x);
}

inline double
PowInteger(const double a, const int b)
{
  if (b==0)
    return 1;

  if (b<0)
    return 1.0/PowInteger(a, -b);

  double result = a;
  for ( int i = 1; i < b; i += 1 ) 
    result *= a;
  return result;
}

inline double
PowHalfInteger(const double a, const double b)
{

  if (b==0)
    return 1;

  if (b<0)
    return 1.0/PowHalfInteger(a, -b);
  
  if (utl::IsInt(b))
    return PowInteger(a, (int)b);

  utlException(!utl::IsInt(2*b), "b need to be an integer or a half integer.");

  int b_int = int(b-0.5);
  return std::sqrt(a)*PowInteger(a, b_int);
}

inline unsigned long 
Factorial(const int n)
{
  utlException(n<0, "n should be non-negative");
  if(n <= 1) // n=0,1 return 1
    return 1;

  utlException(n>20, "n is too big for the table");
  
  return FactorialTable[n];
}

template <typename T>
inline T 
Factorial ( const T v1, const int times )
{
  T result = v1;
  for ( int i = 1; i < times; i += 1 ) 
    {
    result *= v1 - i;
    }
  return result;
}

template <typename T>
inline T 
Binomial ( const T v1, const int times )
{
  utlException(times<0, "negative value is invalid");
  if (times==0)
    return 1;
  T fac_1 = Factorial(v1,times);
  unsigned long fac_2 = Factorial(times);
  return fac_1 / T(fac_2);
}

/* Compute a scaled Dawson integral 
            FADDEEVA(w_im)(x) = 2*Dawson(x)/sqrt(pi)
   equivalent to the imaginary part w(x) for real x.

   Uses methods similar to the erfcx calculation above: continued fractions
   for large |x|, a lookup table of Chebyshev polynomials for smaller |x|,
   and finally a Taylor expansion for |x|<0.01.
   
   Steven G. Johnson, October 2012. */

/* Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x).

   Uses a look-up table of 100 different Chebyshev polynomials
   for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
   with the help of Maple and a little shell script.   This allows
   the Chebyshev polynomials to be of significantly lower degree (about 1/30)
   compared to fitting the whole [0,1] interval with a single polynomial. */
static double w_im_y100(double y100, double x) {
  switch ((int) y100) {
    case 0: {
      double t = 2*y100 - 1;
      return 0.28351593328822191546e-2 + (0.28494783221378400759e-2 + (0.14427470563276734183e-4 + (0.10939723080231588129e-6 + (0.92474307943275042045e-9 + (0.89128907666450075245e-11 + 0.92974121935111111110e-13 * t) * t) * t) * t) * t) * t;
    }
    case 1: {
      double t = 2*y100 - 3;
      return 0.85927161243940350562e-2 + (0.29085312941641339862e-2 + (0.15106783707725582090e-4 + (0.11716709978531327367e-6 + (0.10197387816021040024e-8 + (0.10122678863073360769e-10 + 0.10917479678400000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 2: {
      double t = 2*y100 - 5;
      return 0.14471159831187703054e-1 + (0.29703978970263836210e-2 + (0.15835096760173030976e-4 + (0.12574803383199211596e-6 + (0.11278672159518415848e-8 + (0.11547462300333495797e-10 + 0.12894535335111111111e-12 * t) * t) * t) * t) * t) * t;
    }
    case 3: {
      double t = 2*y100 - 7;
      return 0.20476320420324610618e-1 + (0.30352843012898665856e-2 + (0.16617609387003727409e-4 + (0.13525429711163116103e-6 + (0.12515095552507169013e-8 + (0.13235687543603382345e-10 + 0.15326595042666666667e-12 * t) * t) * t) * t) * t) * t;
    }
    case 4: {
      double t = 2*y100 - 9;
      return 0.26614461952489004566e-1 + (0.31034189276234947088e-2 + (0.17460268109986214274e-4 + (0.14582130824485709573e-6 + (0.13935959083809746345e-8 + (0.15249438072998932900e-10 + 0.18344741882133333333e-12 * t) * t) * t) * t) * t) * t;
    }
    case 5: {
      double t = 2*y100 - 11;
      return 0.32892330248093586215e-1 + (0.31750557067975068584e-2 + (0.18369907582308672632e-4 + (0.15761063702089457882e-6 + (0.15577638230480894382e-8 + (0.17663868462699097951e-10 + (0.22126732680711111111e-12 + 0.30273474177737853668e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 6: {
      double t = 2*y100 - 13;
      return 0.39317207681134336024e-1 + (0.32504779701937539333e-2 + (0.19354426046513400534e-4 + (0.17081646971321290539e-6 + (0.17485733959327106250e-8 + (0.20593687304921961410e-10 + (0.26917401949155555556e-12 + 0.38562123837725712270e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 7: {
      double t = 2*y100 - 15;
      return 0.45896976511367738235e-1 + (0.33300031273110976165e-2 + (0.20423005398039037313e-4 + (0.18567412470376467303e-6 + (0.19718038363586588213e-8 + (0.24175006536781219807e-10 + (0.33059982791466666666e-12 + 0.49756574284439426165e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 8: {
      double t = 2*y100 - 17;
      return 0.52640192524848962855e-1 + (0.34139883358846720806e-2 + (0.21586390240603337337e-4 + (0.20247136501568904646e-6 + (0.22348696948197102935e-8 + (0.28597516301950162548e-10 + (0.41045502119111111110e-12 + 0.65151614515238361946e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 9: {
      double t = 2*y100 - 19;
      return 0.59556171228656770456e-1 + (0.35028374386648914444e-2 + (0.22857246150998562824e-4 + (0.22156372146525190679e-6 + (0.25474171590893813583e-8 + (0.34122390890697400584e-10 + (0.51593189879111111110e-12 + 0.86775076853908006938e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 10: {
      double t = 2*y100 - 21;
      return 0.66655089485108212551e-1 + (0.35970095381271285568e-2 + (0.24250626164318672928e-4 + (0.24339561521785040536e-6 + (0.29221990406518411415e-8 + (0.41117013527967776467e-10 + (0.65786450716444444445e-12 + 0.11791885745450623331e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 11: {
      double t = 2*y100 - 23;
      return 0.73948106345519174661e-1 + (0.36970297216569341748e-2 + (0.25784588137312868792e-4 + (0.26853012002366752770e-6 + (0.33763958861206729592e-8 + (0.50111549981376976397e-10 + (0.85313857496888888890e-12 + 0.16417079927706899860e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 12: {
      double t = 2*y100 - 25;
      return 0.81447508065002963203e-1 + (0.38035026606492705117e-2 + (0.27481027572231851896e-4 + (0.29769200731832331364e-6 + (0.39336816287457655076e-8 + (0.61895471132038157624e-10 + (0.11292303213511111111e-11 + 0.23558532213703884304e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 13: {
      double t = 2*y100 - 27;
      return 0.89166884027582716628e-1 + (0.39171301322438946014e-2 + (0.29366827260422311668e-4 + (0.33183204390350724895e-6 + (0.46276006281647330524e-8 + (0.77692631378169813324e-10 + (0.15335153258844444444e-11 + 0.35183103415916026911e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 14: {
      double t = 2*y100 - 29;
      return 0.97121342888032322019e-1 + (0.40387340353207909514e-2 + (0.31475490395950776930e-4 + (0.37222714227125135042e-6 + (0.55074373178613809996e-8 + (0.99509175283990337944e-10 + (0.21552645758222222222e-11 + 0.55728651431872687605e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 15: {
      double t = 2*y100 - 31;
      return 0.10532778218603311137e0 + (0.41692873614065380607e-2 + (0.33849549774889456984e-4 + (0.42064596193692630143e-6 + (0.66494579697622432987e-8 + (0.13094103581931802337e-9 + (0.31896187409777777778e-11 + 0.97271974184476560742e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 16: {
      double t = 2*y100 - 33;
      return 0.11380523107427108222e0 + (0.43099572287871821013e-2 + (0.36544324341565929930e-4 + (0.47965044028581857764e-6 + (0.81819034238463698796e-8 + (0.17934133239549647357e-9 + (0.50956666166186293627e-11 + (0.18850487318190638010e-12 + 0.79697813173519853340e-14 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 17: {
      double t = 2*y100 - 35;
      return 0.12257529703447467345e0 + (0.44621675710026986366e-2 + (0.39634304721292440285e-4 + (0.55321553769873381819e-6 + (0.10343619428848520870e-7 + (0.26033830170470368088e-9 + (0.87743837749108025357e-11 + (0.34427092430230063401e-12 + 0.10205506615709843189e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 18: {
      double t = 2*y100 - 37;
      return 0.13166276955656699478e0 + (0.46276970481783001803e-2 + (0.43225026380496399310e-4 + (0.64799164020016902656e-6 + (0.13580082794704641782e-7 + (0.39839800853954313927e-9 + (0.14431142411840000000e-10 + 0.42193457308830027541e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 19: {
      double t = 2*y100 - 39;
      return 0.14109647869803356475e0 + (0.48088424418545347758e-2 + (0.47474504753352150205e-4 + (0.77509866468724360352e-6 + (0.18536851570794291724e-7 + (0.60146623257887570439e-9 + (0.18533978397305276318e-10 + (0.41033845938901048380e-13 - 0.46160680279304825485e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 20: {
      double t = 2*y100 - 41;
      return 0.15091057940548936603e0 + (0.50086864672004685703e-2 + (0.52622482832192230762e-4 + (0.95034664722040355212e-6 + (0.25614261331144718769e-7 + (0.80183196716888606252e-9 + (0.12282524750534352272e-10 + (-0.10531774117332273617e-11 - 0.86157181395039646412e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 21: {
      double t = 2*y100 - 43;
      return 0.16114648116017010770e0 + (0.52314661581655369795e-2 + (0.59005534545908331315e-4 + (0.11885518333915387760e-5 + (0.33975801443239949256e-7 + (0.82111547144080388610e-9 + (-0.12357674017312854138e-10 + (-0.24355112256914479176e-11 - 0.75155506863572930844e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 22: {
      double t = 2*y100 - 45;
      return 0.17185551279680451144e0 + (0.54829002967599420860e-2 + (0.67013226658738082118e-4 + (0.14897400671425088807e-5 + (0.40690283917126153701e-7 + (0.44060872913473778318e-9 + (-0.52641873433280000000e-10 - 0.30940587864543343124e-11 * t) * t) * t) * t) * t) * t) * t;
    }
    case 23: {
      double t = 2*y100 - 47;
      return 0.18310194559815257381e0 + (0.57701559375966953174e-2 + (0.76948789401735193483e-4 + (0.18227569842290822512e-5 + (0.41092208344387212276e-7 + (-0.44009499965694442143e-9 + (-0.92195414685628803451e-10 + (-0.22657389705721753299e-11 + 0.10004784908106839254e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 24: {
      double t = 2*y100 - 49;
      return 0.19496527191546630345e0 + (0.61010853144364724856e-2 + (0.88812881056342004864e-4 + (0.21180686746360261031e-5 + (0.30652145555130049203e-7 + (-0.16841328574105890409e-8 + (-0.11008129460612823934e-9 + (-0.12180794204544515779e-12 + 0.15703325634590334097e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 25: {
      double t = 2*y100 - 51;
      return 0.20754006813966575720e0 + (0.64825787724922073908e-2 + (0.10209599627522311893e-3 + (0.22785233392557600468e-5 + (0.73495224449907568402e-8 + (-0.29442705974150112783e-8 + (-0.94082603434315016546e-10 + (0.23609990400179321267e-11 + 0.14141908654269023788e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 26: {
      double t = 2*y100 - 53;
      return 0.22093185554845172146e0 + (0.69182878150187964499e-2 + (0.11568723331156335712e-3 + (0.22060577946323627739e-5 + (-0.26929730679360840096e-7 + (-0.38176506152362058013e-8 + (-0.47399503861054459243e-10 + (0.40953700187172127264e-11 + 0.69157730376118511127e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 27: {
      double t = 2*y100 - 55;
      return 0.23524827304057813918e0 + (0.74063350762008734520e-2 + (0.12796333874615790348e-3 + (0.18327267316171054273e-5 + (-0.66742910737957100098e-7 + (-0.40204740975496797870e-8 + (0.14515984139495745330e-10 + (0.44921608954536047975e-11 - 0.18583341338983776219e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 28: {
      double t = 2*y100 - 57;
      return 0.25058626331812744775e0 + (0.79377285151602061328e-2 + (0.13704268650417478346e-3 + (0.11427511739544695861e-5 + (-0.10485442447768377485e-6 + (-0.34850364756499369763e-8 + (0.72656453829502179208e-10 + (0.36195460197779299406e-11 - 0.84882136022200714710e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 29: {
      double t = 2*y100 - 59;
      return 0.26701724900280689785e0 + (0.84959936119625864274e-2 + (0.14112359443938883232e-3 + (0.17800427288596909634e-6 + (-0.13443492107643109071e-6 + (-0.23512456315677680293e-8 + (0.11245846264695936769e-9 + (0.19850501334649565404e-11 - 0.11284666134635050832e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 30: {
      double t = 2*y100 - 61;
      return 0.28457293586253654144e0 + (0.90581563892650431899e-2 + (0.13880520331140646738e-3 + (-0.97262302362522896157e-6 + (-0.15077100040254187366e-6 + (-0.88574317464577116689e-9 + (0.12760311125637474581e-9 + (0.20155151018282695055e-12 - 0.10514169375181734921e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 31: {
      double t = 2*y100 - 63;
      return 0.30323425595617385705e0 + (0.95968346790597422934e-2 + (0.12931067776725883939e-3 + (-0.21938741702795543986e-5 + (-0.15202888584907373963e-6 + (0.61788350541116331411e-9 + (0.11957835742791248256e-9 + (-0.12598179834007710908e-11 - 0.75151817129574614194e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 32: {
      double t = 2*y100 - 65;
      return 0.32292521181517384379e0 + (0.10082957727001199408e-1 + (0.11257589426154962226e-3 + (-0.33670890319327881129e-5 + (-0.13910529040004008158e-6 + (0.19170714373047512945e-8 + (0.94840222377720494290e-10 + (-0.21650018351795353201e-11 - 0.37875211678024922689e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 33: {
      double t = 2*y100 - 67;
      return 0.34351233557911753862e0 + (0.10488575435572745309e-1 + (0.89209444197248726614e-4 + (-0.43893459576483345364e-5 + (-0.11488595830450424419e-6 + (0.28599494117122464806e-8 + (0.61537542799857777779e-10 - 0.24935749227658002212e-11 * t) * t) * t) * t) * t) * t) * t;
    }
    case 34: {
      double t = 2*y100 - 69;
      return 0.36480946642143669093e0 + (0.10789304203431861366e-1 + (0.60357993745283076834e-4 + (-0.51855862174130669389e-5 + (-0.83291664087289801313e-7 + (0.33898011178582671546e-8 + (0.27082948188277716482e-10 + (-0.23603379397408694974e-11 + 0.19328087692252869842e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 35: {
      double t = 2*y100 - 71;
      return 0.38658679935694939199e0 + (0.10966119158288804999e-1 + (0.27521612041849561426e-4 + (-0.57132774537670953638e-5 + (-0.48404772799207914899e-7 + (0.35268354132474570493e-8 + (-0.32383477652514618094e-11 + (-0.19334202915190442501e-11 + 0.32333189861286460270e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 36: {
      double t = 2*y100 - 73;
      return 0.40858275583808707870e0 + (0.11006378016848466550e-1 + (-0.76396376685213286033e-5 + (-0.59609835484245791439e-5 + (-0.13834610033859313213e-7 + (0.33406952974861448790e-8 + (-0.26474915974296612559e-10 + (-0.13750229270354351983e-11 + 0.36169366979417390637e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 37: {
      double t = 2*y100 - 75;
      return 0.43051714914006682977e0 + (0.10904106549500816155e-1 + (-0.43477527256787216909e-4 + (-0.59429739547798343948e-5 + (0.17639200194091885949e-7 + (0.29235991689639918688e-8 + (-0.41718791216277812879e-10 + (-0.81023337739508049606e-12 + 0.33618915934461994428e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 38: {
      double t = 2*y100 - 77;
      return 0.45210428135559607406e0 + (0.10659670756384400554e-1 + (-0.78488639913256978087e-4 + (-0.56919860886214735936e-5 + (0.44181850467477733407e-7 + (0.23694306174312688151e-8 + (-0.49492621596685443247e-10 + (-0.31827275712126287222e-12 + 0.27494438742721623654e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 39: {
      double t = 2*y100 - 79;
      return 0.47306491195005224077e0 + (0.10279006119745977570e-1 + (-0.11140268171830478306e-3 + (-0.52518035247451432069e-5 + (0.64846898158889479518e-7 + (0.17603624837787337662e-8 + (-0.51129481592926104316e-10 + (0.62674584974141049511e-13 + 0.20055478560829935356e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 40: {
      double t = 2*y100 - 81;
      return 0.49313638965719857647e0 + (0.97725799114772017662e-2 + (-0.14122854267291533334e-3 + (-0.46707252568834951907e-5 + (0.79421347979319449524e-7 + (0.11603027184324708643e-8 + (-0.48269605844397175946e-10 + (0.32477251431748571219e-12 + 0.12831052634143527985e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 41: {
      double t = 2*y100 - 83;
      return 0.51208057433416004042e0 + (0.91542422354009224951e-2 + (-0.16726530230228647275e-3 + (-0.39964621752527649409e-5 + (0.88232252903213171454e-7 + (0.61343113364949928501e-9 + (-0.42516755603130443051e-10 + (0.47910437172240209262e-12 + 0.66784341874437478953e-14 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 42: {
      double t = 2*y100 - 85;
      return 0.52968945458607484524e0 + (0.84400880445116786088e-2 + (-0.18908729783854258774e-3 + (-0.32725905467782951931e-5 + (0.91956190588652090659e-7 + (0.14593989152420122909e-9 + (-0.35239490687644444445e-10 + 0.54613829888448694898e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 43: {
      double t = 2*y100 - 87;
      return 0.54578857454330070965e0 + (0.76474155195880295311e-2 + (-0.20651230590808213884e-3 + (-0.25364339140543131706e-5 + (0.91455367999510681979e-7 + (-0.23061359005297528898e-9 + (-0.27512928625244444444e-10 + 0.54895806008493285579e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 44: {
      double t = 2*y100 - 89;
      return 0.56023851910298493910e0 + (0.67938321739997196804e-2 + (-0.21956066613331411760e-3 + (-0.18181127670443266395e-5 + (0.87650335075416845987e-7 + (-0.51548062050366615977e-9 + (-0.20068462174044444444e-10 + 0.50912654909758187264e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 45: {
      double t = 2*y100 - 91;
      return 0.57293478057455721150e0 + (0.58965321010394044087e-2 + (-0.22841145229276575597e-3 + (-0.11404605562013443659e-5 + (0.81430290992322326296e-7 + (-0.71512447242755357629e-9 + (-0.13372664928000000000e-10 + 0.44461498336689298148e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 46: {
      double t = 2*y100 - 93;
      return 0.58380635448407827360e0 + (0.49717469530842831182e-2 + (-0.23336001540009645365e-3 + (-0.51952064448608850822e-6 + (0.73596577815411080511e-7 + (-0.84020916763091566035e-9 + (-0.76700972702222222221e-11 + 0.36914462807972467044e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 47: {
      double t = 2*y100 - 95;
      return 0.59281340237769489597e0 + (0.40343592069379730568e-2 + (-0.23477963738658326185e-3 + (0.34615944987790224234e-7 + (0.64832803248395814574e-7 + (-0.90329163587627007971e-9 + (-0.30421940400000000000e-11 + 0.29237386653743536669e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 48: {
      double t = 2*y100 - 97;
      return 0.59994428743114271918e0 + (0.30976579788271744329e-2 + (-0.23308875765700082835e-3 + (0.51681681023846925160e-6 + (0.55694594264948268169e-7 + (-0.91719117313243464652e-9 + (0.53982743680000000000e-12 + 0.22050829296187771142e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 49: {
      double t = 2*y100 - 99;
      return 0.60521224471819875444e0 + (0.21732138012345456060e-2 + (-0.22872428969625997456e-3 + (0.92588959922653404233e-6 + (0.46612665806531930684e-7 + (-0.89393722514414153351e-9 + (0.31718550353777777778e-11 + 0.15705458816080549117e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 50: {
      double t = 2*y100 - 101;
      return 0.60865189969791123620e0 + (0.12708480848877451719e-2 + (-0.22212090111534847166e-3 + (0.12636236031532793467e-5 + (0.37904037100232937574e-7 + (-0.84417089968101223519e-9 + (0.49843180828444444445e-11 + 0.10355439441049048273e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 51: {
      double t = 2*y100 - 103;
      return 0.61031580103499200191e0 + (0.39867436055861038223e-3 + (-0.21369573439579869291e-3 + (0.15339402129026183670e-5 + (0.29787479206646594442e-7 + (-0.77687792914228632974e-9 + (0.61192452741333333334e-11 + 0.60216691829459295780e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 52: {
      double t = 2*y100 - 105;
      return 0.61027109047879835868e0 + (-0.43680904508059878254e-3 + (-0.20383783788303894442e-3 + (0.17421743090883439959e-5 + (0.22400425572175715576e-7 + (-0.69934719320045128997e-9 + (0.67152759655111111110e-11 + 0.26419960042578359995e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 53: {
      double t = 2*y100 - 107;
      return 0.60859639489217430521e0 + (-0.12305921390962936873e-2 + (-0.19290150253894682629e-3 + (0.18944904654478310128e-5 + (0.15815530398618149110e-7 + (-0.61726850580964876070e-9 + 0.68987888999111111110e-11 * t) * t) * t) * t) * t) * t;
    }
    case 54: {
      double t = 2*y100 - 109;
      return 0.60537899426486075181e0 + (-0.19790062241395705751e-2 + (-0.18120271393047062253e-3 + (0.19974264162313241405e-5 + (0.10055795094298172492e-7 + (-0.53491997919318263593e-9 + (0.67794550295111111110e-11 - 0.17059208095741511603e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 55: {
      double t = 2*y100 - 111;
      return 0.60071229457904110537e0 + (-0.26795676776166354354e-2 + (-0.16901799553627508781e-3 + (0.20575498324332621581e-5 + (0.51077165074461745053e-8 + (-0.45536079828057221858e-9 + (0.64488005516444444445e-11 - 0.29311677573152766338e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 56: {
      double t = 2*y100 - 113;
      return 0.59469361520112714738e0 + (-0.33308208190600993470e-2 + (-0.15658501295912405679e-3 + (0.20812116912895417272e-5 + (0.93227468760614182021e-9 + (-0.38066673740116080415e-9 + (0.59806790359111111110e-11 - 0.36887077278950440597e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 57: {
      double t = 2*y100 - 115;
      return 0.58742228631775388268e0 + (-0.39321858196059227251e-2 + (-0.14410441141450122535e-3 + (0.20743790018404020716e-5 + (-0.25261903811221913762e-8 + (-0.31212416519526924318e-9 + (0.54328422462222222221e-11 - 0.40864152484979815972e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 58: {
      double t = 2*y100 - 117;
      return 0.57899804200033018447e0 + (-0.44838157005618913447e-2 + (-0.13174245966501437965e-3 + (0.20425306888294362674e-5 + (-0.53330296023875447782e-8 + (-0.25041289435539821014e-9 + (0.48490437205333333334e-11 - 0.42162206939169045177e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 59: {
      double t = 2*y100 - 119;
      return 0.56951968796931245974e0 + (-0.49864649488074868952e-2 + (-0.11963416583477567125e-3 + (0.19906021780991036425e-5 + (-0.75580140299436494248e-8 + (-0.19576060961919820491e-9 + (0.42613011928888888890e-11 - 0.41539443304115604377e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 60: {
      double t = 2*y100 - 121;
      return 0.55908401930063918964e0 + (-0.54413711036826877753e-2 + (-0.10788661102511914628e-3 + (0.19229663322982839331e-5 + (-0.92714731195118129616e-8 + (-0.14807038677197394186e-9 + (0.36920870298666666666e-11 - 0.39603726688419162617e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 61: {
      double t = 2*y100 - 123;
      return 0.54778496152925675315e0 + (-0.58501497933213396670e-2 + (-0.96582314317855227421e-4 + (0.18434405235069270228e-5 + (-0.10541580254317078711e-7 + (-0.10702303407788943498e-9 + (0.31563175582222222222e-11 - 0.36829748079110481422e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 62: {
      double t = 2*y100 - 125;
      return 0.53571290831682823999e0 + (-0.62147030670760791791e-2 + (-0.85782497917111760790e-4 + (0.17553116363443470478e-5 + (-0.11432547349815541084e-7 + (-0.72157091369041330520e-10 + (0.26630811607111111111e-11 - 0.33578660425893164084e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 63: {
      double t = 2*y100 - 127;
      return 0.52295422962048434978e0 + (-0.65371404367776320720e-2 + (-0.75530164941473343780e-4 + (0.16613725797181276790e-5 + (-0.12003521296598910761e-7 + (-0.42929753689181106171e-10 + (0.22170894940444444444e-11 - 0.30117697501065110505e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 64: {
      double t = 2*y100 - 129;
      return 0.50959092577577886140e0 + (-0.68197117603118591766e-2 + (-0.65852936198953623307e-4 + (0.15639654113906716939e-5 + (-0.12308007991056524902e-7 + (-0.18761997536910939570e-10 + (0.18198628922666666667e-11 - 0.26638355362285200932e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 65: {
      double t = 2*y100 - 131;
      return 0.49570040481823167970e0 + (-0.70647509397614398066e-2 + (-0.56765617728962588218e-4 + (0.14650274449141448497e-5 + (-0.12393681471984051132e-7 + (0.92904351801168955424e-12 + (0.14706755960177777778e-11 - 0.23272455351266325318e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 66: {
      double t = 2*y100 - 133;
      return 0.48135536250935238066e0 + (-0.72746293327402359783e-2 + (-0.48272489495730030780e-4 + (0.13661377309113939689e-5 + (-0.12302464447599382189e-7 + (0.16707760028737074907e-10 + (0.11672928324444444444e-11 - 0.20105801424709924499e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 67: {
      double t = 2*y100 - 135;
      return 0.46662374675511439448e0 + (-0.74517177649528487002e-2 + (-0.40369318744279128718e-4 + (0.12685621118898535407e-5 + (-0.12070791463315156250e-7 + (0.29105507892605823871e-10 + (0.90653314645333333334e-12 - 0.17189503312102982646e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 68: {
      double t = 2*y100 - 137;
      return 0.45156879030168268778e0 + (-0.75983560650033817497e-2 + (-0.33045110380705139759e-4 + (0.11732956732035040896e-5 + (-0.11729986947158201869e-7 + (0.38611905704166441308e-10 + (0.68468768305777777779e-12 - 0.14549134330396754575e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 69: {
      double t = 2*y100 - 139;
      return 0.43624909769330896904e0 + (-0.77168291040309554679e-2 + (-0.26283612321339907756e-4 + (0.10811018836893550820e-5 + (-0.11306707563739851552e-7 + (0.45670446788529607380e-10 + (0.49782492549333333334e-12 - 0.12191983967561779442e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 70: {
      double t = 2*y100 - 141;
      return 0.42071877443548481181e0 + (-0.78093484015052730097e-2 + (-0.20064596897224934705e-4 + (0.99254806680671890766e-6 + (-0.10823412088884741451e-7 + (0.50677203326904716247e-10 + (0.34200547594666666666e-12 - 0.10112698698356194618e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 71: {
      double t = 2*y100 - 143;
      return 0.40502758809710844280e0 + (-0.78780384460872937555e-2 + (-0.14364940764532853112e-4 + (0.90803709228265217384e-6 + (-0.10298832847014466907e-7 + (0.53981671221969478551e-10 + (0.21342751381333333333e-12 - 0.82975901848387729274e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 72: {
      double t = 2*y100 - 145;
      return 0.38922115269731446690e0 + (-0.79249269708242064120e-2 + (-0.91595258799106970453e-5 + (0.82783535102217576495e-6 + (-0.97484311059617744437e-8 + (0.55889029041660225629e-10 + (0.10851981336888888889e-12 - 0.67278553237853459757e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 73: {
      double t = 2*y100 - 147;
      return 0.37334112915460307335e0 + (-0.79519385109223148791e-2 + (-0.44219833548840469752e-5 + (0.75209719038240314732e-6 + (-0.91848251458553190451e-8 + (0.56663266668051433844e-10 + (0.23995894257777777778e-13 - 0.53819475285389344313e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 74: {
      double t = 2*y100 - 149;
      return 0.35742543583374223085e0 + (-0.79608906571527956177e-2 + (-0.12530071050975781198e-6 + (0.68088605744900552505e-6 + (-0.86181844090844164075e-8 + (0.56530784203816176153e-10 + (-0.43120012248888888890e-13 - 0.42372603392496813810e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 75: {
      double t = 2*y100 - 151;
      return 0.34150846431979618536e0 + (-0.79534924968773806029e-2 + (0.37576885610891515813e-5 + (0.61419263633090524326e-6 + (-0.80565865409945960125e-8 + (0.55684175248749269411e-10 + (-0.95486860764444444445e-13 - 0.32712946432984510595e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 76: {
      double t = 2*y100 - 153;
      return 0.32562129649136346824e0 + (-0.79313448067948884309e-2 + (0.72539159933545300034e-5 + (0.55195028297415503083e-6 + (-0.75063365335570475258e-8 + (0.54281686749699595941e-10 - 0.13545424295111111111e-12 * t) * t) * t) * t) * t) * t;
    }
    case 77: {
      double t = 2*y100 - 155;
      return 0.30979191977078391864e0 + (-0.78959416264207333695e-2 + (0.10389774377677210794e-4 + (0.49404804463196316464e-6 + (-0.69722488229411164685e-8 + (0.52469254655951393842e-10 - 0.16507860650666666667e-12 * t) * t) * t) * t) * t) * t;
    }
    case 78: {
      double t = 2*y100 - 157;
      return 0.29404543811214459904e0 + (-0.78486728990364155356e-2 + (0.13190885683106990459e-4 + (0.44034158861387909694e-6 + (-0.64578942561562616481e-8 + (0.50354306498006928984e-10 - 0.18614473550222222222e-12 * t) * t) * t) * t) * t) * t;
    }
    case 79: {
      double t = 2*y100 - 159;
      return 0.27840427686253660515e0 + (-0.77908279176252742013e-2 + (0.15681928798708548349e-4 + (0.39066226205099807573e-6 + (-0.59658144820660420814e-8 + (0.48030086420373141763e-10 - 0.20018995173333333333e-12 * t) * t) * t) * t) * t) * t;
    }
    case 80: {
      double t = 2*y100 - 161;
      return 0.26288838011163800908e0 + (-0.77235993576119469018e-2 + (0.17886516796198660969e-4 + (0.34482457073472497720e-6 + (-0.54977066551955420066e-8 + (0.45572749379147269213e-10 - 0.20852924954666666667e-12 * t) * t) * t) * t) * t) * t;
    }
    case 81: {
      double t = 2*y100 - 163;
      return 0.24751539954181029717e0 + (-0.76480877165290370975e-2 + (0.19827114835033977049e-4 + (0.30263228619976332110e-6 + (-0.50545814570120129947e-8 + (0.43043879374212005966e-10 - 0.21228012028444444444e-12 * t) * t) * t) * t) * t) * t;
    }
    case 82: {
      double t = 2*y100 - 165;
      return 0.23230087411688914593e0 + (-0.75653060136384041587e-2 + (0.21524991113020016415e-4 + (0.26388338542539382413e-6 + (-0.46368974069671446622e-8 + (0.40492715758206515307e-10 - 0.21238627815111111111e-12 * t) * t) * t) * t) * t) * t;
    }
    case 83: {
      double t = 2*y100 - 167;
      return 0.21725840021297341931e0 + (-0.74761846305979730439e-2 + (0.23000194404129495243e-4 + (0.22837400135642906796e-6 + (-0.42446743058417541277e-8 + (0.37958104071765923728e-10 - 0.20963978568888888889e-12 * t) * t) * t) * t) * t) * t;
    }
    case 84: {
      double t = 2*y100 - 169;
      return 0.20239979200788191491e0 + (-0.73815761980493466516e-2 + (0.24271552727631854013e-4 + (0.19590154043390012843e-6 + (-0.38775884642456551753e-8 + (0.35470192372162901168e-10 - 0.20470131678222222222e-12 * t) * t) * t) * t) * t) * t;
    }
    case 85: {
      double t = 2*y100 - 171;
      return 0.18773523211558098962e0 + (-0.72822604530339834448e-2 + (0.25356688567841293697e-4 + (0.16626710297744290016e-6 + (-0.35350521468015310830e-8 + (0.33051896213898864306e-10 - 0.19811844544000000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 86: {
      double t = 2*y100 - 173;
      return 0.17327341258479649442e0 + (-0.71789490089142761950e-2 + (0.26272046822383820476e-4 + (0.13927732375657362345e-6 + (-0.32162794266956859603e-8 + (0.30720156036105652035e-10 - 0.19034196304000000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 87: {
      double t = 2*y100 - 175;
      return 0.15902166648328672043e0 + (-0.70722899934245504034e-2 + (0.27032932310132226025e-4 + (0.11474573347816568279e-6 + (-0.29203404091754665063e-8 + (0.28487010262547971859e-10 - 0.18174029063111111111e-12 * t) * t) * t) * t) * t) * t;
    }
    case 88: {
      double t = 2*y100 - 177;
      return 0.14498609036610283865e0 + (-0.69628725220045029273e-2 + (0.27653554229160596221e-4 + (0.92493727167393036470e-7 + (-0.26462055548683583849e-8 + (0.26360506250989943739e-10 - 0.17261211260444444444e-12 * t) * t) * t) * t) * t) * t;
    }
    case 89: {
      double t = 2*y100 - 179;
      return 0.13117165798208050667e0 + (-0.68512309830281084723e-2 + (0.28147075431133863774e-4 + (0.72351212437979583441e-7 + (-0.23927816200314358570e-8 + (0.24345469651209833155e-10 - 0.16319736960000000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 90: {
      double t = 2*y100 - 181;
      return 0.11758232561160626306e0 + (-0.67378491192463392927e-2 + (0.28525664781722907847e-4 + (0.54156999310046790024e-7 + (-0.21589405340123827823e-8 + (0.22444150951727334619e-10 - 0.15368675584000000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 91: {
      double t = 2*y100 - 183;
      return 0.10422112945361673560e0 + (-0.66231638959845581564e-2 + (0.28800551216363918088e-4 + (0.37758983397952149613e-7 + (-0.19435423557038933431e-8 + (0.20656766125421362458e-10 - 0.14422990012444444444e-12 * t) * t) * t) * t) * t) * t;
    }
    case 92: {
      double t = 2*y100 - 185;
      return 0.91090275493541084785e-1 + (-0.65075691516115160062e-2 + (0.28982078385527224867e-4 + (0.23014165807643012781e-7 + (-0.17454532910249875958e-8 + (0.18981946442680092373e-10 - 0.13494234691555555556e-12 * t) * t) * t) * t) * t) * t;
    }
    case 93: {
      double t = 2*y100 - 187;
      return 0.78191222288771379358e-1 + (-0.63914190297303976434e-2 + (0.29079759021299682675e-4 + (0.97885458059415717014e-8 + (-0.15635596116134296819e-8 + (0.17417110744051331974e-10 - 0.12591151763555555556e-12 * t) * t) * t) * t) * t) * t;
    }
    case 94: {
      double t = 2*y100 - 189;
      return 0.65524757106147402224e-1 + (-0.62750311956082444159e-2 + (0.29102328354323449795e-4 + (-0.20430838882727954582e-8 + (-0.13967781903855367270e-8 + (0.15958771833747057569e-10 - 0.11720175765333333333e-12 * t) * t) * t) * t) * t) * t;
    }
    case 95: {
      double t = 2*y100 - 191;
      return 0.53091065838453612773e-1 + (-0.61586898417077043662e-2 + (0.29057796072960100710e-4 + (-0.12597414620517987536e-7 + (-0.12440642607426861943e-8 + (0.14602787128447932137e-10 - 0.10885859114666666667e-12 * t) * t) * t) * t) * t) * t;
    }
    case 96: {
      double t = 2*y100 - 193;
      return 0.40889797115352738582e-1 + (-0.60426484889413678200e-2 + (0.28953496450191694606e-4 + (-0.21982952021823718400e-7 + (-0.11044169117553026211e-8 + (0.13344562332430552171e-10 - 0.10091231402844444444e-12 * t) * t) * t) * t) * t) * t;
    }
  case 97: case 98:
  case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.0309...)
      //  (2/sqrt(pi)) * (x - 2/3 x^3  + 4/15 x^5  - 8/105 x^7 + 16/945 x^9) 
      double x2 = x*x;
      return x * (1.1283791670955125739
                  - x2 * (0.75225277806367504925
                          - x2 * (0.30090111122547001970
                                  - x2 * (0.085971746064420005629
                                          - x2 * 0.016931216931216931217))));
    }
  }
  /* Since 0 <= y100 < 101, this is only reached if x is NaN,
     in which case we should return NaN. */
  return std::numeric_limits<double>::quiet_NaN();
}

inline double 
w_im(double x)
{
  if (x >= 0) {
    if (x > 45) { // continued-fraction expansion is faster
      const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
      if (x > 5e7) // 1-term expansion, important to avoid overflow
        return ispi / x;
      /* 5-term expansion (rely on compiler for CSE), simplified from:
                ispi / (x-0.5/(x-1/(x-1.5/(x-2/x))))  */
      return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
    }
    return w_im_y100(100/(1+x), x);
  }
  else { // = -FADDEEVA(w_im)(-x)
    if (x < -45) { // continued-fraction expansion is faster
      const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
      if (x < -5e7) // 1-term expansion, important to avoid overflow
        return ispi / x;
      /* 5-term expansion (rely on compiler for CSE), simplified from:
                ispi / (x-0.5/(x-1/(x-1.5/(x-2/x))))  */
      return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
    }
    return -w_im_y100(100/(1-x), -x);
  }
}

inline double 
Erfi(double x, Func1To1 expF=&std::exp)
{
  return x*x > 720 ? (x > 0 ? std::numeric_limits<double>::infinity() : -std::numeric_limits<double>::infinity())
    : expF(x*x) * w_im(x);
}

inline double 
Erf(double x)
{
  return std::tr1::erf(x);
}

template <class VectorType>
inline double
Entropy( const VectorType& pdfVec, const int N )
{
  double entropy = 0;
  for ( int j = 0; j < N; ++j ) 
    entropy += -pdfVec[j] * std::log(pdfVec[j]);
  return entropy;
}

inline double
DawsonF(double x)
{
  return utl::w_im_y100(100./(1+x), x)*utl::SQRTPI*0.5;
}

inline double 
LegendrePolynomialAt0(const int order)
{
  if(order == 0)
    return 1.0;
  if(order % 2 != 0)
    return 0.0;
  else 
    {
    double odd = 1;
    double even = 2;

    for(int i = 3; i <= order - 1; i+=2) 
      odd = odd * i;
    for(int i = 4; i <= order; i+=2) 
      even = even * i;

    if((order / 2) % 2 == 0)
      return odd / even;
    else
      return -1 * odd / even;
    }
}

inline std::vector<double>
GetCoefLaguerre (const int n, const double a=0.5)
{
  utlException(n<0, "n should > 0");
  std::vector<double> coef(n+1);
  if (n==0)
    {
    coef[0] = 1;
    }
  else if (n==1)
    {
    coef[0]=1+a; coef[1]=-1.0;
    }
  else if (n==2)
    {
    double a2= a*a;
    coef[0]=1+(3*a)/2.0+a2/2.0; coef[1]=-2-a; coef[2]=1.0/2.0;
    }
  else if (n==3)
    {
    double a2= a*a;
    double a3= a2*a;
    coef[0]=1+(11*a)/6.0+a2+a3/6.0; coef[1]=-3-(5*a)/2.0-a2/2.0; 
    coef[2]=3.0/2.0+a/2.0; 
    coef[3]=-1.0/6.0; 
    }
  else if (n==4)
    {
    double a2= a*a;
    double a3= a2*a;
    double a4= a3*a;
    coef[0]=1+ (25*a)/12.0+(35*a2)/24.0+(5*a3)/12.0+a4/24.0; 
    coef[1]=-4-(13*a)/3.0-(3*a2)/2.0-a3/6.0;
    coef[2]=3+(7*a)/4.0+a2/4.0; 
    coef[3]=-(2.0/3.0)-a/6.0; 
    coef[4]=1/24.0;
    }
  else if (n==5)
    {
    double a2= a*a;
    double a3= a2*a;
    double a4= a3*a;
    double a5= a4*a;
    coef[0]=1+(137*a)/60.0+(15*a2)/8.0+(17*a3)/24.0+a4/8.0+a5/120.0; 
    coef[1]=-5-(77*a)/12.0-(71*a2)/24.0-(7*a3)/12.0-a4/24.0; 
    coef[2]=5+(47*a)/12.0+a2+a3/12.0;
    coef[3]=-(5.0/3.0)-(3*a)/4.0-a2/12.0; 
    coef[4]=5.0/24.0+a/24.0;
    coef[5]=-(1/120.0);
    }
  else if (n==6)
    {
    double a2= a*a;
    double a3= a2*a;
    double a4= a3*a;
    double a5= a4*a;
    double a6= a5*a;
    coef[0]=1+(49*a)/20.0+(203*a2)/90.0+(49*a3)/48.0+(35*a4)/144.0+(7*a5)/240.0+a6/720.0;
    coef[1]=-6-(87*a)/10.0-(29*a2)/6.0-(31*a3)/24.0-a4/6.0-a5/120.0;
    coef[2]=15/2.0+(57*a)/8.0+(119*a2)/48.0+(3*a3)/8.0+a4/48.0;
    coef[3]=-(10/3.0)-(37*a)/18.0-(5*a2)/12.0-a3/36.0;
    coef[4]=5/8.0+(11*a)/48.0+a2/48.0;
    coef[5]=-(1/20.0)-a/120.0;
    coef[6]=1/720.0;
    }
  else if (n==7)
    {
    double a2= a*a;
    double a3= a2*a;
    double a4= a3*a;
    double a5= a4*a;
    double a6= a5*a;
    double a7= a6*a;
    coef[0]=1.+(363.*a)/140.0+(469.*a2)/180.0+(967.*a3)/720.0+(7.*a4)/18.0+(23.*a5)/360.0+a6/180.0+a7/5040.0;
    coef[1]=-7-(223.*a)/20.0-(319*a2)/45.0-(37*a3)/16.0-(59*a4)/144.0-(3*a5)/80.0-a6/720.0;
    coef[2]=21/2.0+(459*a)/40.0+(235*a2)/48.0+(49*a3)/48.0+(5*a4)/48.0+a5/240.0;
    coef[3]=-(35.0/6.0)-(319*a)/72.0-(179*a2)/144.0-(11*a3)/72.0-a4/144.0;
    coef[4]=35.0/24.0+(107*a)/144.0+a2/8.0+a3/144.0;
    coef[5]=-(7.0/40.0)-(13*a)/240.0-a2/240.0;
    coef[6]=7.0/720.0+a/720.0;
    coef[7]=-(1.0/5040);
    }
  else if (n==8)
    {
    double a2= a*a;
    double a3= a2*a;
    double a4= a3*a;
    double a5= a4*a;
    double a6= a5*a;
    double a7= a6*a;
    double a8= a7*a;
    coef[0]=1+(761*a)/280.0+(29531*a2)/10080.0+(267*a3)/160.0+(1069*a4)/1920.0+(9*a5)/80.0+(13*a6)/960.0+a7/1120.0+a8/40320.0;
    coef[1]=-8-(481*a)/35.0-(349*a2)/36.0-(329*a3)/90.0-(115*a4)/144.0-(73*a5)/720.0-a6/144.0-a7/5040.0;
    coef[2]=14+(341*a)/20.0+(6077*a2)/720.0+(209*a3)/96.0+(89*a4)/288.0+(11*a5)/480.0+a6/1440.0;
    coef[3]=-(28/3.0)-(743*a)/90.0-(23*a2)/8.0-(71*a3)/144.0-a4/24.0-a5/720.0;
    coef[4]=35/12.0+(533*a)/288.0+(251*a2)/576.0+(13*a3)/288.0+a4/576.0;
    coef[5]=-(7/15.0)-(73*a)/360.0-(7*a2)/240.0-a3/720.0;
    coef[6]=7/180.0+a/96.0+a2/1440.0;
    coef[7]=-(1/630.0)-a/5040.0;
    coef[8]=1./40320.0;
    }
  else
    {
    std::vector<double> coef_1 = GetCoefLaguerre(n-1,a);
    std::vector<double> coef_2 = GetCoefLaguerre(n-2,a);
    for ( int i = 1; i < n-1; i += 1 ) 
      coef[i] = coef_1[i]* (2.0+(a-1)/double(n)) + coef_1[i-1]*(-1.0/double(n)) - coef_2[i]*(1+(a-1)/n);
    coef[0] = coef_1[0]*(2.0+(a-1)/double(n)) - coef_2[0]*(1+(a-1)/n);
    coef[n-1] = coef_1[n-1]* (2.0+(a-1)/double(n)) + coef_1[n-2]*(-1.0/double(n));
    coef[n] = coef_1[n-1]*(-1.0/double(n));
    }
  return coef;
}    // -----  end of method getCoefLaguerre  -----

inline std::vector<double>
GetCoefLaguerreProduct ( const int n1, const double a1, const int n2, const double a2)
{
  utlException(n1<0 || n2<0, "n1 and n2 should > 0");
  std::vector<double> coef(n1+n2+1,0);
  std::vector<double> coef1 = GetCoefLaguerre(n1,a1);
  std::vector<double> coef2 = GetCoefLaguerre(n2,a2);
  for ( int nn = 0; nn < n1+n2+1; nn += 1 ) 
    {
    for ( int j = 0; j<n1+1 && j<=nn; j += 1 ) 
      {
      int i = nn-j;
      if (i<n2+1)
        coef[nn] += coef1[j]*coef2[i];
      }
    }
  return coef;
}    // -----  end of method getCoefLaguerreProduct  -----

inline double 
GammaHalfInteger(const double x)
{
  utlException(x<=0, "x should be half of a positive integer");
  if (utl::IsInt(x))
    return Factorial(x-1);

  utlException(!utl::IsInt(2*x), "x should be half of a positive integer");
  
  const static int TABLE_LENGTH = 30; 
  
  if (x>TABLE_LENGTH)
    {
    double result=1, tmp=x-1;
    for ( long x_floor = std::floor(x); x_floor>=TABLE_LENGTH; x_floor-- ) 
      {
      result *= tmp;
      tmp -= 1;
      }
    return result*GammaHalfIntegerTable[TABLE_LENGTH-1];
    }
  else 
    {
    int x_floor = std::floor(x);
    return GammaHalfIntegerTable[x_floor];
    }
}

inline double
GetExpLegendreCoef(const double a, const int l, Func1To1 expF=&std::exp )
{
  // coefficient for odd order is 0;
  if (!utl::IsInt(0.5*l))
    return 0;

  if (std::fabs(a)<1e-12)
    return l==0 ? 1.0 : 0.0;
  
  double b = std::fabs(a);

  double result = -1000.0;
  switch ( l )
    {
    case 0 :
      if (a>0)
       result = (utl::SQRTPI*std::tr1::erf(std::sqrt(a)))/(2.*std::sqrt(a));
      else
       result = (utl::SQRTPI*utl::Erfi(std::sqrt(b)))/(2.*std::sqrt(b)); 
       break;
    case 2 :
      if (a>0)
       result = (5*(-3/(2.*a*expF(a)) + ((3 - 2*a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a)))/(4.*std::pow(a,1.5))))/2. ;
      else
       result = (5*(6*std::sqrt(b)*expF(b) - (3 + 2*b)*utl::SQRTPI*Erfi(std::sqrt(b))))/(8.*std::pow(b,1.5)); 
       break;
    case 4 :
      if (a>0)
       result = (9*((-5*(21 + 2*a))/(16.*std::pow(a,2)*expF(a)) + (3*(35 + 4*(-5 + a)*a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a)))/(32.*std::pow(a,2.5))))/2. ;
      else
       result = (9*(10*std::sqrt(b)*(-21 + 2*b)*expF(b) + 3*(35 + 20*b + 4*std::pow(b,2))*utl::SQRTPI*Erfi(std::sqrt(b))))/(64.*std::pow(b,2.5)); 
       break;
    case 6 : 
      if (a>0)
       result = (13*(-42*std::sqrt(a)*(165 + 20*a + 4*std::pow(a,2)) - 5*(-693 + 378*a - 84*std::pow(a,2) + 8*std::pow(a,3))*expF(a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a))))/(256.*std::pow(a,3.5)*expF(a));
      else
       result = (13*(42*std::sqrt(b)*(165 - 20*b + 4*std::pow(b,2))*expF(b) - 5*(693 + 378*b + 84*std::pow(b,2) + 8*std::pow(b,3))*utl::SQRTPI*Erfi(std::sqrt(b))))/(256.*std::pow(b,3.5));
       break;
    case 8 :
      if (a>0)
       result = (17*(-6*std::sqrt(a)*(225225 + 2*a*(15015 + 2*a*(1925 + 62*a))) + 35*(19305 + 8*a*(-1287 + a*(297 + 2*(-18 + a)*a)))*expF(a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a))))/(4096.*std::pow(a,4.5)*expF(a)) ;
      else
       result = (17*(6*std::sqrt(b)*(-225225 + 30030*b - 7700*std::pow(b,2) + 248*std::pow(b,3))*expF(b) + 35*(19305 + 10296*b + 2376*std::pow(b,2) + 288*std::pow(b,3) + 16*std::pow(b,4))*utl::SQRTPI* utl::Erfi(std::sqrt(b)))) / (4096.*std::pow(b,4.5)); 
       break;
    case 10 : 
      if (a>0)
       result = (21*(-22*std::sqrt(a)*(3968055 + 556920*a + 157248*std::pow(a,2) + 7488*std::pow(a,3) + 464*std::pow(a,4)) - 63*(-692835 + 364650*a - 85800*std::pow(a,2) + 11440*std::pow(a,3) - 880*std::pow(a,4) + 32*std::pow(a,5))*expF(a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a))))/(16384.*std::pow(a,5.5)*expF(a)) ; 
      else
       result = (21*(22*std::sqrt(b)*(3968055 - 556920*b + 157248*std::pow(b,2) - 7488*std::pow(b,3) + 464*std::pow(b,4))*expF(b) - 63*(692835 + 364650*b + 85800*std::pow(b,2) + 11440*std::pow(b,3) + 880*std::pow(b,4) + 32*std::pow(b,5))*utl::SQRTPI*Erfi(std::sqrt(b))))/(16384.*std::pow(b,5.5)); 
      break;
    case 12 : 
      if (a>0)
      result = (25*(-26*std::sqrt(a)*(540571185 + 2*a*(39171825 + 4*a*(2909907 + 2*a*(82467 + a*(7469 + 122*a))))) + 231*(30421755 + 4*a*(-3968055 + a*(944775 + 4*a*(-33150 + a*(2925 + 4*(-39 + a)*a)))))*expF(a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a))))/(131072.*std::pow(a,6.5)*expF(a)) ; 
      else
      result = (25*(26*std::sqrt(b)*(-540571185 + 78343650*b - 23279256*std::pow(b,2) + 1319472*std::pow(b,3) - 119504*std::pow(b,4) + 1952*std::pow(b,5))*expF(b) + 231*(30421755 + 15872220*b + 3779100*std::pow(b,2) + 530400*std::pow(b,3) + 46800*std::pow(b,4) + 2496*std::pow(b,5) + 64*std::pow(b,6))*utl::SQRTPI*Erfi(std::sqrt(b))))/(131072.*std::pow(b,6.5)); 
      break;
    case 14 : 
      if (a>0)
      result = (29*(-2*std::sqrt(a)*(677644592625 + 100391791500*a + 30786816060*std::pow(a,2) + 1928852640*std::pow(a,3) + 206187696*std::pow(a,4) + 5360576*std::pow(a,5) + 158528*std::pow(a,6)) - 429*(-1579591125 + 819047250*a - 196571340*std::pow(a,2) + 28488600*std::pow(a,3) - 2713200*std::pow(a,4) + 171360*std::pow(a,5) - 6720*std::pow(a,6) + 128*std::pow(a,7))*expF(a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a))))/(524288.*std::pow(a,7.5)*expF(a)) ; 
      else
      result = (29*(2*std::sqrt(b)*(677644592625 - 100391791500*b + 30786816060*std::pow(b,2) - 1928852640*std::pow(b,3) + 206187696*std::pow(b,4) - 5360576*std::pow(b,5) + 158528*std::pow(b,6))*expF(b) - 429*(1579591125 + 819047250*b + 196571340*std::pow(b,2) + 28488600*std::pow(b,3) + 2713200*std::pow(b,4) + 171360*std::pow(b,5) + 6720*std::pow(b,6) + 128*std::pow(b,7))*utl::SQRTPI*Erfi(std::sqrt(b))))/(524288.*std::pow(b,7.5));
      break;
    case 16 :
      if (a>0)
      result = (33*(-34*std::sqrt(a)*(3583544051587.5e1 + 2*a*(269729122162.5e1 + 2*a*(42253133422.5e1 + 2*a*(1413077737.5e1 + 2.0*a*(82956802.5e1 + 2*a*(1343803.5e1 + 62415.0e1*a + 6196*std::pow(a,2))))))) + 6435*(94670161425 + 16*a*(-305387617.5e1 + a*(73714252.5e1 + 2*a*(-5460315.0e1 + a*(5460315 + 8*a*(-47481 + a*(2261 + (-68 + a)*a)))))))*expF(a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a))))/(1.6777216e7*std::pow(a,8.5)*expF(a)) ; 
      else
      result = (33*(34*std::sqrt(b)*(-35835440515875 + 5394582443250*b - 1690125336900*std::pow(b,2) + 113046219000*std::pow(b,3) - 13273088400*std::pow(b,4) + 430017120*std::pow(b,5) - 19972800*std::pow(b,6) + 198272*std::pow(b,7))*expF(b) + 6435*(94670161425 + 48862018800*b + 11794280400*std::pow(b,2) + 1747300800*std::pow(b,3) + 174730080*std::pow(b,4) + 12155136*std::pow(b,5) + 578816*std::pow(b,6) + 17408*std::pow(b,7) + 256*std::pow(b,8))*utl::SQRTPI*Erfi(std::sqrt(b))))/(1.6777216e7*std::pow(b,8.5));
      break;
    case 18 : 
      if (a>0)
      result = (3.7e1*(-11.4e1*std::sqrt(a)*(137159624664562.5e1 + 8*a*(2612564279325.0e1 + a*(831270452512.5e1 + 4*a*(14556809767.5e1 + 2*a*(907887337.5e1 + 2*a*(16773900.0e1 + a*(961273.5e1 + 2*a*(7939.0e1 + 136.3e1*a)))))))) - 1215.5e1*(-643200214387.5e1 + 2*a*(165394340842.5e1 + 8*a*(-5011949722.5e1 + 2*a*(377243527.5e1 + a*(-39025192.5e1 + 2*a*(1445377.5e1 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a)))))))))*expF(a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a))))/(6.7108864e7*std::pow(a,9.5)*expF(a)) ; 
      else
      result = (37*(114*std::sqrt(b)*(1371596246645625 - 209005142346000*b + 66501636201000*std::pow(b,2) - 4658179125600*std::pow(b,3) + 581047896000*std::pow(b,4) - 21470592000*std::pow(b,5) + 1230430080*std::pow(b,6) - 20323840*std::pow(b,7) + 348928*std::pow(b,8))*expF(b) - 12155*(6432002143875 + 3307886816850*b + 801911955600*std::pow(b,2) + 120717928800*std::pow(b,3) + 12488061600*std::pow(b,4) + 925041600*std::pow(b,5) + 49335552*std::pow(b,6) + 1838592*std::pow(b,7) + 43776*std::pow(b,8) + 512*std::pow(b,9))*utl::SQRTPI*Erfi(std::sqrt(b))))/(6.7108864e7*std::pow(b,9.5)); 
      break;
    case 20 : 
      if (a>0)
      result = (41.0*(-30.0*std::sqrt(a)*(150420217177132402.5e1 + 2*a*(11570785936702492.5e1 + 8*a*(465611205861301.5e1 + 2*a*(16877165244439.5e1 + a*(2196137366923.5e1 + 2*a*(44317398625.5e1 + 4*a*(718677745.5e1 + a*(15248765.4e1 + a*(430226.5e1 + 2879.4e1*a))))))))) + 46189*(48849363650587.5e1 + 4*a*(-6262738929562.5e1 + a*(1523368928812.5e1 + 8*a*(-29016551025.0e1 + a*(3077512987.5e1 + 4*a*(-59564767.5e1 + a*(3423262.5e1 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)))))))))*expF(a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a))))/(5.36870912e8*std::pow(a,10.5)*expF(a)) ; 
      else
      result = (41*(30*std::sqrt(b)*(-1504202171771324025 + 231415718734049850*b - 74497792937808240*std::pow(b,2) + 5400692878220640*std::pow(b,3) - 702763957415520*std::pow(b,4) + 28363135120320*std::pow(b,5) - 1839815028480*std::pow(b,6) + 39036839424*std::pow(b,7) - 1101379840*std::pow(b,8) + 7371264*std::pow(b,9))*expF(b) + 46189*(488493636505875 + 250509557182500*b + 60934757152500*std::pow(b,2) + 9285296328000*std::pow(b,3) + 984804156000*std::pow(b,4) + 76242902400*std::pow(b,5) + 4381776000*std::pow(b,6) + 185472000*std::pow(b,7) + 5564160*std::pow(b,8) + 107520*std::pow(b,9) + 1024*std::pow(b,10))*utl::SQRTPI*Erfi(std::sqrt(b))))/(5.36870912e8*std::pow(b,10.5)); 
   //    result =(41*expF(b)*(15*std::sqrt(b)*(-1504202171771324025 + 2*b*(115707859367024925 + 8*b*(-4656112058613015 + 
   //              2*b*(168771652444395 + b*(-21961373669235 + 2*b*(443173986255 + 4*b*(-7186777455 + b*(152487654 + b*(-4302265 + 28794*b))))))))) + 
   //     46189*(488493636505875 + 4*b*(62627389295625 + b*(15233689288125 + 8*b*(290165510250 + b*(30775129875 + 4*b*(595647675 + b*(34232625 + 2*b*(724500 + b*(21735 + 4*b*(105 + b))))))))))*DawsonF(std::sqrt(b))))/
   // (2.68435456e8*std::pow(b,10.5)); 
      break;
    case 22 : 
      if (a>0)
        {
        result = 45*(-46* std::sqrt(a)*(15722777246044531162.5e1 + 4*a*(609409970776919812.5e1 + a*(197983922213379802.5e1 + 8*a*(1842440222557935.0e1 + a*(247349944474132.5e1 + 4*a*(2659408317247.5e1 + a*(187282275862.5e1 + 2.0*a*(2316815046.0e1 + a*(83165764.5e1 + 4*a*(238528.5e1 + 2662.9e1*a)))))))))) -  8817.9e1*(-4101020386475512.5e1 + 2*a*(1049098238400712.5e1 + 2*a*(-127938809561062.5e1 +    2*a*(9841446889312.5e1 + 4*a*(-265985051062.5e1 + 2*a*(10639402042.5e1 + 2*a*(-322406122.5e1 + 2*a*(7428712.5e1 + a*(-256162.5e1 + 2*a*(3162.5e1 - 50.6e1*a + 4.0* std::pow(a,2)))))))))))*   expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)))  ; 
      result /= (2.147483648e9* std::pow(a,11.5)* expF(a));
        }
      else
        result =(45*(46*std::sqrt(b)*(157227772460445311625. - 24376398831076792500.*b + 7919356888535192100*std::pow(b,2) - 589580871218539200.*std::pow(b,3) + 79151982231722400.*std::pow(b,4) - 3404042646076800.*std::pow(b,5) + 239721313104000.*std::pow(b,6) - 5931046517760.*std::pow(b,7) + 212904357120.*std::pow(b,8) - 2442531840.*std::pow(b,9) + 27268096.*std::pow(b,10))*expF(b) - 88179*(41010203864755125. + 20981964768014250.*b + 5117552382442500.*std::pow(b,2) + 787315751145000.*std::pow(b,3) + 85115216340000.*std::pow(b,4) + 6809217307200.*std::pow(b,5) + 412679836800.*std::pow(b,6) + 19017504000.*std::pow(b,7) + 655776000.*std::pow(b,8) + 16192000.*std::pow(b,9) + 259072.*std::pow(b,10) + 2048*std::pow(b,11))*utl::SQRTPI*Erfi(std::sqrt(b))))/(2.147483648e9*std::pow(b,11.5)); 
      break;
    case 24 :
      if (a>0)
      result =(49.0*(-30.0* std::sqrt(a)*(16996322202974138186662.5e1 + 2*a*(1325954214416422128037.5e1 + 
             2*a*(216974325995414530042.5e1 + 2*a*(8249669832974417347.5e1 + 4*a*
                    (283606141073075320.5e1 + 2*a*(6398463995484151.5e1 + 2*a*(238903097649898.5e1 + 2*a*(3282810567991.5e1 + a*(136370408128.5e1 + 2*a*(1026865776.3e1 + 2*a*(9717230.9e1 + 46910.2e1*a)))))))
                   )))) + 67603.9e1*(377115570321552562.5e1 + 8.0*a*(-24071206616269312.5e1 + 
             a*(5884072728421387.5e1 + 2*a*(-456129668869875.0e1 + a*(50063012436937.5e1 + 
                      16.0*a*(-256733397112.5e1 + a*(16190394412.5e1 + a*(-792998910.0e1 + a*(30037837.5e1 + 8.0*a*(-107662.5e1 + a*(2227.5e1 + 2*(-15.0e1 + a)*a)))))))))))*expF(a)*utl::SQRTPI*std::tr1::erf(std::sqrt(a))))/
   (3.4359738368e10*std::pow(a,12.5)*expF(a)) ;
      else
      // result = (49*(30*std::sqrt(b)*(-169963222029741381866625. + 26519084288328442560750.*b - 8678973039816581201700.*std::pow(b,2) + 659973586637953387800.*std::pow(b,3) - 90753965143384102560.*std::pow(b,4) + 4095016957109856960.*std::pow(b,5) - 305795964991870080.*std::pow(b,6) + 8403995054058240.*std::pow(b,7) - 349108244808960.*std::pow(b,8) + 5257552774656.*std::pow(b,9) - 99504444416.*std::pow(b,10) + 480360448.*std::pow(b,11))*expF(b) + 676039.*(3771155703215525625. + 1925696529301545000.*b + 470725818273711000.*std::pow(b,2) + 72980747019180000.*std::pow(b,3) + 8010081989910000.*std::pow(b,4) + 657237496608000.*std::pow(b,5) + 41447409696000.*std::pow(b,6) + 2030077209600.*std::pow(b,7) + 76896864000.*std::pow(b,8) + 2204928000.*std::pow(b,9) + 45619200.*std::pow(b,10) + 614400.*std::pow(b,11) + 4096.*std::pow(b,12))*utl::SQRTPI*Erfi(std::sqrt(b))))/(3.4359738368e10*std::pow(b,12.5));
      result =(49*expF(b)*(15*std::sqrt(b)*(-169963222029741381866625. + 2*b*(13259542144164221280375. + 2*b*
              (-2169743259954145300425. + 2*b*(82496698329744173475. + 4*b*(-2836061410730753205. + 
                      2*b*(63984639954841515. + 2*b*(-2389030976498985. + 2*b*(32828105679915. + b*(-1363704081285. + 2*b*(10268657763. + 2*b*(-97172309. + 469102.*b))))))))))) + 
       676039*(3771155703215525625. + 8*b*(240712066162693125. + b*(58840727284213875 + 2*b*
                 (4561296688698750. + b*(500630124369375. + 16*b*(2567333971125 + b*(161903944125 + b*(7929989100. + b*(300378375. + 8*b*(1076625. + b*(22275. + 2*b*(150. + b))))))))))))*DawsonF(std::sqrt(b))))/
   (1.7179869184e10*std::pow(b,12.5)); 
      break;
    default :
      utlException(true, "l is too big. l=" << l);
       break;
    }
  return result;
}

inline double
GetExpLegendreCoefDerivative(const double a, const int l, Func1To1 expF=&std::exp )
{
  // coefficient for odd order is 0;
  if (!utl::IsInt(0.5*l))
    return 0;

  double result = -1000.0;
  switch ( l )
    {
    case 0 :
       result = 1/(2.*a* expF(a)) - ( utl::SQRTPI* std::tr1::erf( std::sqrt(a)))/(4.* std::pow(a,1.5)) ; 
       break;
    case 2 :
       result = (5*(3/(2.* std::pow(a,2)* expF(a)) + (3 - 2*a)/(4.* std::pow(a,2)* expF(a)) + 3/(2.*a* expF(a)) - (3*(3 - 2*a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)))/(8.* std::pow(a,2.5)) - ( utl::SQRTPI* std::tr1::erf( std::sqrt(a)))/(2.* std::pow(a,1.5))))/2. ; 
       break;
    case 4 :
       result = (9*(-5/(8.* std::pow(a,2)* expF(a)) + (5*(21 + 2*a))/(8.* std::pow(a,3)* expF(a)) + (5*(21 + 2*a))/(16.* std::pow(a,2)* expF(a)) + (3*(35 + 4*(-5 + a)*a))/(32.* std::pow(a,3)* expF(a)) + (3*(4*(-5 + a) + 4*a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)))/(32.* std::pow(a,2.5)) - (15*(35 + 4*(-5 + a)*a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)))/(64.* std::pow(a,3.5))))/2. ; 
       break;
    case 6 : 
       result = (-91*(-42* std::sqrt(a)*(165 + 20*a + 4* std::pow(a,2)) - 5*(-693 + 378*a - 84* std::pow(a,2) + 8* std::pow(a,3))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(512.* std::pow(a,4.5)* expF(a)) - (13*(-42* std::sqrt(a)*(165 + 20*a + 4* std::pow(a,2)) - 5*(-693 + 378*a - 84* std::pow(a,2) + 8* std::pow(a,3))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(256.* std::pow(a,3.5)* expF(a)) + (13*(-42* std::sqrt(a)*(20 + 8*a) - (21*(165 + 20*a + 4* std::pow(a,2)))/ std::sqrt(a) - (5*(-693 + 378*a - 84* std::pow(a,2) + 8* std::pow(a,3)))/ std::sqrt(a) - 5*(378 - 168*a + 24* std::pow(a,2))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)) - 5*(-693 + 378*a - 84* std::pow(a,2) + 8* std::pow(a,3))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(256.* std::pow(a,3.5)* expF(a)) ;
       break;
    case 8 :
       result = (-153*(-6* std::sqrt(a)*(225225 + 2*a*(15015 + 2*a*(1925 + 62*a))) + 35*(19305 + 8*a*(-1287 + a*(297 + 2*(-18 + a)*a)))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(8192.* std::pow(a,5.5)* expF(a)) - (17*(-6* std::sqrt(a)*(225225 + 2*a*(15015 + 2*a*(1925 + 62*a))) + 35*(19305 + 8*a*(-1287 + a*(297 + 2*(-18 + a)*a)))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(4096.* std::pow(a,4.5)* expF(a)) + (17*(-6* std::sqrt(a)*(2*a*(124*a + 2*(1925 + 62*a)) + 2*(15015 + 2*a*(1925 + 62*a))) - (3*(225225 + 2*a*(15015 + 2*a*(1925 + 62*a))))/ std::sqrt(a) + (35*(19305 + 8*a*(-1287 + a*(297 + 2*(-18 + a)*a))))/ std::sqrt(a) + 35*(8*a*(297 + 2*(-18 + a)*a + a*(2*(-18 + a) + 2*a)) + 8*(-1287 + a*(297 + 2*(-18 + a)*a)))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)) + 35*(19305 + 8*a*(-1287 + a*(297 + 2*(-18 + a)*a)))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(4096.* std::pow(a,4.5)* expF(a)) ;
       break;
    case 10 : 
       result = (-231*(-22* std::sqrt(a)*(3968055 + 556920*a + 157248* std::pow(a,2) + 7488* std::pow(a,3) + 464* std::pow(a,4)) - 63*(-692835 + 364650*a - 85800* std::pow(a,2) + 11440* std::pow(a,3) - 880* std::pow(a,4) + 32* std::pow(a,5))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(32768.* std::pow(a,6.5)* expF(a)) - (21*(-22* std::sqrt(a)*(3968055 + 556920*a + 157248* std::pow(a,2) + 7488* std::pow(a,3) + 464* std::pow(a,4)) - 63*(-692835 + 364650*a - 85800* std::pow(a,2) + 11440* std::pow(a,3) - 880* std::pow(a,4) + 32* std::pow(a,5))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(16384.* std::pow(a,5.5)* expF(a)) + (21*(-22* std::sqrt(a)*(556920 + 314496*a + 22464* std::pow(a,2) + 1856* std::pow(a,3)) - (11*(3968055 + 556920*a + 157248* std::pow(a,2) + 7488* std::pow(a,3) + 464* std::pow(a,4)))/ std::sqrt(a) - (63*(-692835 + 364650*a - 85800* std::pow(a,2) + 11440* std::pow(a,3) - 880* std::pow(a,4) + 32* std::pow(a,5)))/ std::sqrt(a) - 63*(364650 - 171600*a + 34320* std::pow(a,2) - 3520* std::pow(a,3) + 160* std::pow(a,4))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)) - 63*(-692835 + 364650*a - 85800* std::pow(a,2) + 11440* std::pow(a,3) - 880* std::pow(a,4) + 32* std::pow(a,5))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(16384.* std::pow(a,5.5)* expF(a)) ;
      break;
    case 12 : 
      result = (-325*(-26* std::sqrt(a)*(540571185 + 2*a*(39171825 + 4*a*(2909907 + 2*a*(82467 + a*(7469 + 122*a))))) + 231*(30421755 + 4*a*(-3968055 + a*(944775 + 4*a*(-33150 + a*(2925 + 4*(-39 + a)*a)))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(262144.* std::pow(a,7.5)* expF(a)) - (25*(-26* std::sqrt(a)*(540571185 + 2*a*(39171825 + 4*a*(2909907 + 2*a*(82467 + a*(7469 + 122*a))))) + 231*(30421755 + 4*a*(-3968055 + a*(944775 + 4*a*(-33150 + a*(2925 + 4*(-39 + a)*a)))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(131072.* std::pow(a,6.5)* expF(a)) + (25*(-26* std::sqrt(a)*(2*a*(4*a*(2*a*(7469 + 244*a) + 2*(82467 + a*(7469 + 122*a))) + 4*(2909907 + 2*a*(82467 + a*(7469 + 122*a)))) + 2*(39171825 + 4*a*(2909907 + 2*a*(82467 + a*(7469 + 122*a))))) - (13*(540571185 + 2*a*(39171825 + 4*a*(2909907 + 2*a*(82467 + a*(7469 + 122*a))))))/ std::sqrt(a) + (231*(30421755 + 4*a*(-3968055 + a*(944775 + 4*a*(-33150 + a*(2925 + 4*(-39 + a)*a))))))/ std::sqrt(a) + 231*(4*a*(944775 + 4*a*(-33150 + a*(2925 + 4*(-39 + a)*a)) + a*(4*a*(2925 + 4*(-39 + a)*a + a*(4*(-39 + a) + 4*a)) + 4*(-33150 + a*(2925 + 4*(-39 + a)*a)))) + 4*(-3968055 + a*(944775 + 4*a*(-33150 + a*(2925 + 4*(-39 + a)*a)))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)) + 231*(30421755 + 4*a*(-3968055 + a*(944775 + 4*a*(-33150 + a*(2925 + 4*(-39 + a)*a)))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(131072.* std::pow(a,6.5)* expF(a)) ;
      break;
    case 14 : 
      result = (-435*(-2* std::sqrt(a)*(677644592625 + 100391791500*a + 30786816060* std::pow(a,2) + 1928852640* std::pow(a,3) + 206187696* std::pow(a,4) + 5360576* std::pow(a,5) + 158528* std::pow(a,6)) - 429*(-1579591125 + 819047250*a - 196571340* std::pow(a,2) + 28488600* std::pow(a,3) - 2713200* std::pow(a,4) + 171360* std::pow(a,5) - 6720* std::pow(a,6) + 128* std::pow(a,7))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(1.048576e6* std::pow(a,8.5)* expF(a)) - (29*(-2* std::sqrt(a)*(677644592625 + 100391791500*a + 30786816060* std::pow(a,2) + 1928852640* std::pow(a,3) + 206187696* std::pow(a,4) + 5360576* std::pow(a,5) + 158528* std::pow(a,6)) - 429*(-1579591125 + 819047250*a - 196571340* std::pow(a,2) + 28488600* std::pow(a,3) - 2713200* std::pow(a,4) + 171360* std::pow(a,5) - 6720* std::pow(a,6) + 128* std::pow(a,7))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(524288.* std::pow(a,7.5)* expF(a)) + (29*(-2* std::sqrt(a)*(100391791500 + 61573632120*a + 5786557920* std::pow(a,2) + 824750784* std::pow(a,3) + 26802880* std::pow(a,4) + 951168* std::pow(a,5)) - (677644592625 + 100391791500*a + 30786816060* std::pow(a,2) + 1928852640* std::pow(a,3) + 206187696* std::pow(a,4) + 5360576* std::pow(a,5) + 158528* std::pow(a,6))/ std::sqrt(a) - (429*(-1579591125 + 819047250*a - 196571340* std::pow(a,2) + 28488600* std::pow(a,3) - 2713200* std::pow(a,4) + 171360* std::pow(a,5) - 6720* std::pow(a,6) + 128* std::pow(a,7)))/ std::sqrt(a) - 429*(819047250 - 393142680*a + 85465800* std::pow(a,2) - 10852800* std::pow(a,3) + 856800* std::pow(a,4) - 40320* std::pow(a,5) + 896* std::pow(a,6))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)) - 429*(-1579591125 + 819047250*a - 196571340* std::pow(a,2) + 28488600* std::pow(a,3) - 2713200* std::pow(a,4) + 171360* std::pow(a,5) - 6720* std::pow(a,6) + 128* std::pow(a,7))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(524288.* std::pow(a,7.5)* expF(a)) ;
      break;
    case 16 :
      result = (-0.00001671910285949707*(-34.* std::sqrt(a)*(3.5835440515875e13 + 2.*a*(2.697291221625e12 + 2.*a*(4.22531334225e11 + 2.*a*(1.4130777375e10 + 2.*a*(8.29568025e8 + 2.*a*(1.3438035e7 + 624150.*a + 6196.* std::pow(a,2))))))) + 11405.740530576995* std::pow(2.718281828459045,a)*(9.4670161425e10 + 16.*a*(-3.053876175e9 + a*(7.37142525e8 + 2.*a*(-5.460315e7 + a*(5.460315e6 + 8.*a*(-47481. + a*(2261. + (-68. + a)*a)))))))* std::tr1::erf( std::sqrt(a))))/( std::pow(2.718281828459045,1.*a)* std::pow(a,9.5)) - (1.9669532775878906e-6*(-34.* std::sqrt(a)*(3.5835440515875e13 + 2.*a*(2.697291221625e12 + 2.*a*(4.22531334225e11 + 2.*a*(1.4130777375e10 + 2.*a*(8.29568025e8 + 2.*a*(1.3438035e7 + 624150.*a + 6196.* std::pow(a,2))))))) + 11405.740530576995* std::pow(2.718281828459045,a)*(9.4670161425e10 + 16.*a*(-3.053876175e9 + a*(7.37142525e8 + 2.*a*(-5.460315e7 + a*(5.460315e6 + 8.*a*(-47481. + a*(2261. + (-68. + a)*a)))))))* std::tr1::erf( std::sqrt(a))))/( std::pow(2.718281828459045,1.*a)* std::pow(a,8.5)) + (1.9669532775878906e-6*(-34.* std::sqrt(a)*(2.*a*(2.*a*(2.*a*(2.*a*(2.*a*(624150. + 12392.*a) + 2.*(1.3438035e7 + 624150.*a + 6196.* std::pow(a,2))) + 2.*(8.29568025e8 + 2.*a*(1.3438035e7 + 624150.*a + 6196.* std::pow(a,2)))) + 2.*(1.4130777375e10 + 2.*a*(8.29568025e8 + 2.*a*(1.3438035e7 + 624150.*a + 6196.* std::pow(a,2))))) + 2.*(4.22531334225e11 + 2.*a*(1.4130777375e10 + 2.*a*(8.29568025e8 + 2.*a*(1.3438035e7 + 624150.*a + 6196.* std::pow(a,2)))))) + 2.*(2.697291221625e12 + 2.*a*(4.22531334225e11 + 2.*a*(1.4130777375e10 + 2.*a*(8.29568025e8 + 2.*a*(1.3438035e7 + 624150.*a + 6196.* std::pow(a,2))))))) - (17.*(3.5835440515875e13 + 2.*a*(2.697291221625e12 + 2.*a*(4.22531334225e11 + 2.*a*(1.4130777375e10 + 2.*a*(8.29568025e8 + 2.*a*(1.3438035e7 + 624150.*a + 6196.* std::pow(a,2))))))))/ std::sqrt(a) + (6435.*(9.4670161425e10 + 16.*a*(-3.053876175e9 + a*(7.37142525e8 + 2.*a*(-5.460315e7 + a*(5.460315e6 + 8.*a*(-47481. + a*(2261. + (-68. + a)*a))))))))/ std::sqrt(a) + 11405.740530576995* std::pow(2.718281828459045,a)*(16.*a*(7.37142525e8 + 2.*a*(-5.460315e7 + a*(5.460315e6 + 8.*a*(-47481. + a*(2261. + (-68. + a)*a)))) + a*(2.*a*(5.460315e6 + 8.*a*(-47481. + a*(2261. + (-68. + a)*a)) + a*(8.*a*(2261. + (-68. + a)*a + a*(-68. + 2.*a)) + 8.*(-47481. + a*(2261. + (-68. + a)*a)))) + 2.*(-5.460315e7 + a*(5.460315e6 + 8.*a*(-47481. + a*(2261. + (-68. + a)*a)))))) + 16.*(-3.053876175e9 + a*(7.37142525e8 + 2.*a*(-5.460315e7 + a*(5.460315e6 + 8.*a*(-47481. + a*(2261. + (-68. + a)*a)))))))* std::tr1::erf( std::sqrt(a)) + 11405.740530576995* std::pow(2.718281828459045,a)*(9.4670161425e10 + 16.*a*(-3.053876175e9 + a*(7.37142525e8 + 2.*a*(-5.460315e7 + a*(5.460315e6 + 8.*a*(-47481. + a*(2261. + (-68. + a)*a)))))))* std::tr1::erf( std::sqrt(a))))/( std::pow(2.718281828459045,1.*a)* std::pow(a,8.5)); 
      // result = (-561.*(-34.* std::sqrt(a)*(3583544051587.5e1 + 2*a*(269729122162.5e1 + 2*a*(42253133422.5e1 + 2*a*(1413077737.5e1 + 2*a*(82956802.5e1 + 2*a*(1343803.5e1 + 62415.0e1*a + 619.6e1* std::pow(a,2.))))))) + 643.5e1*(9467016142.5e1 + 1.6e1*a*(-305387617.5e1 + a*(73714252.5e1 + 2*a*(-5460315.0e1 + a*(546031.5e1 + 8*a*(-4748.1e1 + a*(226.1e1 + (-68 + a)*a)))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(3.3554432e7* std::pow(a,9.5)* expF(a)) - (33*(-34* std::sqrt(a)*(3583544051587.5e1 + 2*a*(269729122162.5e1 + 2*a*(42253133422.5e1 + 2*a*(1413077737.5e1 + 2*a*(82956802.5e1 + 2*a*(1343803.5e1 + 62415.0e1*a + 619.6e1* std::pow(a,2))))))) + 643.5e1*(9467016142.5e1 + 16*a*(-305387617.5e1 + a*(73714252.5e1 + 2*a*(-5460315.0e1 + a*(546031.5e1 + 8*a*(-4748.1e1 + a*(226.1e1 + (-6.8e1 + a)*a)))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(1.6777216e7* std::pow(a,8.5)* expF(a)) + (33*(-34* std::sqrt(a)*(2*a*(2*a*(2*a*(2*a*(2*a*(62415.0e1 + 1239.2e1*a) + 2*(1343803.5e1 + 62415.0e1*a + 619.6e1* std::pow(a,2))) + 2*(82956802.5e1 + 2*a*(1343803.5e1 + 62415.0e1*a + 6196.* std::pow(a,2)))) + 2*(14130777375. + 2*a*(829568025. + 2*a*(13438035. + 624150.*a + 6196.* std::pow(a,2))))) + 2*(422531334225. + 2*a*(14130777375. + 2*a*(829568025. + 2*a*(13438035. + 624150.*a + 6196.* std::pow(a,2)))))) + 2*(2697291221625. + 2*a*(422531334225. + 2*a*(14130777375. + 2*a*(829568025. + 2*a*(13438035. + 624150.*a + 6196.* std::pow(a,2))))))) - (17*(35835440515875. + 2*a*(2697291221625. + 2*a*(422531334225. + 2*a*(14130777375. + 2*a*(829568025. + 2*a*(13438035. + 624150.*a + 6196.* std::pow(a,2))))))))/ std::sqrt(a) + (6435.*(94670161425. + 16.*a*(-3053876175. + a*(737142525. + 2*a*(-54603150. + a*(5460315. + 8*a*(-47481. + a*(2261. + (-68. + a)*a))))))))/ std::sqrt(a) + 6435.*(16*a*(737142525. + 2*a*(-54603150. + a*(5460315. + 8*a*(-47481. + a*(2261. + (-68. + a)*a)))) + a*(2*a*(5460315. + 8*a*(-47481. + a*(2261. + (-68 + a)*a)) + a*(8*a*(2261. + (-68. + a)*a + a*(-68 + 2*a)) + 8*(-47481. + a*(2261. + (-68 + a)*a)))) + 2*(-54603150. + a*(5460315. + 8*a*(-47481. + a*(2261. + (-68. + a)*a)))))) + 16*(-3053876175. + a*(737142525. + 2*a*(-54603150 + a*(5460315. + 8*a*(-47481. + a*(2261. + (-68 + a)*a)))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)) + 6435.*(94670161425. + 16.*a*(-3053876175 + a*(737142525. + 2.*a*(-54603150. + a*(5460315. + 8.*a*(-47481. + a*(2261. + (-68 + a)*a)))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(1.6777216e7* std::pow(a,8.5)* expF(a)) ;
      break;
    case 18 : 
      result = (-5.237758159637451e-6*(-114.* std::sqrt(a)*(1.371596246645625e15 + 8.*a*(2.612564279325e13 + a*(8.312704525125e12 + 4.*a*(1.45568097675e11 + 2.*a*(9.078873375e9 + 2.*a*(1.67739e8 + a*(9.612735e6 + 2.*a*(79390. + 1363.*a)))))))) - 21544.176557756546* std::pow(2.718281828459045,a)*(-6.432002143875e12 + 2.*a*(1.653943408425e12 + 8.*a*(-5.0119497225e10 + 2.*a*(3.772435275e9 + a*(-3.90251925e8 + 2.*a*(1.4453775e7 + 4.*a*(-192717. + a*(7182. + a*(-171. + 2.*a)))))))))* std::tr1::erf( std::sqrt(a))))/( std::pow(2.718281828459045,1.*a)* std::pow(a,10.5)) - (5.513429641723633e-7*(-114.* std::sqrt(a)*(1.371596246645625e15 + 8.*a*(2.612564279325e13 + a*(8.312704525125e12 + 4.*a*(1.45568097675e11 + 2.*a*(9.078873375e9 + 2.*a*(1.67739e8 + a*(9.612735e6 + 2.*a*(79390. + 1363.*a)))))))) - 21544.176557756546* std::pow(2.718281828459045,a)*(-6.432002143875e12 + 2.*a*(1.653943408425e12 + 8.*a*(-5.0119497225e10 + 2.*a*(3.772435275e9 + a*(-3.90251925e8 + 2.*a*(1.4453775e7 + 4.*a*(-192717. + a*(7182. + a*(-171. + 2.*a)))))))))* std::tr1::erf( std::sqrt(a))))/( std::pow(2.718281828459045,1.*a)* std::pow(a,9.5)) + (5.513429641723633e-7*(-114.* std::sqrt(a)*(8.*a*(8.312704525125e12 + 4.*a*(1.45568097675e11 + 2.*a*(9.078873375e9 + 2.*a*(1.67739e8 + a*(9.612735e6 + 2.*a*(79390. + 1363.*a))))) + a*(4.*a*(2.*a*(2.*a*(9.612735e6 + 2.*a*(79390. + 1363.*a) + a*(2726.*a + 2.*(79390. + 1363.*a))) + 2.*(1.67739e8 + a*(9.612735e6 + 2.*a*(79390. + 1363.*a)))) + 2.*(9.078873375e9 + 2.*a*(1.67739e8 + a*(9.612735e6 + 2.*a*(79390. + 1363.*a))))) + 4.*(1.45568097675e11 + 2.*a*(9.078873375e9 + 2.*a*(1.67739e8 + a*(9.612735e6 + 2.*a*(79390. + 1363.*a))))))) + 8.*(2.612564279325e13 + a*(8.312704525125e12 + 4.*a*(1.45568097675e11 + 2.*a*(9.078873375e9 + 2.*a*(1.67739e8 + a*(9.612735e6 + 2.*a*(79390. + 1363.*a)))))))) - (57.*(1.371596246645625e15 + 8.*a*(2.612564279325e13 + a*(8.312704525125e12 + 4.*a*(1.45568097675e11 + 2.*a*(9.078873375e9 + 2.*a*(1.67739e8 + a*(9.612735e6 + 2.*a*(79390. + 1363.*a)))))))))/ std::sqrt(a) - (12155.*(-6.432002143875e12 + 2.*a*(1.653943408425e12 + 8.*a*(-5.0119497225e10 + 2.*a*(3.772435275e9 + a*(-3.90251925e8 + 2.*a*(1.4453775e7 + 4.*a*(-192717. + a*(7182. + a*(-171. + 2.*a))))))))))/ std::sqrt(a) - 21544.176557756546* std::pow(2.718281828459045,a)*(2.*a*(8.*a*(2.*a*(-3.90251925e8 + 2.*a*(1.4453775e7 + 4.*a*(-192717. + a*(7182. + a*(-171. + 2.*a)))) + a*(2.*a*(4.*a*(7182. + a*(-171. + 2.*a) + a*(-171. + 4.*a)) + 4.*(-192717. + a*(7182. + a*(-171. + 2.*a)))) + 2.*(1.4453775e7 + 4.*a*(-192717. + a*(7182. + a*(-171. + 2.*a)))))) + 2.*(3.772435275e9 + a*(-3.90251925e8 + 2.*a*(1.4453775e7 + 4.*a*(-192717. + a*(7182. + a*(-171. + 2.*a))))))) + 8.*(-5.0119497225e10 + 2.*a*(3.772435275e9 + a*(-3.90251925e8 + 2.*a*(1.4453775e7 + 4.*a*(-192717. + a*(7182. + a*(-171. + 2.*a)))))))) + 2.*(1.653943408425e12 + 8.*a*(-5.0119497225e10 + 2.*a*(3.772435275e9 + a*(-3.90251925e8 + 2.*a*(1.4453775e7 + 4.*a*(-192717. + a*(7182. + a*(-171. + 2.*a)))))))))* std::tr1::erf( std::sqrt(a)) - 21544.176557756546* std::pow(2.718281828459045,a)*(-6.432002143875e12 + 2.*a*(1.653943408425e12 + 8.*a*(-5.0119497225e10 + 2.*a*(3.772435275e9 + a*(-3.90251925e8 + 2.*a*(1.4453775e7 + 4.*a*(-192717. + a*(7182. + a*(-171. + 2.*a)))))))))* std::tr1::erf( std::sqrt(a))))/( std::pow(2.718281828459045,1.*a)* std::pow(a,9.5)); 
      // result = (-703*(-114* std::sqrt(a)*(1371596246645625 + 8*a*(26125642793250 + a*(8312704525125 + 4*a*(145568097675 + 2*a*(9078873375 + 2*a*(167739000 + a*(9612735 + 2*a*(79390 + 1363*a)))))))) - 12155*(-6432002143875 + 2*a*(1653943408425 + 8*a*(-50119497225 + 2*a*(3772435275 + a*(-390251925 + 2*a*(14453775 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a)))))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(1.34217728e8* std::pow(a,10.5)* expF(a)) - (37*(-114* std::sqrt(a)*(1371596246645625 + 8*a*(26125642793250 + a*(8312704525125 + 4*a*(145568097675 + 2*a*(9078873375 + 2*a*(167739000 + a*(9612735 + 2*a*(79390 + 1363*a)))))))) - 12155*(-6432002143875 + 2*a*(1653943408425 + 8*a*(-50119497225 + 2*a*(3772435275 + a*(-390251925 + 2*a*(14453775 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a)))))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(6.7108864e7* std::pow(a,9.5)* expF(a)) + (37*(-114* std::sqrt(a)*(8*a*(8312704525125 + 4*a*(145568097675 + 2*a*(9078873375 + 2*a*(167739000 + a*(9612735 + 2*a*(79390 + 1363*a))))) + a*(4*a*(2*a*(2*a*(9612735 + 2*a*(79390 + 1363*a) + a*(2726*a + 2*(79390 + 1363*a))) + 2*(167739000 + a*(9612735 + 2*a*(79390 + 1363*a)))) + 2*(9078873375 + 2*a*(167739000 + a*(9612735 + 2*a*(79390 + 1363*a))))) + 4*(145568097675 + 2*a*(9078873375 + 2*a*(167739000 + a*(9612735 + 2*a*(79390 + 1363*a))))))) + 8*(26125642793250 + a*(8312704525125 + 4*a*(145568097675 + 2*a*(9078873375 + 2*a*(167739000 + a*(9612735 + 2*a*(79390 + 1363*a)))))))) - (57*(1371596246645625 + 8*a*(26125642793250 + a*(8312704525125 + 4*a*(145568097675 + 2*a*(9078873375 + 2*a*(167739000 + a*(9612735 + 2*a*(79390 + 1363*a)))))))))/ std::sqrt(a) - (12155*(-6432002143875 + 2*a*(1653943408425 + 8*a*(-50119497225 + 2*a*(3772435275 + a*(-390251925 + 2*a*(14453775 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a))))))))))/ std::sqrt(a) - 12155*(2*a*(8*a*(2*a*(-390251925 + 2*a*(14453775 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a)))) + a*(2*a*(4*a*(7182 + a*(-171 + 2*a) + a*(-171 + 4*a)) + 4*(-192717 + a*(7182 + a*(-171 + 2*a)))) + 2*(14453775 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a)))))) + 2*(3772435275 + a*(-390251925 + 2*a*(14453775 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a))))))) + 8*(-50119497225 + 2*a*(3772435275 + a*(-390251925 + 2*a*(14453775 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a)))))))) + 2*(1653943408425 + 8*a*(-50119497225 + 2*a*(3772435275 + a*(-390251925 + 2*a*(14453775 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a)))))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)) - 12155*(-6432002143875 + 2*a*(1653943408425 + 8*a*(-50119497225 + 2*a*(3772435275 + a*(-390251925 + 2*a*(14453775 + 4*a*(-192717 + a*(7182 + a*(-171 + 2*a)))))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(6.7108864e7* std::pow(a,9.5)* expF(a)) ; 
      break;
    case 20 : 
      result = (-8.01868736743927e-7*(-30.* std::sqrt(a)*(1.504202171771324e18 + 2.*a*(1.1570785936702493e17 + 8.*a*(4.656112058613015e15 + 2.*a*(1.68771652444395e14 + a*(2.1961373669235e13 + 2.*a*(4.43173986255e11 + 4.*a*(7.186777455e9 + a*(1.52487654e8 + a*(4.302265e6 + 28794.*a))))))))) + 81867.87091947487* std::pow(2.718281828459045,a)*(4.88493636505875e14 + 4.*a*(-6.2627389295625e13 + a*(1.5233689288125e13 + 8.*a*(-2.9016551025e11 + a*(3.0775129875e10 + 4.*a*(-5.95647675e8 + a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a)))))))))* std::tr1::erf( std::sqrt(a))))/( std::pow(2.718281828459045,1.*a)* std::pow(a,11.5)) - (7.636845111846924e-8*(-30.* std::sqrt(a)*(1.504202171771324e18 + 2.*a*(1.1570785936702493e17 + 8.*a*(4.656112058613015e15 + 2.*a*(1.68771652444395e14 + a*(2.1961373669235e13 + 2.*a*(4.43173986255e11 + 4.*a*(7.186777455e9 + a*(1.52487654e8 + a*(4.302265e6 + 28794.*a))))))))) + 81867.87091947487* std::pow(2.718281828459045,a)*(4.88493636505875e14 + 4.*a*(-6.2627389295625e13 + a*(1.5233689288125e13 + 8.*a*(-2.9016551025e11 + a*(3.0775129875e10 + 4.*a*(-5.95647675e8 + a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a)))))))))* std::tr1::erf( std::sqrt(a))))/( std::pow(2.718281828459045,1.*a)* std::pow(a,10.5)) + (7.636845111846924e-8*(-30.* std::sqrt(a)*(2.*a*(8.*a*(2.*a*(2.1961373669235e13 + 2.*a*(4.43173986255e11 + 4.*a*(7.186777455e9 + a*(1.52487654e8 + a*(4.302265e6 + 28794.*a)))) + a*(2.*a*(4.*a*(1.52487654e8 + a*(4.302265e6 + 28794.*a) + a*(4.302265e6 + 57588.*a)) + 4.*(7.186777455e9 + a*(1.52487654e8 + a*(4.302265e6 + 28794.*a)))) + 2.*(4.43173986255e11 + 4.*a*(7.186777455e9 + a*(1.52487654e8 + a*(4.302265e6 + 28794.*a)))))) + 2.*(1.68771652444395e14 + a*(2.1961373669235e13 + 2.*a*(4.43173986255e11 + 4.*a*(7.186777455e9 + a*(1.52487654e8 + a*(4.302265e6 + 28794.*a))))))) + 8.*(4.656112058613015e15 + 2.*a*(1.68771652444395e14 + a*(2.1961373669235e13 + 2.*a*(4.43173986255e11 + 4.*a*(7.186777455e9 + a*(1.52487654e8 + a*(4.302265e6 + 28794.*a)))))))) + 2.*(1.1570785936702493e17 + 8.*a*(4.656112058613015e15 + 2.*a*(1.68771652444395e14 + a*(2.1961373669235e13 + 2.*a*(4.43173986255e11 + 4.*a*(7.186777455e9 + a*(1.52487654e8 + a*(4.302265e6 + 28794.*a))))))))) - (15.*(1.504202171771324e18 + 2.*a*(1.1570785936702493e17 + 8.*a*(4.656112058613015e15 + 2.*a*(1.68771652444395e14 + a*(2.1961373669235e13 + 2.*a*(4.43173986255e11 + 4.*a*(7.186777455e9 + a*(1.52487654e8 + a*(4.302265e6 + 28794.*a))))))))))/ std::sqrt(a) + (46189.*(4.88493636505875e14 + 4.*a*(-6.2627389295625e13 + a*(1.5233689288125e13 + 8.*a*(-2.9016551025e11 + a*(3.0775129875e10 + 4.*a*(-5.95647675e8 + a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a))))))))))/ std::sqrt(a) + 81867.87091947487* std::pow(2.718281828459045,a)*(4.*a*(1.5233689288125e13 + 8.*a*(-2.9016551025e11 + a*(3.0775129875e10 + 4.*a*(-5.95647675e8 + a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a)))))) + a*(8.*a*(3.0775129875e10 + 4.*a*(-5.95647675e8 + a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a)))) + a*(4.*a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a)) + a*(2.*a*(21735. + 4.*(-105. + a)*a + a*(4.*(-105. + a) + 4.*a)) + 2.*(-724500. + a*(21735. + 4.*(-105. + a)*a)))) + 4.*(-5.95647675e8 + a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a)))))) + 8.*(-2.9016551025e11 + a*(3.0775129875e10 + 4.*a*(-5.95647675e8 + a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a)))))))) + 4.*(-6.2627389295625e13 + a*(1.5233689288125e13 + 8.*a*(-2.9016551025e11 + a*(3.0775129875e10 + 4.*a*(-5.95647675e8 + a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a)))))))))* std::tr1::erf( std::sqrt(a)) + 81867.87091947487* std::pow(2.718281828459045,a)*(4.88493636505875e14 + 4.*a*(-6.2627389295625e13 + a*(1.5233689288125e13 + 8.*a*(-2.9016551025e11 + a*(3.0775129875e10 + 4.*a*(-5.95647675e8 + a*(3.4232625e7 + 2.*a*(-724500. + a*(21735. + 4.*(-105. + a)*a)))))))))* std::tr1::erf( std::sqrt(a))))/( std::pow(2.718281828459045,1.*a)* std::pow(a,10.5)); 
      // result = (-861*(-30* std::sqrt(a)*(1504202171771324025 + 2*a*(115707859367024925 + 8*a*(4656112058613015 + 2*a*(168771652444395 + a*(21961373669235 + 2*a*(443173986255 + 4*a*(7186777455 + a*(152487654 + a*(4302265 + 28794*a))))))))) + 46189*(488493636505875 + 4*a*(-62627389295625 + a*(15233689288125 + 8*a*(-290165510250 + a*(30775129875 + 4*a*(-595647675 + a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)))))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(1.073741824e9* std::pow(a,11.5)* expF(a)) - (41*(-30* std::sqrt(a)*(1504202171771324025 + 2*a*(115707859367024925 + 8*a*(4656112058613015 + 2*a*(168771652444395 + a*(21961373669235 + 2*a*(443173986255 + 4*a*(7186777455 + a*(152487654 + a*(4302265 + 28794*a))))))))) + 46189*(488493636505875 + 4*a*(-62627389295625 + a*(15233689288125 + 8*a*(-290165510250 + a*(30775129875 + 4*a*(-595647675 + a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)))))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(5.36870912e8* std::pow(a,10.5)* expF(a)) + (41*(-30* std::sqrt(a)*(2*a*(8*a*(2*a*(21961373669235 + 2*a*(443173986255 + 4*a*(7186777455 + a*(152487654 + a*(4302265 + 28794*a)))) + a*(2*a*(4*a*(152487654 + a*(4302265 + 28794*a) + a*(4302265 + 57588*a)) + 4*(7186777455 + a*(152487654 + a*(4302265 + 28794*a)))) + 2*(443173986255 + 4*a*(7186777455 + a*(152487654 + a*(4302265 + 28794*a)))))) + 2*(168771652444395 + a*(21961373669235 + 2*a*(443173986255 + 4*a*(7186777455 + a*(152487654 + a*(4302265 + 28794*a))))))) + 8*(4656112058613015 + 2*a*(168771652444395 + a*(21961373669235 + 2*a*(443173986255 + 4*a*(7186777455 + a*(152487654 + a*(4302265 + 28794*a)))))))) + 2*(115707859367024925 + 8*a*(4656112058613015 + 2*a*(168771652444395 + a*(21961373669235 + 2*a*(443173986255 + 4*a*(7186777455 + a*(152487654 + a*(4302265 + 28794*a))))))))) - (15*(1504202171771324025 + 2*a*(115707859367024925 + 8*a*(4656112058613015 + 2*a*(168771652444395 + a*(21961373669235 + 2*a*(443173986255 + 4*a*(7186777455 + a*(152487654 + a*(4302265 + 28794*a))))))))))/ std::sqrt(a) + (46189*(488493636505875 + 4*a*(-62627389295625 + a*(15233689288125 + 8*a*(-290165510250 + a*(30775129875 + 4*a*(-595647675 + a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a))))))))))/ std::sqrt(a) + 46189*(4*a*(15233689288125 + 8*a*(-290165510250 + a*(30775129875 + 4*a*(-595647675 + a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)))))) + a*(8*a*(30775129875 + 4*a*(-595647675 + a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)))) + a*(4*a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)) + a*(2*a*(21735 + 4*(-105 + a)*a + a*(4*(-105 + a) + 4*a)) + 2*(-724500 + a*(21735 + 4*(-105 + a)*a)))) + 4*(-595647675 + a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)))))) + 8*(-290165510250 + a*(30775129875 + 4*a*(-595647675 + a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)))))))) + 4*(-62627389295625 + a*(15233689288125 + 8*a*(-290165510250 + a*(30775129875 + 4*a*(-595647675 + a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)))))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a)) + 46189*(488493636505875 + 4*a*(-62627389295625 + a*(15233689288125 + 8*a*(-290165510250 + a*(30775129875 + 4*a*(-595647675 + a*(34232625 + 2*a*(-724500 + a*(21735 + 4*(-105 + a)*a)))))))))* expF(a)* utl::SQRTPI* std::tr1::erf( std::sqrt(a))))/(5.36870912e8* std::pow(a,10.5)* expF(a)) ; 
      break;
    default :
      utlException(true, "l is too big. l=" << l);
       break;
    }
  return result;

}

inline std::vector<std::complex<double> >
PolynomialRoot(const std::vector<double>& coef)
{
  int order = coef.size()-1;
  utlSAException(order==0 || order>=5)(order).msg("polynomial order should be 1,2,3,4");
  std::vector<std::complex<double> > result(order);
  utlException(coef[order]==0, "the first coef cannot be zero");
  switch ( order )
    {
    case 1 :
       {
       result[0] = std::complex<double>(-coef[0]/coef[1], 0);
       return result;
       }
    case 2 :
       {
       double tmp = coef[1]*coef[1] - 4.0*coef[0]*coef[2];
       std::complex<double> tmp1(tmp, 0);
       tmp1 = std::sqrt(tmp1);
       result[0] = (-coef[1] - tmp1)/(2.0*coef[2]);
       result[1] = (-coef[1] + tmp1)/(2.0*coef[2]);
       return result;
       }
    case 3 :
       {
       double a0 = coef[0];
       double a1 = coef[1];
       double a2 = coef[2];
       double a3 = coef[3];

       double temp_2_3 = std::pow(2.0,1.0/3.0);
       std::complex<double> temp0 (4.0* PowInteger(-a2*a2+3.0*a1*a3, 3) + PowInteger(-2.0*a2*a2*a2+9.0*a1*a2*a3-27.0*a0*a3*a3, 2), 0 );

       std::complex<double> temp1 = -2.0*a2*a2*a2 + 9.0*a1*a2*a3 - 27.0*a0*a3*a3 + std::sqrt(temp0);
       temp1 = std::pow(temp1, 1.0/3.0);
       double temp2 = a2/(3.0*a3);
       double temp3 = -a2*a2+3.0*a1*a3;

       result[0] = -temp2 - temp_2_3*temp3 / (3.0*a3*temp1) + temp1/(3.0*temp_2_3*a3);
       result[1] = -temp2 + (std::complex<double>(1, utl::SQRT3) *temp3) / (3.0*std::pow(2,2.0/3.0)*a3*temp1) - (std::complex<double>(1,-utl::SQRT3) * temp1) / (6.0*temp_2_3*a3);
       result[2] = -temp2 + (std::complex<double>(1, -utl::SQRT3) *temp3) / (3.0*std::pow(2,2.0/3.0)*a3*temp1) - (std::complex<double>(1,utl::SQRT3) * temp1) / (6.0*temp_2_3*a3);

       return result;
       }
    case 4 :
       {
       utlGlobalException(true, "TODO: for order 4");
       // return;
       }
    default :
       {
       utlSAException(true)(order).msg("polynomial order should be 1,2,3,4");
       // return;
       }
    }
  return std::vector<std::complex<double> >();
}

template <class TMatrixType>
inline auto
DeterminantSmallMatrix ( const TMatrixType& mat, const int row  ) ->utl::remove_reference_t<decltype(mat(0,0))>
{
  switch ( row )
    {
    case 1 : return mat(0,0);
    case 2 : return mat(0,0)*mat(1,1) - mat(0,1)*mat(1,0);
    case 3 : return 
             -mat(0,2)*mat(1,1)*mat(2,0) 
             +mat(0,1)*mat(1,2)*mat(2,0)
             +mat(0,2)*mat(1,0)*mat(2,1) 
             -mat(0,0)*mat(1,2)*mat(2,1)
             -mat(0,1)*mat(1,0)*mat(2,2) 
             +mat(0,0)*mat(1,1)*mat(2,2);
    case 4:
            return 
            + mat(0,0)*mat(1,1)*mat(2,2)*mat(3,3)
            - mat(0,0)*mat(1,1)*mat(3,2)*mat(2,3)
            - mat(0,0)*mat(2,1)*mat(1,2)*mat(3,3)
            + mat(0,0)*mat(2,1)*mat(3,2)*mat(1,3)
            + mat(0,0)*mat(3,1)*mat(1,2)*mat(2,3)
            - mat(0,0)*mat(3,1)*mat(2,2)*mat(1,3)
            - mat(1,0)*mat(0,1)*mat(2,2)*mat(3,3)
            + mat(1,0)*mat(0,1)*mat(3,2)*mat(2,3)
            + mat(1,0)*mat(2,1)*mat(0,2)*mat(3,3)
            - mat(1,0)*mat(2,1)*mat(3,2)*mat(0,3)
            - mat(1,0)*mat(3,1)*mat(0,2)*mat(2,3)
            + mat(1,0)*mat(3,1)*mat(2,2)*mat(0,3)
            + mat(2,0)*mat(0,1)*mat(1,2)*mat(3,3)
            - mat(2,0)*mat(0,1)*mat(3,2)*mat(1,3)
            - mat(2,0)*mat(1,1)*mat(0,2)*mat(3,3)
            + mat(2,0)*mat(1,1)*mat(3,2)*mat(0,3)
            + mat(2,0)*mat(3,1)*mat(0,2)*mat(1,3)
            - mat(2,0)*mat(3,1)*mat(1,2)*mat(0,3)
            - mat(3,0)*mat(0,1)*mat(1,2)*mat(2,3)
            + mat(3,0)*mat(0,1)*mat(2,2)*mat(1,3)
            + mat(3,0)*mat(1,1)*mat(0,2)*mat(2,3)
            - mat(3,0)*mat(1,1)*mat(2,2)*mat(0,3)
            - mat(3,0)*mat(2,1)*mat(0,2)*mat(1,3)
            + mat(3,0)*mat(2,1)*mat(1,2)*mat(0,3);
    default :
         utlGlobalException(true, "size too large");
    }
}

template <class TMatrixType>
inline void
InverseSmallMatrix ( const TMatrixType& mat, TMatrixType& result, const int row  )
{
  typedef utl::remove_reference_t<decltype(mat(0,0))> ValueType;
  switch ( row )
    {
    case 1 :
       {
       result(0,0)=1.0/mat(0,0);
       return;
       }
    case 2 :
       {
       ValueType det = DeterminantSmallMatrix(mat,row);
       utlException(std::abs(det)==0, "det==0, cannot invert the matrix");
       ValueType detInv = 1.0/det;
       result(0,0) =  mat(1,1)*detInv;
       result(0,1) = -mat(0,1)*detInv;
       result(1,0) = -mat(1,0)*detInv;
       result(1,1) =  mat(0,0)*detInv;
       return;
       }
    case 3 :
       {
       ValueType det = DeterminantSmallMatrix(mat,row);
       utlException(std::abs(det)==0, "det==0, cannot invert the matrix");
       ValueType detInv = 1.0/det;
       result(0,0) = (-mat(1,2)*mat(2,1) + mat(1,1)*mat(2,2)) * detInv;
       result(0,1) = ( mat(0,2)*mat(2,1) - mat(0,1)*mat(2,2)) * detInv;
       result(0,2) = (-mat(0,2)*mat(1,1) + mat(0,1)*mat(1,2)) * detInv;
       result(1,0) = ( mat(1,2)*mat(2,0) - mat(1,0)*mat(2,2)) * detInv;
       result(1,1) = (-mat(0,2)*mat(2,0) + mat(0,0)*mat(2,2)) * detInv;
       result(1,2) = ( mat(0,2)*mat(1,0) - mat(0,0)*mat(1,2)) * detInv;
       result(2,0) = (-mat(1,1)*mat(2,0) + mat(1,0)*mat(2,1)) * detInv;
       result(2,1) = ( mat(0,1)*mat(2,0) - mat(0,0)*mat(2,1)) * detInv;
       result(2,2) = (-mat(0,1)*mat(1,0) + mat(0,0)*mat(1,1)) * detInv;
       return;
       }
    case 4 :
       {
       ValueType det = DeterminantSmallMatrix(mat,row);
       utlException(std::abs(det)==0, "det==0, cannot invert the matrix");
       ValueType detInv = 1.0/det;
       result(0,0) = ( mat(1,1)*mat(2,2)*mat(3,3) - mat(1,1)*mat(2,3)*mat(3,2) - mat(2,1)*mat(1,2)*mat(3,3)
                       +mat(2,1)*mat(1,3)*mat(3,2) + mat(3,1)*mat(1,2)*mat(2,3) - mat(3,1)*mat(1,3)*mat(2,2))*detInv;
       result(0,1) = (-mat(0,1)*mat(2,2)*mat(3,3) + mat(0,1)*mat(2,3)*mat(3,2) + mat(2,1)*mat(0,2)*mat(3,3)
                      - mat(2,1)*mat(0,3)*mat(3,2) - mat(3,1)*mat(0,2)*mat(2,3) + mat(3,1)*mat(0,3)*mat(2,2))*detInv;
       result(0,2) = ( mat(0,1)*mat(1,2)*mat(3,3) - mat(0,1)*mat(1,3)*mat(3,2) - mat(1,1)*mat(0,2)*mat(3,3)
                       +mat(1,1)*mat(0,3)*mat(3,2) + mat(3,1)*mat(0,2)*mat(1,3) - mat(3,1)*mat(0,3)*mat(1,2))*detInv;
       result(0,3) = (-mat(0,1)*mat(1,2)*mat(2,3) + mat(0,1)*mat(1,3)*mat(2,2) + mat(1,1)*mat(0,2)*mat(2,3)
                      - mat(1,1)*mat(0,3)*mat(2,2) - mat(2,1)*mat(0,2)*mat(1,3) + mat(2,1)*mat(0,3)*mat(1,2))*detInv;
       result(1,0) = (-mat(1,0)*mat(2,2)*mat(3,3) + mat(1,0)*mat(2,3)*mat(3,2) + mat(2,0)*mat(1,2)*mat(3,3)
                      - mat(2,0)*mat(1,3)*mat(3,2) - mat(3,0)*mat(1,2)*mat(2,3) + mat(3,0)*mat(1,3)*mat(2,2))*detInv;
       result(1,1) = ( mat(0,0)*mat(2,2)*mat(3,3) - mat(0,0)*mat(2,3)*mat(3,2) - mat(2,0)*mat(0,2)*mat(3,3)
                       +mat(2,0)*mat(0,3)*mat(3,2) + mat(3,0)*mat(0,2)*mat(2,3) - mat(3,0)*mat(0,3)*mat(2,2))*detInv;
       result(1,2) = (-mat(0,0)*mat(1,2)*mat(3,3) + mat(0,0)*mat(1,3)*mat(3,2) + mat(1,0)*mat(0,2)*mat(3,3)
                      - mat(1,0)*mat(0,3)*mat(3,2) - mat(3,0)*mat(0,2)*mat(1,3) + mat(3,0)*mat(0,3)*mat(1,2))*detInv;
       result(1,3) = ( mat(0,0)*mat(1,2)*mat(2,3) - mat(0,0)*mat(1,3)*mat(2,2) - mat(1,0)*mat(0,2)*mat(2,3)
                       +mat(1,0)*mat(0,3)*mat(2,2) + mat(2,0)*mat(0,2)*mat(1,3) - mat(2,0)*mat(0,3)*mat(1,2))*detInv;
       result(2,0) = ( mat(1,0)*mat(2,1)*mat(3,3) - mat(1,0)*mat(2,3)*mat(3,1) - mat(2,0)*mat(1,1)*mat(3,3)
                       +mat(2,0)*mat(1,3)*mat(3,1) + mat(3,0)*mat(1,1)*mat(2,3) - mat(3,0)*mat(1,3)*mat(2,1))*detInv;
       result(2,1) = (-mat(0,0)*mat(2,1)*mat(3,3) + mat(0,0)*mat(2,3)*mat(3,1) + mat(2,0)*mat(0,1)*mat(3,3)
                      - mat(2,0)*mat(0,3)*mat(3,1) - mat(3,0)*mat(0,1)*mat(2,3) + mat(3,0)*mat(0,3)*mat(2,1))*detInv;
       result(2,2) = ( mat(0,0)*mat(1,1)*mat(3,3) - mat(0,0)*mat(1,3)*mat(3,1) - mat(1,0)*mat(0,1)*mat(3,3)
                       +mat(1,0)*mat(0,3)*mat(3,1) + mat(3,0)*mat(0,1)*mat(1,3) - mat(3,0)*mat(0,3)*mat(1,1))*detInv;
       result(2,3) = (-mat(0,0)*mat(1,1)*mat(2,3) + mat(0,0)*mat(1,3)*mat(2,1) + mat(1,0)*mat(0,1)*mat(2,3)
                      - mat(1,0)*mat(0,3)*mat(2,1) - mat(2,0)*mat(0,1)*mat(1,3) + mat(2,0)*mat(0,3)*mat(1,1))*detInv;
       result(3,0) = (-mat(1,0)*mat(2,1)*mat(3,2) + mat(1,0)*mat(2,2)*mat(3,1) + mat(2,0)*mat(1,1)*mat(3,2)
                      - mat(2,0)*mat(1,2)*mat(3,1) - mat(3,0)*mat(1,1)*mat(2,2) + mat(3,0)*mat(1,2)*mat(2,1))*detInv;
       result(3,1) = ( mat(0,0)*mat(2,1)*mat(3,2) - mat(0,0)*mat(2,2)*mat(3,1) - mat(2,0)*mat(0,1)*mat(3,2)
                       +mat(2,0)*mat(0,2)*mat(3,1) + mat(3,0)*mat(0,1)*mat(2,2) - mat(3,0)*mat(0,2)*mat(2,1))*detInv;
       result(3,2) = (-mat(0,0)*mat(1,1)*mat(3,2) + mat(0,0)*mat(1,2)*mat(3,1) + mat(1,0)*mat(0,1)*mat(3,2)
                      - mat(1,0)*mat(0,2)*mat(3,1) - mat(3,0)*mat(0,1)*mat(1,2) + mat(3,0)*mat(0,2)*mat(1,1))*detInv;
       result(3,3) = ( mat(0,0)*mat(1,1)*mat(2,2) - mat(0,0)*mat(1,2)*mat(2,1) - mat(1,0)*mat(0,1)*mat(2,2)
                       +mat(1,0)*mat(0,2)*mat(2,1) + mat(2,0)*mat(0,1)*mat(1,2) - mat(2,0)*mat(0,2)*mat(1,1))*detInv;
       return;
       }
    default :
     utlGlobalException(true, "size too big");
     return;
    }
}


template <class TVector1, class TVector2, class TMatrix>
inline void 
OuterProduct ( const TVector1& v1, const int N1, const TVector2& v2, const int N2, TMatrix& mat )
{
  for ( int i = 0; i < N1; ++i ) 
    for ( int j = 0; j < N2; ++j ) 
      mat(i,j) = v1[i]*v2[j];
}

template <class TVector1, class TMatrix>
inline void 
OuterProduct ( const TVector1& v1, const int N1, TMatrix& mat )
{
  for ( int i = 0; i < N1; ++i ) 
    for ( int j = i; j < N1; ++j ) 
      {
      double tmp = v1[i]*v1[j];
      mat(i,j) = tmp;
      mat(j,i) = tmp;
      }
}

template <class TVector1, class TVector2>
inline double
InnerProduct ( const TVector1& v1, const TVector2& v2, const int N1 )
{
  double result=0.0;
  for ( int i = 0; i < N1; ++i ) 
    result += v1[i]*v2[i];
  return result;
}

template <class TVector1, class TVector2>
inline double
DotProduct ( const TVector1& v1, const TVector2& v2, const int N1 )
{
  return InnerProduct(v1,v2,N1);
}

template <class TVector1>
inline double
SquaredTwoNorm ( const TVector1& v1, const int N1 )
{
  double result=0.0;
  for ( int i = 0; i < N1; ++i ) 
    result += v1[i]*v1[i];
  return result;
}

template <class TVector1, class TVector2, class TVector3>
inline double
CrossProduct ( const TVector1& v1, const TVector2& v2, TVector3& v3 )
{
  v3[0]=v1[1]*v2[2] - v1[2]*v2[1];
  v3[1]=v1[2]*v2[0] - v1[0]*v2[2];
  v3[2]=v1[0]*v2[1] - v1[1]*v2[0];
}

template <class TMatrix, class TVector1, class TVector2>
inline void
ProductMv ( const TMatrix& mat, int rows, int cols, const TVector1& v1, TVector2& v2 )
{
  for ( int i = 0; i < rows; ++i ) 
    {
    v2[i]=0;
    for ( int j = 0; j < cols; ++j ) 
      v2[i] += mat(i,j)*v1[j];
    }
}

template <class TVector1, class TMatrix, class TVector2>
inline void
ProductvM ( const TVector1& v1, int rows, const TMatrix& mat, int cols, TVector2& v2 )
{
  for ( int j = 0; j < cols; ++j ) 
    {
    v2[j]=0;
    for ( int i = 0; i < rows; ++i ) 
      v2[j] += mat(i,j)*v1[i];
    }
}

template <class TMatrix1, class TMatrix2, class TMatrix3>
inline void
ProductMM ( const TMatrix1& mat1, int rows, int cols, const TMatrix2& mat2, int cols2, TMatrix3& mat3)
{
  for ( int i = 0; i < rows; ++i ) 
    {
    for ( int j = 0; j < cols2; ++j ) 
      {
      mat3(i,j)=0;
      for ( int k = 0; k < cols; ++k ) 
        mat3(i,j) += mat1(i,k)*mat2(k,j);
      }
    }
}

}


#endif