From d286182d253c99a286e044bdc10189f7e3a662e0 Mon Sep 17 00:00:00 2001
From: RD42 <42702181+dashr9230@users.noreply.github.com>
Date: Mon, 13 Nov 2023 20:32:38 +0800
Subject: [PATCH] Add ttmath 0.9.1 library
---
saco/saco.vcproj | 10 +
server/server.vcproj | 10 +
ttmath/ttmath.h | 2835 ++++++++++++++++
ttmath/ttmathbig.h | 5236 +++++++++++++++++++++++++++++
ttmath/ttmathint.h | 1547 +++++++++
ttmath/ttmathmisc.h | 243 ++
ttmath/ttmathobjects.h | 766 +++++
ttmath/ttmathparser.h | 2754 +++++++++++++++
ttmath/ttmaththreads.h | 250 ++
ttmath/ttmathtypes.h | 646 ++++
ttmath/ttmathuint.h | 3520 +++++++++++++++++++
ttmath/ttmathuint_noasm.h | 1013 ++++++
ttmath/ttmathuint_x86.h | 1602 +++++++++
ttmath/ttmathuint_x86_64.h | 1222 +++++++
ttmath/ttmathuint_x86_64_msvc.asm | 548 +++
15 files changed, 22202 insertions(+)
create mode 100644 ttmath/ttmath.h
create mode 100644 ttmath/ttmathbig.h
create mode 100644 ttmath/ttmathint.h
create mode 100644 ttmath/ttmathmisc.h
create mode 100644 ttmath/ttmathobjects.h
create mode 100644 ttmath/ttmathparser.h
create mode 100644 ttmath/ttmaththreads.h
create mode 100644 ttmath/ttmathtypes.h
create mode 100644 ttmath/ttmathuint.h
create mode 100644 ttmath/ttmathuint_noasm.h
create mode 100644 ttmath/ttmathuint_x86.h
create mode 100644 ttmath/ttmathuint_x86_64.h
create mode 100644 ttmath/ttmathuint_x86_64_msvc.asm
diff --git a/saco/saco.vcproj b/saco/saco.vcproj
index 1bfd1ee..615add3 100644
--- a/saco/saco.vcproj
+++ b/saco/saco.vcproj
@@ -230,6 +230,16 @@
RelativePath=".\archive\TinyEncrypt.h">
+
+
+
+
+
+
diff --git a/server/server.vcproj b/server/server.vcproj
index b4daa51..0c1b618 100644
--- a/server/server.vcproj
+++ b/server/server.vcproj
@@ -108,6 +108,16 @@
+
+
+
+
+
+
diff --git a/ttmath/ttmath.h b/ttmath/ttmath.h
new file mode 100644
index 0000000..8bd24f5
--- /dev/null
+++ b/ttmath/ttmath.h
@@ -0,0 +1,2835 @@
+/*
+ * This file is a part of TTMath Bignum Library
+ * and is distributed under the (new) BSD licence.
+ * Author: Tomasz Sowa
+ */
+
+/*
+ * Copyright (c) 2006-2009, Tomasz Sowa
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * * Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ *
+ * * Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * * Neither the name Tomasz Sowa nor the names of contributors to this
+ * project may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+ * THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+
+
+#ifndef headerfilettmathmathtt
+#define headerfilettmathmathtt
+
+/*!
+ \file ttmath.h
+ \brief Mathematics functions.
+*/
+
+#ifdef _MSC_VER
+//warning C4127: conditional expression is constant
+#pragma warning( disable: 4127 )
+//warning C4702: unreachable code
+#pragma warning( disable: 4702 )
+//warning C4800: forcing value to bool 'true' or 'false' (performance warning)
+#pragma warning( disable: 4800 )
+#endif
+
+
+#include "ttmathbig.h"
+#include "ttmathobjects.h"
+
+
+namespace ttmath
+{
+ /*
+ *
+ * functions defined here are used only with Big<> types
+ *
+ *
+ */
+
+
+ /*
+ *
+ * functions for rounding
+ *
+ *
+ */
+
+
+ /*!
+ this function skips the fraction from x
+ e.g 2.2 = 2
+ 2.7 = 2
+ -2.2 = 2
+ -2.7 = 2
+ */
+ template
+ ValueType SkipFraction(const ValueType & x)
+ {
+ ValueType result( x );
+ result.SkipFraction();
+
+ return result;
+ }
+
+
+ /*!
+ this function rounds to the nearest integer value
+ e.g 2.2 = 2
+ 2.7 = 3
+ -2.2 = -2
+ -2.7 = -3
+ */
+ template
+ ValueType Round(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType result( x );
+ uint c = result.Round();
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+
+ /*!
+ this function returns a value representing the smallest integer
+ that is greater than or equal to x
+
+ Ceil(-3.7) = -3
+ Ceil(-3.1) = -3
+ Ceil(-3.0) = -3
+ Ceil(4.0) = 4
+ Ceil(4.2) = 5
+ Ceil(4.8) = 5
+ */
+ template
+ ValueType Ceil(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType result(x);
+ uint c = 0;
+
+ result.SkipFraction();
+
+ if( result != x )
+ {
+ // x is with fraction
+ // if x is negative we don't have to do anything
+ if( !x.IsSign() )
+ {
+ ValueType one;
+ one.SetOne();
+
+ c += result.Add(one);
+ }
+ }
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ this function returns a value representing the largest integer
+ that is less than or equal to x
+
+ Floor(-3.6) = -4
+ Floor(-3.1) = -4
+ Floor(-3) = -3
+ Floor(2) = 2
+ Floor(2.3) = 2
+ Floor(2.8) = 2
+ */
+ template
+ ValueType Floor(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType result(x);
+ uint c = 0;
+
+ result.SkipFraction();
+
+ if( result != x )
+ {
+ // x is with fraction
+ // if x is positive we don't have to do anything
+ if( x.IsSign() )
+ {
+ ValueType one;
+ one.SetOne();
+
+ c += result.Sub(one);
+ }
+ }
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+
+ /*
+ *
+ * logarithms and the exponent
+ *
+ *
+ */
+
+
+ /*!
+ this function calculates the natural logarithm (logarithm with the base 'e')
+ */
+ template
+ ValueType Ln(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType result;
+ uint state = result.Ln(x);
+
+ if( err )
+ {
+ switch( state )
+ {
+ case 0:
+ *err = err_ok;
+ break;
+ case 1:
+ *err = err_overflow;
+ break;
+ case 2:
+ *err = err_improper_argument;
+ break;
+ default:
+ *err = err_internal_error;
+ break;
+ }
+ }
+
+
+ return result;
+ }
+
+
+ /*!
+ this function calculates the logarithm
+ */
+ template
+ ValueType Log(const ValueType & x, const ValueType & base, ErrorCode * err = 0)
+ {
+ if( x.IsNan() || base.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return ValueType(); // default NaN
+ }
+
+ ValueType result;
+ uint state = result.Log(x, base);
+
+ if( err )
+ {
+ switch( state )
+ {
+ case 0:
+ *err = err_ok;
+ break;
+ case 1:
+ *err = err_overflow;
+ break;
+ case 2:
+ case 3:
+ *err = err_improper_argument;
+ break;
+ default:
+ *err = err_internal_error;
+ break;
+ }
+ }
+
+ return result;
+ }
+
+
+ /*!
+ this function calculates the expression e^x
+ */
+ template
+ ValueType Exp(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType result;
+ uint c = result.Exp(x);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ *
+ * trigonometric functions
+ *
+ */
+
+
+ /*
+ this namespace consists of auxiliary functions
+ (something like 'private' in a class)
+ */
+ namespace auxiliaryfunctions
+ {
+
+ /*!
+ an auxiliary function for calculating the Sine
+ (you don't have to call this function)
+ */
+ template
+ uint PrepareSin(ValueType & x, bool & change_sign)
+ {
+ ValueType temp;
+
+ change_sign = false;
+
+ if( x.IsSign() )
+ {
+ // we're using the formula 'sin(-x) = -sin(x)'
+ change_sign = !change_sign;
+ x.ChangeSign();
+ }
+
+ // we're reducing the period 2*PI
+ // (for big values there'll always be zero)
+ temp.Set2Pi();
+
+ if( x.Mod(temp) )
+ return 1;
+
+
+ // we're setting 'x' as being in the range of <0, 0.5PI>
+
+ temp.SetPi();
+
+ if( x > temp )
+ {
+ // x is in (pi, 2*pi>
+ x.Sub( temp );
+ change_sign = !change_sign;
+ }
+
+ temp.Set05Pi();
+
+ if( x > temp )
+ {
+ // x is in (0.5pi, pi>
+ x.Sub( temp );
+ x = temp - x;
+ }
+
+ return 0;
+ }
+
+
+ /*!
+ an auxiliary function for calculating the Sine
+ (you don't have to call this function)
+
+ it returns Sin(x) where 'x' is from <0, PI/2>
+ we're calculating the Sin with using Taylor series in zero or PI/2
+ (depending on which point of these two points is nearer to the 'x')
+
+ Taylor series:
+ sin(x) = sin(a) + cos(a)*(x-a)/(1!)
+ - sin(a)*((x-a)^2)/(2!) - cos(a)*((x-a)^3)/(3!)
+ + sin(a)*((x-a)^4)/(4!) + ...
+
+ when a=0 it'll be:
+ sin(x) = (x)/(1!) - (x^3)/(3!) + (x^5)/(5!) - (x^7)/(7!) + (x^9)/(9!) ...
+
+ and when a=PI/2:
+ sin(x) = 1 - ((x-PI/2)^2)/(2!) + ((x-PI/2)^4)/(4!) - ((x-PI/2)^6)/(6!) ...
+ */
+ template
+ ValueType Sin0pi05(const ValueType & x)
+ {
+ ValueType result;
+ ValueType numerator, denominator;
+ ValueType d_numerator, d_denominator;
+ ValueType one, temp, old_result;
+
+ // temp = pi/4
+ temp.Set05Pi();
+ temp.exponent.SubOne();
+
+ one.SetOne();
+
+ if( x < temp )
+ {
+ // we're using the Taylor series with a=0
+ result = x;
+ numerator = x;
+ denominator = one;
+
+ // d_numerator = x^2
+ d_numerator = x;
+ d_numerator.Mul(x);
+
+ d_denominator = 2;
+ }
+ else
+ {
+ // we're using the Taylor series with a=PI/2
+ result = one;
+ numerator = one;
+ denominator = one;
+
+ // d_numerator = (x-pi/2)^2
+ ValueType pi05;
+ pi05.Set05Pi();
+
+ temp = x;
+ temp.Sub( pi05 );
+ d_numerator = temp;
+ d_numerator.Mul( temp );
+
+ d_denominator = one;
+ }
+
+ uint c = 0;
+ bool addition = false;
+
+ old_result = result;
+ for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
+ {
+ // we're starting from a second part of the formula
+ c += numerator. Mul( d_numerator );
+ c += denominator. Mul( d_denominator );
+ c += d_denominator.Add( one );
+ c += denominator. Mul( d_denominator );
+ c += d_denominator.Add( one );
+ temp = numerator;
+ c += temp.Div(denominator);
+
+ if( c )
+ // Sin is from <-1,1> and cannot make an overflow
+ // but the carry can be from the Taylor series
+ // (then we only break our calculations)
+ break;
+
+ if( addition )
+ result.Add( temp );
+ else
+ result.Sub( temp );
+
+
+ addition = !addition;
+
+ // we're testing whether the result has changed after adding
+ // the next part of the Taylor formula, if not we end the loop
+ // (it means 'x' is zero or 'x' is PI/2 or this part of the formula
+ // is too small)
+ if( result == old_result )
+ break;
+
+ old_result = result;
+ }
+
+ return result;
+ }
+
+ } // namespace auxiliaryfunctions
+
+
+
+ /*!
+ this function calculates the Sine
+ */
+ template
+ ValueType Sin(ValueType x, ErrorCode * err = 0)
+ {
+ using namespace auxiliaryfunctions;
+
+ ValueType one, result;
+ bool change_sign;
+
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ if( err )
+ *err = err_ok;
+
+ if( PrepareSin( x, change_sign ) )
+ {
+ // x is too big, we cannnot reduce the 2*PI period
+ // prior to version 0.8.5 the result was zero
+
+ // result has NaN flag set by default
+
+ if( err )
+ *err = err_overflow; // maybe another error code? err_improper_argument?
+
+ return result; // NaN is set by default
+ }
+
+ result = Sin0pi05( x );
+
+ one.SetOne();
+
+ // after calculations there can be small distortions in the result
+ if( result > one )
+ result = one;
+ else
+ if( result.IsSign() )
+ // we've calculated the sin from <0, pi/2> and the result
+ // should be positive
+ result.SetZero();
+
+ if( change_sign )
+ result.ChangeSign();
+
+ return result;
+ }
+
+
+ /*!
+ this function calulates the Cosine
+ we're using the formula cos(x) = sin(x + PI/2)
+ */
+ template
+ ValueType Cos(ValueType x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType pi05;
+ pi05.Set05Pi();
+
+ uint c = x.Add( pi05 );
+
+ if( c )
+ {
+ if( err )
+ *err = err_overflow;
+
+ return ValueType(); // result is undefined (NaN is set by default)
+ }
+
+ return Sin(x, err);
+ }
+
+
+ /*!
+ this function calulates the Tangent
+ we're using the formula tan(x) = sin(x) / cos(x)
+
+ it takes more time than calculating the Tan directly
+ from for example Taylor series but should be a bit preciser
+ because Tan receives its values from -infinity to +infinity
+ and when we calculate it from any series then we can make
+ a greater mistake than calculating 'sin/cos'
+ */
+ template
+ ValueType Tan(const ValueType & x, ErrorCode * err = 0)
+ {
+ ValueType result = Cos(x, err);
+
+ if( err && *err != err_ok )
+ return result;
+
+ if( result.IsZero() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ result.SetNan();
+
+ return result;
+ }
+
+ return Sin(x, err) / result;
+ }
+
+
+ /*!
+ this function calulates the Tangent
+ look at the description of Tan(...)
+
+ (the abbreviation of Tangent can be 'tg' as well)
+ */
+ template
+ ValueType Tg(const ValueType & x, ErrorCode * err = 0)
+ {
+ return Tan(x, err);
+ }
+
+
+ /*!
+ this function calulates the Cotangent
+ we're using the formula tan(x) = cos(x) / sin(x)
+
+ (why do we make it in this way?
+ look at information in Tan() function)
+ */
+ template
+ ValueType Cot(const ValueType & x, ErrorCode * err = 0)
+ {
+ ValueType result = Sin(x, err);
+
+ if( err && *err != err_ok )
+ return result;
+
+ if( result.IsZero() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ result.SetNan();
+
+ return result;
+ }
+
+ return Cos(x, err) / result;
+ }
+
+
+ /*!
+ this function calulates the Cotangent
+ look at the description of Cot(...)
+
+ (the abbreviation of Cotangent can be 'ctg' as well)
+ */
+ template
+ ValueType Ctg(const ValueType & x, ErrorCode * err = 0)
+ {
+ return Cot(x, err);
+ }
+
+
+ /*
+ *
+ * inverse trigonometric functions
+ *
+ *
+ */
+
+ namespace auxiliaryfunctions
+ {
+
+ /*!
+ an auxiliary function for calculating the Arc Sine
+
+ we're calculating asin from the following formula:
+ asin(x) = x + (1*x^3)/(2*3) + (1*3*x^5)/(2*4*5) + (1*3*5*x^7)/(2*4*6*7) + ...
+ where abs(x) <= 1
+
+ we're using this formula when x is from <0, 1/2>
+ */
+ template
+ ValueType ASin_0(const ValueType & x)
+ {
+ ValueType nominator, denominator, nominator_add, nominator_x, denominator_add, denominator_x;
+ ValueType two, result(x), x2(x);
+ ValueType nominator_temp, denominator_temp, old_result = result;
+ uint c = 0;
+
+ x2.Mul(x);
+ two = 2;
+
+ nominator.SetOne();
+ denominator = two;
+ nominator_add = nominator;
+ denominator_add = denominator;
+ nominator_x = x;
+ denominator_x = 3;
+
+ for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
+ {
+ c += nominator_x.Mul(x2);
+ nominator_temp = nominator_x;
+ c += nominator_temp.Mul(nominator);
+ denominator_temp = denominator;
+ c += denominator_temp.Mul(denominator_x);
+ c += nominator_temp.Div(denominator_temp);
+
+ // if there is a carry somewhere we only break the calculating
+ // the result should be ok -- it's from <-pi/2, pi/2>
+ if( c )
+ break;
+
+ result.Add(nominator_temp);
+
+ if( result == old_result )
+ // there's no sense to calculate more
+ break;
+
+ old_result = result;
+
+
+ c += nominator_add.Add(two);
+ c += denominator_add.Add(two);
+ c += nominator.Mul(nominator_add);
+ c += denominator.Mul(denominator_add);
+ c += denominator_x.Add(two);
+ }
+
+ return result;
+ }
+
+
+
+ /*!
+ an auxiliary function for calculating the Arc Sine
+
+ we're calculating asin from the following formula:
+ asin(x) = pi/2 - sqrt(2)*sqrt(1-x) * asin_temp
+ asin_temp = 1 + (1*(1-x))/((2*3)*(2)) + (1*3*(1-x)^2)/((2*4*5)*(4)) + (1*3*5*(1-x)^3)/((2*4*6*7)*(8)) + ...
+
+ where abs(x) <= 1
+
+ we're using this formula when x is from (1/2, 1>
+ */
+ template
+ ValueType ASin_1(const ValueType & x)
+ {
+ ValueType nominator, denominator, nominator_add, nominator_x, nominator_x_add, denominator_add, denominator_x;
+ ValueType denominator2;
+ ValueType one, two, result;
+ ValueType nominator_temp, denominator_temp, old_result;
+ uint c = 0;
+
+ two = 2;
+
+ one.SetOne();
+ nominator = one;
+ result = one;
+ old_result = result;
+ denominator = two;
+ nominator_add = nominator;
+ denominator_add = denominator;
+ nominator_x = one;
+ nominator_x.Sub(x);
+ nominator_x_add = nominator_x;
+ denominator_x = 3;
+ denominator2 = two;
+
+
+ for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
+ {
+ nominator_temp = nominator_x;
+ c += nominator_temp.Mul(nominator);
+ denominator_temp = denominator;
+ c += denominator_temp.Mul(denominator_x);
+ c += denominator_temp.Mul(denominator2);
+ c += nominator_temp.Div(denominator_temp);
+
+ // if there is a carry somewhere we only break the calculating
+ // the result should be ok -- it's from <-pi/2, pi/2>
+ if( c )
+ break;
+
+ result.Add(nominator_temp);
+
+ if( result == old_result )
+ // there's no sense to calculate more
+ break;
+
+ old_result = result;
+
+ c += nominator_x.Mul(nominator_x_add);
+ c += nominator_add.Add(two);
+ c += denominator_add.Add(two);
+ c += nominator.Mul(nominator_add);
+ c += denominator.Mul(denominator_add);
+ c += denominator_x.Add(two);
+ c += denominator2.Mul(two);
+ }
+
+
+ nominator_x_add.exponent.AddOne(); // *2
+ one.exponent.SubOne(); // =0.5
+ nominator_x_add.Pow(one); // =sqrt(nominator_x_add)
+ result.Mul(nominator_x_add);
+
+ one.Set05Pi();
+ one.Sub(result);
+
+ return one;
+ }
+
+
+ } // namespace auxiliaryfunctions
+
+
+ /*!
+ this function calculates the Arc Sine
+ x is from <-1,1>
+ */
+ template
+ ValueType ASin(ValueType x, ErrorCode * err = 0)
+ {
+ using namespace auxiliaryfunctions;
+
+ ValueType result, one;
+ one.SetOne();
+ bool change_sign = false;
+
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ if( x.GreaterWithoutSignThan(one) )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ if( x.IsSign() )
+ {
+ change_sign = true;
+ x.Abs();
+ }
+
+ one.exponent.SubOne(); // =0.5
+
+ // asin(-x) = -asin(x)
+ if( x.GreaterWithoutSignThan(one) )
+ result = ASin_1(x);
+ else
+ result = ASin_0(x);
+
+ if( change_sign )
+ result.ChangeSign();
+
+ if( err )
+ *err = err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ this function calculates the Arc Cosine
+
+ we're using the formula:
+ acos(x) = pi/2 - asin(x)
+ */
+ template
+ ValueType ACos(const ValueType & x, ErrorCode * err = 0)
+ {
+ ValueType temp;
+
+ temp.Set05Pi();
+ temp.Sub(ASin(x, err));
+
+ return temp;
+ }
+
+
+
+ namespace auxiliaryfunctions
+ {
+
+ /*!
+ an auxiliary function for calculating the Arc Tangent
+
+ arc tan (x) where x is in <0; 0.5)
+ (x can be in (-0.5 ; 0.5) too)
+
+ we're using the Taylor series expanded in zero:
+ atan(x) = x - (x^3)/3 + (x^5)/5 - (x^7)/7 + ...
+ */
+ template
+ ValueType ATan0(const ValueType & x)
+ {
+ ValueType nominator, denominator, nominator_add, denominator_add, temp;
+ ValueType result, old_result;
+ bool adding = false;
+ uint c = 0;
+
+ result = x;
+ old_result = result;
+ nominator = x;
+ nominator_add = x;
+ nominator_add.Mul(x);
+
+ denominator.SetOne();
+ denominator_add = 2;
+
+ for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
+ {
+ c += nominator.Mul(nominator_add);
+ c += denominator.Add(denominator_add);
+
+ temp = nominator;
+ c += temp.Div(denominator);
+
+ if( c )
+ // the result should be ok
+ break;
+
+ if( adding )
+ result.Add(temp);
+ else
+ result.Sub(temp);
+
+ if( result == old_result )
+ // there's no sense to calculate more
+ break;
+
+ old_result = result;
+ adding = !adding;
+ }
+
+ return result;
+ }
+
+
+ /*!
+ an auxiliary function for calculating the Arc Tangent
+
+ where x is in <0 ; 1>
+ */
+ template
+ ValueType ATan01(const ValueType & x)
+ {
+ ValueType half;
+ half.Set05();
+
+ /*
+ it would be better if we chose about sqrt(2)-1=0.41... instead of 0.5 here
+
+ because as you can see below:
+ when x = sqrt(2)-1
+ abs(x) = abs( (x-1)/(1+x) )
+ so when we're calculating values around x
+ then they will be better converged to each other
+
+ for example if we have x=0.4999 then during calculating ATan0(0.4999)
+ we have to make about 141 iterations but when we have x=0.5
+ then during calculating ATan0( (x-1)/(1+x) ) we have to make
+ only about 89 iterations (both for Big<3,9>)
+
+ in the future this 0.5 can be changed
+ */
+ if( x.SmallerWithoutSignThan(half) )
+ return ATan0(x);
+
+
+ /*
+ x>=0.5 and x<=1
+ (x can be even smaller than 0.5)
+
+ y = atac(x)
+ x = tan(y)
+
+ tan(y-b) = (tan(y)-tab(b)) / (1+tan(y)*tan(b))
+ y-b = atan( (tan(y)-tab(b)) / (1+tan(y)*tan(b)) )
+ y = b + atan( (x-tab(b)) / (1+x*tan(b)) )
+
+ let b = pi/4
+ tan(b) = tan(pi/4) = 1
+ y = pi/4 + atan( (x-1)/(1+x) )
+
+ so
+ atac(x) = pi/4 + atan( (x-1)/(1+x) )
+ when x->1 (x converges to 1) the (x-1)/(1+x) -> 0
+ and we can use ATan0() function here
+ */
+
+ ValueType n(x),d(x),one,result;
+
+ one.SetOne();
+ n.Sub(one);
+ d.Add(one);
+ n.Div(d);
+
+ result = ATan0(n);
+
+ n.Set05Pi();
+ n.exponent.SubOne(); // =pi/4
+ result.Add(n);
+
+ return result;
+ }
+
+
+ /*!
+ an auxiliary function for calculating the Arc Tangent
+ where x > 1
+
+ we're using the formula:
+ atan(x) = pi/2 - atan(1/x) for x>0
+ */
+ template
+ ValueType ATanGreaterThanPlusOne(const ValueType & x)
+ {
+ ValueType temp, atan;
+
+ temp.SetOne();
+
+ if( temp.Div(x) )
+ {
+ // if there was a carry here that means x is very big
+ // and atan(1/x) fast converged to 0
+ atan.SetZero();
+ }
+ else
+ atan = ATan01(temp);
+
+ temp.Set05Pi();
+ temp.Sub(atan);
+
+ return temp;
+ }
+
+ } // namespace auxiliaryfunctions
+
+
+ /*!
+ this function calculates the Arc Tangent
+ */
+ template
+ ValueType ATan(ValueType x)
+ {
+ using namespace auxiliaryfunctions;
+
+ ValueType one, result;
+ one.SetOne();
+ bool change_sign = false;
+
+ if( x.IsNan() )
+ return result; // NaN is set by default
+
+ // if x is negative we're using the formula:
+ // atan(-x) = -atan(x)
+ if( x.IsSign() )
+ {
+ change_sign = true;
+ x.Abs();
+ }
+
+ if( x.GreaterWithoutSignThan(one) )
+ result = ATanGreaterThanPlusOne(x);
+ else
+ result = ATan01(x);
+
+ if( change_sign )
+ result.ChangeSign();
+
+ return result;
+ }
+
+
+ /*!
+ this function calculates the Arc Tangent
+ look at the description of ATan(...)
+
+ (the abbreviation of Arc Tangent can be 'atg' as well)
+ */
+ template
+ ValueType ATg(const ValueType & x)
+ {
+ return ATan(x);
+ }
+
+
+ /*!
+ this function calculates the Arc Cotangent
+
+ we're using the formula:
+ actan(x) = pi/2 - atan(x)
+ */
+ template
+ ValueType ACot(const ValueType & x)
+ {
+ ValueType result;
+
+ result.Set05Pi();
+ result.Sub(ATan(x));
+
+ return result;
+ }
+
+
+ /*!
+ this function calculates the Arc Cotangent
+ look at the description of ACot(...)
+
+ (the abbreviation of Arc Cotangent can be 'actg' as well)
+ */
+ template
+ ValueType ACtg(const ValueType & x)
+ {
+ return ACot(x);
+ }
+
+
+ /*
+ *
+ * hyperbolic functions
+ *
+ *
+ */
+
+
+ /*!
+ this function calculates the Hyperbolic Sine
+
+ we're using the formula sinh(x)= ( e^x - e^(-x) ) / 2
+ */
+ template
+ ValueType Sinh(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType ex, emx;
+ uint c = 0;
+
+ c += ex.Exp(x);
+ c += emx.Exp(-x);
+
+ c += ex.Sub(emx);
+ c += ex.exponent.SubOne();
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return ex;
+ }
+
+
+ /*!
+ this function calculates the Hyperbolic Cosine
+
+ we're using the formula cosh(x)= ( e^x + e^(-x) ) / 2
+ */
+ template
+ ValueType Cosh(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType ex, emx;
+ uint c = 0;
+
+ c += ex.Exp(x);
+ c += emx.Exp(-x);
+
+ c += ex.Add(emx);
+ c += ex.exponent.SubOne();
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return ex;
+ }
+
+
+ /*!
+ this function calculates the Hyperbolic Tangent
+
+ we're using the formula tanh(x)= ( e^x - e^(-x) ) / ( e^x + e^(-x) )
+ */
+ template
+ ValueType Tanh(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType ex, emx, nominator, denominator;
+ uint c = 0;
+
+ c += ex.Exp(x);
+ c += emx.Exp(-x);
+
+ nominator = ex;
+ c += nominator.Sub(emx);
+ denominator = ex;
+ c += denominator.Add(emx);
+
+ c += nominator.Div(denominator);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return nominator;
+ }
+
+
+ /*!
+ this function calculates the Hyperbolic Tangent
+ look at the description of Tanh(...)
+
+ (the abbreviation of Hyperbolic Tangent can be 'tgh' as well)
+ */
+ template
+ ValueType Tgh(const ValueType & x, ErrorCode * err = 0)
+ {
+ return Tanh(x, err);
+ }
+
+ /*!
+ this function calculates the Hyperbolic Cotangent
+
+ we're using the formula coth(x)= ( e^x + e^(-x) ) / ( e^x - e^(-x) )
+ */
+ template
+ ValueType Coth(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ if( x.IsZero() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return ValueType(); // NaN is set by default
+ }
+
+ ValueType ex, emx, nominator, denominator;
+ uint c = 0;
+
+ c += ex.Exp(x);
+ c += emx.Exp(-x);
+
+ nominator = ex;
+ c += nominator.Add(emx);
+ denominator = ex;
+ c += denominator.Sub(emx);
+
+ c += nominator.Div(denominator);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return nominator;
+ }
+
+
+ /*!
+ this function calculates the Hyperbolic Cotangent
+ look at the description of Coth(...)
+
+ (the abbreviation of Hyperbolic Cotangent can be 'ctgh' as well)
+ */
+ template
+ ValueType Ctgh(const ValueType & x, ErrorCode * err = 0)
+ {
+ return Coth(x, err);
+ }
+
+
+ /*
+ *
+ * inverse hyperbolic functions
+ *
+ *
+ */
+
+
+ /*!
+ inverse hyperbolic sine
+
+ asinh(x) = ln( x + sqrt(x^2 + 1) )
+ */
+ template
+ ValueType ASinh(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType xx(x), one, result;
+ uint c = 0;
+ one.SetOne();
+
+ c += xx.Mul(x);
+ c += xx.Add(one);
+ one.exponent.SubOne(); // one=0.5
+ // xx is >= 1
+ c += xx.PowFrac(one); // xx=sqrt(xx)
+ c += xx.Add(x);
+ c += result.Ln(xx); // xx > 0
+
+ // here can only be a carry
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ inverse hyperbolic cosine
+
+ acosh(x) = ln( x + sqrt(x^2 - 1) ) x in <1, infinity)
+ */
+ template
+ ValueType ACosh(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType xx(x), one, result;
+ uint c = 0;
+ one.SetOne();
+
+ if( x < one )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ c += xx.Mul(x);
+ c += xx.Sub(one);
+ // xx is >= 0
+ // we can't call a PowFrac when the 'x' is zero
+ // if x is 0 the sqrt(0) is 0
+ if( !xx.IsZero() )
+ {
+ one.exponent.SubOne(); // one=0.5
+ c += xx.PowFrac(one); // xx=sqrt(xx)
+ }
+ c += xx.Add(x);
+ c += result.Ln(xx); // xx >= 1
+
+ // here can only be a carry
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ inverse hyperbolic tangent
+
+ atanh(x) = 0.5 * ln( (1+x) / (1-x) ) x in (-1, 1)
+ */
+ template
+ ValueType ATanh(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType nominator(x), denominator, one, result;
+ uint c = 0;
+ one.SetOne();
+
+ if( !x.SmallerWithoutSignThan(one) )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ c += nominator.Add(one);
+ denominator = one;
+ c += denominator.Sub(x);
+ c += nominator.Div(denominator);
+ c += result.Ln(nominator);
+ c += result.exponent.SubOne();
+
+ // here can only be a carry
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ inverse hyperbolic tantent
+ */
+ template
+ ValueType ATgh(const ValueType & x, ErrorCode * err = 0)
+ {
+ return ATanh(x, err);
+ }
+
+
+ /*!
+ inverse hyperbolic cotangent
+
+ acoth(x) = 0.5 * ln( (x+1) / (x-1) ) x in (-infinity, -1) or (1, infinity)
+ */
+ template
+ ValueType ACoth(const ValueType & x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return x; // NaN
+ }
+
+ ValueType nominator(x), denominator(x), one, result;
+ uint c = 0;
+ one.SetOne();
+
+ if( !x.GreaterWithoutSignThan(one) )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ c += nominator.Add(one);
+ c += denominator.Sub(one);
+ c += nominator.Div(denominator);
+ c += result.Ln(nominator);
+ c += result.exponent.SubOne();
+
+ // here can only be a carry
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ inverse hyperbolic cotantent
+ */
+ template
+ ValueType ACtgh(const ValueType & x, ErrorCode * err = 0)
+ {
+ return ACoth(x, err);
+ }
+
+
+
+
+
+ /*
+ *
+ * functions for converting between degrees, radians and gradians
+ *
+ *
+ */
+
+
+ /*!
+ this function converts degrees to radians
+
+ it returns: x * pi / 180
+ */
+ template
+ ValueType DegToRad(const ValueType & x, ErrorCode * err = 0)
+ {
+ ValueType result, temp;
+ uint c = 0;
+
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ result = x;
+
+ // it is better to make division first and then multiplication
+ // the result is more accurate especially when x is: 90,180,270 or 360
+ temp = 180;
+ c += result.Div(temp);
+
+ temp.SetPi();
+ c += result.Mul(temp);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ this function converts radians to degrees
+
+ it returns: x * 180 / pi
+ */
+ template
+ ValueType RadToDeg(const ValueType & x, ErrorCode * err = 0)
+ {
+ ValueType result, delimiter;
+ uint c = 0;
+
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ result = 180;
+ c += result.Mul(x);
+
+ delimiter.SetPi();
+ c += result.Div(delimiter);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ this function converts degrees in the long format into one value
+
+ long format: (degrees, minutes, seconds)
+ minutes and seconds must be greater than or equal zero
+
+ result:
+ if d>=0 : result= d + ((s/60)+m)/60
+ if d<0 : result= d - ((s/60)+m)/60
+
+ ((s/60)+m)/60 = (s+60*m)/3600 (second version is faster because
+ there's only one division)
+
+ for example:
+ DegToDeg(10, 30, 0) = 10.5
+ DegToDeg(10, 24, 35.6)=10.4098(8)
+ */
+ template
+ ValueType DegToDeg( const ValueType & d, const ValueType & m, const ValueType & s,
+ ErrorCode * err = 0)
+ {
+ ValueType delimiter, multipler;
+ uint c = 0;
+
+ if( d.IsNan() || m.IsNan() || s.IsNan() || m.IsSign() || s.IsSign() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return delimiter ; // NaN is set by default
+ }
+
+ multipler = 60;
+ delimiter = 3600;
+
+ c += multipler.Mul(m);
+ c += multipler.Add(s);
+ c += multipler.Div(delimiter);
+
+ if( d.IsSign() )
+ multipler.ChangeSign();
+
+ c += multipler.Add(d);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return multipler;
+ }
+
+
+ /*!
+ this function converts degrees in the long format to radians
+ */
+ template
+ ValueType DegToRad( const ValueType & d, const ValueType & m, const ValueType & s,
+ ErrorCode * err = 0)
+ {
+ ValueType temp_deg = DegToDeg(d,m,s,err);
+
+ if( err && *err!=err_ok )
+ return temp_deg;
+
+ return DegToRad(temp_deg, err);
+ }
+
+
+ /*!
+ this function converts gradians to radians
+
+ it returns: x * pi / 200
+ */
+ template
+ ValueType GradToRad(const ValueType & x, ErrorCode * err = 0)
+ {
+ ValueType result, temp;
+ uint c = 0;
+
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ result = x;
+
+ // it is better to make division first and then multiplication
+ // the result is more accurate especially when x is: 100,200,300 or 400
+ temp = 200;
+ c += result.Div(temp);
+
+ temp.SetPi();
+ c += result.Mul(temp);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ this function converts radians to gradians
+
+ it returns: x * 200 / pi
+ */
+ template
+ ValueType RadToGrad(const ValueType & x, ErrorCode * err = 0)
+ {
+ ValueType result, delimiter;
+ uint c = 0;
+
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ result = 200;
+ c += result.Mul(x);
+
+ delimiter.SetPi();
+ c += result.Div(delimiter);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ this function converts degrees to gradians
+
+ it returns: x * 200 / 180
+ */
+ template
+ ValueType DegToGrad(const ValueType & x, ErrorCode * err = 0)
+ {
+ ValueType result, temp;
+ uint c = 0;
+
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ result = x;
+
+ temp = 200;
+ c += result.Mul(temp);
+
+ temp = 180;
+ c += result.Div(temp);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+ /*!
+ this function converts degrees in the long format to gradians
+ */
+ template
+ ValueType DegToGrad( const ValueType & d, const ValueType & m, const ValueType & s,
+ ErrorCode * err = 0)
+ {
+ ValueType temp_deg = DegToDeg(d,m,s,err);
+
+ if( err && *err!=err_ok )
+ return temp_deg;
+
+ return DegToGrad(temp_deg, err);
+ }
+
+
+ /*!
+ this function converts degrees to gradians
+
+ it returns: x * 180 / 200
+ */
+ template
+ ValueType GradToDeg(const ValueType & x, ErrorCode * err = 0)
+ {
+ ValueType result, temp;
+ uint c = 0;
+
+ if( x.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ result = x;
+
+ temp = 180;
+ c += result.Mul(temp);
+
+ temp = 200;
+ c += result.Div(temp);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return result;
+ }
+
+
+
+
+ /*
+ *
+ * another functions
+ *
+ *
+ */
+
+
+ /*!
+ this function calculates the square root
+
+ Sqrt(9) = 3
+ */
+ template
+ ValueType Sqrt(ValueType x, ErrorCode * err = 0)
+ {
+ if( x.IsNan() || x.IsSign() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return ValueType(); // NaN is set by default
+ }
+
+ uint c = x.Sqrt();
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return x;
+ }
+
+
+
+ namespace auxiliaryfunctions
+ {
+
+ template
+ bool RootCheckIndexSign(ValueType & x, const ValueType & index, ErrorCode * err)
+ {
+ if( index.IsSign() )
+ {
+ // index cannot be negative
+ if( err )
+ *err = err_improper_argument;
+
+ x.SetNan();
+
+ return true;
+ }
+
+ return false;
+ }
+
+
+ template
+ bool RootCheckIndexZero(ValueType & x, const ValueType & index, ErrorCode * err)
+ {
+ if( index.IsZero() )
+ {
+ if( x.IsZero() )
+ {
+ // there isn't root(0;0) - we assume it's not defined
+ if( err )
+ *err = err_improper_argument;
+
+ x.SetNan();
+
+ return true;
+ }
+
+ // root(x;0) is 1 (if x!=0)
+ x.SetOne();
+
+ if( err )
+ *err = err_ok;
+
+ return true;
+ }
+
+ return false;
+ }
+
+
+ template
+ bool RootCheckIndexOne(const ValueType & index, ErrorCode * err)
+ {
+ ValueType one;
+ one.SetOne();
+
+ if( index == one )
+ {
+ //root(x;1) is x
+ // we do it because if we used the PowFrac function
+ // we would lose the precision
+ if( err )
+ *err = err_ok;
+
+ return true;
+ }
+
+ return false;
+ }
+
+
+ template
+ bool RootCheckIndexTwo(ValueType & x, const ValueType & index, ErrorCode * err)
+ {
+ if( index == 2 )
+ {
+ x = Sqrt(x, err);
+
+ return true;
+ }
+
+ return false;
+ }
+
+
+ template
+ bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err)
+ {
+ if( !index.IsInteger() )
+ {
+ // index must be integer
+ if( err )
+ *err = err_improper_argument;
+
+ x.SetNan();
+
+ return true;
+ }
+
+ return false;
+ }
+
+
+ template
+ bool RootCheckXZero(ValueType & x, ErrorCode * err)
+ {
+ if( x.IsZero() )
+ {
+ // root(0;index) is zero (if index!=0)
+ // RootCheckIndexZero() must be called beforehand
+ x.SetZero();
+
+ if( err )
+ *err = err_ok;
+
+ return true;
+ }
+
+ return false;
+ }
+
+
+ template
+ bool RootCheckIndex(ValueType & x, const ValueType & index, ErrorCode * err, bool * change_sign)
+ {
+ *change_sign = false;
+
+ if( index.Mod2() )
+ {
+ // index is odd (1,3,5...)
+ if( x.IsSign() )
+ {
+ *change_sign = true;
+ x.Abs();
+ }
+ }
+ else
+ {
+ // index is even
+ // x cannot be negative
+ if( x.IsSign() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ x.SetNan();
+
+ return true;
+ }
+ }
+
+ return false;
+ }
+
+
+ template
+ uint RootCorrectInteger(ValueType & old_x, ValueType & x, const ValueType & index)
+ {
+ if( !old_x.IsInteger() || x.IsInteger() || !index.exponent.IsSign() )
+ return 0;
+
+ // old_x is integer,
+ // x is not integer,
+ // index is relatively small (index.exponent<0 or index.exponent<=0)
+ // (because we're using a special powering algorithm Big::PowUInt())
+
+ uint c = 0;
+
+ ValueType temp(x);
+ c += temp.Round();
+
+ ValueType temp_round(temp);
+ c += temp.PowUInt(index);
+
+ if( temp == old_x )
+ x = temp_round;
+
+ return (c==0)? 0 : 1;
+ }
+
+
+
+ } // namespace auxiliaryfunctions
+
+
+
+ /*!
+ indexth Root of x
+ index must be integer and not negative <0;1;2;3....)
+
+ if index==0 the result is one
+ if x==0 the result is zero and we assume root(0;0) is not defined
+
+ if index is even (2;4;6...) the result is x^(1/index) and x>0
+ if index is odd (1;2;3;...) the result is either
+ -(abs(x)^(1/index)) if x<0 or
+ x^(1/index)) if x>0
+
+ (for index==1 the result is equal x)
+ */
+ template
+ ValueType Root(ValueType x, const ValueType & index, ErrorCode * err = 0)
+ {
+ using namespace auxiliaryfunctions;
+
+ if( x.IsNan() || index.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return ValueType(); // NaN is set by default
+ }
+
+ if( RootCheckIndexSign(x, index, err) ) return x;
+ if( RootCheckIndexZero(x, index, err) ) return x;
+ if( RootCheckIndexOne ( index, err) ) return x;
+ if( RootCheckIndexTwo (x, index, err) ) return x;
+ if( RootCheckIndexFrac(x, index, err) ) return x;
+ if( RootCheckXZero (x, err) ) return x;
+
+ // index integer and index!=0
+ // x!=0
+
+ ValueType old_x(x);
+ bool change_sign;
+
+ if( RootCheckIndex(x, index, err, &change_sign ) ) return x;
+
+ ValueType temp;
+ uint c = 0;
+
+ // we're using the formula: root(x ; n) = exp( ln(x) / n )
+ c += temp.Ln(x);
+ c += temp.Div(index);
+ c += x.Exp(temp);
+
+ if( change_sign )
+ {
+ // x is different from zero
+ x.SetSign();
+ }
+
+ c += RootCorrectInteger(old_x, x, index);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return x;
+ }
+
+
+
+ /*!
+ absolute value of x
+ e.g. -2 = 2
+ 2 = 2
+ */
+ template
+ ValueType Abs(const ValueType & x)
+ {
+ ValueType result( x );
+ result.Abs();
+
+ return result;
+ }
+
+
+ /*!
+ it returns the sign of the value
+ e.g. -2 = -1
+ 0 = 0
+ 10 = 1
+ */
+ template
+ ValueType Sgn(ValueType x)
+ {
+ x.Sgn();
+
+ return x;
+ }
+
+
+ /*!
+ the remainder from a division
+
+ e.g.
+ mod( 12.6 ; 3) = 0.6 because 12.6 = 3*4 + 0.6
+ mod(-12.6 ; 3) = -0.6 bacause -12.6 = 3*(-4) + (-0.6)
+ mod( 12.6 ; -3) = 0.6
+ mod(-12.6 ; -3) = -0.6
+ */
+ template
+ ValueType Mod(ValueType a, const ValueType & b, ErrorCode * err = 0)
+ {
+ if( a.IsNan() || b.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return ValueType(); // NaN is set by default
+ }
+
+ uint c = a.Mod(b);
+
+ if( err )
+ *err = c ? err_overflow : err_ok;
+
+ return a;
+ }
+
+
+
+ namespace auxiliaryfunctions
+ {
+
+ /*!
+ this function is used to store factorials in a given container
+ 'more' means how many values should be added at the end
+
+ e.g.
+ std::vector fact;
+ SetFactorialSequence(fact, 3);
+ // now the container has three values: 1 1 2
+
+ SetFactorialSequence(fact, 2);
+ // now the container has five values: 1 1 2 6 24
+ */
+ template
+ void SetFactorialSequence(std::vector & fact, uint more = 20)
+ {
+ if( more == 0 )
+ more = 1;
+
+ uint start = static_cast(fact.size());
+ fact.resize(fact.size() + more);
+
+ if( start == 0 )
+ {
+ fact[0] = 1;
+ ++start;
+ }
+
+ for(uint i=start ; i
+ ValueType SetBernoulliNumbersSum(CGamma & cgamma, const ValueType & n_, uint m,
+ const volatile StopCalculating * stop = 0)
+ {
+ ValueType k_, temp, temp2, temp3, sum;
+
+ sum.SetZero();
+
+ for(uint k=0 ; kWasStopSignal() )
+ return ValueType(); // NaN
+
+ if( k>1 && (k & 1) == 1 ) // for that k the Bernoulli number is zero
+ continue;
+
+ k_ = k;
+
+ temp = n_; // n_ is equal 2
+ temp.Pow(k_);
+ // temp = 2^k
+
+ temp2 = cgamma.fact[m];
+ temp3 = cgamma.fact[k];
+ temp3.Mul(cgamma.fact[m-k]);
+ temp2.Div(temp3);
+ // temp2 = (m k) = m! / ( k! * (m-k)! )
+
+ temp.Mul(temp2);
+ temp.Mul(cgamma.bern[k]);
+
+ sum.Add(temp);
+ // sum += 2^k * (m k) * B(k)
+
+ if( sum.IsNan() )
+ break;
+ }
+
+ return sum;
+ }
+
+
+ /*!
+ an auxiliary function used to calculate Bernoulli numbers
+ start is >= 2
+
+ we use the recurrence formula:
+ B(m) = 1 / (2*(1 - 2^m)) * sum(m)
+ where sum(m) is calculated by SetBernoulliNumbersSum()
+ */
+ template
+ bool SetBernoulliNumbersMore(CGamma & cgamma, uint start, const volatile StopCalculating * stop = 0)
+ {
+ ValueType denominator, temp, temp2, temp3, m_, sum, sum2, n_, k_;
+
+ const uint n = 2;
+ n_ = n;
+
+ // start is >= 2
+ for(uint m=start ; mWasStopSignal() )
+ {
+ cgamma.bern.resize(m); // valid numbers are in [0, m-1]
+ return false;
+ }
+
+ cgamma.bern[m].Div(denominator);
+ }
+ }
+
+ return true;
+ }
+
+
+ /*!
+ this function is used to calculate Bernoulli numbers,
+ returns false if there was a stop signal,
+ 'more' means how many values should be added at the end
+
+ e.g.
+ typedef Big<1,2> MyBig;
+ CGamma cgamma;
+ SetBernoulliNumbers(cgamma, 3);
+ // now we have three first Bernoulli numbers: 1 -0.5 0.16667
+
+ SetBernoulliNumbers(cgamma, 4);
+ // now we have 7 Bernoulli numbers: 1 -0.5 0.16667 0 -0.0333 0 0.0238
+ */
+ template
+ bool SetBernoulliNumbers(CGamma & cgamma, uint more = 20, const volatile StopCalculating * stop = 0)
+ {
+ if( more == 0 )
+ more = 1;
+
+ uint start = static_cast(cgamma.bern.size());
+ cgamma.bern.resize(cgamma.bern.size() + more);
+
+ if( start == 0 )
+ {
+ cgamma.bern[0].SetOne();
+ ++start;
+ }
+
+ if( cgamma.bern.size() == 1 )
+ return true;
+
+ if( start == 1 )
+ {
+ cgamma.bern[1].Set05();
+ cgamma.bern[1].ChangeSign();
+ ++start;
+ }
+
+ // we should have sufficient factorials in cgamma.fact
+ if( cgamma.fact.size() < cgamma.bern.size() )
+ SetFactorialSequence(cgamma.fact, static_cast(cgamma.bern.size() - cgamma.fact.size()));
+
+
+ return SetBernoulliNumbersMore(cgamma, start, stop);
+ }
+
+
+ /*!
+ an auxiliary function used to calculate the Gamma() function
+
+ we calculate a sum:
+ sum(n) = sum_{m=2} { B(m) / ( (m^2 - m) * n^(m-1) ) } = 1/(12*n) - 1/(360*n^3) + 1/(1260*n^5) + ...
+ B(m) means a mth Bernoulli number
+ the sum starts from m=2, we calculate as long as the value will not change after adding a next part
+ */
+ template
+ ValueType GammaFactorialHighSum(const ValueType & n, CGamma & cgamma, ErrorCode & err,
+ const volatile StopCalculating * stop)
+ {
+ ValueType temp, temp2, denominator, sum, oldsum;
+
+ sum.SetZero();
+
+ for(uint m=2 ; mWasStopSignal() )
+ {
+ err = err_interrupt;
+ return ValueType(); // NaN
+ }
+
+ temp = (m-1);
+ denominator = n;
+ denominator.Pow(temp);
+ // denominator = n ^ (m-1)
+
+ temp = m;
+ temp2 = temp;
+ temp.Mul(temp2);
+ temp.Sub(temp2);
+ // temp = m^2 - m
+
+ denominator.Mul(temp);
+ // denominator = (m^2 - m) * n ^ (m-1)
+
+ if( m >= cgamma.bern.size() )
+ {
+ if( !SetBernoulliNumbers(cgamma, m - cgamma.bern.size() + 1 + 3, stop) ) // 3 more than needed
+ {
+ // there was the stop signal
+ err = err_interrupt;
+ return ValueType(); // NaN
+ }
+ }
+
+ temp = cgamma.bern[m];
+ temp.Div(denominator);
+
+ oldsum = sum;
+ sum.Add(temp);
+
+ if( sum.IsNan() || oldsum==sum )
+ break;
+ }
+
+ return sum;
+ }
+
+
+ /*!
+ an auxiliary function used to calculate the Gamma() function
+
+ we calculate a helper function GammaFactorialHigh() by using Stirling's series:
+ n! = (n/e)^n * sqrt(2*pi*n) * exp( sum(n) )
+ where n is a real number (not only an integer) and is sufficient large (greater than TTMATH_GAMMA_BOUNDARY)
+ and sum(n) is calculated by GammaFactorialHighSum()
+ */
+ template
+ ValueType GammaFactorialHigh(const ValueType & n, CGamma & cgamma, ErrorCode & err,
+ const volatile StopCalculating * stop)
+ {
+ ValueType temp, temp2, temp3, denominator, sum;
+
+ temp.Set2Pi();
+ temp.Mul(n);
+ temp2 = Sqrt(temp);
+ // temp2 = sqrt(2*pi*n)
+
+ temp = n;
+ temp3.SetE();
+ temp.Div(temp3);
+ temp.Pow(n);
+ // temp = (n/e)^n
+
+ sum = GammaFactorialHighSum(n, cgamma, err, stop);
+ temp3.Exp(sum);
+ // temp3 = exp(sum)
+
+ temp.Mul(temp2);
+ temp.Mul(temp3);
+
+ return temp;
+ }
+
+
+ /*!
+ an auxiliary function used to calculate the Gamma() function
+
+ Gamma(x) = GammaFactorialHigh(x-1)
+ */
+ template
+ ValueType GammaPlusHigh(ValueType n, CGamma & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
+ {
+ ValueType one;
+
+ one.SetOne();
+ n.Sub(one);
+
+ return GammaFactorialHigh(n, cgamma, err, stop);
+ }
+
+
+ /*!
+ an auxiliary function used to calculate the Gamma() function
+
+ we use this function when n is integer and a small value (from 0 to TTMATH_GAMMA_BOUNDARY]
+ we use the formula:
+ gamma(n) = (n-1)! = 1 * 2 * 3 * ... * (n-1)
+ */
+ template
+ ValueType GammaPlusLowIntegerInt(uint n, CGamma & cgamma)
+ {
+ TTMATH_ASSERT( n > 0 )
+
+ if( n - 1 < static_cast(cgamma.fact.size()) )
+ return cgamma.fact[n - 1];
+
+ ValueType res;
+ uint start = 2;
+
+ if( cgamma.fact.size() < 2 )
+ {
+ res.SetOne();
+ }
+ else
+ {
+ start = static_cast(cgamma.fact.size());
+ res = cgamma.fact[start-1];
+ }
+
+ for(uint i=start ; i
+ ValueType GammaPlusLowInteger(const ValueType & n, CGamma & cgamma)
+ {
+ sint n_;
+
+ n.ToInt(n_);
+
+ return GammaPlusLowIntegerInt(n_, cgamma);
+ }
+
+
+ /*!
+ an auxiliary function used to calculate the Gamma() function
+
+ we use this function when n is a small value (from 0 to TTMATH_GAMMA_BOUNDARY]
+ we use a recurrence formula:
+ gamma(z+1) = z * gamma(z)
+ then: gamma(z) = gamma(z+1) / z
+
+ e.g.
+ gamma(3.89) = gamma(2001.89) / ( 3.89 * 4.89 * 5.89 * ... * 1999.89 * 2000.89 )
+ */
+ template
+ ValueType GammaPlusLow(ValueType n, CGamma & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
+ {
+ ValueType one, denominator, temp, boundary;
+
+ if( n.IsInteger() )
+ return GammaPlusLowInteger(n, cgamma);
+
+ one.SetOne();
+ denominator = n;
+ boundary = TTMATH_GAMMA_BOUNDARY;
+
+ while( n < boundary )
+ {
+ n.Add(one);
+ denominator.Mul(n);
+ }
+
+ n.Add(one);
+
+ // now n is sufficient big
+ temp = GammaPlusHigh(n, cgamma, err, stop);
+ temp.Div(denominator);
+
+ return temp;
+ }
+
+
+ /*!
+ an auxiliary function used to calculate the Gamma() function
+ */
+ template
+ ValueType GammaPlus(const ValueType & n, CGamma & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
+ {
+ if( n > TTMATH_GAMMA_BOUNDARY )
+ return GammaPlusHigh(n, cgamma, err, stop);
+
+ return GammaPlusLow(n, cgamma, err, stop);
+ }
+
+
+ /*!
+ an auxiliary function used to calculate the Gamma() function
+
+ this function is used when n is negative
+ we use the reflection formula:
+ gamma(1-z) * gamma(z) = pi / sin(pi*z)
+ then: gamma(z) = pi / (sin(pi*z) * gamma(1-z))
+
+ */
+ template
+ ValueType GammaMinus(const ValueType & n, CGamma & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
+ {
+ ValueType pi, denominator, temp, temp2;
+
+ if( n.IsInteger() )
+ {
+ // gamma function is not defined when n is negative and integer
+ err = err_improper_argument;
+ return temp; // NaN
+ }
+
+ pi.SetPi();
+
+ temp = pi;
+ temp.Mul(n);
+ temp2 = Sin(temp);
+ // temp2 = sin(pi * n)
+
+ temp.SetOne();
+ temp.Sub(n);
+ temp = GammaPlus(temp, cgamma, err, stop);
+ // temp = gamma(1 - n)
+
+ temp.Mul(temp2);
+ pi.Div(temp);
+
+ return pi;
+ }
+
+ } // namespace auxiliaryfunctions
+
+
+
+ /*!
+ this function calculates the Gamma function
+
+ it's multithread safe, you should create a CGamma<> object and use it whenever you call the Gamma()
+ e.g.
+ typedef Big<1,2> MyBig;
+ MyBig x=234, y=345.53;
+ CGamma cgamma;
+ std::cout << Gamma(x, cgamma) << std::endl;
+ std::cout << Gamma(y, cgamma) << std::endl;
+ in the CGamma<> object the function stores some coefficients (factorials, Bernoulli numbers),
+ and they will be reused in next calls to the function
+
+ each thread should have its own CGamma<> object, and you can use these objects with Factorial() function too
+ */
+ template
+ ValueType Gamma(const ValueType & n, CGamma & cgamma, ErrorCode * err = 0,
+ const volatile StopCalculating * stop = 0)
+ {
+ using namespace auxiliaryfunctions;
+
+ ValueType result;
+ ErrorCode err_tmp;
+
+ if( n.IsNan() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN is set by default
+ }
+
+ if( cgamma.history.Get(n, result, err_tmp) )
+ {
+ if( err )
+ *err = err_tmp;
+
+ return result;
+ }
+
+ err_tmp = err_ok;
+
+ if( n.IsSign() )
+ {
+ result = GammaMinus(n, cgamma, err_tmp, stop);
+ }
+ else
+ if( n.IsZero() )
+ {
+ err_tmp = err_improper_argument;
+ result.SetNan();
+ }
+ else
+ {
+ result = GammaPlus(n, cgamma, err_tmp, stop);
+ }
+
+ if( result.IsNan() && err_tmp==err_ok )
+ err_tmp = err_overflow;
+
+ if( err )
+ *err = err_tmp;
+
+ if( stop && !stop->WasStopSignal() )
+ cgamma.history.Add(n, result, err_tmp);
+
+ return result;
+ }
+
+
+ /*!
+ this function calculates the Gamma function
+
+ note: this function should be used only in a single-thread environment
+ */
+ template
+ ValueType Gamma(const ValueType & n, ErrorCode * err = 0)
+ {
+ // warning: this static object is not thread safe
+ static CGamma cgamma;
+
+ return Gamma(n, cgamma, err);
+ }
+
+
+
+ namespace auxiliaryfunctions
+ {
+
+ /*!
+ an auxiliary function for calculating the factorial function
+
+ we use the formula:
+ x! = gamma(x+1)
+ */
+ template
+ ValueType Factorial2(ValueType x,
+ CGamma * cgamma = 0,
+ ErrorCode * err = 0,
+ const volatile StopCalculating * stop = 0)
+ {
+ ValueType result, one;
+
+ if( x.IsNan() || x.IsSign() || !x.IsInteger() )
+ {
+ if( err )
+ *err = err_improper_argument;
+
+ return result; // NaN set by default
+ }
+
+ one.SetOne();
+ x.Add(one);
+
+ if( cgamma )
+ return Gamma(x, *cgamma, err, stop);
+
+ return Gamma(x, err);
+ }
+
+ } // namespace auxiliaryfunctions
+
+
+
+ /*!
+ the factorial from given 'x'
+ e.g.
+ Factorial(4) = 4! = 1*2*3*4
+
+ it's multithread safe, you should create a CGamma<> object and use it whenever you call the Factorial()
+ e.g.
+ typedef Big<1,2> MyBig;
+ MyBig x=234, y=54345;
+ CGamma cgamma;
+ std::cout << Factorial(x, cgamma) << std::endl;
+ std::cout << Factorial(y, cgamma) << std::endl;
+ in the CGamma<> object the function stores some coefficients (factorials, Bernoulli numbers),
+ and they will be reused in next calls to the function
+
+ each thread should have its own CGamma<> object, and you can use these objects with Gamma() function too
+ */
+ template
+ ValueType Factorial(const ValueType & x, CGamma & cgamma, ErrorCode * err = 0,
+ const volatile StopCalculating * stop = 0)
+ {
+ return auxiliaryfunctions::Factorial2(x, &cgamma, err, stop);
+ }
+
+
+ /*!
+ the factorial from given 'x'
+ e.g.
+ Factorial(4) = 4! = 1*2*3*4
+
+ note: this function should be used only in a single-thread environment
+ */
+ template
+ ValueType Factorial(const ValueType & x, ErrorCode * err = 0)
+ {
+ return auxiliaryfunctions::Factorial2(x, (CGamma*)0, err, 0);
+ }
+
+
+ /*!
+ this method prepares some coefficients: factorials and Bernoulli numbers
+ stored in 'fact' and 'bern' objects
+
+ we're defining the method here because we're using Gamma() function which
+ is not available in ttmathobjects.h
+
+ read the doc info in ttmathobjects.h file where CGamma<> struct is declared
+ */
+ template
+ void CGamma::InitAll()
+ {
+ ValueType x = TTMATH_GAMMA_BOUNDARY + 1;
+
+ // history.Remove(x) removes only one object
+ // we must be sure that there are not others objects with the key 'x'
+ while( history.Remove(x) )
+ {
+ }
+
+ // the simplest way to initialize is to call the Gamma function with (TTMATH_GAMMA_BOUNDARY + 1)
+ // when x is larger then fewer coefficients we need
+ Gamma(x, *this);
+ }
+
+
+
+} // namespace
+
+
+/*!
+ this is for convenience for the user
+ he can only use '#include ' even if he uses the parser
+*/
+#include "ttmathparser.h"
+
+
+#ifdef _MSC_VER
+//warning C4127: conditional expression is constant
+#pragma warning( default: 4127 )
+//warning C4702: unreachable code
+#pragma warning( default: 4702 )
+//warning C4800: forcing value to bool 'true' or 'false' (performance warning)
+#pragma warning( default: 4800 )
+#endif
+
+#endif
diff --git a/ttmath/ttmathbig.h b/ttmath/ttmathbig.h
new file mode 100644
index 0000000..771f9c8
--- /dev/null
+++ b/ttmath/ttmathbig.h
@@ -0,0 +1,5236 @@
+/*
+ * This file is a part of TTMath Bignum Library
+ * and is distributed under the (new) BSD licence.
+ * Author: Tomasz Sowa
+ */
+
+/*
+ * Copyright (c) 2006-2010, Tomasz Sowa
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * * Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ *
+ * * Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * * Neither the name Tomasz Sowa nor the names of contributors to this
+ * project may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+ * THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef headerfilettmathbig
+#define headerfilettmathbig
+
+/*!
+ \file ttmathbig.h
+ \brief A Class for representing floating point numbers
+*/
+
+#include "ttmathint.h"
+#include "ttmaththreads.h"
+
+#include
+
+#ifdef TTMATH_MULTITHREADS
+#include
+#endif
+
+namespace ttmath
+{
+
+
+/*!
+ \brief Big implements the floating point numbers
+*/
+template
+class Big
+{
+
+/*
+ value = mantissa * 2^exponent
+
+ exponent - an integer value with a sign
+ mantissa - an integer value without a sing
+
+ mantissa must be pushed into the left side that is the highest bit from
+ mantissa must be one (of course if there's another value than zero) -- this job
+ (pushing bits into the left side) making Standardizing() method
+
+ for example:
+ if we want to store value one (1) into our Big object we must:
+ set mantissa to 1
+ set exponent to 0
+ set info to 0
+ and call method Standardizing()
+*/
+
+
+public:
+
+Int exponent;
+UInt mantissa;
+unsigned char info;
+
+
+/*!
+ Sign
+ the mask of a bit from 'info' which means that there is a sign
+ (when the bit is set)
+*/
+#define TTMATH_BIG_SIGN 128
+
+
+/*!
+ Not a number
+ if this bit is set that there is not a valid number
+*/
+#define TTMATH_BIG_NAN 64
+
+
+/*!
+ Zero
+ if this bit is set that there is value zero
+ mantissa should be zero and exponent should be zero too
+ (the Standardizing() method does this)
+*/
+#define TTMATH_BIG_ZERO 32
+
+
+ /*!
+ this method sets NaN if there was a carry (and returns 1 in such a case)
+
+ c can be 0, 1 or other value different from zero
+ */
+ uint CheckCarry(uint c)
+ {
+ if( c != 0 )
+ {
+ SetNan();
+ return 1;
+ }
+
+ return 0;
+ }
+
+public:
+
+
+ /*!
+ returning the string represents the currect type of the library
+ we have following types:
+ asm_vc_32 - with asm code designed for Microsoft Visual C++ (32 bits)
+ asm_gcc_32 - with asm code designed for GCC (32 bits)
+ asm_vc_64 - with asm for VC (64 bit)
+ asm_gcc_64 - with asm for GCC (64 bit)
+ no_asm_32 - pure C++ version (32 bit) - without any asm code
+ no_asm_64 - pure C++ version (64 bit) - without any asm code
+ */
+ static const char * LibTypeStr()
+ {
+ return UInt::LibTypeStr();
+ }
+
+
+ /*!
+ returning the currect type of the library
+ */
+ static LibTypeCode LibType()
+ {
+ return UInt::LibType();
+ }
+
+
+
+ /*!
+ this method moves all bits from mantissa into its left side
+ (suitably changes the exponent) or if the mantissa is zero
+ it sets the exponent to zero as well
+ (and clears the sign bit and sets the zero bit)
+
+ it can return a carry
+ the carry will be when we don't have enough space in the exponent
+
+ you don't have to use this method if you don't change the mantissa
+ and exponent directly
+ */
+ uint Standardizing()
+ {
+ if( mantissa.IsTheHighestBitSet() )
+ {
+ ClearInfoBit(TTMATH_BIG_ZERO);
+ return 0;
+ }
+
+ if( CorrectZero() )
+ return 0;
+
+ uint comp = mantissa.CompensationToLeft();
+
+ return exponent.Sub( comp );
+ }
+
+
+private:
+
+ /*!
+ if the mantissa is equal zero this method sets exponent to zero and
+ info without the sign
+
+ it returns true if there was the correction
+ */
+ bool CorrectZero()
+ {
+ if( mantissa.IsZero() )
+ {
+ SetInfoBit(TTMATH_BIG_ZERO);
+ ClearInfoBit(TTMATH_BIG_SIGN);
+ exponent.SetZero();
+
+ return true;
+ }
+ else
+ {
+ ClearInfoBit(TTMATH_BIG_ZERO);
+ }
+
+ return false;
+ }
+
+
+public:
+
+ /*!
+ this method clears a specific bit in the 'info' variable
+
+ bit is one of: TTMATH_BIG_SIGN, TTMATH_BIG_NAN etc.
+ */
+ void ClearInfoBit(unsigned char bit)
+ {
+ info = info & (~bit);
+ }
+
+
+ /*!
+ this method sets a specific bit in the 'info' variable
+
+ bit is one of: TTMATH_BIG_SIGN, TTMATH_BIG_NAN etc.
+
+ */
+ void SetInfoBit(unsigned char bit)
+ {
+ info = info | bit;
+ }
+
+
+ /*!
+ this method returns true if a specific bit in the 'info' variable is set
+
+ bit is one of: TTMATH_BIG_SIGN, TTMATH_BIG_NAN etc.
+ */
+ bool IsInfoBit(unsigned char bit) const
+ {
+ return (info & bit) != 0;
+ }
+
+
+ /*!
+ this method sets zero
+ */
+ void SetZero()
+ {
+ info = TTMATH_BIG_ZERO;
+ exponent.SetZero();
+ mantissa.SetZero();
+
+ /*
+ we don't have to compensate zero
+ */
+ }
+
+
+ /*!
+ this method sets one
+ */
+ void SetOne()
+ {
+ FromUInt(1);
+ }
+
+
+ /*!
+ this method sets value 0.5
+ */
+ void Set05()
+ {
+ FromUInt(1);
+ exponent.SubOne();
+ }
+
+
+ /*!
+ this method sets NaN flag (Not a Number)
+ when this flag is set that means there is no a valid number
+ */
+ void SetNan()
+ {
+ SetInfoBit(TTMATH_BIG_NAN);
+ }
+
+
+private:
+
+ /*!
+ this method sets the mantissa of the value of pi
+ */
+ void SetMantissaPi()
+ {
+ // this is a static table which represents the value of Pi (mantissa of it)
+ // (first is the highest word)
+ // we must define this table as 'unsigned int' because
+ // both on 32bit and 64bit platforms this table is 32bit
+ static const unsigned int temp_table[] = {
+ 0xc90fdaa2, 0x2168c234, 0xc4c6628b, 0x80dc1cd1, 0x29024e08, 0x8a67cc74, 0x020bbea6, 0x3b139b22,
+ 0x514a0879, 0x8e3404dd, 0xef9519b3, 0xcd3a431b, 0x302b0a6d, 0xf25f1437, 0x4fe1356d, 0x6d51c245,
+ 0xe485b576, 0x625e7ec6, 0xf44c42e9, 0xa637ed6b, 0x0bff5cb6, 0xf406b7ed, 0xee386bfb, 0x5a899fa5,
+ 0xae9f2411, 0x7c4b1fe6, 0x49286651, 0xece45b3d, 0xc2007cb8, 0xa163bf05, 0x98da4836, 0x1c55d39a,
+ 0x69163fa8, 0xfd24cf5f, 0x83655d23, 0xdca3ad96, 0x1c62f356, 0x208552bb, 0x9ed52907, 0x7096966d,
+ 0x670c354e, 0x4abc9804, 0xf1746c08, 0xca18217c, 0x32905e46, 0x2e36ce3b, 0xe39e772c, 0x180e8603,
+ 0x9b2783a2, 0xec07a28f, 0xb5c55df0, 0x6f4c52c9, 0xde2bcbf6, 0x95581718, 0x3995497c, 0xea956ae5,
+ 0x15d22618, 0x98fa0510, 0x15728e5a, 0x8aaac42d, 0xad33170d, 0x04507a33, 0xa85521ab, 0xdf1cba64,
+ 0xecfb8504, 0x58dbef0a, 0x8aea7157, 0x5d060c7d, 0xb3970f85, 0xa6e1e4c7, 0xabf5ae8c, 0xdb0933d7,
+ 0x1e8c94e0, 0x4a25619d, 0xcee3d226, 0x1ad2ee6b, 0xf12ffa06, 0xd98a0864, 0xd8760273, 0x3ec86a64,
+ 0x521f2b18, 0x177b200c, 0xbbe11757, 0x7a615d6c, 0x770988c0, 0xbad946e2, 0x08e24fa0, 0x74e5ab31,
+ 0x43db5bfc, 0xe0fd108e, 0x4b82d120, 0xa9210801, 0x1a723c12, 0xa787e6d7, 0x88719a10, 0xbdba5b26,
+ 0x99c32718, 0x6af4e23c, 0x1a946834, 0xb6150bda, 0x2583e9ca, 0x2ad44ce8, 0xdbbbc2db, 0x04de8ef9,
+ 0x2e8efc14, 0x1fbecaa6, 0x287c5947, 0x4e6bc05d, 0x99b2964f, 0xa090c3a2, 0x233ba186, 0x515be7ed,
+ 0x1f612970, 0xcee2d7af, 0xb81bdd76, 0x2170481c, 0xd0069127, 0xd5b05aa9, 0x93b4ea98, 0x8d8fddc1,
+ 0x86ffb7dc, 0x90a6c08f, 0x4df435c9, 0x34028492, 0x36c3fab4, 0xd27c7026, 0xc1d4dcb2, 0x602646de,
+ 0xc9751e76, 0x3dba37bd, 0xf8ff9406, 0xad9e530e, 0xe5db382f, 0x413001ae, 0xb06a53ed, 0x9027d831,
+ 0x179727b0, 0x865a8918, 0xda3edbeb, 0xcf9b14ed, 0x44ce6cba, 0xced4bb1b, 0xdb7f1447, 0xe6cc254b,
+ 0x33205151, 0x2bd7af42, 0x6fb8f401, 0x378cd2bf, 0x5983ca01, 0xc64b92ec, 0xf032ea15, 0xd1721d03,
+ 0xf482d7ce, 0x6e74fef6, 0xd55e702f, 0x46980c82, 0xb5a84031, 0x900b1c9e, 0x59e7c97f, 0xbec7e8f3,
+ 0x23a97a7e, 0x36cc88be, 0x0f1d45b7, 0xff585ac5, 0x4bd407b2, 0x2b4154aa, 0xcc8f6d7e, 0xbf48e1d8,
+ 0x14cc5ed2, 0x0f8037e0, 0xa79715ee, 0xf29be328, 0x06a1d58b, 0xb7c5da76, 0xf550aa3d, 0x8a1fbff0,
+ 0xeb19ccb1, 0xa313d55c, 0xda56c9ec, 0x2ef29632, 0x387fe8d7, 0x6e3c0468, 0x043e8f66, 0x3f4860ee,
+ 0x12bf2d5b, 0x0b7474d6, 0xe694f91e, 0x6dbe1159, 0x74a3926f, 0x12fee5e4, 0x38777cb6, 0xa932df8c,
+ 0xd8bec4d0, 0x73b931ba, 0x3bc832b6, 0x8d9dd300, 0x741fa7bf, 0x8afc47ed, 0x2576f693, 0x6ba42466,
+ 0x3aab639c, 0x5ae4f568, 0x3423b474, 0x2bf1c978, 0x238f16cb, 0xe39d652d, 0xe3fdb8be, 0xfc848ad9,
+ 0x22222e04, 0xa4037c07, 0x13eb57a8, 0x1a23f0c7, 0x3473fc64, 0x6cea306b, 0x4bcbc886, 0x2f8385dd,
+ 0xfa9d4b7f, 0xa2c087e8, 0x79683303, 0xed5bdd3a, 0x062b3cf5, 0xb3a278a6, 0x6d2a13f8, 0x3f44f82d,
+ 0xdf310ee0, 0x74ab6a36, 0x4597e899, 0xa0255dc1, 0x64f31cc5, 0x0846851d, 0xf9ab4819, 0x5ded7ea1,
+ 0xb1d510bd, 0x7ee74d73, 0xfaf36bc3, 0x1ecfa268, 0x359046f4, 0xeb879f92, 0x4009438b, 0x481c6cd7,
+ 0x889a002e, 0xd5ee382b, 0xc9190da6, 0xfc026e47, 0x9558e447, 0x5677e9aa, 0x9e3050e2, 0x765694df,
+ 0xc81f56e8, 0x80b96e71, 0x60c980dd, 0x98a573ea, 0x4472065a, 0x139cd290, 0x6cd1cb72, 0x9ec52a53 // last one was: 0x9ec52a52
+ //0x86d44014, ...
+ // (the last word 0x9ec52a52 was rounded up because the next one is 0x86d44014 -- first bit is one 0x8..)
+ // 256 32bit words for the mantissa -- about 2464 valid decimal digits
+ };
+ // the value of PI is comming from the website http://zenwerx.com/pi.php
+ // 3101 digits were taken from this website
+ // (later the digits were compared with:
+ // http://www.eveandersson.com/pi/digits/1000000 and http://www.geom.uiuc.edu/~huberty/math5337/groupe/digits.html )
+ // and they were set into Big<1,400> type (using operator=(const char*) on a 32bit platform)
+ // and then the first 256 words were taken into this table
+ // (TTMATH_BUILTIN_VARIABLES_SIZE on 32bit platform should have the value 256,
+ // and on 64bit platform value 128 (256/2=128))
+
+ mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(int));
+ }
+
+public:
+
+
+ /*!
+ this method sets the value of pi
+ */
+ void SetPi()
+ {
+ SetMantissaPi();
+ info = 0;
+ exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT) + 2;
+ }
+
+
+ /*!
+ this method sets the value of 0.5 * pi
+ */
+ void Set05Pi()
+ {
+ SetMantissaPi();
+ info = 0;
+ exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT) + 1;
+ }
+
+
+ /*!
+ this method sets the value of 2 * pi
+ */
+ void Set2Pi()
+ {
+ SetMantissaPi();
+ info = 0;
+ exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT) + 3;
+ }
+
+
+ /*!
+ this method sets the value of e
+ (the base of the natural logarithm)
+ */
+ void SetE()
+ {
+ static const unsigned int temp_table[] = {
+ 0xadf85458, 0xa2bb4a9a, 0xafdc5620, 0x273d3cf1, 0xd8b9c583, 0xce2d3695, 0xa9e13641, 0x146433fb,
+ 0xcc939dce, 0x249b3ef9, 0x7d2fe363, 0x630c75d8, 0xf681b202, 0xaec4617a, 0xd3df1ed5, 0xd5fd6561,
+ 0x2433f51f, 0x5f066ed0, 0x85636555, 0x3ded1af3, 0xb557135e, 0x7f57c935, 0x984f0c70, 0xe0e68b77,
+ 0xe2a689da, 0xf3efe872, 0x1df158a1, 0x36ade735, 0x30acca4f, 0x483a797a, 0xbc0ab182, 0xb324fb61,
+ 0xd108a94b, 0xb2c8e3fb, 0xb96adab7, 0x60d7f468, 0x1d4f42a3, 0xde394df4, 0xae56ede7, 0x6372bb19,
+ 0x0b07a7c8, 0xee0a6d70, 0x9e02fce1, 0xcdf7e2ec, 0xc03404cd, 0x28342f61, 0x9172fe9c, 0xe98583ff,
+ 0x8e4f1232, 0xeef28183, 0xc3fe3b1b, 0x4c6fad73, 0x3bb5fcbc, 0x2ec22005, 0xc58ef183, 0x7d1683b2,
+ 0xc6f34a26, 0xc1b2effa, 0x886b4238, 0x611fcfdc, 0xde355b3b, 0x6519035b, 0xbc34f4de, 0xf99c0238,
+ 0x61b46fc9, 0xd6e6c907, 0x7ad91d26, 0x91f7f7ee, 0x598cb0fa, 0xc186d91c, 0xaefe1309, 0x85139270,
+ 0xb4130c93, 0xbc437944, 0xf4fd4452, 0xe2d74dd3, 0x64f2e21e, 0x71f54bff, 0x5cae82ab, 0x9c9df69e,
+ 0xe86d2bc5, 0x22363a0d, 0xabc52197, 0x9b0deada, 0x1dbf9a42, 0xd5c4484e, 0x0abcd06b, 0xfa53ddef,
+ 0x3c1b20ee, 0x3fd59d7c, 0x25e41d2b, 0x669e1ef1, 0x6e6f52c3, 0x164df4fb, 0x7930e9e4, 0xe58857b6,
+ 0xac7d5f42, 0xd69f6d18, 0x7763cf1d, 0x55034004, 0x87f55ba5, 0x7e31cc7a, 0x7135c886, 0xefb4318a,
+ 0xed6a1e01, 0x2d9e6832, 0xa907600a, 0x918130c4, 0x6dc778f9, 0x71ad0038, 0x092999a3, 0x33cb8b7a,
+ 0x1a1db93d, 0x7140003c, 0x2a4ecea9, 0xf98d0acc, 0x0a8291cd, 0xcec97dcf, 0x8ec9b55a, 0x7f88a46b,
+ 0x4db5a851, 0xf44182e1, 0xc68a007e, 0x5e0dd902, 0x0bfd64b6, 0x45036c7a, 0x4e677d2c, 0x38532a3a,
+ 0x23ba4442, 0xcaf53ea6, 0x3bb45432, 0x9b7624c8, 0x917bdd64, 0xb1c0fd4c, 0xb38e8c33, 0x4c701c3a,
+ 0xcdad0657, 0xfccfec71, 0x9b1f5c3e, 0x4e46041f, 0x388147fb, 0x4cfdb477, 0xa52471f7, 0xa9a96910,
+ 0xb855322e, 0xdb6340d8, 0xa00ef092, 0x350511e3, 0x0abec1ff, 0xf9e3a26e, 0x7fb29f8c, 0x183023c3,
+ 0x587e38da, 0x0077d9b4, 0x763e4e4b, 0x94b2bbc1, 0x94c6651e, 0x77caf992, 0xeeaac023, 0x2a281bf6,
+ 0xb3a739c1, 0x22611682, 0x0ae8db58, 0x47a67cbe, 0xf9c9091b, 0x462d538c, 0xd72b0374, 0x6ae77f5e,
+ 0x62292c31, 0x1562a846, 0x505dc82d, 0xb854338a, 0xe49f5235, 0xc95b9117, 0x8ccf2dd5, 0xcacef403,
+ 0xec9d1810, 0xc6272b04, 0x5b3b71f9, 0xdc6b80d6, 0x3fdd4a8e, 0x9adb1e69, 0x62a69526, 0xd43161c1,
+ 0xa41d570d, 0x7938dad4, 0xa40e329c, 0xcff46aaa, 0x36ad004c, 0xf600c838, 0x1e425a31, 0xd951ae64,
+ 0xfdb23fce, 0xc9509d43, 0x687feb69, 0xedd1cc5e, 0x0b8cc3bd, 0xf64b10ef, 0x86b63142, 0xa3ab8829,
+ 0x555b2f74, 0x7c932665, 0xcb2c0f1c, 0xc01bd702, 0x29388839, 0xd2af05e4, 0x54504ac7, 0x8b758282,
+ 0x2846c0ba, 0x35c35f5c, 0x59160cc0, 0x46fd8251, 0x541fc68c, 0x9c86b022, 0xbb709987, 0x6a460e74,
+ 0x51a8a931, 0x09703fee, 0x1c217e6c, 0x3826e52c, 0x51aa691e, 0x0e423cfc, 0x99e9e316, 0x50c1217b,
+ 0x624816cd, 0xad9a95f9, 0xd5b80194, 0x88d9c0a0, 0xa1fe3075, 0xa577e231, 0x83f81d4a, 0x3f2fa457,
+ 0x1efc8ce0, 0xba8a4fe8, 0xb6855dfe, 0x72b0a66e, 0xded2fbab, 0xfbe58a30, 0xfafabe1c, 0x5d71a87e,
+ 0x2f741ef8, 0xc1fe86fe, 0xa6bbfde5, 0x30677f0d, 0x97d11d49, 0xf7a8443d, 0x0822e506, 0xa9f4614e,
+ 0x011e2a94, 0x838ff88c, 0xd68c8bb7, 0xc51eef6d, 0x49ea8ab4, 0xf2c3df5b, 0xb4e0735a, 0xb0d68749
+ // 0x2fe26dd4, ...
+ // 256 32bit words for the mantissa -- about 2464 valid decimal digits
+ };
+
+ // above value was calculated using Big<1,400> type on a 32bit platform
+ // and then the first 256 words were taken,
+ // the calculating was made by using ExpSurrounding0(1) method
+ // which took 1420 iterations
+ // (the result was compared with e taken from http://antwrp.gsfc.nasa.gov/htmltest/gifcity/e.2mil)
+ // (TTMATH_BUILTIN_VARIABLES_SIZE on 32bit platform should have the value 256,
+ // and on 64bit platform value 128 (256/2=128))
+
+ mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(int));
+ exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT) + 2;
+ info = 0;
+ }
+
+
+ /*!
+ this method sets the value of ln(2)
+ the natural logarithm from 2
+ */
+ void SetLn2()
+ {
+ static const unsigned int temp_table[] = {
+ 0xb17217f7, 0xd1cf79ab, 0xc9e3b398, 0x03f2f6af, 0x40f34326, 0x7298b62d, 0x8a0d175b, 0x8baafa2b,
+ 0xe7b87620, 0x6debac98, 0x559552fb, 0x4afa1b10, 0xed2eae35, 0xc1382144, 0x27573b29, 0x1169b825,
+ 0x3e96ca16, 0x224ae8c5, 0x1acbda11, 0x317c387e, 0xb9ea9bc3, 0xb136603b, 0x256fa0ec, 0x7657f74b,
+ 0x72ce87b1, 0x9d6548ca, 0xf5dfa6bd, 0x38303248, 0x655fa187, 0x2f20e3a2, 0xda2d97c5, 0x0f3fd5c6,
+ 0x07f4ca11, 0xfb5bfb90, 0x610d30f8, 0x8fe551a2, 0xee569d6d, 0xfc1efa15, 0x7d2e23de, 0x1400b396,
+ 0x17460775, 0xdb8990e5, 0xc943e732, 0xb479cd33, 0xcccc4e65, 0x9393514c, 0x4c1a1e0b, 0xd1d6095d,
+ 0x25669b33, 0x3564a337, 0x6a9c7f8a, 0x5e148e82, 0x074db601, 0x5cfe7aa3, 0x0c480a54, 0x17350d2c,
+ 0x955d5179, 0xb1e17b9d, 0xae313cdb, 0x6c606cb1, 0x078f735d, 0x1b2db31b, 0x5f50b518, 0x5064c18b,
+ 0x4d162db3, 0xb365853d, 0x7598a195, 0x1ae273ee, 0x5570b6c6, 0x8f969834, 0x96d4e6d3, 0x30af889b,
+ 0x44a02554, 0x731cdc8e, 0xa17293d1, 0x228a4ef9, 0x8d6f5177, 0xfbcf0755, 0x268a5c1f, 0x9538b982,
+ 0x61affd44, 0x6b1ca3cf, 0x5e9222b8, 0x8c66d3c5, 0x422183ed, 0xc9942109, 0x0bbb16fa, 0xf3d949f2,
+ 0x36e02b20, 0xcee886b9, 0x05c128d5, 0x3d0bd2f9, 0x62136319, 0x6af50302, 0x0060e499, 0x08391a0c,
+ 0x57339ba2, 0xbeba7d05, 0x2ac5b61c, 0xc4e9207c, 0xef2f0ce2, 0xd7373958, 0xd7622658, 0x901e646a,
+ 0x95184460, 0xdc4e7487, 0x156e0c29, 0x2413d5e3, 0x61c1696d, 0xd24aaebd, 0x473826fd, 0xa0c238b9,
+ 0x0ab111bb, 0xbd67c724, 0x972cd18b, 0xfbbd9d42, 0x6c472096, 0xe76115c0, 0x5f6f7ceb, 0xac9f45ae,
+ 0xcecb72f1, 0x9c38339d, 0x8f682625, 0x0dea891e, 0xf07afff3, 0xa892374e, 0x175eb4af, 0xc8daadd8,
+ 0x85db6ab0, 0x3a49bd0d, 0xc0b1b31d, 0x8a0e23fa, 0xc5e5767d, 0xf95884e0, 0x6425a415, 0x26fac51c,
+ 0x3ea8449f, 0xe8f70edd, 0x062b1a63, 0xa6c4c60c, 0x52ab3316, 0x1e238438, 0x897a39ce, 0x78b63c9f,
+ 0x364f5b8a, 0xef22ec2f, 0xee6e0850, 0xeca42d06, 0xfb0c75df, 0x5497e00c, 0x554b03d7, 0xd2874a00,
+ 0x0ca8f58d, 0x94f0341c, 0xbe2ec921, 0x56c9f949, 0xdb4a9316, 0xf281501e, 0x53daec3f, 0x64f1b783,
+ 0x154c6032, 0x0e2ff793, 0x33ce3573, 0xfacc5fdc, 0xf1178590, 0x3155bbd9, 0x0f023b22, 0x0224fcd8,
+ 0x471bf4f4, 0x45f0a88a, 0x14f0cd97, 0x6ea354bb, 0x20cdb5cc, 0xb3db2392, 0x88d58655, 0x4e2a0e8a,
+ 0x6fe51a8c, 0xfaa72ef2, 0xad8a43dc, 0x4212b210, 0xb779dfe4, 0x9d7307cc, 0x846532e4, 0xb9694eda,
+ 0xd162af05, 0x3b1751f3, 0xa3d091f6, 0x56658154, 0x12b5e8c2, 0x02461069, 0xac14b958, 0x784934b8,
+ 0xd6cce1da, 0xa5053701, 0x1aa4fb42, 0xb9a3def4, 0x1bda1f85, 0xef6fdbf2, 0xf2d89d2a, 0x4b183527,
+ 0x8fd94057, 0x89f45681, 0x2b552879, 0xa6168695, 0xc12963b0, 0xff01eaab, 0x73e5b5c1, 0x585318e7,
+ 0x624f14a5, 0x1a4a026b, 0x68082920, 0x57fd99b6, 0x6dc085a9, 0x8ac8d8ca, 0xf9eeeea9, 0x8a2400ca,
+ 0xc95f260f, 0xd10036f9, 0xf91096ac, 0x3195220a, 0x1a356b2a, 0x73b7eaad, 0xaf6d6058, 0x71ef7afb,
+ 0x80bc4234, 0x33562e94, 0xb12dfab4, 0x14451579, 0xdf59eae0, 0x51707062, 0x4012a829, 0x62c59cab,
+ 0x347f8304, 0xd889659e, 0x5a9139db, 0x14efcc30, 0x852be3e8, 0xfc99f14d, 0x1d822dd6, 0xe2f76797,
+ 0xe30219c8, 0xaa9ce884, 0x8a886eb3, 0xc87b7295, 0x988012e8, 0x314186ed, 0xbaf86856, 0xccd3c3b6,
+ 0xee94e62f, 0x110a6783, 0xd2aae89c, 0xcc3b76fc, 0x435a0ce1, 0x34c2838f, 0xd571ec6c, 0x1366a993 // last one was: 0x1366a992
+ //0xcbb9ac40, ...
+ // (the last word 0x1366a992 was rounded up because the next one is 0xcbb9ac40 -- first bit is one 0xc..)
+ // 256 32bit words for the mantissa -- about 2464 valid decimal digits
+ };
+
+ // above value was calculated using Big<1,400> type on a 32bit platform
+ // and then the first 256 words were taken,
+ // the calculating was made by using LnSurrounding1(2) method
+ // which took 4035 iterations
+ // (the result was compared with ln(2) taken from http://ja0hxv.calico.jp/pai/estart.html)
+ // (TTMATH_BUILTIN_VARIABLES_SIZE on 32bit platform should have the value 256,
+ // and on 64bit platform value 128 (256/2=128))
+
+ mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(int));
+ exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT);
+ info = 0;
+ }
+
+
+ /*!
+ this method sets the value of ln(10)
+ the natural logarithm from 10
+
+ I introduced this constant especially to make the conversion ToString()
+ being faster. In fact the method ToString() is keeping values of logarithms
+ it has calculated but it must calculate the logarithm at least once.
+ If a program, which uses this library, is running for a long time this
+ would be ok, but for programs which are running shorter, for example for
+ CGI applications which only once are printing values, this would be much
+ inconvenience. Then if we're printing with base (radix) 10 and the mantissa
+ of our value is smaller than or equal to TTMATH_BUILTIN_VARIABLES_SIZE
+ we don't calculate the logarithm but take it from this constant.
+ */
+ void SetLn10()
+ {
+ static const unsigned int temp_table[] = {
+ 0x935d8ddd, 0xaaa8ac16, 0xea56d62b, 0x82d30a28, 0xe28fecf9, 0xda5df90e, 0x83c61e82, 0x01f02d72,
+ 0x962f02d7, 0xb1a8105c, 0xcc70cbc0, 0x2c5f0d68, 0x2c622418, 0x410be2da, 0xfb8f7884, 0x02e516d6,
+ 0x782cf8a2, 0x8a8c911e, 0x765aa6c3, 0xb0d831fb, 0xef66ceb0, 0x4ab3c6fa, 0x5161bb49, 0xd219c7bb,
+ 0xca67b35b, 0x23605085, 0x8e93368d, 0x44789c4f, 0x5b08b057, 0xd5ede20f, 0x469ea58e, 0x9305e981,
+ 0xe2478fca, 0xad3aee98, 0x9cd5b42e, 0x6a271619, 0xa47ecb26, 0x978c5d4f, 0xdb1d28ea, 0x57d4fdc0,
+ 0xe40bf3cc, 0x1e14126a, 0x45765cde, 0x268339db, 0xf47fa96d, 0xeb271060, 0xaf88486e, 0xa9b7401e,
+ 0x3dfd3c51, 0x748e6d6e, 0x3848c8d2, 0x5faf1bca, 0xe88047f1, 0x7b0d9b50, 0xa949eaaa, 0xdf69e8a5,
+ 0xf77e3760, 0x4e943960, 0xe38a5700, 0xffde2db1, 0xad6bfbff, 0xd821ba0a, 0x4cb0466d, 0x61ba648e,
+ 0xef99c8e5, 0xf6974f36, 0x3982a78c, 0xa45ddfc8, 0x09426178, 0x19127a6e, 0x3b70fcda, 0x2d732d47,
+ 0xb5e4b1c8, 0xc0e5a10a, 0xaa6604a5, 0x324ec3dc, 0xbc64ea80, 0x6e198566, 0x1f1d366c, 0x20663834,
+ 0x4d5e843f, 0x20642b97, 0x0a62d18e, 0x478f7bd5, 0x8fcd0832, 0x4a7b32a6, 0xdef85a05, 0xeb56323a,
+ 0x421ef5e0, 0xb00410a0, 0xa0d9c260, 0x794a976f, 0xf6ff363d, 0xb00b6b33, 0xf42c58de, 0xf8a3c52d,
+ 0xed69b13d, 0xc1a03730, 0xb6524dc1, 0x8c167e86, 0x99d6d20e, 0xa2defd2b, 0xd006f8b4, 0xbe145a2a,
+ 0xdf3ccbb3, 0x189da49d, 0xbc1261c8, 0xb3e4daad, 0x6a36cecc, 0xb2d5ae5b, 0x89bf752f, 0xb5dfb353,
+ 0xff3065c4, 0x0cfceec8, 0x1be5a9a9, 0x67fddc57, 0xc4b83301, 0x006bf062, 0x4b40ed7a, 0x56c6cdcd,
+ 0xa2d6fe91, 0x388e9e3e, 0x48a93f5f, 0x5e3b6eb4, 0xb81c4a5b, 0x53d49ea6, 0x8e668aea, 0xba83c7f8,
+ 0xfb5f06c3, 0x58ac8f70, 0xfa9d8c59, 0x8c574502, 0xbaf54c96, 0xc84911f0, 0x0482d095, 0x1a0af022,
+ 0xabbab080, 0xec97efd3, 0x671e4e0e, 0x52f166b6, 0xcd5cd226, 0x0dc67795, 0x2e1e34a3, 0xf799677f,
+ 0x2c1d48f1, 0x2944b6c5, 0x2ba1307e, 0x704d67f9, 0x1c1035e4, 0x4e927c63, 0x03cf12bf, 0xe2cd2e31,
+ 0xf8ee4843, 0x344d51b0, 0xf37da42b, 0x9f0b0fd9, 0x134fb2d9, 0xf815e490, 0xd966283f, 0x23962766,
+ 0xeceab1e4, 0xf3b5fc86, 0x468127e2, 0xb606d10d, 0x3a45f4b6, 0xb776102d, 0x2fdbb420, 0x80c8fa84,
+ 0xd0ff9f45, 0xc58aef38, 0xdb2410fd, 0x1f1cebad, 0x733b2281, 0x52ca5f36, 0xddf29daa, 0x544334b8,
+ 0xdeeaf659, 0x4e462713, 0x1ed485b4, 0x6a0822e1, 0x28db471c, 0xa53938a8, 0x44c3bef7, 0xf35215c8,
+ 0xb382bc4e, 0x3e4c6f15, 0x6285f54c, 0x17ab408e, 0xccbf7f5e, 0xd16ab3f6, 0xced2846d, 0xf457e14f,
+ 0xbb45d9c5, 0x646ad497, 0xac697494, 0x145de32e, 0x93907128, 0xd263d521, 0x79efb424, 0xd64651d6,
+ 0xebc0c9f0, 0xbb583a44, 0xc6412c84, 0x85bb29a6, 0x4d31a2cd, 0x92954469, 0xa32b1abd, 0xf7f5202c,
+ 0xa4aa6c93, 0x2e9b53cf, 0x385ab136, 0x2741f356, 0x5de9c065, 0x6009901c, 0x88abbdd8, 0x74efcf73,
+ 0x3f761ad4, 0x35f3c083, 0xfd6b8ee0, 0x0bef11c7, 0xc552a89d, 0x58ce4a21, 0xd71e54f2, 0x4157f6c7,
+ 0xd4622316, 0xe98956d7, 0x450027de, 0xcbd398d8, 0x4b98b36a, 0x0724c25c, 0xdb237760, 0xe9324b68,
+ 0x7523e506, 0x8edad933, 0x92197f00, 0xb853a326, 0xb330c444, 0x65129296, 0x34bc0670, 0xe177806d,
+ 0xe338dac4, 0x5537492a, 0xe19add83, 0xcf45000f, 0x5b423bce, 0x6497d209, 0xe30e18a1, 0x3cbf0687,
+ 0x67973103, 0xd9485366, 0x81506bba, 0x2e93a9a4, 0x7dd59d3f, 0xf17cd746, 0x8c2075be, 0x552a4348 // last one was: 0x552a4347
+ // 0xb4a638ef, ...
+ //(the last word 0x552a4347 was rounded up because the next one is 0xb4a638ef -- first bit is one 0xb..)
+ // 256 32bit words for the mantissa -- about 2464 valid digits (decimal)
+ };
+
+ // above value was calculated using Big<1,400> type on a 32bit platform
+ // and then the first 256 32bit words were taken,
+ // the calculating was made by using LnSurrounding1(10) method
+ // which took 22080 iterations
+ // (the result was compared with ln(10) taken from http://ja0hxv.calico.jp/pai/estart.html)
+ // (the formula used in LnSurrounding1(x) converges badly when
+ // the x is greater than one but in fact we can use it, only the
+ // number of iterations will be greater)
+ // (TTMATH_BUILTIN_VARIABLES_SIZE on 32bit platform should have the value 256,
+ // and on 64bit platform value 128 (256/2=128))
+
+ mantissa.SetFromTable(temp_table, sizeof(temp_table) / sizeof(int));
+ exponent = -sint(man)*sint(TTMATH_BITS_PER_UINT) + 2;
+ info = 0;
+ }
+
+
+ /*!
+ this method sets the maximum value which can be held in this type
+ */
+ void SetMax()
+ {
+ info = 0;
+ mantissa.SetMax();
+ exponent.SetMax();
+
+ // we don't have to use 'Standardizing()' because the last bit from
+ // the mantissa is set
+ }
+
+
+ /*!
+ this method sets the minimum value which can be held in this type
+ */
+ void SetMin()
+ {
+ info = 0;
+
+ mantissa.SetMax();
+ exponent.SetMax();
+ SetSign();
+
+ // we don't have to use 'Standardizing()' because the last bit from
+ // the mantissa is set
+ }
+
+
+ /*!
+ testing whether there is a value zero or not
+ */
+ bool IsZero() const
+ {
+ return IsInfoBit(TTMATH_BIG_ZERO);
+ }
+
+
+ /*!
+ this method returns true when there's the sign set
+ also we don't check the NaN flag
+ */
+ bool IsSign() const
+ {
+ return IsInfoBit(TTMATH_BIG_SIGN);
+ }
+
+
+ /*!
+ this method returns true when there is not a valid number
+ */
+ bool IsNan() const
+ {
+ return IsInfoBit(TTMATH_BIG_NAN);
+ }
+
+
+
+ /*!
+ this method clears the sign
+ (there'll be an absolute value)
+
+ e.g.
+ -1 -> 1
+ 2 -> 2
+ */
+ void Abs()
+ {
+ ClearInfoBit(TTMATH_BIG_SIGN);
+ }
+
+
+ /*!
+ this method remains the 'sign' of the value
+ e.g. -2 = -1
+ 0 = 0
+ 10 = 1
+ */
+ void Sgn()
+ {
+ // we have to check the NaN flag, because the next SetOne() method would clear it
+ if( IsNan() )
+ return;
+
+ if( IsSign() )
+ {
+ SetOne();
+ SetSign();
+ }
+ else
+ if( IsZero() )
+ SetZero();
+ else
+ SetOne();
+ }
+
+
+
+ /*!
+ this method sets the sign
+
+ e.g.
+ -1 -> -1
+ 2 -> -2
+
+ we do not check whether there is a zero or not, if you're using this method
+ you must be sure that the value is (or will be afterwards) different from zero
+ */
+ void SetSign()
+ {
+ SetInfoBit(TTMATH_BIG_SIGN);
+ }
+
+
+ /*!
+ this method changes the sign
+ when there is a value of zero then the sign is not changed
+
+ e.g.
+ -1 -> 1
+ 2 -> -2
+ */
+ void ChangeSign()
+ {
+ // we don't have to check the NaN flag here
+
+ if( IsZero() )
+ return;
+
+ if( IsSign() )
+ ClearInfoBit(TTMATH_BIG_SIGN);
+ else
+ SetInfoBit(TTMATH_BIG_SIGN);
+ }
+
+
+
+private:
+
+ /*!
+ this method does the half-to-even rounding (banker's rounding)
+
+ if is_half is:
+ true - that means the rest was equal the half (0.5 decimal)
+ false - that means the rest was greater than a half (greater than 0.5 decimal)
+
+ if the rest was less than a half then don't call this method
+ (the rounding should does nothing then)
+ */
+ uint RoundHalfToEven(bool is_half, bool rounding_up = true)
+ {
+ uint c = 0;
+
+ if( !is_half || mantissa.IsTheLowestBitSet() )
+ {
+ if( rounding_up )
+ {
+ if( mantissa.AddOne() )
+ {
+ mantissa.Rcr(1, 1);
+ c = exponent.AddOne();
+ }
+ }
+ else
+ {
+ #ifdef TTMATH_DEBUG
+ uint c_from_zero =
+ #endif
+ mantissa.SubOne();
+
+ // we're using rounding_up=false in Add() when the mantissas have different signs
+ // mantissa can be zero only when previous mantissa was equal to ss2.mantissa
+ // but in such a case 'last_bit_set' will not be set and consequently 'do_rounding' will be false
+ TTMATH_ASSERT( c_from_zero == 0 )
+ }
+ }
+
+ return c;
+ }
+
+
+
+
+
+ /*!
+ *
+ * basic mathematic functions
+ *
+ */
+
+private:
+
+
+ /*!
+ an auxiliary method for adding
+ */
+ void AddCheckExponents( Big & ss2,
+ Int & exp_offset,
+ bool & last_bit_set,
+ bool & rest_zero,
+ bool & do_adding,
+ bool & do_rounding)
+ {
+ Int mantissa_size_in_bits( man * TTMATH_BITS_PER_UINT );
+
+ if( exp_offset == mantissa_size_in_bits )
+ {
+ last_bit_set = ss2.mantissa.IsTheHighestBitSet();
+ rest_zero = ss2.mantissa.AreFirstBitsZero(man*TTMATH_BITS_PER_UINT - 1);
+ do_rounding = true; // we'are only rounding
+ }
+ else
+ if( exp_offset < mantissa_size_in_bits )
+ {
+ uint moved = exp_offset.ToInt(); // how many times we must move ss2.mantissa
+ rest_zero = true;
+
+ if( moved > 0 )
+ {
+ last_bit_set = static_cast( ss2.mantissa.GetBit(moved-1) );
+
+ if( moved > 1 )
+ rest_zero = ss2.mantissa.AreFirstBitsZero(moved - 1);
+
+ // (2) moving 'exp_offset' times
+ ss2.mantissa.Rcr(moved, 0);
+ }
+
+ do_adding = true;
+ do_rounding = true;
+ }
+
+ // if exp_offset is greater than mantissa_size_in_bits then we do nothing
+ // ss2 is too small for taking into consideration in the sum
+ }
+
+
+ /*!
+ an auxiliary method for adding
+ */
+ uint AddMantissas( Big & ss2,
+ bool & last_bit_set,
+ bool & rest_zero,
+ bool & rounding_up)
+ {
+ uint c = 0;
+
+ if( IsSign() == ss2.IsSign() )
+ {
+ // values have the same signs
+ if( mantissa.Add(ss2.mantissa) )
+ {
+ // we have one bit more from addition (carry)
+ // now rest_zero means the old rest_zero with the old last_bit_set
+ rest_zero = (!last_bit_set && rest_zero);
+ last_bit_set = mantissa.Rcr(1,1);
+ c += exponent.AddOne();
+ }
+
+ rounding_up = true;
+ }
+ else
+ {
+ // values have different signs
+ // there shouldn't be a carry here because
+ // (1) (2) guarantee that the mantissa of this
+ // is greater than or equal to the mantissa of the ss2
+
+ #ifdef TTMATH_DEBUG
+ uint c_temp =
+ #endif
+ mantissa.Sub(ss2.mantissa);
+
+ TTMATH_ASSERT( c_temp == 0 )
+ }
+
+ return c;
+ }
+
+
+public:
+
+
+ /*!
+ Addition this = this + ss2
+
+ it returns carry if the sum is too big
+ */
+ uint Add(Big ss2, bool round = true)
+ {
+ bool last_bit_set, rest_zero, do_adding, do_rounding, rounding_up;
+ Int exp_offset( exponent );
+ uint c = 0;
+
+ if( IsNan() || ss2.IsNan() )
+ return CheckCarry(1);
+
+ exp_offset.Sub( ss2.exponent );
+ exp_offset.Abs();
+
+ // (1) abs(this) will be >= abs(ss2)
+ if( SmallerWithoutSignThan(ss2) )
+ {
+ // !! use Swap here (not implemented yet)
+ Big temp(ss2);
+
+ ss2 = *this;
+ *this = temp;
+ }
+
+ if( ss2.IsZero() )
+ return 0;
+
+ last_bit_set = rest_zero = do_adding = do_rounding = rounding_up = false;
+ AddCheckExponents(ss2, exp_offset, last_bit_set, rest_zero, do_adding, do_rounding);
+
+ if( do_adding )
+ c += AddMantissas(ss2, last_bit_set, rest_zero, rounding_up);
+
+ if( !round || !last_bit_set )
+ do_rounding = false;
+
+ if( do_rounding )
+ c += RoundHalfToEven(rest_zero, rounding_up);
+
+ if( do_adding || do_rounding )
+ c += Standardizing();
+
+ return CheckCarry(c);
+ }
+
+
+
+ /*!
+ Subtraction this = this - ss2
+
+ it returns carry if the result is too big
+ */
+ uint Sub(Big ss2, bool round = true)
+ {
+ ss2.ChangeSign();
+
+ return Add(ss2, round);
+ }
+
+
+ /*!
+ bitwise AND
+
+ this and ss2 must be >= 0
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - this or ss2 was negative
+ */
+ uint BitAnd(Big ss2)
+ {
+ if( IsNan() || ss2.IsNan() )
+ return CheckCarry(1);
+
+ if( IsSign() || ss2.IsSign() )
+ return 2;
+
+ if( IsZero() )
+ return 0;
+
+ if( ss2.IsZero() )
+ {
+ SetZero();
+ return 0;
+ }
+
+ Int exp_offset( exponent );
+ Int mantissa_size_in_bits( man * TTMATH_BITS_PER_UINT );
+
+ uint c = 0;
+
+ exp_offset.Sub( ss2.exponent );
+ exp_offset.Abs();
+
+ // abs(this) will be >= abs(ss2)
+ if( SmallerWithoutSignThan(ss2) )
+ {
+ Big temp(ss2);
+
+ ss2 = *this;
+ *this = temp;
+ }
+
+ if( exp_offset >= mantissa_size_in_bits )
+ {
+ // the second value is too small
+ SetZero();
+ return 0;
+ }
+
+ // exp_offset < mantissa_size_in_bits, moving 'exp_offset' times
+ ss2.mantissa.Rcr( exp_offset.ToInt(), 0 );
+ mantissa.BitAnd(ss2.mantissa);
+
+ c += Standardizing();
+
+ return CheckCarry(c);
+ }
+
+
+ /*!
+ bitwise OR
+
+ this and ss2 must be >= 0
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - this or ss2 was negative
+ */
+ uint BitOr(Big ss2)
+ {
+ if( IsNan() || ss2.IsNan() )
+ return CheckCarry(1);
+
+ if( IsSign() || ss2.IsSign() )
+ return 2;
+
+ if( IsZero() )
+ {
+ *this = ss2;
+ return 0;
+ }
+
+ if( ss2.IsZero() )
+ return 0;
+
+ Int exp_offset( exponent );
+ Int mantissa_size_in_bits( man * TTMATH_BITS_PER_UINT );
+
+ uint c = 0;
+
+ exp_offset.Sub( ss2.exponent );
+ exp_offset.Abs();
+
+ // abs(this) will be >= abs(ss2)
+ if( SmallerWithoutSignThan(ss2) )
+ {
+ Big temp(ss2);
+
+ ss2 = *this;
+ *this = temp;
+ }
+
+ if( exp_offset >= mantissa_size_in_bits )
+ // the second value is too small
+ return 0;
+
+ // exp_offset < mantissa_size_in_bits, moving 'exp_offset' times
+ ss2.mantissa.Rcr( exp_offset.ToInt(), 0 );
+ mantissa.BitOr(ss2.mantissa);
+
+ c += Standardizing();
+
+ return CheckCarry(c);
+ }
+
+
+ /*!
+ bitwise XOR
+
+ this and ss2 must be >= 0
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - this or ss2 was negative
+ */
+ uint BitXor(Big ss2)
+ {
+ if( IsNan() || ss2.IsNan() )
+ return CheckCarry(1);
+
+ if( IsSign() || ss2.IsSign() )
+ return 2;
+
+ if( ss2.IsZero() )
+ return 0;
+
+ if( IsZero() )
+ {
+ *this = ss2;
+ return 0;
+ }
+
+ Int exp_offset( exponent );
+ Int mantissa_size_in_bits( man * TTMATH_BITS_PER_UINT );
+
+ uint c = 0;
+
+ exp_offset.Sub( ss2.exponent );
+ exp_offset.Abs();
+
+ // abs(this) will be >= abs(ss2)
+ if( SmallerWithoutSignThan(ss2) )
+ {
+ Big temp(ss2);
+
+ ss2 = *this;
+ *this = temp;
+ }
+
+ if( exp_offset >= mantissa_size_in_bits )
+ // the second value is too small
+ return 0;
+
+ // exp_offset < mantissa_size_in_bits, moving 'exp_offset' times
+ ss2.mantissa.Rcr( exp_offset.ToInt(), 0 );
+ mantissa.BitXor(ss2.mantissa);
+
+ c += Standardizing();
+
+ return CheckCarry(c);
+ }
+
+
+
+ /*!
+ Multiplication this = this * ss2 (ss2 is uint)
+
+ ss2 without a sign
+ */
+ uint MulUInt(uint ss2)
+ {
+ UInt man_result;
+ uint i,c = 0;
+
+ if( IsNan() )
+ return 1;
+
+ if( IsZero() )
+ return 0;
+
+ if( ss2 == 0 )
+ {
+ SetZero();
+ return 0;
+ }
+
+ // man_result = mantissa * ss2.mantissa
+ mantissa.MulInt(ss2, man_result);
+
+ sint bit = UInt::FindLeadingBitInWord(man_result.table[man]); // man - last word
+
+ if( bit!=-1 && uint(bit) > (TTMATH_BITS_PER_UINT/2) )
+ {
+ // 'i' will be from 0 to TTMATH_BITS_PER_UINT
+ i = man_result.CompensationToLeft();
+ c = exponent.Add( TTMATH_BITS_PER_UINT - i );
+
+ for(i=0 ; i0 && (tab[len-1] & TTMATH_UINT_HIGHEST_BIT)!=0 )
+
+ for(i=0 ; i & ss2, bool round = true)
+ {
+ TTMATH_REFERENCE_ASSERT( ss2 )
+
+ UInt man_result;
+ uint c = 0;
+ uint i;
+
+ if( IsNan() || ss2.IsNan() )
+ return CheckCarry(1);
+
+ if( IsZero() )
+ return 0;
+
+ if( ss2.IsZero() )
+ {
+ SetZero();
+ return 0;
+ }
+
+ // man_result = mantissa * ss2.mantissa
+ mantissa.MulBig(ss2.mantissa, man_result);
+
+ // 'i' will be from 0 to man*TTMATH_BITS_PER_UINT
+ // because mantissa and ss2.mantissa are standardized
+ // (the highest bit in man_result is set to 1 or
+ // if there is a zero value in man_result the method CompensationToLeft()
+ // returns 0 but we'll correct this at the end in Standardizing() method)
+ i = man_result.CompensationToLeft();
+ uint exp_add = man * TTMATH_BITS_PER_UINT - i;
+
+ if( exp_add )
+ c += exponent.Add( exp_add );
+
+ c += exponent.Add( ss2.exponent );
+
+ for(i=0 ; i & ss2, bool round = true)
+ {
+ TTMATH_REFERENCE_ASSERT( ss2 )
+
+ UInt man1;
+ UInt man2;
+ uint i,c = 0;
+
+ if( IsNan() || ss2.IsNan() )
+ return CheckCarry(1);
+
+ if( ss2.IsZero() )
+ {
+ SetNan();
+ return 2;
+ }
+
+ if( IsZero() )
+ return 0;
+
+ for(i=0 ; i & ss2)
+ {
+ TTMATH_REFERENCE_ASSERT( ss2 )
+
+ uint c = 0;
+
+ if( IsNan() || ss2.IsNan() )
+ return CheckCarry(1);
+
+ if( ss2.IsZero() )
+ {
+ SetNan();
+ return 2;
+ }
+
+ if( !SmallerWithoutSignThan(ss2) )
+ {
+ Big temp(*this);
+
+ c = temp.Div(ss2);
+ temp.SkipFraction();
+ c += temp.Mul(ss2);
+ c += Sub(temp);
+
+ if( !SmallerWithoutSignThan( ss2 ) )
+ c += 1;
+ }
+
+ return CheckCarry(c);
+ }
+
+
+
+
+ /*!
+ power this = this ^ pow
+ (pow without a sign)
+
+ binary algorithm (r-to-l)
+
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - incorrect arguments (0^0)
+ */
+ template
+ uint Pow(UInt pow)
+ {
+ if( IsNan() )
+ return 1;
+
+ if( IsZero() )
+ {
+ if( pow.IsZero() )
+ {
+ // we don't define zero^zero
+ SetNan();
+ return 2;
+ }
+
+ // 0^(+something) is zero
+ return 0;
+ }
+
+ Big start(*this), start_temp;
+ Big result;
+ result.SetOne();
+ uint c = 0;
+
+ while( !c )
+ {
+ if( pow.table[0] & 1 )
+ c += result.Mul(start);
+
+ pow.Rcr(1);
+
+ if( pow.IsZero() )
+ break;
+
+ start_temp = start;
+ c += start.Mul(start_temp);
+ }
+
+ *this = result;
+
+ return CheckCarry(c);
+ }
+
+
+ /*!
+ power this = this ^ pow
+ p can be negative
+
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - incorrect arguments 0^0 or 0^(-something)
+ */
+ template
+ uint Pow(Int pow)
+ {
+ if( IsNan() )
+ return 1;
+
+ if( !pow.IsSign() )
+ return Pow( UInt(pow) );
+
+ if( IsZero() )
+ {
+ // if 'p' is negative then
+ // 'this' must be different from zero
+ SetNan();
+ return 2;
+ }
+
+ uint c = pow.ChangeSign();
+
+ Big t(*this);
+ c += t.Pow( UInt(pow) ); // here can only be a carry (return:1)
+
+ SetOne();
+ c += Div(t);
+
+ return CheckCarry(c);
+ }
+
+
+ /*!
+ this method returns: 'this' mod 2
+ (either zero or one)
+
+ this method is much faster than using Mod( object_with_value_two )
+ */
+ uint Mod2() const
+ {
+ if( exponent>sint(0) || exponent<=-sint(man*TTMATH_BITS_PER_UINT) )
+ return 0;
+
+ sint exp_int = exponent.ToInt();
+ // 'exp_int' is negative (or zero), we set it as positive
+ exp_int = -exp_int;
+
+ return mantissa.GetBit(exp_int);
+ }
+
+
+
+ /*!
+ power this = this ^ abs([pow])
+ pow is treated as a value without a sign and without a fraction
+ if pow has a sign then the method pow.Abs() is used
+ if pow has a fraction the fraction is skipped (not used in calculation)
+
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - incorrect arguments (0^0)
+ */
+ uint PowUInt(Big pow)
+ {
+ if( IsNan() || pow.IsNan() )
+ return CheckCarry(1);
+
+ if( IsZero() )
+ {
+ if( pow.IsZero() )
+ {
+ SetNan();
+ return 2;
+ }
+
+ // 0^(+something) is zero
+ return 0;
+ }
+
+ if( pow.IsSign() )
+ pow.Abs();
+
+ Big start(*this), start_temp;
+ Big result;
+ Big one;
+ Int e_one;
+ uint c = 0;
+
+ e_one.SetOne();
+ one.SetOne();
+ result = one;
+
+ while( !c )
+ {
+ if( pow.Mod2() )
+ c += result.Mul(start);
+
+ c += pow.exponent.Sub( e_one ); // !! may use SubOne() here?
+
+ if( pow < one )
+ break;
+
+ start_temp = start;
+ c += start.Mul(start_temp);
+ }
+
+ *this = result;
+
+ return CheckCarry(c);
+ }
+
+
+ /*!
+ power this = this ^ [pow]
+ pow is treated as a value without a fraction
+ pow can be negative
+
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - incorrect arguments 0^0 or 0^(-something)
+ */
+ uint PowInt(const Big & pow)
+ {
+ TTMATH_REFERENCE_ASSERT( pow )
+
+ if( IsNan() || pow.IsNan() )
+ return CheckCarry(1);
+
+ if( !pow.IsSign() )
+ return PowUInt(pow);
+
+ if( IsZero() )
+ {
+ // if 'pow' is negative then
+ // 'this' must be different from zero
+ SetNan();
+ return 2;
+ }
+
+ Big temp(*this);
+ uint c = temp.PowUInt(pow); // here can only be a carry (result:1)
+
+ SetOne();
+ c += Div(temp);
+
+ return CheckCarry(c);
+ }
+
+
+ /*!
+ power this = this ^ pow
+ this must be greater than zero (this > 0)
+ pow can be negative and with fraction
+
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - incorrect argument ('this' <= 0)
+ */
+ uint PowFrac(const Big & pow)
+ {
+ TTMATH_REFERENCE_ASSERT( pow )
+
+ if( IsNan() || pow.IsNan() )
+ return CheckCarry(1);
+
+ Big temp;
+ uint c = temp.Ln(*this);
+
+ if( c != 0 ) // can be 2 from Ln()
+ {
+ SetNan();
+ return c;
+ }
+
+ c += temp.Mul(pow);
+ c += Exp(temp);
+
+ return CheckCarry(c);
+ }
+
+
+
+ /*!
+ power this = this ^ pow
+ pow can be negative and with fraction
+
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - incorrect argument ('this' or 'pow')
+ */
+ uint Pow(const Big & pow)
+ {
+ TTMATH_REFERENCE_ASSERT( pow )
+
+ if( IsNan() || pow.IsNan() )
+ return CheckCarry(1);
+
+ if( IsZero() )
+ {
+ // 0^pow will be 0 only for pow>0
+ if( pow.IsSign() || pow.IsZero() )
+ {
+ SetNan();
+ return 2;
+ }
+
+ SetZero();
+
+ return 0;
+ }
+
+ if( pow.exponent>-int(man*TTMATH_BITS_PER_UINT) && pow.exponent<=0 )
+ {
+ if( pow.IsInteger() )
+ return PowInt( pow );
+ }
+
+ return PowFrac(pow);
+ }
+
+
+ /*!
+ this function calculates the square root
+ e.g. let this=9 then this.Sqrt() gives 3
+
+ return: 0 - ok
+ 1 - carry
+ 2 - improper argument (this<0 or NaN)
+ */
+ uint Sqrt()
+ {
+ if( IsNan() || IsSign() )
+ {
+ SetNan();
+ return 2;
+ }
+
+ if( IsZero() )
+ return 0;
+
+ Big old(*this);
+ Big ln;
+ uint c = 0;
+
+ // we're using the formula: sqrt(x) = e ^ (ln(x) / 2)
+ c += ln.Ln(*this);
+ c += ln.exponent.SubOne(); // ln = ln / 2
+ c += Exp(ln);
+
+ // above formula doesn't give accurate results for some integers
+ // e.g. Sqrt(81) would not be 9 but a value very closed to 9
+ // we're rounding the result, calculating result*result and comparing
+ // with the old value, if they are equal then the result is an integer too
+
+ if( !c && old.IsInteger() && !IsInteger() )
+ {
+ Big temp(*this);
+ c += temp.Round();
+
+ Big temp2(temp);
+ c += temp.Mul(temp2);
+
+ if( temp == old )
+ *this = temp2;
+ }
+
+ return CheckCarry(c);
+ }
+
+
+private:
+
+#ifdef TTMATH_CONSTANTSGENERATOR
+public:
+#endif
+
+ /*!
+ Exponent this = exp(x) = e^x where x is in (-1,1)
+
+ we're using the formula exp(x) = 1 + (x)/(1!) + (x^2)/(2!) + (x^3)/(3!) + ...
+ */
+ void ExpSurrounding0(const Big & x, uint * steps = 0)
+ {
+ TTMATH_REFERENCE_ASSERT( x )
+
+ Big denominator, denominator_i;
+ Big one, old_value, next_part;
+ Big numerator = x;
+
+ SetOne();
+ one.SetOne();
+ denominator.SetOne();
+ denominator_i.SetOne();
+
+ uint i;
+ old_value = *this;
+
+ // we begin from 1 in order to not test at the beginning
+ #ifdef TTMATH_CONSTANTSGENERATOR
+ for(i=1 ; true ; ++i)
+ #else
+ for(i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
+ #endif
+ {
+ bool testing = ((i & 3) == 0); // it means '(i % 4) == 0'
+
+ next_part = numerator;
+
+ if( next_part.Div( denominator ) )
+ // if there is a carry here we only break the loop
+ // however the result we return as good
+ // it means there are too many parts of the formula
+ break;
+
+ // there shouldn't be a carry here
+ Add( next_part );
+
+ if( testing )
+ {
+ if( old_value == *this )
+ // we've added next few parts of the formula but the result
+ // is still the same then we break the loop
+ break;
+ else
+ old_value = *this;
+ }
+
+ // we set the denominator and the numerator for a next part of the formula
+ if( denominator_i.Add(one) )
+ // if there is a carry here the result we return as good
+ break;
+
+ if( denominator.Mul(denominator_i) )
+ break;
+
+ if( numerator.Mul(x) )
+ break;
+ }
+
+ if( steps )
+ *steps = i;
+ }
+
+public:
+
+
+ /*!
+ Exponent this = exp(x) = e^x
+
+ we're using the fact that our value is stored in form of:
+ x = mantissa * 2^exponent
+ then
+ e^x = e^(mantissa* 2^exponent) or
+ e^x = (e^mantissa)^(2^exponent)
+
+ 'Exp' returns a carry if we can't count the result ('x' is too big)
+ */
+ uint Exp(const Big & x)
+ {
+ uint c = 0;
+
+ if( x.IsNan() )
+ return CheckCarry(1);
+
+ if( x.IsZero() )
+ {
+ SetOne();
+ return 0;
+ }
+
+ // m will be the value of the mantissa in range (-1,1)
+ Big m(x);
+ m.exponent = -sint(man*TTMATH_BITS_PER_UINT);
+
+ // 'e_' will be the value of '2^exponent'
+ // e_.mantissa.table[man-1] = TTMATH_UINT_HIGHEST_BIT; and
+ // e_.exponent.Add(1) mean:
+ // e_.mantissa.table[0] = 1;
+ // e_.Standardizing();
+ // e_.exponent.Add(man*TTMATH_BITS_PER_UINT)
+ // (we must add 'man*TTMATH_BITS_PER_UINT' because we've taken it from the mantissa)
+ Big e_(x);
+ e_.mantissa.SetZero();
+ e_.mantissa.table[man-1] = TTMATH_UINT_HIGHEST_BIT;
+ c += e_.exponent.Add(1);
+ e_.Abs();
+
+ /*
+ now we've got:
+ m - the value of the mantissa in range (-1,1)
+ e_ - 2^exponent
+
+ e_ can be as:
+ ...2^-2, 2^-1, 2^0, 2^1, 2^2 ...
+ ...1/4 , 1/2 , 1 , 2 , 4 ...
+
+ above one e_ is integer
+
+ if e_ is greater than 1 we calculate the exponent as:
+ e^(m * e_) = ExpSurrounding0(m) ^ e_
+ and if e_ is smaller or equal one we calculate the exponent in this way:
+ e^(m * e_) = ExpSurrounding0(m* e_)
+ because if e_ is smaller or equal 1 then the product of m*e_ is smaller or equal m
+ */
+
+ if( e_ <= 1 )
+ {
+ m.Mul(e_);
+ ExpSurrounding0(m);
+ }
+ else
+ {
+ ExpSurrounding0(m);
+ c += PowUInt(e_);
+ }
+
+ return CheckCarry(c);
+ }
+
+
+
+
+private:
+
+#ifdef TTMATH_CONSTANTSGENERATOR
+public:
+#endif
+
+ /*!
+ Natural logarithm this = ln(x) where x in range <1,2)
+
+ we're using the formula:
+ ln x = 2 * [ (x-1)/(x+1) + (1/3)((x-1)/(x+1))^3 + (1/5)((x-1)/(x+1))^5 + ... ]
+ */
+ void LnSurrounding1(const Big & x, uint * steps = 0)
+ {
+ Big old_value, next_part, denominator, one, two, x1(x), x2(x);
+
+ one.SetOne();
+
+ if( x == one )
+ {
+ // LnSurrounding1(1) is 0
+ SetZero();
+ return;
+ }
+
+ two = 2;
+
+ x1.Sub(one);
+ x2.Add(one);
+
+ x1.Div(x2);
+ x2 = x1;
+ x2.Mul(x1);
+
+ denominator.SetOne();
+ SetZero();
+
+ old_value = *this;
+ uint i;
+
+
+ #ifdef TTMATH_CONSTANTSGENERATOR
+ for(i=1 ; true ; ++i)
+ #else
+ // we begin from 1 in order to not test at the beginning
+ for(i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
+ #endif
+ {
+ bool testing = ((i & 3) == 0); // it means '(i % 4) == 0'
+
+ next_part = x1;
+
+ if( next_part.Div(denominator) )
+ // if there is a carry here we only break the loop
+ // however the result we return as good
+ // it means there are too many parts of the formula
+ break;
+
+ // there shouldn't be a carry here
+ Add(next_part);
+
+ if( testing )
+ {
+ if( old_value == *this )
+ // we've added next (step_test) parts of the formula but the result
+ // is still the same then we break the loop
+ break;
+ else
+ old_value = *this;
+ }
+
+ if( x1.Mul(x2) )
+ // if there is a carry here the result we return as good
+ break;
+
+ if( denominator.Add(two) )
+ break;
+ }
+
+ // this = this * 2
+ // ( there can't be a carry here because we calculate the logarithm between <1,2) )
+ exponent.AddOne();
+
+ if( steps )
+ *steps = i;
+ }
+
+
+
+
+public:
+
+
+ /*!
+ Natural logarithm this = ln(x)
+ (a logarithm with the base equal 'e')
+
+ we're using the fact that our value is stored in form of:
+ x = mantissa * 2^exponent
+ then
+ ln(x) = ln (mantissa * 2^exponent) = ln (mantissa) + (exponent * ln (2))
+
+ the mantissa we'll show as a value from range <1,2) because the logarithm
+ is decreasing too fast when 'x' is going to 0
+
+ return values:
+ 0 - ok
+ 1 - overflow (carry)
+ 2 - incorrect argument (x<=0)
+ */
+ uint Ln(const Big & x)
+ {
+ TTMATH_REFERENCE_ASSERT( x )
+
+ if( x.IsNan() )
+ return CheckCarry(1);
+
+ if( x.IsSign() || x.IsZero() )
+ {
+ SetNan();
+ return 2;
+ }
+
+ // m will be the value of the mantissa in range <1,2)
+ Big m(x);
+ m.exponent = -sint(man*TTMATH_BITS_PER_UINT - 1);
+
+ LnSurrounding1(m);
+
+ Big exponent_temp;
+ exponent_temp.FromInt( x.exponent );
+
+ // we must add 'man*TTMATH_BITS_PER_UINT-1' because we've taken it from the mantissa
+ uint c = exponent_temp.Add(man*TTMATH_BITS_PER_UINT-1);
+
+ Big ln2;
+ ln2.SetLn2();
+ c += exponent_temp.Mul(ln2);
+ c += Add(exponent_temp);
+
+ return CheckCarry(c);
+ }
+
+
+
+ /*!
+ Logarithm from 'x' with a 'base'
+
+ we're using the formula:
+ Log(x) with 'base' = ln(x) / ln(base)
+
+ return values:
+ 0 - ok
+ 1 - overflow
+ 2 - incorrect argument (x<=0)
+ 3 - incorrect base (a<=0 lub a=1)
+ */
+ uint Log(const Big & x, const Big & base)
+ {
+ TTMATH_REFERENCE_ASSERT( base )
+ TTMATH_REFERENCE_ASSERT( x )
+
+ if( x.IsNan() || base.IsNan() )
+ return CheckCarry(1);
+
+ if( x.IsSign() || x.IsZero() )
+ {
+ SetNan();
+ return 2;
+ }
+
+ Big denominator;;
+ denominator.SetOne();
+
+ if( base.IsSign() || base.IsZero() || base==denominator )
+ {
+ SetNan();
+ return 3;
+ }
+
+ if( x == denominator ) // (this is: if x == 1)
+ {
+ // log(1) is 0
+ SetZero();
+ return 0;
+ }
+
+ // another error values we've tested at the beginning
+ // there can only be a carry
+ uint c = Ln(x);
+
+ c += denominator.Ln(base);
+ c += Div(denominator);
+
+ return CheckCarry(c);
+ }
+
+
+
+
+ /*!
+ *
+ * converting methods
+ *
+ */
+
+
+ /*!
+ converting from another type of a Big object
+ */
+ template
+ uint FromBig(const Big & another)
+ {
+ info = another.info;
+
+ if( IsNan() )
+ return 1;
+
+ if( exponent.FromInt(another.exponent) )
+ {
+ SetNan();
+ return 1;
+ }
+
+ uint man_len_min = (man < another_man)? man : another_man;
+ uint i;
+ uint c = 0;
+
+ for( i = 0 ; i another_man )' and 'if( man < another_man )' and there'll be no such situation here
+ #ifdef _MSC_VER
+ #pragma warning( disable: 4307 )
+ #endif
+
+ if( man > another_man )
+ {
+ uint man_diff = (man - another_man) * TTMATH_BITS_PER_UINT;
+ c += exponent.SubInt(man_diff, 0);
+ }
+ else
+ if( man < another_man )
+ {
+ uint man_diff = (another_man - man) * TTMATH_BITS_PER_UINT;
+ c += exponent.AddInt(man_diff, 0);
+ }
+
+ #ifdef _MSC_VER
+ #pragma warning( default: 4307 )
+ #endif
+
+ // mantissa doesn't have to be standardized (either the highest bit is set or all bits are equal zero)
+ CorrectZero();
+
+ return CheckCarry(c);
+ }
+
+
+ /*!
+ this method converts 'this' into 'result'
+
+ if the value is too big this method returns a carry (1)
+ */
+ uint ToUInt(uint & result, bool test_sign = true) const
+ {
+ result = 0;
+
+ if( IsZero() )
+ return 0;
+
+ if( test_sign && IsSign() )
+ // the result should be positive
+ return 1;
+
+ sint maxbit = -sint(man*TTMATH_BITS_PER_UINT);
+
+ if( exponent > maxbit + sint(TTMATH_BITS_PER_UINT) )
+ // if exponent > (maxbit + sint(TTMATH_BITS_PER_UINT)) the value can't be passed
+ // into the 'sint' type (it's too big)
+ return 1;
+
+ if( exponent <= maxbit )
+ // our value is from the range of (-1,1) and we return zero
+ return 0;
+
+ UInt mantissa_temp(mantissa);
+ // exponent is from a range of (maxbit, maxbit + sint(TTMATH_BITS_PER_UINT) >
+ sint how_many_bits = exponent.ToInt();
+
+ // how_many_bits is negative, we'll make it positive
+ how_many_bits = -how_many_bits;
+
+ // we're taking into account only the last word in a mantissa table
+ mantissa_temp.Rcr( how_many_bits % TTMATH_BITS_PER_UINT, 0 );
+ result = mantissa_temp.table[ man-1 ];
+
+ return 0;
+ }
+
+
+
+ /*!
+ this method converts 'this' into 'result'
+
+ if the value is too big this method returns a carry (1)
+ */
+ uint ToInt(sint & result) const
+ {
+ result = 0;
+ uint result_uint;
+
+ if( ToUInt(result_uint, false) )
+ return 1;
+
+ result = static_cast( result_uint );
+
+ // the exception for the minimal value
+ if( IsSign() && result_uint == TTMATH_UINT_HIGHEST_BIT )
+ return 0;
+
+ if( (result_uint & TTMATH_UINT_HIGHEST_BIT) != 0 )
+ // the value is too big
+ return 1;
+
+ if( IsSign() )
+ result = -result;
+
+ return 0;
+ }
+
+
+ /*!
+ this method converts 'this' into 'result'
+
+ if the value is too big this method returns a carry (1)
+ */
+ template
+ uint ToInt(Int & result) const
+ {
+ result.SetZero();
+
+ if( IsZero() )
+ return 0;
+
+ sint maxbit = -sint(man*TTMATH_BITS_PER_UINT);
+
+ if( exponent > maxbit + sint(int_size*TTMATH_BITS_PER_UINT) )
+ // if exponent > (maxbit + sint(int_size*TTMATH_BITS_PER_UINT)) the value can't be passed
+ // into the 'Int' type (it's too big)
+ return 1;
+
+ if( exponent <= maxbit )
+ // our value is from range (-1,1) and we return zero
+ return 0;
+
+ sint how_many_bits = exponent.ToInt();
+
+ if( how_many_bits < 0 )
+ {
+ how_many_bits = -how_many_bits;
+ uint index = how_many_bits / TTMATH_BITS_PER_UINT;
+
+ UInt mantissa_temp(mantissa);
+ mantissa_temp.Rcr( how_many_bits % TTMATH_BITS_PER_UINT, 0 );
+
+ for(uint i=index, a=0 ; i min;
+ min.SetMin();
+
+ if( result == min )
+ return 0;
+ }
+
+ if( (result.table[int_size-1] & TTMATH_UINT_HIGHEST_BIT) != 0 )
+ // the value is too big
+ return 1;
+
+ if( IsSign() )
+ result.ChangeSign();
+
+ return 0;
+ }
+
+
+ /*!
+ a method for converting 'uint' to this class
+ */
+ void FromUInt(uint value)
+ {
+ info = 0;
+
+ for(uint i=0 ; i> 20;
+ uint m1 = ((temp.u[1] & 0xFFFFFu) << 11) | (temp.u[0] >> 21);
+ uint m2 = temp.u[0] << 11;
+
+ if( e == 2047 )
+ {
+ // If E=2047 and F is nonzero, then V=NaN ("Not a number")
+ // If E=2047 and F is zero and S is 1, then V=-Infinity
+ // If E=2047 and F is zero and S is 0, then V=Infinity
+
+ // we do not support -Infinity and +Infinity
+ // we assume that there is always NaN
+
+ SetNan();
+ }
+ else
+ if( e > 0 )
+ {
+ // If 0 m;
+ m.table[1] = m1;
+ m.table[0] = m2;
+ uint moved = m.CompensationToLeft();
+
+ FromDouble_SetExpAndMan((temp.u[1] & 0x80000000u) != 0,
+ e - 1022 - man*TTMATH_BITS_PER_UINT + 1 - moved, 0,
+ m.table[1], m.table[2]);
+ }
+ else
+ {
+ // If E=0 and F is zero and S is 1, then V=-0
+ // If E=0 and F is zero and S is 0, then V=0
+
+ // we do not support -0 or 0, only is one 0
+ SetZero();
+ }
+ }
+ }
+
+
+private:
+
+ void FromDouble_SetExpAndMan(bool is_sign, int e, uint mhighest, uint m1, uint m2)
+ {
+ exponent = e;
+
+ if( man > 1 )
+ {
+ mantissa.table[man-1] = m1 | mhighest;
+ mantissa.table[sint(man-2)] = m2;
+ // although man>1 we're using casting into sint
+ // to get rid from a warning which generates Microsoft Visual:
+ // warning C4307: '*' : integral constant overflow
+
+ for(uint i=0 ; i> 52;
+ uint m = (temp.u & 0xFFFFFFFFFFFFFul) << 11;
+
+ if( e == 2047 )
+ {
+ // If E=2047 and F is nonzero, then V=NaN ("Not a number")
+ // If E=2047 and F is zero and S is 1, then V=-Infinity
+ // If E=2047 and F is zero and S is 0, then V=Infinity
+
+ // we do not support -Infinity and +Infinity
+ // we assume that there is always NaN
+
+ SetNan();
+ }
+ else
+ if( e > 0 )
+ {
+ // If 0= 1024 - e_correction )
+ {
+ // +/- infinity
+ result = ToDouble_SetDouble( 0, 2047, 0, true);
+
+ return 1;
+ }
+ else
+ if( exponent <= -1023 - 52 - e_correction )
+ {
+ // too small value - we assume that there'll be a zero
+ result = 0;
+
+ // and return a carry
+ return 1;
+ }
+
+ sint e = exponent.ToInt() + e_correction;
+
+ if( e <= -1023 )
+ {
+ // -1023-52 < e <= -1023 (unnormalized value)
+ result = ToDouble_SetDouble( IsSign(), 0, -(e + 1023));
+ }
+ else
+ {
+ // -1023 < e < 1024
+ result = ToDouble_SetDouble( IsSign(), e + 1023, -1);
+ }
+
+ return 0;
+ }
+
+private:
+
+#ifdef TTMATH_PLATFORM32
+
+ // 32bit platforms
+ double ToDouble_SetDouble(bool is_sign, uint e, sint move, bool infinity = false, bool nan = false) const
+ {
+ union
+ {
+ double d;
+ uint u[2]; // two 32bit words
+ } temp;
+
+ temp.u[0] = temp.u[1] = 0;
+
+ if( is_sign )
+ temp.u[1] |= 0x80000000u;
+
+ temp.u[1] |= (e << 20) & 0x7FF00000u;
+
+ if( nan )
+ {
+ temp.u[0] |= 1;
+ return temp.d;
+ }
+
+ if( infinity )
+ return temp.d;
+
+ UInt<2> m;
+ m.table[1] = mantissa.table[man-1];
+ m.table[0] = ( man > 1 ) ? mantissa.table[sint(man-2)] : 0;
+ // although man>1 we're using casting into sint
+ // to get rid from a warning which generates Microsoft Visual:
+ // warning C4307: '*' : integral constant overflow
+
+ m.Rcr( 12 + move );
+ m.table[1] &= 0xFFFFFu; // cutting the 20 bit (when 'move' was -1)
+
+ temp.u[1] |= m.table[1];
+ temp.u[0] |= m.table[0];
+
+ return temp.d;
+ }
+
+#else
+
+ // 64bit platforms
+ double ToDouble_SetDouble(bool is_sign, uint e, sint move, bool infinity = false, bool nan = false) const
+ {
+ union
+ {
+ double d;
+ uint u; // 64bit word
+ } temp;
+
+ temp.u = 0;
+
+ if( is_sign )
+ temp.u |= 0x8000000000000000ul;
+
+ temp.u |= (e << 52) & 0x7FF0000000000000ul;
+
+ if( nan )
+ {
+ temp.u |= 1;
+ return temp.d;
+ }
+
+ if( infinity )
+ return temp.d;
+
+ uint m = mantissa.table[man-1];
+
+ m >>= ( 12 + move );
+ m &= 0xFFFFFFFFFFFFFul; // cutting the 20 bit (when 'move' was -1)
+ temp.u |= m;
+
+ return temp.d;
+ }
+
+#endif
+
+
+public:
+
+
+ /*!
+ an operator= for converting 'sint' to this class
+ */
+ Big & operator=(sint value)
+ {
+ FromInt(value);
+
+ return *this;
+ }
+
+
+ /*!
+ an operator= for converting 'uint' to this class
+ */
+ Big & operator=(uint value)
+ {
+ FromUInt(value);
+
+ return *this;
+ }
+
+
+ /*!
+ an operator= for converting 'double' to this class
+ */
+ Big & operator=(double value)
+ {
+ FromDouble(value);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting 'sint' to this class
+ */
+ Big(sint value)
+ {
+ FromInt(value);
+ }
+
+ /*!
+ a constructor for converting 'uint' to this class
+ */
+ Big(uint value)
+ {
+ FromUInt(value);
+ }
+
+
+ /*!
+ a constructor for converting 'double' to this class
+ */
+ Big(double value)
+ {
+ FromDouble(value);
+ }
+
+
+#ifdef TTMATH_PLATFORM64
+
+ /*!
+ in 64bit platforms we must define additional operators and contructors
+ in order to allow a user initializing the objects in this way:
+ Big<...> type = 20;
+ or
+ Big<...> type;
+ type = 30;
+
+ decimal constants such as 20, 30 etc. are integer literal of type int,
+ if the value is greater it can even be long int,
+ 0 is an octal integer of type int
+ (ISO 14882 p2.13.1 Integer literals)
+ */
+
+ /*!
+ an operator= for converting 'signed int' to this class
+ ***this operator is created only on a 64bit platform***
+ it takes one argument of 32bit
+
+
+ */
+ Big & operator=(signed int value)
+ {
+ FromInt(sint(value));
+
+ return *this;
+ }
+
+
+ /*!
+ an operator= for converting 'unsigned int' to this class
+ ***this operator is created only on a 64bit platform***
+ it takes one argument of 32bit
+ */
+ Big & operator=(unsigned int value)
+ {
+ FromUInt(uint(value));
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting 'signed int' to this class
+ ***this constructor is created only on a 64bit platform***
+ it takes one argument of 32bit
+ */
+ Big(signed int value)
+ {
+ FromInt(sint(value));
+ }
+
+ /*!
+ a constructor for converting 'unsigned int' to this class
+ ***this constructor is created only on a 64bit platform***
+ it takes one argument of 32bit
+ */
+ Big(unsigned int value)
+ {
+ FromUInt(uint(value));
+ }
+
+#endif
+
+private:
+
+ /*!
+ an auxiliary method for converting from UInt and Int
+
+ we assume that there'll never be a carry here
+ (we have an exponent and the value in Big can be bigger than
+ that one from the UInt)
+ */
+ template
+ void FromUIntOrInt(const UInt & value, sint compensation)
+ {
+ uint minimum_size = (int_size < man)? int_size : man;
+ exponent = (sint(int_size)-sint(man)) * sint(TTMATH_BITS_PER_UINT) - compensation;
+
+ // copying the highest words
+ uint i;
+ for(i=1 ; i<=minimum_size ; ++i)
+ mantissa.table[man-i] = value.table[int_size-i];
+
+ // setting the rest of mantissa.table into zero (if some has left)
+ for( ; i<=man ; ++i)
+ mantissa.table[man-i] = 0;
+
+ // the highest bit is either one or zero (when the whole mantissa is zero)
+ // we can only call CorrectZero()
+ CorrectZero();
+ }
+
+
+public:
+
+
+ /*!
+ a method for converting from 'UInt' to this class
+ */
+ template
+ void FromUInt(UInt value)
+ {
+ info = 0;
+ sint compensation = (sint)value.CompensationToLeft();
+
+ return FromUIntOrInt(value, compensation);
+ }
+
+
+ /*!
+ a method for converting from 'Int' to this class
+ */
+ template
+ void FromInt(Int value)
+ {
+ info = 0;
+ bool is_sign = false;
+
+ if( value.IsSign() )
+ {
+ value.ChangeSign();
+ is_sign = true;
+ }
+
+ sint compensation = (sint)value.CompensationToLeft();
+ FromUIntOrInt(value, compensation);
+
+ if( is_sign )
+ SetSign();
+ }
+
+
+ /*!
+ an operator= for converting from 'Int' to this class
+ */
+ template
+ Big & operator=(const Int & value)
+ {
+ FromInt(value);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting from 'Int' to this class
+ */
+ template
+ Big(const Int & value)
+ {
+ FromInt(value);
+ }
+
+
+ /*!
+ an operator= for converting from 'UInt' to this class
+ */
+ template
+ Big & operator=(const UInt & value)
+ {
+ FromUInt(value);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting from 'UInt' to this class
+ */
+ template
+ Big(const UInt & value)
+ {
+ FromUInt(value);
+ }
+
+
+ /*!
+ an operator= for converting from 'Big' to this class
+ */
+ template
+ Big & operator=(const Big & value)
+ {
+ FromBig(value);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting from 'Big' to this class
+ */
+ template
+ Big(const Big & value)
+ {
+ FromBig(value);
+ }
+
+
+ /*!
+ a default constructor
+
+ we don't set any of the members to zero
+ only NaN flag is set
+ */
+ Big()
+ {
+ info = TTMATH_BIG_NAN;
+
+ /*
+ we're directly setting 'info' (instead of calling SetNan())
+ in order to get rid of a warning saying that 'info' is uninitialized
+ */
+ }
+
+
+ /*!
+ a destructor
+ */
+ ~Big()
+ {
+ }
+
+
+ /*!
+ the default assignment operator
+ */
+ Big & operator=(const Big & value)
+ {
+ info = value.info;
+ exponent = value.exponent;
+ mantissa = value.mantissa;
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for copying from another object of this class
+ */
+
+ Big(const Big & value)
+ {
+ operator=(value);
+ }
+
+
+
+ /*!
+ a method for converting into a string
+ struct Conv is defined in ttmathtypes.h, look there for more information about parameters
+
+ output:
+ return value:
+ 0 - ok and 'result' will be an object of type std::string (or std::wstring) which holds the value
+ 1 - if there is a carry (it shoudn't be in a normal situation - if it is that means there
+ is somewhere an error in the library)
+ */
+ uint ToString( std::string & result,
+ uint base = 10,
+ bool scient = false,
+ sint scient_from = 15,
+ sint round = -1,
+ bool trim_zeroes = true,
+ wchar_t comma = '.' ) const
+ {
+ Conv conv;
+
+ conv.base = base;
+ conv.scient = scient;
+ conv.scient_from = scient_from;
+ conv.round = round;
+ conv.trim_zeroes = trim_zeroes;
+ conv.comma = static_cast(comma);
+
+ return ToStringBase(result, conv);
+ }
+
+
+ /*!
+ a method for converting into a string
+ struct Conv is defined in ttmathtypes.h, look there for more information about parameters
+ */
+ uint ToString( std::wstring & result,
+ uint base = 10,
+ bool scient = false,
+ sint scient_from = 15,
+ sint round = -1,
+ bool trim_zeroes = true,
+ wchar_t comma = '.' ) const
+ {
+ Conv conv;
+
+ conv.base = base;
+ conv.scient = scient;
+ conv.scient_from = scient_from;
+ conv.round = round;
+ conv.trim_zeroes = trim_zeroes;
+ conv.comma = static_cast(comma);
+
+ return ToStringBase(result, conv);
+ }
+
+
+ /*!
+ a method for converting into a string
+ struct Conv is defined in ttmathtypes.h, look there for more information about parameters
+ */
+ uint ToString(std::string & result, const Conv & conv) const
+ {
+ return ToStringBase(result, conv);
+ }
+
+
+ /*!
+ a method for converting into a string
+ struct Conv is defined in ttmathtypes.h, look there for more information about parameters
+ */
+ uint ToString(std::wstring & result, const Conv & conv) const
+ {
+ return ToStringBase(result, conv);
+ }
+
+
+ /*!
+ a method for converting into a string
+ struct Conv is defined in ttmathtypes.h, look there for more information about parameters
+ */
+ std::string ToString(const Conv & conv) const
+ {
+ std::string result;
+ ToStringBase(result, conv);
+
+ return result;
+ }
+
+
+ /*!
+ a method for converting into a string
+ struct Conv is defined in ttmathtypes.h, look there for more information about parameters
+ */
+ std::string ToString() const
+ {
+ Conv conv;
+
+ return ToString(conv);
+ }
+
+
+ /*!
+ a method for converting into a string
+ struct Conv is defined in ttmathtypes.h, look there for more information about parameters
+ */
+ std::wstring ToWString(const Conv & conv) const
+ {
+ std::wstring result;
+ ToStringBase(result, conv);
+
+ return result;
+ }
+
+
+ /*!
+ a method for converting into a string
+ struct Conv is defined in ttmathtypes.h, look there for more information about parameters
+ */
+ std::wstring ToWString() const
+ {
+ Conv conv;
+
+ return ToWString(conv);
+ }
+
+
+private:
+
+
+ /*!
+ an auxiliary method for converting into the string
+ */
+ template
+ uint ToStringBase(string_type & result, const Conv & conv) const
+ {
+ static char error_overflow_msg[] = "overflow";
+ static char error_nan_msg[] = "NaN";
+ result.erase();
+
+ if( IsNan() )
+ {
+ Misc::AssignString(result, error_nan_msg);
+ return 0;
+ }
+
+ if( conv.base<2 || conv.base>16 )
+ {
+ Misc::AssignString(result, error_overflow_msg);
+ return 1;
+ }
+
+ if( IsZero() )
+ {
+ result = '0';
+
+ return 0;
+ }
+
+ /*
+ since 'base' is greater or equal 2 that 'new_exp' of type 'Int' should
+ hold the new value of exponent but we're using 'Int' because
+ if the value for example would be 'max()' then we couldn't show it
+
+ max() -> 11111111 * 2 ^ 11111111111 (bin)(the mantissa and exponent have all bits set)
+ if we were using 'Int' we couldn't show it in this format:
+ 1,1111111 * 2 ^ 11111111111 (bin)
+ because we have to add something to the mantissa and because
+ mantissa is full we can't do it and it'll be a carry
+ (look at ToString_SetCommaAndExponent(...))
+
+ when the base would be greater than two (for example 10)
+ we could use 'Int' here
+ */
+ Int