File indexing completed on 2024-05-12 17:22:39

 // HMath: C++ high precision math routines
 // Copyright (C) 2004 Ariya Hidayat <ariya.hidayat@gmail.com>
 // Copyright (C) 2007-2008, 2014 @heldercorreia
 // Copyright (C) 2008, 2009 Wolf Lammen
 //
 // This program is free software; you can redistribute it and/or
 // modify it under the terms of the GNU General Public License
 // as published by the Free Software Foundation; either version 2
 // of the License, or (at your option) any later version.
 //
 // This program is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 // GNU General Public License for more details.
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; see the file COPYING.  If not, write to
 // the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 // Boston, MA 02110-1301, USA.
 
 #include "hmath.h"
 
 #include "floatcommon.h"
 #include "floatconst.h"
 #include "floatconvert.h"
 #include "floathmath.h"
 #include "rational.h"
 
 #include <QMap>
 #include <QString>
 #include <QStringList>
 
 #include <cstdio>
 #include <cstdlib>
 #include <cstring>
 #include <sstream>
 
 #define RATIONAL_TOL HNumber("1e-20")
 
 //TODO make this configurable
 #define HMATH_WORKING_PREC (DECPRECISION + 3)
 #define HMATH_EVAL_PREC (HMATH_WORKING_PREC + 2)
 
 //TODO should go into a separate format file
 // default scale for fall back in formatting
 #define HMATH_MAX_SHOWN 20
 #define HMATH_BIN_MAX_SHOWN ((33219*HMATH_MAX_SHOWN)/10000 + 1)
 #define HMATH_OCT_MAX_SHOWN ((11073*HMATH_MAX_SHOWN)/10000 + 1)
 #define HMATH_HEX_MAX_SHOWN ((8305*HMATH_MAX_SHOWN)/10000 + 1)
 
 /*------------------------   Helper routines  -------------------------*/
 
 static void checkfullcancellation(cfloatnum op1, cfloatnum op2,
                                    floatnum r)
 {
     if (float_getlength(op1) != 0
       && float_getlength(op2) != 0
       && float_getlength(r) != 0) {
         // NaN or zero not involved in computation.
         int expr = float_getexponent(r);
         if (float_getexponent(op1) - expr >= HMATH_WORKING_PREC - 1
             || float_getexponent(op2) - expr >= HMATH_WORKING_PREC - 1)
             float_setzero(r);
     }
 }
 
 static char checkAdd(floatnum dest, cfloatnum s1, cfloatnum s2, int digits)
 {
     if (float_add(dest, s1, s2, digits)
       && float_getsign(s1) + float_getsign(s2) == 0)
         checkfullcancellation(s1, s2, dest);
     return float_isnan(dest);
 }
 
 static char checkSub(floatnum dest, cfloatnum s1, cfloatnum s2, int digits)
 {
     if (float_sub(dest, s1, s2, digits)
       && float_getsign(s1) - float_getsign(s2) == 0)
         checkfullcancellation(s1, s2, dest);
     return float_isnan(dest);
 }
 
 static void h_init()
 {
     static bool h_initialized = false;
     if (!h_initialized) {
         h_initialized = true;
         //TODO related to formats, get rid of it.
         float_stdconvert();
     }
 }
 
 static void checkpoleorzero(floatnum result, floatnum x)
 {
     if (float_getlength(result) == 0 || float_getlength(x) == 0)
         return;
     int expx = float_getexponent(x)-HMATH_WORKING_PREC+1;
     int expr = float_getexponent(result);
     if (expr <= expx)
         float_setzero(result);
     else if (expr >= -expx) {
         float_setnan(result);
         float_seterror(EvalUnstable);
     }
 }
 
 static char roundResult(floatnum x)
 {
     if (float_isnan(x)) // Avoids setting float_error.
         return 0;
     return float_round(x, x, HMATH_WORKING_PREC, TONEAREST);
 }
 
 /*---------------------------   HNumberPrivate   --------------------*/
 
 class HNumberPrivate
 {
 public:
     HNumberPrivate();
     ~HNumberPrivate();
     //TODO make this a variant
     floatstruct fnum;
     Error error;
 };
 
 HNumberPrivate::HNumberPrivate()
   : error(Success)
 {
     h_init();
     float_create(&fnum);
 }
 
 HNumberPrivate::~HNumberPrivate()
 {
     float_free(&fnum);
 }
 
 typedef char (*Float1ArgND)(floatnum x);
 typedef char (*Float1Arg)(floatnum x, int digits);
 typedef char (*Float2ArgsND)(floatnum result, cfloatnum p1, cfloatnum p2);
 typedef char (*Float2Args)(floatnum result, cfloatnum p1, cfloatnum p2, int digits);
 
 static Error checkNaNParam(const HNumberPrivate& v1,
                             const HNumberPrivate* v2 = 0)
 {
     if (!float_isnan(&v1.fnum) && (!v2 || !float_isnan(&v2->fnum)))
         return Success;
     Error error = v1.error;
     if (error == Success && v2)
         error = v2->error;
     return error == 0? NoOperand : error;
 }
 
 void roundSetError(HNumberPrivate* dest)
 {
     dest->error = float_geterror();
     floatnum dfnum = &dest->fnum;
     if (dest->error != Success)
         float_setnan(dfnum);
     if (!float_isnan(dfnum))
         float_round(dfnum, dfnum, HMATH_WORKING_PREC, TONEAREST);
 }
 
 void call2Args(HNumberPrivate* dest, HNumberPrivate* n1, HNumberPrivate* n2, Float2Args func)
 {
     dest->error = checkNaNParam(*n1, n2);
     if (dest->error == Success) {
         floatnum dfnum = &dest->fnum;
         func(dfnum, &n1->fnum, &n2->fnum, HMATH_EVAL_PREC);
         roundSetError(dest);
     }
 }
 
 void call2ArgsND(HNumberPrivate* dest, HNumberPrivate* n1, HNumberPrivate* n2, Float2ArgsND func)
 {
     dest->error = checkNaNParam(*n1, n2);
     if (dest->error == Success) {
         floatnum dfnum = &dest->fnum;
         func(dfnum, &n1->fnum, &n2->fnum);
         roundSetError(dest);
     }
 }
 
 void call1Arg(HNumberPrivate* dest, HNumberPrivate* n, Float1Arg func)
 {
     dest->error = checkNaNParam(*n);
     if (dest->error == Success) {
         floatnum dfnum = &dest->fnum;
         float_copy(dfnum, &n->fnum, HMATH_EVAL_PREC);
         func(dfnum, HMATH_EVAL_PREC);
         roundSetError(dest);
     }
 }
 
 void call1ArgPoleCheck(HNumberPrivate* dest, HNumberPrivate* n, Float1Arg func)
 {
     dest->error = checkNaNParam(*n);
     if (dest->error == Success) {
         floatnum dfnum = &dest->fnum;
         float_copy(dfnum, &n->fnum, HMATH_EVAL_PREC);
         if (func(dfnum, HMATH_EVAL_PREC))
           checkpoleorzero(dfnum, &n->fnum);
         roundSetError(dest);
     }
 }
 
 void call1ArgND(HNumberPrivate* dest, HNumberPrivate* n, Float1ArgND func)
 {
     dest->error = checkNaNParam(*n);
     if (dest->error == Success)
     {
         floatnum dfnum = &dest->fnum;
         float_copy(dfnum, &n->fnum, HMATH_EVAL_PREC);
         func(dfnum);
         roundSetError(dest);
     }
 }
 
 static char modwrap(floatnum result, cfloatnum p1, cfloatnum p2, int /*digits*/)
 {
     floatstruct tmp;
     float_create(&tmp);
     char ok = float_divmod(&tmp, result, p1, p2, INTQUOT);
     float_free(&tmp);
     return ok;
 }
 
 char idivwrap(floatnum result, cfloatnum p1, cfloatnum p2)
 {
     floatstruct tmp;
     int save = float_setprecision(DECPRECISION);
     float_create(&tmp);
     char ok = float_divmod(result, &tmp, p1, p2, INTQUOT);
     float_free(&tmp);
     float_setprecision(save);
     return ok;
 }
 
 /*--------------------------   HNumber   -------------------*/
 
 /**
  * Creates a new number.
  */
 HNumber::HNumber() : d(new HNumberPrivate)
 {
 }
 
 /**
  * Copies from another number.
  */
 HNumber::HNumber(const HNumber& hn) : d(new HNumberPrivate)
 {
     operator=(hn);
 }
 
 /**
  * Creates a new number from an integer value.
  */
 HNumber::HNumber(int i) : d(new HNumberPrivate)
 {
     float_setinteger(&d->fnum, i);
 }
 
 /**
  * Creates a new number from a string.
  */
 HNumber::HNumber(const char* str) : d(new HNumberPrivate)
 {
     t_itokens tokens;
 
     if ((d->error = parse(&tokens, &str)) == Success && *str == 0)
         d->error = float_in(&d->fnum, &tokens);
     float_geterror();
 }
 
 HNumber::HNumber(const QJsonObject& json)
     : d(new HNumberPrivate)
 {
     *this = deSerialize(json);
 }
 
 /**
  * Destroys this number.
  */
 HNumber::~HNumber()
 {
     delete d;
 }
 
 /**
  * Returns the error code kept with a NaN.
  */
 Error HNumber::error() const
 {
     return d->error;
 }
 
 /**
  * Returns true if this number is Not a Number (NaN).
  */
 bool HNumber::isNan() const
 {
     return float_isnan(&d->fnum) != 0;
 }
 
 /**
  * Returns true if this number is zero.
  */
 bool HNumber::isZero() const
 {
     return float_iszero(&d->fnum) != 0;
 }
 
 /**
  * Returns true if this number is more than zero.
  */
 bool HNumber::isPositive() const
 {
     return float_getsign(&d->fnum) > 0;
 }
 
 /**
  * Returns true if this number is less than zero.
  */
 bool HNumber::isNegative() const
 {
     return float_getsign(&d->fnum) < 0;
 }
 
 /**
  * Returns true if this number is integer.
  */
 bool HNumber::isInteger() const
 {
     return float_isinteger(&d->fnum) != 0;
 }
 
 void HNumber::serialize(QJsonObject& json) const
 {
     json["value"] = HMath::format(*this, Format::Fixed() + Format::Precision(DECPRECISION));
 }
 
 HNumber HNumber::deSerialize(const QJsonObject& json)
 {
     HNumber result;
     if (json.contains("value")) {
         QString str = json["value"].toString();
         str.replace(",", ".");
         result = HNumber(str.toLatin1().constData());
     }
     return result;
 }
 
 /**
  * Returns the number as an int.
  * It is meant to convert small (integer) numbers only and no
  * checking is done whatsoever.
  */
 int HNumber::toInt() const
 {
     return float_asinteger(&d->fnum);
 }
 
 /**
  * Assigns from another number.
  */
 HNumber& HNumber::operator=(const HNumber& hn)
 {
     d->error = hn.error();
 
     float_copy(&d->fnum, &hn.d->fnum, EXACT);
 
     return *this;
 }
 
 
 /**
  * Adds another number.
  */
 HNumber HNumber::operator+(const HNumber& num) const
 {
     if (this->isZero())
         return num;
     if (num.isZero())
         return *this;
     HNumber result;
     call2Args(result.d, d, num.d, checkAdd);
     return result;
 }
 
 /**
  * Adds another number.
  */
 HNumber& HNumber::operator+=(const HNumber& num)
 {
     return operator=(*this + num);
 }
 
 /**
  * Subtract from another number.
  */
 HNumber operator-(const HNumber& n1, const HNumber& n2)
 {
     HNumber result;
     call2Args(result.d, n1.d, n2.d, checkSub);
     return result;
 }
 
 /**
  * Subtract from another number.
  */
 HNumber& HNumber::operator-=(const HNumber& num)
 {
     return operator=(*this - num);
 }
 
 /**
  * Multiplies with another number.
  */
 HNumber HNumber::operator*(const HNumber& num) const
 {
     HNumber result;
     call2Args(result.d, d, num.d, float_mul);
     return result;
 }
 
 /**
  * Multiplies with another number.
  */
 HNumber& HNumber::operator*=(const HNumber& num)
 {
     return operator=(*this * num);
 }
 
 /**
  * Divides with another number.
  */
 HNumber HNumber::operator/(const HNumber& num) const
 {
     HNumber result;
     call2Args(result.d, d, num.d, float_div);
     return result;
 }
 
 /**
  * Divides with another number.
  */
 HNumber& HNumber::operator/=(const HNumber& num)
 {
     return operator=(*this / num);
 }
 
 /**
  * Modulo (rest of integer division)
  */
 HNumber HNumber::operator%(const HNumber& num) const
 {
     HNumber result;
     call2Args(result.d, d, num.d, modwrap);
     return result;
 }
 
 /**
  * Returns -1, 0, 1 if *this is less than, equal to, or more than other.
  */
 int HNumber::compare(const HNumber& other) const
 {
     int result = float_relcmp(&d->fnum, &other.d->fnum, HMATH_EVAL_PREC-1);
     float_geterror(); // clears error, if one operand was a NaN
     return result;
 }
 
 /**
  * Returns true if l is greater than r.
  */
 bool operator>(const HNumber& l, const HNumber& r)
 {
     return l.compare(r) > 0;
 }
 
 /**
  * Returns true if l is less than r.
  */
 bool operator<(const HNumber& l, const HNumber& r)
 {
     return l.compare(r) < 0;
 }
 
 /**
  * Returns true if l is greater than or equal to r.
  */
 bool operator>=(const HNumber& l, const HNumber& r)
 {
     return l.compare(r) >= 0;
 }
 
 /**
  * Returns true if l is less than or equal to r.
  */
 bool operator<=(const HNumber& l, const HNumber& r)
 {
     return l.compare(r) <= 0;
 }
 
 /**
  * Returns true if l is equal to r.
  */
 bool operator==(const HNumber& l, const HNumber& r)
 {
     return l.compare(r) == 0;
 }
 
 /**
  * Returns true if l is not equal to r.
  */
 bool operator!=(const HNumber& l, const HNumber& r)
 {
     return l.compare(r) != 0;
 }
 
 /**
  * Bitwise ANDs the integral parts of both operands.
  * Yields NaN, if any operand exeeds the logic range
  */
 HNumber HNumber::operator&(const HNumber& num) const
 {
     HNumber result;
     call2ArgsND(result.d, d, num.d, float_and);
     return result;
 }
 
 /**
  * Bitwise ANDs the integral parts of both operands.
  * Yields NaN, if any operand exeeds the logic range
  */
 HNumber& HNumber::operator&=(const HNumber& num)
 {
     return operator=(*this & num);
 }
 
 /**
  * Bitwise ORs the integral parts of both operands.
  * Yields NaN, if any operand exeeds the logic range
  */
 HNumber HNumber::operator|(const HNumber& num) const
 {
     HNumber result;
     call2ArgsND(result.d, d, num.d, float_or);
     return result;
 }
 
 /**
  * Bitwise ORs the integral parts of both operands.
  * Yields NaN, if any operand exeeds the logic range
  */
 HNumber& HNumber::operator|=(const HNumber& num)
 {
     return operator=(*this | num);
 }
 
 /**
  * Bitwise XORs the integral parts of both operands.
  * Yields NaN, if any operand exeeds the logic range
  */
 HNumber HNumber::operator^(const HNumber& num) const
 {
     HNumber result;
     call2ArgsND(result.d, d, num.d, float_xor);
     return result;
 }
 
 /**
  * Bitwise XORs the integral parts of both operands.
  * Yields NaN, if any operand exeeds the logic range
  */
 HNumber& HNumber::operator^=(const HNumber& num)
 {
     return operator=(*this ^ num);
 }
 
 /**
  * Bitwise NOTs the integral part *this.
  * Yields NaN, if *this exeeds the logic range
  */
 HNumber HNumber::operator~() const
 {
     HNumber result;
     call1ArgND(result.d, d, float_not);
     return result;
 }
 
 /**
  * Changes the sign.
  */
 HNumber operator-(const HNumber& x)
 {
     HNumber result;
     call1ArgND(result.d, x.d, float_neg);
     return result;
 }
 
 /**
  * Shifts the integral part of <*this> to the left by
  * the parameters value's bits. Zeros are shifted in
  * to the right, shifted out bits are dropped.
  * Yields NaN, if the operand exeeds the logic range,
  * or the shift count is not a non-negative integer.
  */
 HNumber HNumber::operator<<(const HNumber& num) const
 {
     HNumber result;
     call2ArgsND(result.d, d, num.d, float_shl);
     return result;
 }
 
 /**
  * Shifts the integral part of <*this> to the right by
  * the parameters value's bits. The most significand
  * bit is duplicated to the left, shifted out bits are dropped
  * (signed or arithmethic shift right).
  * Yields NaN, if the operand exeeds the logic range,
  * or the shift count is not a non-negative integer.
  */
 HNumber HNumber::operator>>(const HNumber& num) const
 {
     HNumber result;
     call2ArgsND(result.d, d, num.d, float_shr);
     return result;
 }
 
 HNumber::Format::Format()
     : base(Base::Null)
     , radixChar(RadixChar::Null)
     , mode(Mode::Null)
     , precision(PrecisionNull)
 {
 }
 
 HNumber::Format::Format(const HNumber::Format& other)
     : base(other.base)
     , radixChar(other.radixChar)
     , mode(other.mode)
     , precision(other.precision)
 {
 }
 
 HNumber::Format HNumber::Format::operator+(const HNumber::Format& other) const
 {
     Format result;
     result.base = (this->base != Base::Null) ? this->base : other.base;
     result.radixChar = (this->radixChar != RadixChar::Null) ? this->radixChar : other.radixChar;
     result.mode = (this->mode != Mode::Null) ? this->mode : other.mode;
     result.precision = (this->precision != PrecisionNull) ? this->precision : other.precision;
     return result;
 }
 
 HNumber::Format HNumber::Format::Binary()
 {
     Format result;
     result.base = Base::Binary;
     return result;
 }
 
 HNumber::Format HNumber::Format::Octal()
 {
     Format result;
     result.base = Base::Octal;
     return result;
 }
 
 HNumber::Format HNumber::Format::Decimal()
 {
     Format result;
     result.base = Base::Decimal;
     return result;
 }
 
 HNumber::Format HNumber::Format::Hexadecimal()
 {
     Format result;
     result.base = Base::Hexadecimal;
     return result;
 }
 
 HNumber::Format HNumber::Format::Precision(int prec)
 {
     Format result;
     result.precision = prec;
     return result;
 }
 
 HNumber::Format HNumber::Format::Point()
 {
     Format result;
     result.radixChar = RadixChar::Point;
     return result;
 }
 
 HNumber::Format HNumber::Format::Comma()
 {
     Format result;
     result.radixChar = RadixChar::Comma;
     return result;
 }
 
 HNumber::Format HNumber::Format::General()
 {
     Format result;
     result.mode = Mode::General;
     return result;
 }
 
 HNumber::Format HNumber::Format::Fixed()
 {
     Format result;
     result.mode = Mode::Fixed;
     return result;
 }
 
 HNumber::Format HNumber::Format::Scientific()
 {
     Format result;
     result.mode = Mode::Scientific;
     return result;
 }
 
 HNumber::Format HNumber::Format::Engineering()
 {
     Format result;
     result.mode = Mode::Engineering;
     return result;
 }
 
 namespace {
 
 char* _doFormat(cfloatnum x, signed char base, signed char expbase, char outmode, int prec, unsigned flags)
 {
     t_otokens tokens;
     char intbuf[BINPRECISION+1];
     char fracbuf[BINPRECISION+1];
     int sz = 0;
     char* str = nullptr;
     switch (base) {
     case 2:
         sz = BINPRECISION+1;
         break;
     case 8:
         sz = OCTPRECISION+1;
         break;
     case 10:
         sz = DECPRECISION+1;
         break;
     case 16:
         sz = HEXPRECISION+1;
         break;
     }
     tokens.intpart.sz = sz;
     tokens.intpart.buf = intbuf;
     tokens.fracpart.sz = sz;
     tokens.fracpart.buf = fracbuf;
 
     floatstruct tmp;
     float_create(&tmp);
     float_copy(&tmp, x, DECPRECISION + 2);
     if (float_out(&tokens, &tmp, prec, base, outmode) == Success)
     {
         sz = cattokens(nullptr, -1, &tokens, expbase, flags);
         str = (char*)malloc(sz);
         cattokens(str, sz, &tokens, expbase, flags);
     }
     float_free(&tmp);
     return str;
 }
 
 /**
  * Formats the given number as string, using specified decimal digits.
  * Note that the returned string must be freed.
  */
 char* formatFixed(cfloatnum x, int prec, int base = 10)
 {
     unsigned flags = IO_FLAG_SUPPRESS_PLUS + IO_FLAG_SUPPRESS_DOT + IO_FLAG_SUPPRESS_EXPZERO;
     if (base != 10)
         flags += IO_FLAG_SHOW_BASE + IO_FLAG_SHOW_EXPBASE;
     if (prec < 0) {
         prec = HMATH_MAX_SHOWN;
     }
     flags |= IO_FLAG_SUPPRESS_TRL_ZERO;
     char* result = _doFormat(x, base, base, IO_MODE_FIXPOINT, prec, flags);
     return result ? result : _doFormat(x, base, base, IO_MODE_SCIENTIFIC, HMATH_MAX_SHOWN, flags);
 }
 
 /**
  * Formats the given number as string, in scientific format.
  * Note that the returned string must be freed.
  */
 char* formatScientific(cfloatnum x, int prec, int base = 10)
 {
     unsigned flags = IO_FLAG_SUPPRESS_PLUS + IO_FLAG_SUPPRESS_DOT + IO_FLAG_SUPPRESS_EXPPLUS;
     if (base != 10)
         flags += IO_FLAG_SHOW_BASE + IO_FLAG_SHOW_EXPBASE;
     if (prec < 0) {
         flags |= IO_FLAG_SUPPRESS_TRL_ZERO;
         prec = HMATH_MAX_SHOWN;
     }
     return _doFormat(x, base, base, IO_MODE_SCIENTIFIC, prec, flags);
 }
 
 /**
  * Formats the given number as string, in engineering notation.
  * Note that the returned string must be freed.
  */
 char* formatEngineering(cfloatnum x, int prec, int base = 10)
 {
     unsigned flags = IO_FLAG_SUPPRESS_PLUS + IO_FLAG_SUPPRESS_EXPPLUS;
     if (base != 10)
         flags += IO_FLAG_SHOW_BASE + IO_FLAG_SHOW_EXPBASE;
     if (prec <= 1) {
         flags |= IO_FLAG_SUPPRESS_TRL_ZERO + IO_FLAG_SUPPRESS_DOT;
         prec = HMATH_MAX_SHOWN;
     }
     return _doFormat(x, base, base, IO_MODE_ENG, prec, flags);
 }
 
 /**
  * Formats the given number as string, using specified decimal digits.
  * Note that the returned string must be freed.
  */
 char* formatGeneral(cfloatnum x, int prec, int base = 10)
 {
     // find the exponent and the factor
     int expd = float_getexponent(x);
 
     char* str;
     if (expd > 5)
         str = formatScientific(x, prec, base);
     else if (expd < -4)
         str = formatScientific(x, prec, base);
     else if ((expd < 0) && (prec >= 0) && (expd < -prec))
         str = formatScientific(x, prec, base);
     else
         str = formatFixed(x, prec, base);
 
     return str;
 }
 
 
 } /* unnamed namespace */
 
 /**
  * Formats the given number as string, using specified decimal digits.
  */
 QString HMath::format(const HNumber& hn, HNumber::Format format)
 {
     char* rs = 0;
 
     if (format.precision < 0)  // This includes PrecisionNull
         format.precision = -1;
 
     int base;
     switch (format.base) {
     case HNumber::Format::Base::Binary:
         base = 2;
         break;
     case HNumber::Format::Base::Octal:
         base = 8;
         break;
     case HNumber::Format::Base::Hexadecimal:
         base = 16;
         break;
     case HNumber::Format::Base::Decimal:
     case HNumber::Format::Base::Null:
     default:
         base = 10;
         break;
     }
 
     switch (format.mode) {
     case HNumber::Format::Mode::Fixed:
         rs = formatFixed(&hn.d->fnum, format.precision, base);
         break;
     case HNumber::Format::Mode::Scientific:
         rs = formatScientific(&hn.d->fnum, format.precision, base) ;
         break;
     case HNumber::Format::Mode::Engineering:
         rs = formatEngineering(&hn.d->fnum, format.precision, base);
         break;
     case HNumber::Format::Mode::General:
     default:
         rs = formatGeneral(&hn.d->fnum, format.precision, base);
     }
 
     QString result(rs);
     free(rs);
     return result;
 }
 
 /**
  * Converts radians to degrees.
  */
 HNumber HMath::rad2deg(const HNumber& angle)
 {
     return angle * (HNumber(180) / HMath::pi());
 }
 
 /**
  * Converts degrees to radians.
  */
 HNumber HMath::deg2rad(const HNumber& angle)
 {
     return angle * (HMath::pi() / HNumber(180));
 }
 
 /**
  * Converts radians to gons.
  */
 HNumber HMath::rad2gon(const HNumber& angle)
 {
     return angle * (HNumber(200) / HMath::pi());
 }
 
 /**
  * Converts gons to radians.
  */
 HNumber HMath::gon2rad(const HNumber& angle)
 {
     return angle * (HMath::pi() / HNumber(200));
 }
 
 /**
  * Returns the constant e (Euler's number).
  */
 HNumber HMath::e()
 {
     HNumber value;
     float_copy(&value.d->fnum, &cExp, HMATH_EVAL_PREC);
     return value;
 }
 
 /**
  * Returns the constant Pi.
  */
 HNumber HMath::pi()
 {
     HNumber value;
     float_copy(&value.d->fnum, &cPi, HMATH_EVAL_PREC);
     return value;
 }
 
 /**
  * Returns the constant Phi (golden number).
  */
 HNumber HMath::phi()
 {
     HNumber value;
     float_copy(&value.d->fnum, &cPhi, HMATH_EVAL_PREC);
     return value;
 }
 
 /**
  * Returns a NaN (Not a Number) with error set to
  * passed parameter.
  */
 HNumber HMath::nan(Error error)
 {
     HNumber result;
     result.d->error = error;
     return result;
 }
 
 /**
  * Returns the maximum of two numbers.
  */
 HNumber HMath::max(const HNumber& n1, const HNumber& n2)
 {
     switch (float_cmp(&n1.d->fnum, &n2.d->fnum))
     {
         case 0:
         case 1:  return n1;
         case -1: return n2;
         default: return HMath::nan(checkNaNParam(*n1.d, n2.d));
     }
 }
 
 /**
  * Returns the minimum of two numbers.
  */
 HNumber HMath::min(const HNumber& n1, const HNumber& n2)
 {
     switch (float_cmp(&n1.d->fnum, &n2.d->fnum))
     {
         case 0:
         case 1:  return n2;
         case -1: return n1;
         default: return HMath::nan(checkNaNParam(*n1.d, n2.d));
     }
 }
 
 /**
  * Returns the absolute value of n.
  */
 HNumber HMath::abs(const HNumber& n)
 {
     HNumber result;
     call1ArgND(result.d, n.d, float_abs);
     return result;
 }
 
 /**
  * Rounds n to the specified decimal digits.
  */
 HNumber HMath::round(const HNumber& n, int prec)
 {
     if (n.isNan())
         return HMath::nan(checkNaNParam(*n.d));
     HNumber result(n);
     floatnum rnum = &result.d->fnum;
     int exp = float_getexponent(rnum);
 
     // Avoid exponent overflow later.
     if (prec > HMATH_WORKING_PREC && exp > 0)
     prec = HMATH_WORKING_PREC;
     if (prec < 0 && -exp-1 > prec)
         float_setzero(rnum);
     else if (exp + prec < HMATH_WORKING_PREC) {
         float_addexp(rnum, prec);
         float_roundtoint(rnum, TONEAREST);
         float_addexp(rnum, -prec);
     }
     return result;
 }
 
 /**
  * Truncates n to the specified decimal digits.
  */
 HNumber HMath::trunc(const HNumber& n, int prec)
 {
     if (n.isNan())
         return HMath::nan(checkNaNParam(*n.d));
     HNumber result(n);
     floatnum rnum = &result.d->fnum;
     int exp = float_getexponent(rnum);
     // Avoid exponent overflow later on.
     if (prec > HMATH_WORKING_PREC && exp > 0)
         prec = HMATH_WORKING_PREC;
     if (prec < 0 && -exp-1 > prec)
         float_setzero(rnum);
     else if (exp + prec < HMATH_WORKING_PREC) {
         float_addexp(rnum, prec);
         float_roundtoint(rnum, TOZERO);
         float_addexp(rnum, -prec);
     }
     return result;
 }
 
 /**
  * Returns the integer part of n.
  */
 HNumber HMath::integer(const HNumber& n)
 {
     HNumber result;
     call1ArgND(result.d, n.d, float_int);
     return result;
 }
 
 /**
  * Returns the fraction part of n.
  */
 HNumber HMath::frac(const HNumber& n)
 {
     HNumber result;
     call1ArgND(result.d, n.d, float_frac);
     return result;
 }
 
 #define CHECK_NAN \
     if (n.isNan()) \
         return HMath::nan(checkNaNParam(*n.d));
 
 #define RETURN_IF_NEAR_INT \
     HNumber nearest_int(n); \
     float_roundtoint(&nearest_int.d->fnum, TONEAREST); \
     /* Note: float_relcmp doesn't work here, because it's doesn't check the relative */ \
     /* tolerance if exponents are not the same. */ \
     /* FIXME: Put this value as parameter. */ \
     /* Kudos to the guy who can figure out why we need such a small tolerance here. */ \
     /* 1e-70 does not work, but 1e-10000 does. */ \
     if (HMath::abs(n - nearest_int) < HNumber("1e-100") * HMath::abs(n + nearest_int)) /* FIXME: Make configurable. */ \
         return nearest_int; \
 
 /**
  * Returns the floor of n.
  */
 HNumber HMath::floor(const HNumber& n)
 {
     CHECK_NAN;
     RETURN_IF_NEAR_INT;
     // Actual rounding, if needed.
     HNumber r(n);
     float_roundtoint(&r.d->fnum, TOMINUSINFINITY);
     return r;
 }
 
 /**
  * Returns the ceiling of n.
  */
 HNumber HMath::ceil(const HNumber& n)
 {
     CHECK_NAN;
     RETURN_IF_NEAR_INT;
     // Actual rounding, if needed.
     HNumber r(n);
     float_roundtoint(&r.d->fnum, TOPLUSINFINITY);
     return r;
 }
 
 /**
  * Returns the greatest common divisor of n1 and n2.
  */
 HNumber HMath::gcd(const HNumber& n1, const HNumber& n2)
 {
     if (!n1.isInteger() || !n2.isInteger()) {
         Error error = checkNaNParam(*n1.d, n2.d);
         if (error != Success)
             return HMath::nan(error);
         return HMath::nan(TypeMismatch);
     }
 
     HNumber a = abs(n1);
     HNumber b = abs(n2);
 
     if (a == 0)
         return b;
     if (b == 0)
         return a;
 
     // Run Euclidean algorithm.
     while (true) {
         a = a % b;
 
         if (a == 0)
             return b;
 
         b = b % a;
 
         if (b == 0)
             return a;
     }
 }
 
 /**
  * Performs an integer divide.
  */
 HNumber HMath::idiv(const HNumber& dividend, const HNumber& divisor)
 {
     HNumber result;
     call2ArgsND(result.d, dividend.d, divisor.d, idivwrap);
     if (result.error() == TooExpensive)
         result.d->error = Overflow;
     return result;
 }
 
 /**
  * Returns the square root of n. If n is negative, returns NaN.
  */
 HNumber HMath::sqrt(const HNumber& n)
 {
     HNumber result;
     call1Arg(result.d, n.d, float_sqrt);
     return result;
 }
 
 /**
  * Returns the cube root of n.
  */
 HNumber HMath::cbrt(const HNumber& n)
 {
     if (n.isNan())
         return HMath::nan(checkNaNParam(*n.d));
     if (n.isZero())
         return n;
     HNumber r;
     floatnum rnum = &r.d->fnum;
 
     // iterations to approximate result
     // X[i+1] = (2/3)X[i] + n / (3 * X[i]^2))
     // initial guess = sqrt(n)
     // r = X[i], q = X[i+1], a = n
 
     floatstruct a, q;
     float_create(&a);
     float_create(&q);
     float_copy(&a, &n.d->fnum, HMATH_EVAL_PREC);
     signed char sign = float_getsign(&a);
     float_abs(&a);
     int expn = float_getexponent(&a);
     float_setexponent(&a, expn % 3);
     expn /= 3;
 
     int digits = 0;
     float_copy(&q, &a, 2);
     float_sqrt(&q, 2);
 
     while (digits < HMATH_EVAL_PREC/2 + 3) {
         digits = 4 * digits + 2;
         if (digits > HMATH_EVAL_PREC+2)
             digits = HMATH_EVAL_PREC+2;
         float_move(rnum, &q);
         float_mul(&q, rnum, rnum, digits);
         float_div(&q, &a, &q, digits);
         float_add(&q, &q, rnum, digits);
         float_add(&q, &q, rnum, digits);
         float_div(&q, &q, &c3, digits);
         float_sub(rnum, rnum, &q, 3);
         if (!float_iszero(rnum))
             digits = float_getexponent(&q) - float_getexponent(rnum);
     }
     float_move(rnum, &q);
     float_free(&q);
     float_free(&a);
     float_setsign(rnum, sign);
     float_addexp(rnum, expn);
 
     roundResult(&r.d->fnum);
     return r;
 }
 
 /**
  * Raises n1 to an integer n.
  */
 HNumber HMath::raise(const HNumber& n1, int n)
 {
     HNumber r;
     float_raisei(&r.d->fnum, &n1.d->fnum, n, HMATH_EVAL_PREC);
     roundSetError(r.d);
     return r;
 }
 
 /**
  * Raises n1 to n2.
  */
 HNumber HMath::raise(const HNumber& n1, const HNumber& n2)
 {
     HNumber result;
     HNumber temp = n1;
     Rational exp;
     bool change_sgn=false;
     if (n1.isNegative() && !n2.isInteger()){
         //Try to convert n2 to a Rational. If n2 is not rational, return NaN.
         // For negative bases only allow odd denominators.
         exp = Rational(n2);
         if (abs(exp.toHNumber() - n2) >= RATIONAL_TOL
             || (n1.isNegative() && exp.denominator()%2 == 0))
             return HMath::nan(OutOfDomain);
         if (n1.isNegative() && !n2.isInteger()) {
             temp = -temp;
             change_sgn = true;
         }
     }
 
     call2Args(result.d, temp.d, n2.d, float_raise);
     result *= (change_sgn) ? -1 : 1;
     return result;
 }
 
 /**
  * Returns e raised to x.
  */
 HNumber HMath::exp(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_exp);
     return result;
 };
 
 /**
  * Returns the natural logarithm of x.
  * If x is non positive, returns NaN.
  */
 HNumber HMath::ln(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_ln);
     return result;
 
 }
 
 /**
  * Returns the common logarithm of x.
  * If x is non positive, returns NaN.
  */
 HNumber HMath::lg(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_lg);
     return result;
 }
 
 /**
  * Returns the binary logarithm of x.
  * If x is non positive, returns NaN.
  */
 HNumber HMath::lb(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_lb);
     return result;
 }
 
 /**
  * Returns the logarithm of x to base.
  * If x is non positive, returns NaN.
  */
 HNumber HMath::log(const HNumber& base, const HNumber& x)
 {
     return lg(x) / lg(base);
 }
 
 /**
  * Returns the sine of x. Note that x must be in radians.
  */
 HNumber HMath::sin(const HNumber& x)
 {
     HNumber result;
     call1ArgPoleCheck(result.d, x.d, float_sin);
     return result;
 }
 
 /**
  * Returns the cosine of x. Note that x must be in radians.
  */
 HNumber HMath::cos(const HNumber& x)
 {
     HNumber result;
     call1ArgPoleCheck(result.d, x.d, float_cos);
     return result;
 }
 
 /**
  * Returns the tangent of x. Note that x must be in radians.
  */
 HNumber HMath::tan(const HNumber& x)
 {
     HNumber result;
     call1ArgPoleCheck(result.d, x.d, float_tan);
     return result;
 }
 
 /**
  * Returns the cotangent of x. Note that x must be in radians.
  */
 HNumber HMath::cot(const HNumber& x)
 {
     return cos(x) / sin(x);
 }
 
 /**
  * Returns the secant of x. Note that x must be in radians.
  */
 HNumber HMath::sec(const HNumber& x)
 {
     return HNumber(1) / cos(x);
 }
 
 /**
  * Returns the cosecant of x. Note that x must be in radians.
  */
 HNumber HMath::csc(const HNumber& x)
 {
     return HNumber(1) / sin(x);
 }
 
 /**
  * Returns the arc tangent of x.
  */
 HNumber HMath::arctan(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_arctan);
     return result;
 }
 
 /**
  * Returns the angle formed by the vector (x, y) and the x-axis.
  */
 HNumber HMath::arctan2(const HNumber& x, const HNumber& y)
 {
     if (y.isZero()) {
         if (x.isNegative())
             return HMath::pi();
         if (x.isZero())
             return HMath::nan(OutOfDomain);
         return HNumber(0);
     } else {
         HNumber phi = HMath::arctan(HMath::abs(y / x));
         if (x.isNegative())
             return (HMath::pi() - phi) * HMath::sgn(y);
         if (x.isZero())
             return HMath::pi()/HNumber(2)*HMath::sgn(y);
         return phi * HMath::sgn(y);
     }
 }
 
 /**
  * Returns the arc sine of x.
  */
 HNumber HMath::arcsin(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_arcsin);
     return result;
 };
 
 /**
  * Returns the arc cosine of x.
  */
 HNumber HMath::arccos(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_arccos);
     return result;
 };
 
 /**
  * Returns the hyperbolic sine of x.
  */
 HNumber HMath::sinh(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_sinh);
     return result;
 }
 
 /**
  * Returns the hyperbolic cosine of x.
  */
 HNumber HMath::cosh(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_cosh);
     return result;
 }
 
 /**
  * Returns the hyperbolic tangent of x.
  */
 HNumber HMath::tanh(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_tanh);
     return result;
 }
 
 /**
  * Returns the area hyperbolic sine of x.
  */
 HNumber HMath::arsinh(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_arsinh);
     return result;
 }
 
 /**
  * Returns the area hyperbolic cosine of x.
  */
 HNumber HMath::arcosh(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_arcosh);
     return result;
 }
 
 /**
  * Returns the area hyperbolic tangent of x.
  */
 HNumber HMath::artanh(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_artanh);
     return result;
 }
 
 /**
  * Returns the Gamma function.
  */
 HNumber HMath::gamma(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_gamma);
     return result;
 }
 
 /**
  * Returns ln(abs(Gamma(x))).
  */
 HNumber HMath::lnGamma(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_lngamma);
     return result;
 }
 
 /**
  * Returns signum x.
  */
 HNumber HMath::sgn(const HNumber& x)
 {
     if (x.isNan())
         return HMath::nan(checkNaNParam(*x.d));
     return float_getsign(&x.d->fnum);
 }
 
 /**
  * Returns the binomial coefficient of n and r.
  * Is any of n and r negative or a non-integer,
  * 1/(n+1)*B(r+1, n-r+1)) is returned, where
  * B(x,y) is the complete Beta function
  */
 HNumber HMath::nCr(const HNumber& n, const HNumber& r)
 {
     Error error = checkNaNParam(*n.d, r.d);
     if (error != Success)
         return HMath::nan(error);
 
     if (r.isInteger() && r < 0) return HNumber(0);
 
     // use symmetry nCr(n, r) == nCr(n, n-r) to find r1 such
     // 2*r1 <= n and nCr(n, r) == nCr(n, r1)
     HNumber r1 = (r + r > n) ? n - r : r;
     HNumber r2 = n - r1;
 
     if (r1 >= 0) {
         if (n.isInteger() && r1.isInteger() && n <= 1000 && r1 <= 50)
             return factorial(n, r2+1) / factorial(r1, 1);
         HNumber result(n);
         floatnum rnum = &result.d->fnum;
         floatstruct fn, fr;
         float_create(&fn);
         float_create(&fr);
         float_copy(&fr, &r1.d->fnum, HMATH_EVAL_PREC);
         float_copy(&fn, rnum, EXACT);
         float_sub(rnum, rnum, &fr, HMATH_EVAL_PREC)
             && float_add(&fn, &fn, &c1, HMATH_EVAL_PREC)
             && float_add(&fr, &fr, &c1, HMATH_EVAL_PREC)
             && float_add(rnum, rnum, &c1, HMATH_EVAL_PREC)
             && float_lngamma(&fn, HMATH_EVAL_PREC)
             && float_lngamma(&fr, HMATH_EVAL_PREC)
             && float_lngamma(rnum, HMATH_EVAL_PREC)
             && float_add(rnum, rnum, &fr, HMATH_EVAL_PREC)
             && float_sub(rnum, &fn, rnum, HMATH_EVAL_PREC)
             && float_exp(rnum, HMATH_EVAL_PREC);
         float_free(&fn);
         float_free(&fr);
         roundSetError(result.d);
         return result;
     }
     if (r2 >= 0 || !r2.isInteger())
         return factorial(n, r1+1)/factorial(r2, 1);
     return 0;
 }
 
 /**
  * Returns the permutation of n elements chosen r elements.
  */
 HNumber HMath::nPr(const HNumber& n, const HNumber& r)
 {
     return factorial(n, (n-r+1));
 }
 
 /**
  * Returns the falling Pochhammer symbol x*(x-1)*..*base.
  * For base == 1, this is the usual factorial x!, what this
  * function is named after.
  * This function has been extended using the Gamma function,
  * so that actually Gamma(x+1)/Gamma(base) is computed, a
  * value that equals the falling Pochhammer symbol, when
  * x - base is an integer, but allows other differences as well.
  */
 HNumber HMath::factorial(const HNumber& x, const HNumber& base)
 {
     floatstruct tmp;
     if (float_cmp(&c1, &base.d->fnum) == 0) {
         HNumber result;
         call1Arg(result.d, x.d, float_factorial);
         return result;
     }
     float_create(&tmp);
     HNumber r(base);
     float_sub(&tmp, &x.d->fnum, &base.d->fnum, HMATH_EVAL_PREC)
     && float_add(&tmp, &tmp, &c1, HMATH_EVAL_PREC)
     && float_pochhammer(&r.d->fnum, &tmp, HMATH_EVAL_PREC);
     roundSetError(r.d);
     float_free(&tmp);
     return r;
 }
 
 static bool checkpn(const HNumber& p, const HNumber& n)
 {
     return n.isInteger() && !n.isNegative() && !p.isNan() && !p.isNegative() && p <= 1;
 }
 
 /**
  * Calculates the binomial discrete distribution probability mass function:
  * \f[X{\sim}B(n,p)\f]
  * \f[\Pr(X=k|n,p)={n\choose k}p^{k}(1-p)^{n-k}\f]
  *
  * \param[in] k the number of probed exact successes
  * \param[in] n the number of trials
  * \param[in] p the probability of success in a single trial
  *
  * \return the probability of exactly \p k successes, otherwise \p NaN if the
  * function is not defined for the specified parameters.
  */
 HNumber HMath::binomialPmf(const HNumber& k, const HNumber& n, const HNumber& p)
 {
     if (!k.isInteger() || !checkpn(p, n))
         return HMath::nan(InvalidParam);
 
     HNumber result = nCr(n, k);
     if (result.isZero())
         return result;
 
     // special case: powers of zero, 0^0 == 1 in this context
     if (p.isInteger())
         return (int)(p.isZero()? k.isZero() : n == k);
 
     return result * raise(p, k) * raise(HNumber(1) - p, n - k);
 }
 
 /**
  * Calculates the binomial cumulative distribution function:
  * \f[X{\sim}B(n,p)\f]
  * \f[\Pr(X \leq k|n,p)=\sum_{i=0}^{k}{n\choose i}p^{i}(1-p)^{n-i}\f]
  *
  * \param[in] k the maximum number of probed successes
  * \param[in] n the number of trials
  * \param[in] p the probability of success in a single trial
  *
  * \return the probability of up to \p k successes, otherwise \p NaN if the
  * function is not defined for the specified parameters.
  */
 HNumber HMath::binomialCdf(const HNumber& k, const HNumber& n, const HNumber& p)
 {
     // FIXME: Use the regularized incomplete Beta function to avoid the potentially very expensive loop.
     if (!k.isInteger() || n.isNan())
         return HMath::nan();
 
     // Initiates summation, checks arguments as well.
     HNumber summand = binomialPmf(0, n, p);
     if (summand.isNan())
         return summand;
 
     HNumber one(1);
 
     // Some early out results.
     if (k.isNegative())
         return 0;
     if (k >= n)
         return one;
 
     HNumber pcompl = one - p;
 
     if (p.isInteger())
         return pcompl;
 
     // use reflexion formula to limit summation
     if (k + k > n && k + one >= p * (n + one))
         return one - binomialCdf(n - k - one, n, pcompl);
 
     // loop adding binomialPdf
     HNumber result(summand);
     for (HNumber i(0); i < k;) {
         summand *= p * (n-i);
         i += one;
         summand /= pcompl * i;
         result += summand;
     }
     return result;
 }
 
 /**
  * Calculates the expected value of a binomially distributed random variable:
  * \f[X{\sim}B(n,p)\f]
  * \f[E(X)=np\f]
  *
  * \param[in] n the number of trials
  * \param[in] p the probability of success in a single trial
  *
  * \return the expected value of the variable, otherwise \p NaN if the
  * function is not defined for the specified parameters.
  */
 HNumber HMath::binomialMean(const HNumber& n, const HNumber& p)
 {
     if (!checkpn(p, n))
         return HMath::nan();
     return n * p;
 }
 
 /**
  * Calculates the variance of a binomially distributed random variable:
  * \f[X{\sim}B(n,p)\f]
  * \f[Var(X)=np(1-p)\f]
  *
  * \param[in] n the number of trials
  * \param[in] p the probability of success in a single trial
  *
  * \return the variance of the variable, otherwise \p NaN if the
  * function is not defined for the specified parameters.
  */
 HNumber HMath::binomialVariance(const HNumber& n, const HNumber& p)
 {
     return binomialMean(n, p) * (HNumber(1) - p);
 }
 
 static bool checkNMn(const HNumber& N, const HNumber& M, const HNumber& n)
 {
     return N.isInteger() && !N.isNegative()
         && M.isInteger() && !M.isNegative()
         && n.isInteger() && !n.isNegative()
         && HMath::max(M, n) <= N;
 }
 
 /**
  * Calculates the hypergeometric discrete distribution probability mass
  * function:
  * \f[X{\sim}H(N,M,n)\f]
  * \f[\Pr(X=k|N,M,n)=\frac{{M\choose k}{N-M\choose n-k}}{{N\choose n}}\f]
  *
  * \param[in] k the number of probed exact successes
  * \param[in] N the number of total elements
  * \param[in] M the number of success elements
  * \param[in] n the number of selected elements
  *
  * \return the probability of exactly \p k successes, otherwise \p NaN if the
  * function is not defined for the specified parameters.
  */
 HNumber HMath::hypergeometricPmf(const HNumber& k, const HNumber& N, const HNumber& M, const HNumber& n)
 {
     if (!k.isInteger() || !checkNMn(N, M, n))
         return HMath::nan();
     return HMath::nCr(M, k) * HMath::nCr(N - M, n - k) / HMath::nCr(N, n);
 }
 
 /**
  * Calculates the hypergeometric cumulative distribution function:
  * \f[X{\sim}H(N,M,n)\f]
  * \f[\Pr(X\leq k|N,M,n)=
  *   \sum_{i=0}^{k}\frac{{M\choose k}{N-M\choose n-k}}{{N\choose n}}\f]
  *
  * \param[in] k the maximum number of probed successes
  * \param[in] N the number of total elements
  * \param[in] M the number of success elements
  * \param[in] n the number of selected elements
  *
  * \return the probability of up to \p k successes, otherwise \p NaN if the
  * function is not defined for the specified parameters.
  */
 HNumber HMath::hypergeometricCdf(const HNumber& k, const HNumber& N, const HNumber& M, const HNumber& n)
 {
     // Lowest index of non-zero summand in loop.
     HNumber c = M + n - N;
     HNumber i = max(c, 0);
 
     // First summand in loop, do the parameter checking here.
     HNumber summand = HMath::hypergeometricPmf(i, N, M, n);
     if (!k.isInteger() || summand.isNan())
         return HMath::nan();
 
     // Some early out results.
     HNumber one(1);
     if (k >= M || k >= n)
         return one;
     if (i > k)
         return 0;
 
     // Use reflexion formula to limit summations. Numerically unstable where the result is near 0, disable for now.
     //   if (k + k > n)
     //     return one - hypergeometricCdf(n - k - 1, N, N - M, n);
 
     HNumber result = summand;
     for (; i < k;) {
         summand *= (M - i) * (n - i);
         i += one;
         summand /= i * (i - c);
         result += summand;
     }
     return result;
 }
 
 /**
  * Calculates the expected value of a hypergeometrically distributed random
  * variable:
  * \f[X{\sim}H(N,M,n)\f]
  * \f[E(X)=n\frac{M}{N}\f]
  *
  * \param[in] N the number of total elements
  * \param[in] M the number of success elements
  * \param[in] n the number of selected elements
  *
  * \return the expected value of the variable, otherwise \p NaN if the
  * function is not defined for the specified parameter.
  */
 HNumber HMath::hypergeometricMean(const HNumber& N, const HNumber& M, const HNumber& n)
 {
     if (!checkNMn(N, M, n))
         return HMath::nan();
     return n * M / N;
 }
 
 /**
  *
  * Calculates the variance of a hypergeometrically distributed random variable:
  * \f[X{\sim}H(N,M,n)\f]
  * \f[Var(X)=\frac{n\frac{M}{N}(1-\frac{M}{N})(N-n)}{N-1}\f]
  *
  * \param[in] N the number of total elements
  * \param[in] M the number of success elements
  * \param[in] n the number of selected elements
  *
  * \return the variance of the variable, otherwise \p NaN if the function is
  * not defined for the specified parameter.
  */
 HNumber HMath::hypergeometricVariance(const HNumber& N, const HNumber& M, const HNumber& n)
 {
   return (hypergeometricMean(N, M, n) * (HNumber(1) - M / N) * (N - n)) / (N - HNumber(1));
 }
 
 /**
  *
  * Calculates the poissonian discrete distribution probability mass function:
  * \f[X{\sim}P(\lambda)\f]
  * \f[\Pr(X=k|\lambda)=\frac{e^{-\lambda}\lambda^k}{k!}\f]
  *
  * \param[in] k the number of event occurrences
  * \param[in] l the expected number of occurrences that occur in an interval
  *
  * \return the probability of exactly \p k event occurrences, otherwise \p NaN
  * if the function is not defined for the specified parameters.
  */
 HNumber HMath::poissonPmf(const HNumber& k, const HNumber& l)
 {
     if (!k.isInteger() || l.isNan() || l.isNegative())
         return HMath::nan();
 
     if (k.isNegative())
         return 0;
     if (l.isZero())
         return int(k.isZero());
 
     return exp(-l) * raise(l, k) / factorial(k);
 }
 
 /**
  * Calculates the poissonian cumulative distribution function:
  * \f[X{\sim}P(\lambda)\f]
  * \f[\Pr(X\leq k|\lambda)=\sum_{i=0}^{k}\frac{e^{-\lambda}\lambda^k}{k!}\f]
  *
  * \param[in] k the maximum number of event occurrences
  * \param[in] l the expected number of occurrences that occur in an interval
  *
  * \return the probability of up to \p k event occurrences, otherwise \p NaN
  * if the function is not defined for the specified parameters.
  */
 HNumber HMath::poissonCdf(const HNumber& k, const HNumber& l)
 {
     // FIXME: Use the incomplete gamma function to avoid a potentially expensive loop.
     if (!k.isInteger() || l.isNan() || l.isNegative())
         return HMath::nan();
 
     if (k.isNegative())
         return 0;
 
     HNumber one(1);
     if (l.isZero())
         return one;
 
     HNumber summand = one;
     HNumber result = one;
     for (HNumber i = one; i <= k; i += one) {
         summand *= l / i;
         result += summand;
     }
     result = exp(-l) * result;
     return result;
 }
 
 /**
  * Calculates the expected value of a Poisson distributed random variable:
  * \f[X{\sim}P(\lambda)\f]
  * \f[E(X)=\lambda\f]
  *
  * \param[in] l the expected number of occurrences that occur in an interval
  *
  * \return the expected value of the variable, otherwise \p NaN if the
  * function is not defined for the specified parameter.
  */
 HNumber HMath::poissonMean(const HNumber& l)
 {
     if (l.isNan() || l.isNegative())
         return HMath::nan();
     return l;
 }
 
 /**
  * Calculates the variance of a Poisson distributed random variable:
  * \f[X{\sim}P(\lambda)\f]
  * \f[Var(X)=\lambda\f]
  *
  * \param[in] l the expected number of occurrences that occur in an interval
  *
  * \return the variance of the variable, otherwise \p NaN if the function is
  * not defined for the specified parameter.
  */
 HNumber HMath::poissonVariance(const HNumber& l)
 {
     return poissonMean(l);
 }
 
 /**
  * Returns the erf function (related to normal distribution).
  */
 HNumber HMath::erf(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_erf);
     return result;
 }
 
 /**
  * Returns the complementary erf function (related to normal distribution).
  */
 HNumber HMath::erfc(const HNumber& x)
 {
     HNumber result;
     call1Arg(result.d, x.d, float_erfc);
     return result;
 }
 
 /**
  * Restricts a logic value to a given bit size.
  */
 HNumber HMath::mask(const HNumber& val, const HNumber& bits)
 {
     if (val.isNan() || bits == 0 || bits >= LOGICRANGE || !bits.isInteger())
         return HMath::nan();
     return val & ~(HNumber(-1) << HNumber(bits));
 }
 
 /**
  * sign-extends an unsigned value
  */
 HNumber HMath::sgnext(const HNumber& val, const HNumber& bits)
 {
     if (val.isNan() || bits == 0 || bits >= LOGICRANGE || !bits.isInteger())
         return HMath::nan();
     HNumber ofs = HNumber(LOGICRANGE) - bits;
     return (val << ofs) >> ofs;
 }
 
 /**
  * For bits >= 0 does an arithmetic shift right, for bits < 0 a shift left.
  */
 HNumber HMath::ashr(const HNumber& val, const HNumber& bits)
 {
     if (val.isNan() || bits <= -LOGICRANGE || bits >= LOGICRANGE || !bits.isInteger())
         return HMath::nan();
     if (bits >= 0)
         return val >> bits;
     return val << -bits;
 }
 
 /**
  * Decode an IEEE-754 bit pattern with the default exponent bias
  */
 HNumber HMath::decodeIeee754(const HNumber& val, const HNumber& exp_bits, const HNumber& significand_bits)
 {
     return HMath::decodeIeee754(val, exp_bits, significand_bits, HMath::raise(2, exp_bits - 1) - 1);
 }
 
 /**
  * Decode an IEEE-754 bit pattern with the given parameters
  */
 HNumber HMath::decodeIeee754(const HNumber& val, const HNumber& exp_bits, const HNumber& significand_bits,
                              const HNumber& exp_bias)
 {
     if (val.isNan()
         || exp_bits <= 0 || exp_bits >= LOGICRANGE || !exp_bits.isInteger()
         || significand_bits <= 0 || significand_bits >= LOGICRANGE || !significand_bits.isInteger()
         || !exp_bias.isInteger())
         return HMath::nan();
 
     HNumber sign(HMath::mask(val >> (exp_bits + significand_bits), 1).isZero() ? 1 : -1);
     HNumber exp = HMath::mask(val >> significand_bits, exp_bits);
     // <=> '0.' b_x b_x-1 b_x-2 ... b_0
     HNumber significand = HMath::mask(val, significand_bits) * HMath::raise(2, -significand_bits);
 
     if (exp.isZero()) {
         // Exponent 0, subnormal value or zero.
         return sign * significand * HMath::raise(2, exp - exp_bias + 1);
     }
     if (exp - HMath::raise(2, exp_bits) == -1) {
         // Exponent all 1...
         if (significand.isZero())
             // ...and signficand 0, infinity.
             // TODO: Represent infinity as something other than NaN?
             return HNumber();
         // ...and significand not 0, NaN.
         return HMath::nan();
     }
     return sign * (significand + 1) * HMath::raise(2, exp - exp_bias); // Normalised value.
 }
 
 /**
  * Encode a value in a IEEE-754 binary representation with the default exponent bias
  */
 HNumber HMath::encodeIeee754(const HNumber& val, const HNumber& exp_bits, const HNumber& significand_bits)
 {
     return HMath::encodeIeee754(val, exp_bits, significand_bits, HMath::raise(2, exp_bits - 1) - 1);
 }
 
 /**
  * Encode a value in a IEEE-754 binary representation
  */
 HNumber HMath::encodeIeee754(const HNumber& val, const HNumber& exp_bits, const HNumber& significand_bits,
                              const HNumber& exp_bias)
 {
     if (exp_bits <= 0 || exp_bits >= LOGICRANGE || !exp_bits.isInteger()
         || significand_bits <= 0 || significand_bits >= LOGICRANGE || !significand_bits.isInteger()
         || !exp_bias.isInteger())
         return HMath::nan();
 
     HNumber sign_bit;
     HNumber significand;
     HNumber exponent;
 
     HNumber min_exp(1 - exp_bias);
     HNumber max_exp(HMath::raise(2, exp_bits) - 2 - exp_bias);
 
     if (val.isNan()) {
         // Encode a NaN.
         sign_bit = 0;
         significand = HMath::raise(2, significand_bits) - 1;
         exponent = HMath::raise(2, exp_bits) - 1;
     } else if (val.isZero()) {
         // Encode a basic 0.
         sign_bit = 0;
         significand = 0;
         exponent = 0;
     } else {
         // Regular input value.
         sign_bit = val.isNegative() ? 1 : 0;
         // Determine exponent.
         HNumber search_min = min_exp;
         HNumber search_max = max_exp;
         exponent = HMath::ceil((search_max - search_min) / 2) + search_min;
         significand = HMath::abs(val) * HMath::raise(2, -exponent);
 
         do {
             if (significand >= 1 && significand < 2) {
                 // Integer part is 1, stop here.
                 break;
             } else if (significand >= 2) {
                 // Increase exponent.
                 search_min = exponent + 1;
             } else if (significand < 1) {
                 // Decrease exponent.
                 search_max = exponent - 1;
             }
             exponent = HMath::ceil((search_max - search_min) / 2) + search_min;
             significand = HMath::abs(val) * HMath::raise(2, -exponent);
         } while (exponent != min_exp && exponent != max_exp);
 
         HNumber rounded = HMath::round(significand * HMath::raise(2, significand_bits));
         HNumber intpart = HMath::integer(rounded * HMath::raise(2, -significand_bits));
         significand = HMath::mask(rounded, significand_bits);
         if (intpart.isZero()) {
             // Subnormal value.
             exponent = 0;
         } else if (intpart > 1) {
             // Infinity.
             exponent = HMath::raise(2, exp_bits) - 1;
             significand = 0;
         } else {
             // Normalised value.
             exponent = exponent + exp_bias;
         }
     }
 
     HNumber result = sign_bit << (exp_bits + significand_bits) | exponent << (significand_bits) | significand;
     return result;
 }
 
 std::ostream& operator<<(std::ostream& s, const HNumber& n)
 {
     QString str = HMath::format(n);
     s << str.toLatin1().constData();
     return s;
 }
 
 struct MathInit {
     MathInit(){ floatmath_init(); }
 };
 
 MathInit mathinit;
 
 /**
  * Parses a string containing a real number.
  *
  * Parameters :
  *   str_in : pointer towards the string to parse
  *   str_out : pointer towards a pointer towards the remaining of the string after parsing
  */
 HNumber HMath::parse_str(const char* str_in, const char** str_out) {
 
     // FIXME: Duplicate code.
     // FIXME: Error management.
 
     const char* str = str_in;
     t_itokens tokens;
 
     HNumber x;
     delete x.d;
 
     x.d = new HNumberPrivate;
     if ((x.d->error = parse(&tokens, &str)) == Success)
     x.d->error = float_in(&x.d->fnum, &tokens);
     float_geterror();
 
     /* Store remaining of the string */
     if (str_out)
         *str_out = str;
 
     return x;
 }
 
 bool HNumber::isNearZero() const
 {
     return float_iszero(&(d->fnum)) || float_getexponent(&(d->fnum)) <= -80;
 }