{

MODULE_INIT

initCommon ();

#ifdef _MSC_VER

#pragma warning (disable : 4127)

#endif

QC_ASSERT (isNan (NaN));

QC_ASSERT (isNan (inf - inf));

QC_ASSERT (isNan (inf / inf));

#ifdef _MSC_VER

#pragma warning (default : 4127)

#endif

#ifndef _MSC_VER

const Real zero = 0.0;

QC_ASSERT (isNan (zero / zero));

#endif

log_2 = log (2.0);

log_10 = log (10.0);

sqrt_2 = sqrt (2.0);

// Initialize local static variables

lnFactorial (200 - 1);

return true;

streamsize decimals,

bool scientific_arg)

{

ostringstream os;

os << (scientific_arg ? scientific : fixed); os. precision (decimals); os << x;

return os. str ();

{

string s1 (s);

trim (s1);

if (s1. empty ())

throw runtime_error ("Cannot convert an empty string to number");

strLower (s1);

if (s1 == "inf")

return inf;

if (s1 == "-inf")

return -inf;

if ( s1 == "nan"

|| s1 == "-nan"

|| s1 == "?"

)

return NaN;

#if 0

return str2<double> (s1); // slow

#else

// 0

trimSuffix (s1, "e+00");

while (! s1. empty () && s1 [0] == '0')

s1. erase (0, 1);

if (! contains (s1, 'e') && contains (s1, '.'))

while (! s1. empty () && s1. back () == '0')

s1. erase (s1. size () - 1);

if ( s1. empty ()

|| s1 == "."

|| s1 == "-0."

)

return 0.0;

const double res = atof (s1. c_str ());

if (! res)

throw runtime_error ("Cannot convert " + strQuote (s) + " to number");

return res;

#endif

{

if (isNan (a))

return numeric_limits<long>::min ();

return (long) floor (a + 0.5);

: n (324)

, step (30)

, stepInc (0.9)

, logStep (log (step))

, bucket (nullptr)

, maxAbs (NaN)

, logStep_MaxAbs (NaN)

{

ASSERT (n > 0);

ASSERT (step > 2.0 * 10.0);

ASSERT (stepInc <= 0.9);

ASSERT (step * stepInc > 1.0);

bucket = new Real [n];

reset (initMaxAbs);

{

ASSERT (maxAbs >= 0.0);

ASSERT (! isNan (logStep_MaxAbs));

ASSERT (a > 0.0);

const int i = (int) floor (logStep_MaxAbs - log (a) / logStep);

ASSERT (i <= (int) n);

return i;

{

ASSERT (i < n);

// If i is flagged then bucket [i] == 0.0

if (bucket [i] == 0.0 ||

absAddendum2index (fabs (bucket [i])) == (int) i)

return;

// Should happen rarely

const Real X = bucket [i];

bucket [i] = 0.0;

add (X);

{

maxAbs = initMaxAbs;

if (isNan (maxAbs))

logStep_MaxAbs = NaN;

else

{

ASSERT (maxAbs >= 0.0);

logStep_MaxAbs = log (maxAbs) / logStep;

}

{

FOR (size_t, i, n)

bucket [i] = 0.0;

hasPlusInf = false;

hasMinusInf = false;

setMaxAbs (initMaxAbs);

{

ASSERT (! isNan (x));

if (x == inf)

{

hasPlusInf = true;

return;

}

if (x == -inf)

{

hasMinusInf = true;

return;

}

if (x == 0.0)

return;

const Real a = fabs (x);

ASSERT (a != inf);

if (isNan (maxAbs))

{

setMaxAbs (a * step * stepInc);

ASSERT (absAddendum2index (a) == 0);

ASSERT (bucket [0] == 0.0);

}

else if (a > maxAbs)

{

setMaxAbs (a * step * stepInc);

ASSERT (absAddendum2index (a) == 0);

// Going up may violate the bucket compatibility conditions

FOR_REV (size_t, i, n)

stabilize (i);

}

const int i = absAddendum2index (a);

ASSERT (i >= 0);

if (bucket [i] == 0.0)

bucket [i] = x;

else

{

// bucket compatibility conditions

bucket [i] += x;

stabilize ((size_t) i);

}

{

if (hasPlusInf)

if (hasMinusInf)

return NaN;

else

return inf;

else

if (hasMinusInf)

return -inf;

Real s = 0.0;

FOR_REV (size_t, i, n)

s += bucket [i];

return s;

{

maxLnX = -inf;

sum. reset ();

{

if (lnX == -inf)

return;

if (lnX - maxLnX > 3.0) // PAR

{

const Real s = sum. get ();

sum. reset ();

sum. add (s / exp (lnX - maxLnX));

maxLnX = lnX;

}

sum. add (exp (lnX - maxLnX));

{

ASSERT (maxError > 0);

// integral(x+1) < \sum_{x'=x+1}^\infty f(x') < integral(x)

// 0 < integral(x) - integral(x+1) < f(x)

// Sum ??

Real s = 0;

FOR_START (uint, x, start, numeric_limits<uint>::max ())

{

const Real y = f (x);

s += y;

// error = integral (x) - integral (x + 1)

if (y <= maxError)

return s + (integral (x) + integral (x + 1)) / 2;

}

NEVER_CALL;

Real &result) const

{

ASSERT (maxIter > 1);

setNonZero (result, getB (0));

Real c = result;

Real d = 0.0;

FOR (size_t, i, maxIter)

{

const Real a = getA (i + 1);

const Real b = getB (i + 1);

setNonZero (d, b + a * d);

setNonZero (c, b + a / c);

d = 1.0 / d;

const Real delta = c * d;

if (eqReal (delta, 1.0))

return true;

result *= delta;

}

return false;

{

ASSERT (x > 0.0);

static constexpr Real c [6] = { 76.18009172947146,

-86.50532032941677,

24.01409824083091,

-1.231739572450155,

0.001208650973866179,

-0.5395239384953e-5};

const Real a = x + 4.5;

const Real b = (x - 0.5) * log (a) - a;

Real s = 1.000000000190015;

FOR (size_t, i, 6)

s += c [i] / (x + (Real) i);

return b + log (2.5066282746310005 * s);

{

Real x;

Real a;

ContinuedFraction_LnComIncGamma (Real x_arg,

Real a_arg)

: ContinuedFraction ()

, x (x_arg)

, a (a_arg)

{ ASSERT (x > 0.0);

ASSERT (a > 0.0);

}

private:

Real getA (size_t index) const final

{ ASSERT (index > 0);

if (index == 1)

return 1.0;

const Real j = (Real) (index - 1);

return - j * (j - a);

}

Real getB (size_t index) const final

{ return index ? x - a + 1.0 + 2.0 * (Real) (index - 1) : 0.0; }

Real a)

{

ASSERT (x >= 0.0);

ASSERT (a > 0.0);

if (! x)

return lnGamma (a);

const ContinuedFraction_LnComIncGamma cf (x, a);

Real g;

if (! cf. lentz (1000, g)) // PAR

return NaN;

return - x + a * log (x) + log (g);

{

ASSERT (delta >= 0.0);

if (isNan (x))

return x;

if (isProb (x))

return x;

if (nullReal (x, delta))

return 0.0;

if (eqReal (x, 1, delta))

return 1.0;

throw runtime_error ("Not a probability: " + real2str (x));

{

ASSERT (isProb (x));

const streamsize decimals = 4; // PAR

if (x < 0.9)

return real2str (x, decimals);

return "1-" + real2str (1 - x, decimals);

{

constexpr size_t resMax = 200; // Cf. initNumeric()

static Real res [resMax] = {log(1.), log(1.), log(2.), log(6.), log(24.), log(120.), log(720.)};

// 0 1 2 3 4 5 6

static size_t resN = 6; // res[0..resN] are defined

if (n >= resMax)

return lnGamma (n + 1);

while (resN < n)

{

resN++;

res [resN] = res [resN - 1] + log ((Real) resN);

}

return res [n];

{

const Real z = fabs (x);

const Prob t = 1 / (1 + z/2);

const Real ans = t * exp (- sqr (z) - 1.26551223

+ t * (1.00002368

+ t * (0.37409196

+ t * (0.09678418

+ t * (-0.18628806

+ t * (0.27886807

+ t * (-1.13520398

+ t * (1.48851587

+ t * (-0.82215223

+ t * 0.17087277

)

)

)

)

)

)

)

)

);

return 1 - (x >= 0 ? ans : (2 - ans));

{

Real s;

ZetaSeries (Real s_arg,

uint start_arg)

: s (s_arg)

{ ASSERT (s > 1);

start = start_arg;

ASSERT (start >= 1);

}

Real f (uint x) const final

{ return 1 / pow (x, s); }

Real integral (uint x) const final

{ return 1 / (pow (x, s - 1) * (s - 1)); }

uint from)

{

ASSERT (alpha > 1.0);

constexpr size_t size = 1000; // PAR

static Real res [size];

{

static bool init = false;

if (! init)

{

FOR (size_t, i, size)

res [i] = 0;

init = true;

}

}

size_t index = no_index;

if (from == 1 && alpha > 1.0 && alpha < 2.0)

{

index = (size_t) round ((alpha - 1.0) * size);

if (index < size)

if (const Real res_ = res [index])

return res_;

}

ZetaSeries zs (alpha, from);

const Real res_ = zs. get (1e-4); // PAR

if (from == 1 && index < size)

res [index] = res_;

return res_;

uint from,

uint to)

{

Real s = 0;

FOR_REV_END (uint, i, from, to + 1)

s += pow (i, - alpha);

return s;

{

Real s;

ZetaLnSeries (Real s_arg,

uint start_arg)

: s (s_arg)

{ ASSERT (s > 1);

start = start_arg;

ASSERT (start >= 1);

}

Real f (uint x) const final

{ return 1 / pow (x, s) * log (x); }

Real integral (uint x) const final

{ const Real b = 1 - s;

return pow (x, b) * (1 / sqr (b) - log (x) / b);

}

uint from)

{

ASSERT (alpha > 1.0);

constexpr size_t size = 1000; // PAR

static Real res [size];

{

static bool init = false;

if (! init)

{

FOR (size_t, i, size)

res [i] = 0;

init = true;

}

}

size_t index = no_index;

if (from == 1 && alpha > 1 && alpha < 2)

{

index = (size_t) round ((alpha - 1) * size);

if (index < size)

if (const Real res_ = res [index])

return res_;

}

ZetaLnSeries zs (alpha, from);

const Real res_ = zs. get (1e-4); // PAR

if (from == 1 && index < size)

res [index] = res_;

return res_;

uint from,

uint to)

{

Real s = 0;

FOR_REV_END (uint, i, from, to + 1)

s += pow (i, - alpha) * log (i);

return s;

Real weight)

{

ASSERT (! isNan (weight));

if (! weight)

return;

ASSERT (! isNan (x));

const Real x2 = sqr (x);

weightedSum += x * weight;

weightedSum2 += x2 * weight;

weights += weight;

{

weightedSum += other. weightedSum;

weightedSum2 += other. weightedSum2;

weights += other. weights;

{

ASSERT (weights >= other. weights);

weightedSum -= other. weightedSum;

weightedSum2 -= other. weightedSum2;

weights -= other. weights;

{

if (! qc_on)

return;

QC_ASSERT (n >= 0);

QC_ASSERT (s2 >= 0);

QC_IMPLY (s2 > 0, n > 0);

Real stop_arg,

Real binRange_arg)

: start (start_arg)

, stop (stop_arg)

, binRange (binRange_arg)

, bins ((size_t) round (ceil ((stop - start) / binRange)), 0)

{

ASSERT (start < stop);

ASSERT (binRange > 0.0);

ASSERT (bins. size ());

{

os << "bin\tcount" << endl;

FFOR (size_t, i, bins. size ())

os << start + binRange * (Real) i << '\t' << bins [i] << endl;

{

if (! betweenEqual (x, start, stop))

return no_index;

size_t bin = (size_t) round (floor ((x - start) / binRange));

if (bin == bins. size ())

bin--;

ASSERT (bin < bins. size ());

return bin;