diff --git a/appveyor.yml b/appveyor.yml
index 99f821570d..fa9564dd4c 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -27,9 +27,11 @@ install:
- cmd: for /r %%i in (*.bat) do unix2dos "%%i"
- cmd: if "%PLATFORM%" == "x64" set FB_PROCESSOR_ARCHITECTURE=AMD64
- cmd: if "%PLATFORM%" == "x64" set FB_OUTPUT_SUFFIX=x64
+ - cmd: if "%PLATFORM%" == "x64" set FB_VS_ARCH=amd64
- cmd: if "%PLATFORM%" == "x86" set FB_PROCESSOR_ARCHITECTURE=x86
- cmd: if "%PLATFORM%" == "x86" set FB_OUTPUT_SUFFIX=win32
- - cmd: if "%APPVEYOR_BUILD_WORKER_IMAGE%" == "Visual Studio 2017" call "%ProgramFiles(x86)%\Microsoft Visual Studio\2017\Community\Common7\Tools\VsDevCmd.bat"
+ - cmd: if "%PLATFORM%" == "x86" set FB_VS_ARCH=x86
+ - cmd: if "%APPVEYOR_BUILD_WORKER_IMAGE%" == "Visual Studio 2017" call "%ProgramFiles(x86)%\Microsoft Visual Studio\2017\Community\Common7\Tools\VsDevCmd.bat" -arch=%FB_VS_ARCH%
- cmd: cd builds\win32
- cmd: run_all.bat JUSTBUILD
- cmd: set ARTIFACTS_PATH=output_%FB_OUTPUT_SUFFIX%
diff --git a/builds/win32/make_boot.bat b/builds/win32/make_boot.bat
index b4d62b90dd..882c0829e7 100644
--- a/builds/win32/make_boot.bat
+++ b/builds/win32/make_boot.bat
@@ -41,6 +41,9 @@ if "%ERRLEV%"=="1" goto :END
call :decNumber
if "%ERRLEV%"=="1" goto :END
+if "%FB_TARGET_PLATFORM%"=="x64" call :ttmath
+if "%ERRLEV%"=="1" goto :END
+
call :re2
if "%ERRLEV%"=="1" goto :END
@@ -162,6 +165,23 @@ if errorlevel 1 call :boot2 decNumber_%FB_OBJ_DIR%
@call set_build_target.bat %*
goto :EOF
+::===================
+:: BUILD ttmath
+:ttmath
+@echo.
+@call set_build_target.bat %* RELEASE
+@echo Building ttmath (%FB_OBJ_DIR%)...
+@mkdir %FB_TEMP_DIR%\..\%FB_OBJ_DIR%\common 2>nul
+@ml64.exe /c /Fo %FB_TEMP_DIR%\..\%FB_OBJ_DIR%\common\ttmathuint_x86_64_msvc.obj %FB_ROOT_PATH%\extern\ttmath\ttmathuint_x86_64_msvc.asm
+if errorlevel 1 call :boot2 ttmath_%FB_OBJ_DIR%
+@call set_build_target.bat %* DEBUG
+@echo Building decNumber (%FB_OBJ_DIR%)...
+@mkdir %FB_TEMP_DIR%\..\%FB_OBJ_DIR%\common 2>nul
+@ml64.exe /c /Zi /Fo %FB_TEMP_DIR%\..\%FB_OBJ_DIR%\common\ttmathuint_x86_64_msvc.obj %FB_ROOT_PATH%\extern\ttmath\ttmathuint_x86_64_msvc.asm
+if errorlevel 1 call :boot2 ttmath_%FB_OBJ_DIR%
+@call set_build_target.bat %*
+goto :EOF
+
::===================
:: BUILD re2
:re2
diff --git a/builds/win32/msvc10/common.vcxproj b/builds/win32/msvc10/common.vcxproj
index 99df722d88..f8f9271413 100644
--- a/builds/win32/msvc10/common.vcxproj
+++ b/builds/win32/msvc10/common.vcxproj
@@ -72,6 +72,7 @@
+
@@ -180,6 +181,7 @@
+
@@ -335,7 +337,7 @@
0x041d
- ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;%(AdditionalDependencies)
+ ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;$(TargetDir)\ttmathuint_x86_64_msvc.obj;%(AdditionalDependencies)
@@ -353,7 +355,7 @@
0x041d
- ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;%(AdditionalDependencies)
+ ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;$(TargetDir)\ttmathuint_x86_64_msvc.obj;%(AdditionalDependencies)
diff --git a/builds/win32/msvc10/common.vcxproj.filters b/builds/win32/msvc10/common.vcxproj.filters
index 978f26a315..f6d9bf98c0 100644
--- a/builds/win32/msvc10/common.vcxproj.filters
+++ b/builds/win32/msvc10/common.vcxproj.filters
@@ -240,6 +240,9 @@
common
+
+ common
+
@@ -581,5 +584,8 @@
headers
+
+ headers
+
\ No newline at end of file
diff --git a/builds/win32/msvc12/common.vcxproj b/builds/win32/msvc12/common.vcxproj
index 0492407080..8b2e9805a0 100644
--- a/builds/win32/msvc12/common.vcxproj
+++ b/builds/win32/msvc12/common.vcxproj
@@ -68,6 +68,7 @@
+
@@ -176,6 +177,7 @@
+
@@ -335,7 +337,7 @@
0x041d
- ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;%(AdditionalDependencies)
+ ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;$(TargetDir)\ttmathuint_x86_64_msvc.obj;%(AdditionalDependencies)
@@ -353,7 +355,7 @@
0x041d
- ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;%(AdditionalDependencies)
+ ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;$(TargetDir)\ttmathuint_x86_64_msvc.obj;%(AdditionalDependencies)
diff --git a/builds/win32/msvc12/common.vcxproj.filters b/builds/win32/msvc12/common.vcxproj.filters
index 328eabca1d..ffc6d9c8fe 100644
--- a/builds/win32/msvc12/common.vcxproj.filters
+++ b/builds/win32/msvc12/common.vcxproj.filters
@@ -240,6 +240,9 @@
common
+
+ common
+
@@ -581,5 +584,8 @@
headers
+
+ headers
+
\ No newline at end of file
diff --git a/builds/win32/msvc14/common.vcxproj b/builds/win32/msvc14/common.vcxproj
index 77823eb058..615a314773 100644
--- a/builds/win32/msvc14/common.vcxproj
+++ b/builds/win32/msvc14/common.vcxproj
@@ -68,6 +68,7 @@
+
@@ -177,6 +178,7 @@
+
@@ -337,7 +339,7 @@
0x041d
- ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;../../../extern/re2/builds/$(PlatformName)\$(Configuration)\re2.lib;%(AdditionalDependencies)
+ ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;../../../extern/re2/builds/$(PlatformName)\$(Configuration)\re2.lib;$(TargetDir)\ttmathuint_x86_64_msvc.obj;%(AdditionalDependencies)
@@ -355,7 +357,7 @@
0x041d
- ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;../../../extern/re2/builds/$(PlatformName)\$(Configuration)\re2.lib;%(AdditionalDependencies)
+ ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;../../../extern/re2/builds/$(PlatformName)\$(Configuration)\re2.lib;$(TargetDir)\ttmathuint_x86_64_msvc.obj;%(AdditionalDependencies)
diff --git a/builds/win32/msvc14/common.vcxproj.filters b/builds/win32/msvc14/common.vcxproj.filters
index a0b2d29362..4da18bbba5 100644
--- a/builds/win32/msvc14/common.vcxproj.filters
+++ b/builds/win32/msvc14/common.vcxproj.filters
@@ -243,6 +243,9 @@
common
+
+ common
+
@@ -587,5 +590,8 @@
headers
+
+ headers
+
diff --git a/builds/win32/msvc15/common.vcxproj b/builds/win32/msvc15/common.vcxproj
index c2ccd9f972..e8d6f56042 100644
--- a/builds/win32/msvc15/common.vcxproj
+++ b/builds/win32/msvc15/common.vcxproj
@@ -68,6 +68,7 @@
+
@@ -177,6 +178,7 @@
+
@@ -338,7 +340,7 @@
0x041d
- ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;../../../extern/re2/builds/$(PlatformName)\$(Configuration)\re2.lib;%(AdditionalDependencies)
+ ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;$(TargetDir)\ttmathuint_x86_64_msvc.obj;../../../extern/re2/builds/$(PlatformName)\$(Configuration)\re2.lib;%(AdditionalDependencies)
@@ -356,7 +358,7 @@
0x041d
- ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;../../../extern/re2/builds/$(PlatformName)\$(Configuration)\re2.lib;%(AdditionalDependencies)
+ ws2_32.lib;../../../extern/libtommath/lib/$(PlatformName)\$(Configuration)\tommath.lib;../../../extern/libtomcrypt/lib/$(PlatformName)\$(Configuration)\tomcrypt.lib;../../../extern/decNumber/lib/$(PlatformName)\$(Configuration)\decnumber.lib;../../../extern/re2/builds/$(PlatformName)\$(Configuration)\re2.lib;$(TargetDir)\ttmathuint_x86_64_msvc.obj;%(AdditionalDependencies)
diff --git a/builds/win32/msvc15/common.vcxproj.filters b/builds/win32/msvc15/common.vcxproj.filters
index 0944d0518d..776d4a9db7 100644
--- a/builds/win32/msvc15/common.vcxproj.filters
+++ b/builds/win32/msvc15/common.vcxproj.filters
@@ -243,6 +243,9 @@
common
+
+ common
+
@@ -587,5 +590,8 @@
headers
+
+ headers
+
\ No newline at end of file
diff --git a/doc/sql.extensions/README.data_types b/doc/sql.extensions/README.data_types
index 4df707dc4c..ac6dbc10b0 100644
--- a/doc/sql.extensions/README.data_types
+++ b/doc/sql.extensions/README.data_types
@@ -200,7 +200,7 @@ Enhancement in precision of calculations with NUMERIC/DECIMAL (FB 4.0)
--------------
Function:
- Maximum precision of NUMERIC and DECIMAL data types is increased to 34 digits.
+ Maximum precision of NUMERIC and DECIMAL data types is increased to 38 digits.
Author:
Alex Peshkoff
@@ -208,20 +208,35 @@ Enhancement in precision of calculations with NUMERIC/DECIMAL (FB 4.0)
Syntax rules:
NUMERIC ( P {, N} )
DECIMAL ( P {, N} )
- where P is precision (P <= 34, was limited prior with 18 digits) and N is optional number
+ where P is precision (P <= 38, was limited prior with 18 digits) and N is optional number
of digits after decimal separator (as before).
Storage:
- 128-bit, format according to IEEE 754.
+ 128-bit signed integer.
Example(s):
1. DECLARE VARIABLE VAR1 DECIMAL(25);
- 2. CREATE TABLE TABLE1 (FIELD1 NUMERIC(34, 17));
+ 2. CREATE TABLE TABLE1 (FIELD1 NUMERIC(38, 19));
Note(s):
Numerics with precision less than 19 digits use SMALLINT, INTEGER, BIGINT or DOUBLE PRECISION
as base datatype depending upon number of digits and dialect. When precision is between 19 and
- 34 digits DECFLOAT(34) is used for it. Actual precision is always increased to 34 digits. For
- complex calculations such digits are casted (internally, in trivial way) to DECFLOAT(34) and
- the result of various math (log, exp, etc.) and aggregate functions using high precision
- numeric argument is DECFLOAT(34).
+ 38 digits 128-bit integer is used for it. Actual precision is always increased to 38 digits.
+ For complex calculations such digits are casted (internally) to DECFLOAT(34) and the result of
+ various math (log, exp, etc.) and aggregate functions using high precision numeric argument is
+ DECFLOAT(34).
+
+ SET INT128 BIND - controls how are INT128 values represented in outer
+ world (i.e. in messages or in XSQLDA). Valid binding types are: NATIVE (use 128-bit
+ binary representation), CHAR/CHARACTER (use ASCII string), DOUBLE PRECISION (use
+ 8-byte FP representation - same as used for DOUBLE PRECISION fields) or BIGINT
+ with possible comma-separated SCALE clause (i.e. 'BIGINT, 3'). Various bindings
+ are useful if one plans to use 128-bit integers with some old client not supporting
+ native format. One can choose between strings (ideal precision, but poor support
+ for further processing), floating point values (ideal support for further processing
+ but poor precision) or scaled integers (good support for further processing and
+ required precision but range of values is very limited). When using in a tool like
+ generic purporse GUI client choice of CHAR binding is OK in most cases. By default
+ NATIVE binding is used.
+ The initial configuration may be specified with DPB isc_dpb_int128_bind followed
+ by a string with its value (case does not matter).
diff --git a/extern/ttmath/ttmath.h b/extern/ttmath/ttmath.h
new file mode 100644
index 0000000000..9b9e8ee6e5
--- /dev/null
+++ b/extern/ttmath/ttmath.h
@@ -0,0 +1,2843 @@
+/*
+ * This file is a part of TTMath Bignum Library
+ * and is distributed under the (new) BSD licence.
+ * Author: Tomasz Sowa
+ */
+
+/*
+ * Copyright (c) 2006-2012, 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
+
+#define TTMATH_DONT_USE_WCHAR
+
+#include "ttmathint.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() )
+ {
+ if( err ) *err = err_improper_argument;
+ return x;
+ }
+
+ if( base.IsNan() )
+ {
+ if( err ) *err = err_improper_argument;
+ return base;
+ }
+
+ 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 x;
+ }
+
+ 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 x;
+ }
+
+ 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 x;
+
+ // 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 x;
+ }
+
+ 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 x;
+ }
+
+ 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;
+
+ delimiter.SetZeroNan(); // not needed, only to get rid of GCC warning about an uninitialized variable
+
+ return delimiter;
+ }
+
+ 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 x;
+ }
+
+ 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 x;
+ }
+
+ 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 x;
+ }
+
+ 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 x;
+ }
+
+ 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;
+
+ x.SetNan();
+
+ return x;
+ }
+
+ 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;
+
+ x.SetNan();
+
+ return x;
+ }
+
+ 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;
+
+ a.SetNan();
+
+ return a;
+ }
+
+ 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 n;
+ }
+
+ 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;
+
+ x.SetNan();
+
+ return x;
+ }
+
+ 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
+
+
+#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/extern/ttmath/ttmathint.h b/extern/ttmath/ttmathint.h
new file mode 100644
index 0000000000..8ad0189f93
--- /dev/null
+++ b/extern/ttmath/ttmathint.h
@@ -0,0 +1,1917 @@
+/*
+ * This file is a part of TTMath Bignum Library
+ * and is distributed under the (new) BSD licence.
+ * Author: Tomasz Sowa
+ */
+
+/*
+ * Copyright (c) 2006-2011, 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 headerfilettmathint
+#define headerfilettmathint
+
+/*!
+ \file ttmathint.h
+ \brief template class Int
+*/
+
+#include "ttmathuint.h"
+
+namespace ttmath
+{
+
+
+/*!
+ \brief Int implements a big integer value with a sign
+
+ value_size - how many bytes specify our value
+ on 32bit platforms: value_size=1 -> 4 bytes -> 32 bits
+ on 64bit platforms: value_size=1 -> 8 bytes -> 64 bits
+ value_size = 1,2,3,4,5,6....
+*/
+template
+class Int : public UInt
+{
+public:
+
+ /*!
+ this method sets the max value which this class can hold
+ (all bits will be one besides the last one)
+ */
+ void SetMax()
+ {
+ UInt::SetMax();
+ UInt::table[value_size-1] = ~ TTMATH_UINT_HIGHEST_BIT;
+ }
+
+
+ /*!
+ this method sets the min value which this class can hold
+ (all bits will be zero besides the last one which is one)
+ */
+ void SetMin()
+ {
+ UInt::SetZero();
+ UInt::table[value_size-1] = TTMATH_UINT_HIGHEST_BIT;
+ }
+
+
+ /*!
+ this method sets -1 as the value
+ (-1 is equal the max value in an unsigned type)
+ */
+ void SetSignOne()
+ {
+ UInt::SetMax();
+ }
+
+
+ /*!
+ we change the sign of the value
+
+ if it isn't possible to change the sign this method returns 1
+ else return 0 and changing the sign
+ */
+ uint ChangeSign()
+ {
+ /*
+ if the value is equal that one which has been returned from SetMin
+ (only the highest bit is set) that means we can't change sign
+ because the value is too big (bigger about one)
+
+ e.g. when value_size = 1 and value is -2147483648 we can't change it to the
+ 2147483648 because the max value which can be held is 2147483647
+
+ we don't change the value and we're using this fact somewhere in some methods
+ (if we look on our value without the sign we get the correct value
+ eg. -2147483648 in Int<1> will be 2147483648 on the UInt<1> type)
+ */
+ if( UInt::IsOnlyTheHighestBitSet() )
+ return 1;
+
+ UInt temp(*this);
+ UInt::SetZero();
+ UInt::Sub(temp);
+
+ return 0;
+ }
+
+
+
+ /*!
+ this method sets the sign
+
+ e.g. 1 -> -1
+ -2 -> -2
+
+ from a positive value we make a negative value,
+ if the value is negative we do nothing
+ */
+ void SetSign()
+ {
+ if( IsSign() )
+ return;
+
+ ChangeSign();
+ }
+
+
+
+ /*!
+ this method returns true if there's the sign
+
+ (the highest bit will be converted to the bool)
+ */
+ bool IsSign() const
+ {
+ return UInt::IsTheHighestBitSet();
+ }
+
+
+
+ /*!
+ it sets an absolute value
+
+ it can return carry (1) (look on ChangeSign() for details)
+ */
+ uint Abs()
+ {
+ if( !IsSign() )
+ return 0;
+
+ return ChangeSign();
+ }
+
+
+
+
+ /*!
+ *
+ * basic mathematic functions
+ *
+ */
+
+private:
+
+ uint CorrectCarryAfterAdding(bool p1_is_sign, bool p2_is_sign)
+ {
+ if( !p1_is_sign && !p2_is_sign )
+ {
+ if( UInt::IsTheHighestBitSet() )
+ return 1;
+ }
+
+ if( p1_is_sign && p2_is_sign )
+ {
+ if( ! UInt::IsTheHighestBitSet() )
+ return 1;
+ }
+
+ return 0;
+ }
+
+
+public:
+
+ /*!
+ this method adds two value with a sign and returns a carry
+
+ we're using methods from the base class because values are stored with U2
+ we must only make the carry correction
+
+ this = p1(=this) + p2
+
+ when p1>=0 i p2>=0 carry is set when the highest bit of value is set
+ when p1<0 i p2<0 carry is set when the highest bit of value is clear
+ when p1>=0 i p2<0 carry will never be set
+ when p1<0 i p2>=0 carry will never be set
+ */
+ uint Add(const Int & ss2)
+ {
+ bool p1_is_sign = IsSign();
+ bool p2_is_sign = ss2.IsSign();
+
+ UInt::Add(ss2);
+
+ return CorrectCarryAfterAdding(p1_is_sign, p2_is_sign);
+ }
+
+
+ /*!
+ this method adds one *unsigned* word (at a specific position)
+ and returns a carry (if it was)
+
+ look at a description in UInt<>::AddInt(...)
+ */
+ uint AddInt(uint value, uint index = 0)
+ {
+ bool p1_is_sign = IsSign();
+
+ UInt::AddInt(value, index);
+
+ return CorrectCarryAfterAdding(p1_is_sign, false);
+ }
+
+
+ /*!
+ this method adds two *unsigned* words to the existing value
+ and these words begin on the 'index' position
+
+ index should be equal or smaller than value_size-2 (index <= value_size-2)
+ x1 - lower word, x2 - higher word
+
+ look at a description in UInt<>::AddTwoInts(...)
+ */
+ uint AddTwoInts(uint x2, uint x1, uint index)
+ {
+ bool p1_is_sign = IsSign();
+
+ UInt::AddTwoInts(x2, x1, index);
+
+ return CorrectCarryAfterAdding(p1_is_sign, false);
+ }
+
+private:
+
+ uint CorrectCarryAfterSubtracting(bool p1_is_sign, bool p2_is_sign)
+ {
+ if( !p1_is_sign && p2_is_sign )
+ {
+ if( UInt::IsTheHighestBitSet() )
+ return 1;
+ }
+
+ if( p1_is_sign && !p2_is_sign )
+ {
+ if( ! UInt::IsTheHighestBitSet() )
+ return 1;
+ }
+
+ return 0;
+ }
+
+public:
+
+ /*!
+ this method subtracts two values with a sign
+
+ we don't use the previous Add because the method ChangeSign can
+ sometimes return carry
+
+ this = p1(=this) - p2
+
+ when p1>=0 i p2>=0 carry will never be set
+ when p1<0 i p2<0 carry will never be set
+ when p1>=0 i p2<0 carry is set when the highest bit of value is set
+ when p1<0 i p2>=0 carry is set when the highest bit of value is clear
+ */
+ uint Sub(const Int & ss2)
+ {
+ bool p1_is_sign = IsSign();
+ bool p2_is_sign = ss2.IsSign();
+
+ UInt::Sub(ss2);
+
+ return CorrectCarryAfterSubtracting(p1_is_sign, p2_is_sign);
+ }
+
+
+ /*!
+ this method subtracts one *unsigned* word (at a specific position)
+ and returns a carry (if it was)
+ */
+ uint SubInt(uint value, uint index = 0)
+ {
+ bool p1_is_sign = IsSign();
+
+ UInt::SubInt(value, index);
+
+ return CorrectCarryAfterSubtracting(p1_is_sign, false);
+ }
+
+
+ /*!
+ this method adds one to the value and returns carry
+ */
+ uint AddOne()
+ {
+ bool p1_is_sign = IsSign();
+
+ UInt::AddOne();
+
+ return CorrectCarryAfterAdding(p1_is_sign, false);
+ }
+
+
+ /*!
+ this method subtracts one from the value and returns carry
+ */
+ uint SubOne()
+ {
+ bool p1_is_sign = IsSign();
+
+ UInt::SubOne();
+
+ return CorrectCarryAfterSubtracting(p1_is_sign, false);
+ }
+
+
+private:
+
+
+ uint CheckMinCarry(bool ss1_is_sign, bool ss2_is_sign)
+ {
+ /*
+ we have to examine the sign of the result now
+ but if the result is with the sign then:
+ 1. if the signs were the same that means the result is too big
+ (the result must be without a sign)
+ 2. if the signs were different that means if the result
+ is different from that one which has been returned from SetMin()
+ that is carry (result too big) but if the result is equal SetMin()
+ there'll be ok (and the next SetSign will has no effect because
+ the value is actually negative -- look at description of that case
+ in ChangeSign())
+ */
+ if( IsSign() )
+ {
+ if( ss1_is_sign != ss2_is_sign )
+ {
+ /*
+ there can be one case where signs are different and
+ the result will be equal the value from SetMin() (only the highest bit is set)
+ (this situation is ok)
+ */
+ if( !UInt::IsOnlyTheHighestBitSet() )
+ return 1;
+ }
+ else
+ {
+ // signs were the same
+ return 1;
+ }
+ }
+
+ return 0;
+ }
+
+
+public:
+
+
+ /*!
+ multiplication: this = this * ss2
+
+ it can return a carry
+ */
+ uint MulInt(sint ss2)
+ {
+ bool ss1_is_sign, ss2_is_sign;
+ uint c;
+
+ ss1_is_sign = IsSign();
+
+ /*
+ we don't have to check the carry from Abs (values will be correct
+ because next we're using the method MulInt from the base class UInt
+ which is without a sign)
+ */
+ Abs();
+
+ if( ss2 < 0 )
+ {
+ ss2 = -ss2;
+ ss2_is_sign = true;
+ }
+ else
+ {
+ ss2_is_sign = false;
+ }
+
+ c = UInt::MulInt((uint)ss2);
+ c += CheckMinCarry(ss1_is_sign, ss2_is_sign);
+
+ if( ss1_is_sign != ss2_is_sign )
+ SetSign();
+
+ return c;
+ }
+
+
+
+ /*!
+ multiplication this = this * ss2
+
+ it returns carry if the result is too big
+ (we're using the method from the base class but we have to make
+ one correction in account of signs)
+ */
+ uint Mul(Int ss2)
+ {
+ bool ss1_is_sign, ss2_is_sign;
+ uint c;
+
+ ss1_is_sign = IsSign();
+ ss2_is_sign = ss2.IsSign();
+
+ /*
+ we don't have to check the carry from Abs (values will be correct
+ because next we're using the method Mul from the base class UInt
+ which is without a sign)
+ */
+ Abs();
+ ss2.Abs();
+
+ c = UInt::Mul(ss2);
+ c += CheckMinCarry(ss1_is_sign, ss2_is_sign);
+
+ if( ss1_is_sign != ss2_is_sign )
+ SetSign();
+
+ return c;
+ }
+
+
+ /*!
+ division this = this / ss2
+ returned values:
+ 0 - ok
+ 1 - division by zero
+
+ for example: (result means 'this')
+ 20 / 3 --> result: 6 remainder: 2
+ -20 / 3 --> result: -6 remainder: -2
+ 20 / -3 --> result: -6 remainder: 2
+ -20 / -3 --> result: 6 remainder: -2
+
+ in other words: this(old) = ss2 * this(new)(result) + remainder
+ */
+ uint Div(Int ss2, Int * remainder = 0)
+ {
+ bool ss1_is_sign, ss2_is_sign;
+
+ ss1_is_sign = IsSign();
+ ss2_is_sign = ss2.IsSign();
+
+ /*
+ we don't have to test the carry from Abs as well as in Mul
+ */
+ Abs();
+ ss2.Abs();
+
+ uint c = UInt::Div(ss2, remainder);
+
+ if( ss1_is_sign != ss2_is_sign )
+ SetSign();
+
+ if( ss1_is_sign && remainder )
+ remainder->SetSign();
+
+ return c;
+ }
+
+ uint Div(const Int & ss2, Int & remainder)
+ {
+ return Div(ss2, &remainder);
+ }
+
+
+ /*!
+ division this = this / ss2 (ss2 is int)
+ returned values:
+ 0 - ok
+ 1 - division by zero
+
+ for example: (result means 'this')
+ 20 / 3 --> result: 6 remainder: 2
+ -20 / 3 --> result: -6 remainder: -2
+ 20 / -3 --> result: -6 remainder: 2
+ -20 / -3 --> result: 6 remainder: -2
+
+ in other words: this(old) = ss2 * this(new)(result) + remainder
+ */
+ uint DivInt(sint ss2, sint * remainder = 0)
+ {
+ bool ss1_is_sign, ss2_is_sign;
+
+ ss1_is_sign = IsSign();
+
+ /*
+ we don't have to test the carry from Abs as well as in Mul
+ */
+ Abs();
+
+ if( ss2 < 0 )
+ {
+ ss2 = -ss2;
+ ss2_is_sign = true;
+ }
+ else
+ {
+ ss2_is_sign = false;
+ }
+
+ uint rem;
+ uint c = UInt::DivInt((uint)ss2, &rem);
+
+ if( ss1_is_sign != ss2_is_sign )
+ SetSign();
+
+ if( remainder )
+ {
+ if( ss1_is_sign )
+ *remainder = -sint(rem);
+ else
+ *remainder = sint(rem);
+ }
+
+ return c;
+ }
+
+
+ uint DivInt(sint ss2, sint & remainder)
+ {
+ return DivInt(ss2, &remainder);
+ }
+
+
+private:
+
+
+ /*!
+ power this = this ^ pow
+ this can be negative
+ pow is >= 0
+ */
+ uint Pow2(const Int & pow)
+ {
+ bool was_sign = IsSign();
+ uint c = 0;
+
+ if( was_sign )
+ c += Abs();
+
+ uint c_temp = UInt::Pow(pow);
+ if( c_temp > 0 )
+ return c_temp; // c_temp can be: 0, 1 or 2
+
+ if( was_sign && (pow.table[0] & 1) == 1 )
+ // negative value to the power of odd number is negative
+ c += ChangeSign();
+
+ return (c==0)? 0 : 1;
+ }
+
+
+public:
+
+
+ /*!
+ power this = this ^ pow
+
+ return values:
+ 0 - ok
+ 1 - carry
+ 2 - incorrect arguments 0^0 or 0^(-something)
+ */
+ uint Pow(Int pow)
+ {
+ if( !pow.IsSign() )
+ return Pow2(pow);
+
+ if( UInt::IsZero() )
+ // if 'pow' is negative then
+ // 'this' must be different from zero
+ return 2;
+
+ if( pow.ChangeSign() )
+ return 1;
+
+ Int t(*this);
+ uint c_temp = t.Pow2(pow);
+ if( c_temp > 0 )
+ return c_temp;
+
+ UInt::SetOne();
+ if( Div(t) )
+ return 1;
+
+ return 0;
+ }
+
+
+
+ /*!
+ *
+ * convertion methods
+ *
+ */
+private:
+
+
+ /*!
+ an auxiliary method for converting both from UInt and Int
+ */
+ template
+ uint FromUIntOrInt(const UInt & p, bool UInt_type)
+ {
+ uint min_size = (value_size < argument_size)? value_size : argument_size;
+ uint i;
+
+ for(i=0 ; i::table[i] = p.table[i];
+
+
+ if( value_size > argument_size )
+ {
+ uint fill;
+
+ if( UInt_type )
+ fill = 0;
+ else
+ fill = (p.table[argument_size-1] & TTMATH_UINT_HIGHEST_BIT)?
+ TTMATH_UINT_MAX_VALUE : 0;
+
+ // 'this' is longer than 'p'
+ for( ; i::table[i] = fill;
+ }
+ else
+ {
+ uint test = (UInt::table[value_size-1] & TTMATH_UINT_HIGHEST_BIT)?
+ TTMATH_UINT_MAX_VALUE : 0;
+
+ if( UInt_type && test!=0 )
+ return 1;
+
+ for( ; i type into this class
+
+ this operation has mainly sense if the value from p
+ can be held in this type
+
+ it returns a carry if the value 'p' is too big
+ */
+ template
+ uint FromInt(const Int & p)
+ {
+ return FromUIntOrInt(p, false);
+ }
+
+
+ /*!
+ this method converts the sint type into this class
+ */
+ uint FromInt(sint value)
+ {
+ uint fill = ( value<0 ) ? TTMATH_UINT_MAX_VALUE : 0;
+
+ for(uint i=1 ; i::table[i] = fill;
+
+ UInt::table[0] = uint(value);
+
+ // there'll never be a carry here
+ return 0;
+ }
+
+
+ /*!
+ this method converts UInt into this class
+ */
+ template
+ uint FromUInt(const UInt & p)
+ {
+ return FromUIntOrInt(p, true);
+ }
+
+
+ /*!
+ this method converts UInt into this class
+ */
+ template
+ uint FromInt(const UInt & p)
+ {
+ return FromUIntOrInt(p, true);
+ }
+
+
+ /*!
+ this method converts the uint type into this class
+ */
+ uint FromUInt(uint value)
+ {
+ for(uint i=1 ; i::table[i] = 0;
+
+ UInt::table[0] = value;
+
+ // there can be a carry here when the size of this value is equal one word
+ // and the 'value' has the highest bit set
+ if( value_size==1 && (value & TTMATH_UINT_HIGHEST_BIT)!=0 )
+ return 1;
+
+ return 0;
+ }
+
+
+ /*!
+ this method converts the uint type into this class
+ */
+ uint FromInt(uint value)
+ {
+ return FromUInt(value);
+ }
+
+
+ /*!
+ the default assignment operator
+ */
+/* Int & operator=(const Int & p)
+ {
+ FromInt(p);
+
+ return *this;
+ }
+ */
+
+ /*!
+ this operator converts an Int type to this class
+
+ it doesn't return a carry
+ */
+/* template
+ Int & operator=(const Int & p)
+ {
+ FromInt(p);
+
+ return *this;
+ }
+ */
+
+ /*!
+ this method converts the sint type to this class
+ */
+ Int & operator=(sint i)
+ {
+ FromInt(i);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting the uint to this class
+ */
+/* Int(sint i)
+ {
+ FromInt(i);
+ }
+*/
+
+ /*!
+ a copy constructor
+ */
+/* Int(const Int & u)
+ {
+ FromInt(u);
+ }
+*/
+
+ /*!
+ a constructor for copying from another types
+ */
+/* template
+ Int(const Int & u)
+ {
+ // look that 'size' we still set as 'value_size' and not as u.value_size
+ FromInt(u);
+ }
+*/
+
+
+ /*!
+ this operator converts an UInt type to this class
+
+ it doesn't return a carry
+ */
+ template
+ Int & operator=(const UInt & p)
+ {
+ FromUInt(p);
+
+ return *this;
+ }
+
+
+ /*!
+ this method converts the Uint type to this class
+ */
+ Int & operator=(uint i)
+ {
+ FromUInt(i);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting the uint to this class
+ */
+/* Int(uint i)
+ {
+ FromUInt(i);
+ }
+*/
+
+ /*!
+ a constructor for copying from another types
+ */
+/* template
+ Int(const UInt & u)
+ {
+ // look that 'size' we still set as 'value_size' and not as u.value_size
+ FromUInt(u);
+ }
+*/
+
+
+#ifdef TTMATH_PLATFORM32
+
+
+ /*!
+ this method converts unsigned 64 bit int type to this class
+ ***this method is created only on a 32bit platform***
+ */
+ uint FromUInt(ulint n)
+ {
+ uint c = UInt::FromUInt(n);
+
+ if( c )
+ return 1;
+
+ if( value_size == 1 )
+ return ((UInt::table[0] & TTMATH_UINT_HIGHEST_BIT) == 0) ? 0 : 1;
+
+ if( value_size == 2 )
+ return ((UInt::table[1] & TTMATH_UINT_HIGHEST_BIT) == 0) ? 0 : 1;
+
+ return 0;
+ }
+
+
+ /*!
+ this method converts unsigned 64 bit int type to this class
+ ***this method is created only on a 32bit platform***
+ */
+ uint FromInt(ulint n)
+ {
+ return FromUInt(n);
+ }
+
+
+ /*!
+ this method converts signed 64 bit int type to this class
+ ***this method is created only on a 32bit platform***
+ */
+ uint FromInt(slint n)
+ {
+ uint mask = (n < 0) ? TTMATH_UINT_MAX_VALUE : 0;
+
+ UInt::table[0] = (uint)(ulint)n;
+
+ if( value_size == 1 )
+ {
+ if( uint(ulint(n) >> 32) != mask )
+ return 1;
+
+ return ((UInt::table[0] & TTMATH_UINT_HIGHEST_BIT) == (mask & TTMATH_UINT_HIGHEST_BIT)) ? 0 : 1;
+ }
+
+ UInt::table[1] = (uint)(ulint(n) >> 32);
+
+ for(uint i=2 ; i::table[i] = mask;
+
+ return 0;
+ }
+
+
+ /*!
+ this operator converts unsigned 64 bit int type to this class
+ ***this operator is created only on a 32bit platform***
+ */
+ Int & operator=(ulint n)
+ {
+ FromUInt(n);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting unsigned 64 bit int to this class
+ ***this constructor is created only on a 32bit platform***
+ */
+/* Int(ulint n)
+ {
+ FromUInt(n);
+ }
+*/
+
+ /*!
+ this operator converts signed 64 bit int type to this class
+ ***this operator is created only on a 32bit platform***
+ */
+ Int & operator=(slint n)
+ {
+ FromInt(n);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting signed 64 bit int to this class
+ ***this constructor is created only on a 32bit platform***
+ */
+/* Int(slint n)
+ {
+ FromInt(n);
+ }
+*/
+#endif
+
+
+
+
+#ifdef TTMATH_PLATFORM64
+
+ /*!
+ this method converts 32 bit unsigned int type to this class
+ ***this operator is created only on a 64bit platform***
+ */
+ uint FromUInt(unsigned int i)
+ {
+ return FromUInt(uint(i));
+ }
+
+
+ /*!
+ this method converts 32 bit unsigned int type to this class
+ ***this operator is created only on a 64bit platform***
+ */
+ uint FromInt(unsigned int i)
+ {
+ return FromUInt(i);
+ }
+
+
+ /*!
+ this method converts 32 bit signed int type to this class
+ ***this operator is created only on a 64bit platform***
+ */
+ uint FromInt(signed int i)
+ {
+ return FromInt(sint(i));
+ }
+
+
+ /*!
+ this method converts 32 bit unsigned int type to this class
+ ***this operator is created only on a 64bit platform***
+ */
+ Int & operator=(unsigned int i)
+ {
+ FromUInt(i);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting 32 bit unsigned int to this class
+ ***this constructor is created only on a 64bit platform***
+ */
+/* Int(unsigned int i)
+ {
+ FromUInt(i);
+ }
+*/
+
+ /*!
+ this operator converts 32 bit signed int type to this class
+ ***this operator is created only on a 64bit platform***
+ */
+ Int & operator=(signed int i)
+ {
+ FromInt(i);
+
+ return *this;
+ }
+
+
+ /*!
+ a constructor for converting 32 bit signed int to this class
+ ***this constructor is created only on a 64bit platform***
+ */
+/* Int(signed int i)
+ {
+ FromInt(i);
+ }
+*/
+#endif
+
+
+
+ /*!
+ a constructor for converting string to this class (with the base=10)
+ */
+/* Int(const char * s)
+ {
+ FromString(s);
+ }
+*/
+
+ /*!
+ a constructor for converting a string to this class (with the base=10)
+ */
+/* Int(const std::string & s)
+ {
+ FromString( s.c_str() );
+ }
+*/
+
+#ifndef TTMATH_DONT_USE_WCHAR
+
+ /*!
+ a constructor for converting string to this class (with the base=10)
+ */
+ Int(const wchar_t * s)
+ {
+ FromString(s);
+ }
+
+
+ /*!
+ a constructor for converting a string to this class (with the base=10)
+ */
+ Int(const std::wstring & s)
+ {
+ FromString( s.c_str() );
+ }
+
+#endif
+
+
+ /*!
+ a default constructor
+
+ we don't clear table etc.
+ */
+/* Int()
+ {
+ }
+*/
+
+ /*!
+ the destructor
+ */
+/* ~Int()
+ {
+ }
+*/
+
+ /*!
+ this method returns the lowest value from table with a sign
+
+ we must be sure when we using this method whether the value
+ will hold in an sint type or not (the rest value from table must be zero or -1)
+ */
+ sint ToInt() const
+ {
+ return sint( UInt::table[0] );
+ }
+
+
+ /*!
+ this method converts the value to uint type
+ can return a carry if the value is too long to store it in uint type
+ */
+ uint ToUInt(uint & result) const
+ {
+ uint c = UInt::ToUInt(result);
+
+ if( value_size == 1 )
+ return (result & TTMATH_UINT_HIGHEST_BIT) == 0 ? 0 : 1;
+
+ return c;
+ }
+
+
+ /*!
+ this method converts the value to uint type
+ can return a carry if the value is too long to store it in uint type
+ */
+ uint ToInt(uint & result) const
+ {
+ return ToUInt(result);
+ }
+
+
+ /*!
+ this method converts the value to sint type
+ can return a carry if the value is too long to store it in sint type
+ */
+ uint ToInt(sint & result) const
+ {
+ result = sint( UInt::table[0] );
+ uint mask = IsSign() ? TTMATH_UINT_MAX_VALUE : 0;
+
+ if( (result & TTMATH_UINT_HIGHEST_BIT) != (mask & TTMATH_UINT_HIGHEST_BIT) )
+ return 1;
+
+ for(uint i=1 ; i::table[i] != mask )
+ return 1;
+
+ return 0;
+ }
+
+
+#ifdef TTMATH_PLATFORM32
+
+ /*!
+ this method converts the value to ulint type (64 bit unsigned integer)
+ can return a carry if the value is too long to store it in ulint type
+ *** this method is created only on a 32 bit platform ***
+ */
+ uint ToUInt(ulint & result) const
+ {
+ uint c = UInt::ToUInt(result);
+
+ if( value_size == 1 )
+ return (UInt::table[0] & TTMATH_UINT_HIGHEST_BIT) == 0 ? 0 : 1;
+
+ if( value_size == 2 )
+ return (UInt::table[1] & TTMATH_UINT_HIGHEST_BIT) == 0 ? 0 : 1;
+
+ return c;
+ }
+
+
+ /*!
+ this method converts the value to ulint type (64 bit unsigned integer)
+ can return a carry if the value is too long to store it in ulint type
+ *** this method is created only on a 32 bit platform ***
+ */
+ uint ToInt(ulint & result) const
+ {
+ return ToUInt(result);
+ }
+
+
+ /*!
+ this method converts the value to slint type (64 bit signed integer)
+ can return a carry if the value is too long to store it in slint type
+ *** this method is created only on a 32 bit platform ***
+ */
+ uint ToInt(slint & result) const
+ {
+ if( value_size == 1 )
+ {
+ result = slint(sint(UInt::table[0]));
+ }
+ else
+ {
+ uint low = UInt::table[0];
+ uint high = UInt::table[1];
+
+ result = low;
+ result |= (ulint(high) << TTMATH_BITS_PER_UINT);
+
+ uint mask = IsSign() ? TTMATH_UINT_MAX_VALUE : 0;
+
+ if( (high & TTMATH_UINT_HIGHEST_BIT) != (mask & TTMATH_UINT_HIGHEST_BIT) )
+ return 1;
+
+ for(uint i=2 ; i::table[i] != mask )
+ return 1;
+ }
+
+ return 0;
+ }
+
+#endif
+
+
+
+#ifdef TTMATH_PLATFORM64
+
+ /*!
+ this method converts the value to a 32 bit unsigned integer
+ can return a carry if the value is too long to store it in this type
+ *** this method is created only on a 64 bit platform ***
+ */
+ uint ToUInt(unsigned int & result) const
+ {
+ uint c = UInt::ToUInt(result);
+
+ if( c || IsSign() )
+ return 1;
+
+ return 0;
+ }
+
+
+ /*!
+ this method converts the value to a 32 bit unsigned integer
+ can return a carry if the value is too long to store it in this type
+ *** this method is created only on a 64 bit platform ***
+ */
+ uint ToInt(unsigned int & result) const
+ {
+ return ToUInt(result);
+ }
+
+
+ /*!
+ this method converts the value to a 32 bit signed integer
+ can return a carry if the value is too long to store it in this type
+ *** this method is created only on a 64 bit platform ***
+ */
+ uint ToInt(int & result) const
+ {
+ uint first = UInt::table[0];
+
+ result = int(first);
+ uint mask = IsSign() ? TTMATH_UINT_MAX_VALUE : 0;
+
+ if( (first >> 31) != (mask >> 31) )
+ return 1;
+
+ for(uint i=1 ; i::table[i] != mask )
+ return 1;
+
+ return 0;
+ }
+
+#endif
+
+
+
+ /*!
+ an auxiliary method for converting to a string
+ */
+ template
+ void ToStringBase(string_type & result, uint b = 10) const
+ {
+ if( IsSign() )
+ {
+ Int temp(*this);
+ temp.Abs();
+ temp.UInt::ToStringBase(result, b, true);
+ }
+ else
+ {
+ UInt::ToStringBase(result, b, false);
+ }
+ }
+
+ /*!
+ this method converts the value to a string with a base equal 'b'
+ */
+ void ToString(std::string & result, uint b = 10) const
+ {
+ return ToStringBase(result, b);
+ }
+
+
+ /*!
+ this method converts the value to a string with a base equal 'b'
+ */
+ std::string ToString(uint b = 10) const
+ {
+ std::string result;
+ ToStringBase(result, b);
+
+ return result;
+ }
+
+
+#ifndef TTMATH_DONT_USE_WCHAR
+
+ /*!
+ this method converts the value to a string with a base equal 'b'
+ */
+ void ToString(std::wstring & result, uint b = 10) const
+ {
+ return ToStringBase(result, b);
+ }
+
+
+ /*!
+ this method converts the value to a string with a base equal 'b'
+ */
+ std::wstring ToWString(uint b = 10) const
+ {
+ std::wstring result;
+ ToStringBase(result, b);
+
+ return result;
+ }
+
+#endif
+
+
+
+private:
+
+ /*!
+ an auxiliary method for converting from a string
+ */
+ template
+ uint FromStringBase(const char_type * s, uint b = 10, const char_type ** after_source = 0, bool * value_read = 0)
+ {
+ bool is_sign = false;
+
+ Misc::SkipWhiteCharacters(s);
+
+ if( *s == '-' )
+ {
+ is_sign = true;
+ Misc::SkipWhiteCharacters(++s);
+ }
+ else
+ if( *s == '+' )
+ {
+ Misc::SkipWhiteCharacters(++s);
+ }
+
+ if( UInt::FromString(s,b,after_source,value_read) )
+ return 1;
+
+ if( is_sign )
+ {
+ Int mmin;
+
+ mmin.SetMin();
+
+ /*
+ the reference to mmin will be automatically converted to the reference
+ to UInt type
+ (this value can be equal mmin -- look at a description in ChangeSign())
+ */
+ if( UInt::operator>( mmin ) )
+ return 1;
+
+ /*
+ if the value is equal mmin the method ChangeSign() does nothing (only returns 1 but we ignore it)
+ */
+ ChangeSign();
+ }
+ else
+ {
+ Int mmax;
+
+ mmax.SetMax();
+
+ if( UInt::operator>( mmax ) )
+ return 1;
+ }
+
+ return 0;
+ }
+
+
+public:
+
+ /*!
+ this method converts a string into its value
+ it returns carry=1 if the value will be too big or an incorrect base 'b' is given
+
+ string is ended with a non-digit value, for example:
+ "-12" will be translated to -12
+ as well as:
+ "- 12foo" will be translated to -12 too
+
+ existing first white characters will be ommited
+ (between '-' and a first digit can be white characters too)
+
+ after_source (if exists) is pointing at the end of the parsed string
+
+ value_read (if exists) tells whether something has actually been read (at least one digit)
+ */
+ uint FromString(const char * s, uint b = 10, const char ** after_source = 0, bool * value_read = 0)
+ {
+ return FromStringBase(s, b, after_source, value_read);
+ }
+
+
+ /*!
+ this method converts a string into its value
+ */
+ uint FromString(const wchar_t * s, uint b = 10, const wchar_t ** after_source = 0, bool * value_read = 0)
+ {
+ return FromStringBase(s, b, after_source, value_read);
+ }
+
+
+ /*!
+ this method converts a string into its value
+ it returns carry=1 if the value will be too big or an incorrect base 'b' is given
+ */
+ uint FromString(const std::string & s, uint b = 10)
+ {
+ return FromString( s.c_str(), b );
+ }
+
+
+ /*!
+ this operator converts a string into its value (with base = 10)
+ */
+ Int & operator=(const char * s)
+ {
+ FromString(s);
+
+ return *this;
+ }
+
+
+#ifndef TTMATH_DONT_USE_WCHAR
+
+
+ /*!
+ this method converts a string into its value
+ it returns carry=1 if the value will be too big or an incorrect base 'b' is given
+ */
+ uint FromString(const std::wstring & s, uint b = 10)
+ {
+ return FromString( s.c_str(), b );
+ }
+
+
+ /*!
+ this operator converts a string into its value (with base = 10)
+ */
+ Int & operator=(const wchar_t * s)
+ {
+ FromString(s);
+
+ return *this;
+ }
+
+
+ /*!
+ this operator converts a string into its value (with base = 10)
+ */
+ Int & operator=(const std::wstring & s)
+ {
+ FromString( s.c_str() );
+
+ return *this;
+ }
+
+#endif
+
+
+ /*!
+ this operator converts a string into its value (with base = 10)
+ */
+ Int & operator=(const std::string & s)
+ {
+ FromString( s.c_str() );
+
+ return *this;
+ }
+
+
+
+ /*!
+ *
+ * methods for comparing
+ *
+ *
+ */
+
+ bool operator==(const Int & l) const
+ {
+ return UInt::operator==(l);
+ }
+
+ bool operator!=(const Int & l) const
+ {
+ return UInt::operator!=(l);
+ }
+
+ bool operator<(const Int & l) const
+ {
+ sint i=value_size-1;
+
+ sint a1 = sint(UInt::table[i]);
+ sint a2 = sint(l.table[i]);
+
+ if( a1 != a2 )
+ return a1 < a2;
+
+
+ for(--i ; i>=0 ; --i)
+ {
+ if( UInt::table[i] != l.table[i] )
+ // comparison as unsigned int
+ return UInt::table[i] < l.table[i];
+ }
+
+ // they're equal
+ return false;
+ }
+
+
+ bool operator>(const Int & l) const
+ {
+ sint i=value_size-1;
+
+ sint a1 = sint(UInt::table[i]);
+ sint a2 = sint(l.table[i]);
+
+ if( a1 != a2 )
+ return a1 > a2;
+
+
+ for(--i ; i>=0 ; --i)
+ {
+ if( UInt::table[i] != l.table[i] )
+ // comparison as unsigned int
+ return UInt::table[i] > l.table[i];
+ }
+
+ // they're equal
+ return false;
+ }
+
+
+ bool operator<=(const Int & l) const
+ {
+ sint i=value_size-1;
+
+ sint a1 = sint(UInt::table[i]);
+ sint a2 = sint(l.table[i]);
+
+ if( a1 != a2 )
+ return a1 < a2;
+
+
+ for(--i ; i>=0 ; --i)
+ {
+ if( UInt::table[i] != l.table[i] )
+ // comparison as unsigned int
+ return UInt::table[i] < l.table[i];
+ }
+
+ // they're equal
+ return true;
+ }
+
+
+ bool operator>=(const Int & l) const
+ {
+ sint i=value_size-1;
+
+ sint a1 = sint(UInt::table[i]);
+ sint a2 = sint(l.table[i]);
+
+ if( a1 != a2 )
+ return a1 > a2;
+
+
+ for(--i ; i>=0 ; --i)
+ {
+ if( UInt::table[i] != l.table[i] )
+ // comparison as unsigned int
+ return UInt::table[i] > l.table[i];
+ }
+
+ // they're equal
+ return true;
+ }
+
+
+
+ /*!
+ *
+ * standard mathematical operators
+ *
+ */
+
+
+ /*!
+ an operator for changing the sign
+
+ it's not changing 'this' but the changed value will be returned
+ */
+ Int operator-() const
+ {
+ Int temp(*this);
+
+ temp.ChangeSign();
+
+ return temp;
+ }
+
+
+ Int operator-(const Int & p2) const
+ {
+ Int temp(*this);
+
+ temp.Sub(p2);
+
+ return temp;
+ }
+
+
+ Int & operator-=(const Int