Implemented precision-of-destination for Float calculations where possible

This commit is contained in:
Alex Dyachenko 2015-12-24 10:57:07 -05:00
parent d895d98428
commit 9870b4d812
9 changed files with 302 additions and 23 deletions

View File

@ -118,6 +118,9 @@
<Compile Include="..\..\..\mpir.net\mpir.net-tests\HugeFloatTests\Math.cs">
<Link>HugeFloatTests\Math.cs</Link>
</Compile>
<Compile Include="..\..\..\mpir.net\mpir.net-tests\HugeFloatTests\Precision.cs">
<Link>HugeFloatTests\Precision.cs</Link>
</Compile>
<Compile Include="..\..\..\mpir.net\mpir.net-tests\HugeIntTests\Arithmetic.cs">
<Link>HugeIntTests\Arithmetic.cs</Link>
</Compile>

View File

@ -0,0 +1,244 @@
/*
Copyright 2014 Alex Dyachenko
This file is part of the MPIR Library.
The MPIR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or (at
your option) any later version.
The MPIR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the MPIR Library. If not, see http://www.gnu.org/licenses/.
*/
using System;
using Microsoft.VisualStudio.TestTools.UnitTesting;
namespace MPIR.Tests.HugeFloatTests
{
/// <summary>
/// tests in this class verify that the correct precision is used when calculating floating point numbers
/// </summary>
[TestClass]
public class Precision
{
[ClassInitialize]
public static void Setup(TestContext context)
{
HugeFloat.DefaultPrecision = 128;
}
[ClassCleanup]
public static void Cleanup()
{
HugeFloat.DefaultPrecision = 64;
}
#region Expression arithmetic
[TestMethod]
public void ExpressionsCalculatedToDestinationPrecision()
{
using (var a = new HugeFloat(1))
using (var b = new HugeFloat(3))
using (var c = new HugeFloat(5))
using (var d = new HugeFloat(a / b + c))
{
Assert.AreEqual("0.5333333333333333333333333333333333333332@1", d.ToString());
d.Reallocate(256);
d.Value = a / b + c;
Assert.AreEqual("0.533333333333333333333333333333333333333333333333333333333333333333333333333333@1", d.ToString());
}
}
[TestMethod]
public void ExpressionsCalculatedToSpecificPrecisionForEquals()
{
using (var a = new HugeFloat(1))
using (var b = new HugeFloat(3))
using (var c = new HugeFloat("12345234589234059823475029384572323452034958723049823408955"))
using (var d = HugeFloat.Allocate(256))
{
Assert.IsTrue(c.Equals(c + a / b, 128));
Assert.IsFalse(c.Equals(c + a / b, 256));
d.SetTo("12345234589234059823475029384572323452034958723049823408955.333333333333333333333333333333333");
Assert.IsTrue(d.Equals(c + a / b, 256));
Assert.IsTrue(d.Equals(c + a / b, 128));
}
}
[TestMethod]
public void ExpressionHashCodeCalculatedToDefaultPrecision()
{
using (var a = new HugeFloat(1))
using (var b = new HugeFloat(13))
using (var c = new HugeFloat("12345234589234059823475029384572323"))
{
c.Value *= 0x4000000000000000L;
var cHash = c.GetHashCode();
var expr = a / b + c;
Assert.AreEqual(cHash, expr.GetHashCode());
HugeFloat.DefaultPrecision = 256;
Assert.AreEqual(cHash, c.GetHashCode());
Assert.AreNotEqual(cHash, expr.GetHashCode());
HugeFloat.DefaultPrecision = 128;
}
}
[TestMethod]
public void CompareToCalculatedToDefaultPrecision()
{
using (var a = new HugeFloat(1))
using (var b = new HugeFloat(13))
using (var c = new HugeFloat("12345234589234059823475029384572323"))
using (var d = HugeFloat.Allocate(256))
{
c.Value *= 0x4000000000000000L;
d.Value = c;
var expr = a / b + c;
Assert.AreEqual(0, c.CompareTo(expr)); //to precision of c
Assert.AreEqual(0, expr.CompareTo(c)); //to precision of c
Assert.IsFalse(expr > c); //to precision of c
Assert.IsTrue(c == expr); //to precision of c
Assert.AreEqual(0, (c + 0).CompareTo(expr)); //to default precision
Assert.AreEqual(0, expr.CompareTo(c + 0)); //to default precision
Assert.IsFalse(expr > c + 0); //to default precision
Assert.IsTrue(c + 0 == expr); //to default precision
HugeFloat.DefaultPrecision = 256;
Assert.AreEqual(0, c.CompareTo(expr)); //to precision of c
Assert.AreEqual(0, expr.CompareTo(c)); //to precision of c
Assert.IsTrue(c == expr); //to precision of c
Assert.IsFalse(expr > c); //to precision of c
Assert.AreEqual(-1, d.CompareTo(expr)); //to precision of d
Assert.AreEqual(1, expr.CompareTo(d)); //to precision of d
Assert.IsFalse(d == expr); //to precision of d
Assert.IsTrue(expr > d); //to precision of d
Assert.AreEqual(-1, (c * 1).CompareTo(expr)); //to default precision
Assert.AreEqual(1, expr.CompareTo(c + 0)); //to default precision
Assert.IsFalse(c + 0 == expr); //to default precision
Assert.IsTrue(expr > c + 0); //to default precision
HugeFloat.DefaultPrecision = 128;
}
}
[TestMethod]
public void CompareToPrimitiveCalculatedToDefaultPrecision()
{
using (var a = new HugeFloat(1))
using (var b = new HugeFloat(13))
using (var c = new HugeFloat("12345234589234059823475029384572323"))
using (var d = HugeFloat.Allocate(256))
{
c.Value *= 0x4000000000000000L;
d.Value = c;
var expr = a / b + c - c;
Assert.AreEqual(0, Math.Sign(expr.CompareTo(0L)));
Assert.AreEqual(0, Math.Sign(expr.CompareTo(0UL)));
Assert.AreEqual(0, Math.Sign(expr.CompareTo(0.0)));
HugeFloat.DefaultPrecision = 256;
Assert.AreEqual(1, Math.Sign(expr.CompareTo(0L)));
Assert.AreEqual(1, Math.Sign(expr.CompareTo(0UL)));
Assert.AreEqual(1, Math.Sign(expr.CompareTo(0.0)));
HugeFloat.DefaultPrecision = 128;
}
}
[TestMethod]
public void EqualsToPrimitiveCalculatedToDefaultPrecision()
{
using (var a = new HugeFloat(1))
using (var b = new HugeFloat(13))
using (var c = new HugeFloat("12345234589234059823475029384572323"))
using (var d = HugeFloat.Allocate(256))
{
c.Value *= 0x4000000000000000L;
d.Value = c;
var expr = a / b + c - c;
Assert.IsTrue(expr.Equals(0L));
Assert.IsTrue(expr.Equals(0UL));
Assert.IsTrue(expr.Equals(0.0));
HugeFloat.DefaultPrecision = 256;
Assert.IsFalse(expr.Equals(0L));
Assert.IsFalse(expr.Equals(0UL));
Assert.IsFalse(expr.Equals(0.0));
HugeFloat.DefaultPrecision = 128;
}
}
[TestMethod]
public void SignCalculatedToDefaultPrecision()
{
using (var a = new HugeFloat(1))
using (var b = new HugeFloat(13))
using (var c = new HugeFloat("12345234589234059823475029384572323"))
using (var d = HugeFloat.Allocate(256))
{
c.Value *= 0x4000000000000000L;
var expr = (a / b + c) - c;
d.Value = expr;
Assert.AreEqual(0, expr.Sign());
Assert.AreEqual(1, d.Sign());
HugeFloat.DefaultPrecision = 256;
Assert.AreEqual(1, expr.Sign());
Assert.AreEqual(1, d.Sign());
d.Precision = 128;
Assert.AreEqual(1, d.Sign());
d.Value = expr;
Assert.AreEqual(0, d.Sign());
HugeFloat.DefaultPrecision = 128;
}
}
[TestMethod]
public void HugeIntSetToPerformedToDefaultPrecision()
{
using (var a = new HugeFloat(14))
using (var b = new HugeFloat(13))
using (var c = new HugeFloat("1234523458923405982347445029384572323"))
using (var d = new HugeInt())
{
c.Value *= 0x4000000000000000L;
c.Value *= 0x4000000000000000L;
var expr = a / b + c - c;
d.SetTo(expr);
Assert.IsTrue(d == 0);
HugeFloat.DefaultPrecision = 256;
d.SetTo(expr);
Assert.IsTrue(d == 1);
HugeFloat.DefaultPrecision = 128;
}
}
[TestMethod]
public void HugeRationalSetToPerformedToDefaultPrecision()
{
using (var a = new HugeFloat(14))
using (var b = new HugeFloat(13))
using (var c = new HugeFloat("1234523458923405982347445029384572323"))
using (var d = new HugeRational())
{
c.Value *= 0x4000000000000000L;
c.Value *= 0x4000000000000000L;
var expr = a / b + c - c;
d.SetTo(expr);
Assert.IsTrue(d == 0);
HugeFloat.DefaultPrecision = 256;
d.SetTo(expr);
Assert.IsTrue(d > 1);
Assert.IsTrue(d < 2);
HugeFloat.DefaultPrecision = 128;
}
}
#endregion
}
}

View File

@ -136,6 +136,7 @@ struct EvaluationContext
};
__int64 Zero;
};
mp_bitcnt_t FloatPrecision;
void inline Initialized(EvaluationOptions flag)
{
@ -145,6 +146,7 @@ struct EvaluationContext
EvaluationContext()
{
Zero = 0;
FloatPrecision = 0;
}
#define CTXT_ADD_RATIONAL_SI(numerator, denominator) \

View File

@ -20,16 +20,13 @@ along with the MPIR Library. If not, see http://www.gnu.org/licenses/.
#pragma region Expression macros
#define IN_CONTEXT_1(a) \
EvaluationContext context; \
a->ASSIGN_TO(context)
#define IN_CONTEXT_2(a, b) \
EvaluationContext context; \
a->ASSIGN_TO(context); \
b->ASSIGN_TO(context)
#define IN_CONTEXT_3(a, b, c) \
EvaluationContext context; \
a->ASSIGN_TO(context); \
b->ASSIGN_TO(context); \
c->ASSIGN_TO(context)
@ -41,7 +38,21 @@ along with the MPIR Library. If not, see http://www.gnu.org/licenses/.
#define MACRO_CHOOSE1(prefix, number) MACRO_CHOOSE2(prefix, number)
#define MACRO_CHOOSE(prefix, number) MACRO_CHOOSE1(prefix, number)
#define MACRO_GLUE(x, y) x y
#define IN_CONTEXT(...) MACRO_GLUE(MACRO_CHOOSE(IN_CONTEXT_, COUNT_ARGS(__VA_ARGS__)), (__VA_ARGS__))
#define IN_CONTEXT(...) \
EvaluationContext context; \
SET_CONTEXT_PRECISION \
MACRO_GLUE(MACRO_CHOOSE(IN_CONTEXT_, COUNT_ARGS(__VA_ARGS__)), (__VA_ARGS__))
#define IN_DEFAULT_CONTEXT(...) \
EvaluationContext context; \
context.FloatPrecision = HugeFloat::DefaultPrecision; \
MACRO_GLUE(MACRO_CHOOSE(IN_CONTEXT_, COUNT_ARGS(__VA_ARGS__)), (__VA_ARGS__))
#define IN_SPECIFIC_CONTEXT(precision, ...) \
EvaluationContext context; \
context.FloatPrecision = precision; \
MACRO_GLUE(MACRO_CHOOSE(IN_CONTEXT_, COUNT_ARGS(__VA_ARGS__)), (__VA_ARGS__))
//defines a unary expression class
#define DEFINE_UNARY_EXPRESSION(base, name, type) \

View File

@ -180,7 +180,7 @@ namespace MPIR
int MPEXPR_NAME::GetHashCode()
{
IN_CONTEXT(this);
IN_DEFAULT_CONTEXT(this);
mp_limb_t hash = CTXT(0)->_mp_exp;
mp_limb_t* ptr = CTXT(0)->_mp_d;
@ -208,23 +208,24 @@ namespace MPIR
if(!IS_NULL(expr))
return CompareTo(expr);
EvaluationContext context;
auto f = dynamic_cast<MPTYPE^>(this);
auto precision = IS_NULL(f) ? MPTYPE::DefaultPrecision : f->Precision;
if(a->GetType() == mpir_ui::typeid)
{
ASSIGN_TO(context);
IN_SPECIFIC_CONTEXT(precision, this);
return MP(cmp_ui)(CTXT(0), (mpir_ui)a);
}
if(a->GetType() == mpir_si::typeid)
{
ASSIGN_TO(context);
IN_SPECIFIC_CONTEXT(precision, this);
return MP(cmp_si)(CTXT(0), (mpir_si)a);
}
if(a->GetType() == double::typeid)
{
ASSIGN_TO(context);
IN_SPECIFIC_CONTEXT(precision, this);
return MP(cmp_d)(CTXT(0), (double)a);
}
@ -248,7 +249,11 @@ namespace MPIR
if (IS_NULL(a))
return 1;
IN_CONTEXT(this, a);
auto f = dynamic_cast<MPTYPE^>(this);
if (IS_NULL(f)) f = dynamic_cast<MPTYPE^>(a);
auto precision = IS_NULL(f) ? MPTYPE::DefaultPrecision : f->Precision;
IN_SPECIFIC_CONTEXT(precision, this, a);
return MP(cmp)(CTXT(0), CTXT(1));
}
@ -340,6 +345,12 @@ namespace MPIR
DEFINE_BINARY_ASSIGNMENT_REF_VAL(Power, Flt, Ui, MP(pow_ui))
DEFINE_BINARY_ASSIGNMENT_REF_REF(RelativeDifferenceFrom, Flt, MP(reldiff))
int MPEXPR_NAME::Sign()
{
IN_DEFAULT_CONTEXT(this);
return MP(sgn)(CTXT(0));
}
#pragma endregion
#pragma region IO
@ -437,13 +448,13 @@ namespace MPIR
void HugeInt::SetTo(MPEXPR_NAME^ value)
{
IN_CONTEXT(value);
IN_DEFAULT_CONTEXT(value);
mpz_set_f(_value, CTXT(0));
}
void HugeRational::SetTo(MPEXPR_NAME^ value)
{
IN_CONTEXT(value);
IN_DEFAULT_CONTEXT(value);
mpq_set_f(_value, CTXT(0));
}

View File

@ -36,6 +36,7 @@ using namespace System::Runtime::InteropServices;
#undef CTXTI
#undef ASSIGN_TO
#undef Mpt
#undef SET_CONTEXT_PRECISION
#endif
#define SPECIALIZE_EXPRESSIONS
#define Mpt Flt
@ -50,6 +51,7 @@ using namespace System::Runtime::InteropServices;
#define CTXTI(x) context.IntArgs[x]
#define CTXTR(x) context.RationalArgs[x]
#define ASSIGN_TO CONCAT(AssignTo, LIT(MPTYPE_NAME))
#define SET_CONTEXT_PRECISION context.FloatPrecision = mpf_get_prec(destination);
#include "ExpressionMacros.h"
namespace MPIR
@ -86,7 +88,7 @@ namespace MPIR
context.Initialized(FloatInitialized);
auto ptr = &context.Temp[context.Index].MPTYPE_NAME;
CTXT(context.Index++) = ptr;
MP(init)(ptr);
MP(init2)(ptr, context.FloatPrecision);
AssignTo(ptr);
}
@ -333,14 +335,17 @@ namespace MPIR
/// <summary>Compares two numbers.
/// <para>If any argument is an expression, it is evaluated into a temporary variable before the comparison is performed.
/// </para>When the argument is a double, it may be an infinity, but results are undefined for a NaN.</summary>
/// </para>The precision of the calculation is the precision of this instance if it is a computed number, otherwise the precision of <paramref name="a"/> if that is a computed number,
/// otherwise the current default float precision.
/// <para>When the argument is a double, it may be an infinity, but results are undefined for a NaN.</para></summary>
/// <param name="a">Value to compare the source with</param>
/// <returns>A positive number if the source is greater than <paramref name="a"/>, negative if less, and zero if they are equal.</returns>
virtual int CompareTo(Object^ a) sealed;
/// <summary>Compares two numbers.
/// <para>If any argument is an expression, it is evaluated into a temporary variable before the comparison is performed.
/// </para></summary>
/// </para>The precision of the calculation is the precision of this instance if it is a computed number, otherwise the precision of <paramref name="a"/> if that is a computed number,
/// otherwise the current default float precision if both are expressions.</summary>
/// <param name="a">Value to compare the source with</param>
/// <returns>A positive number if the source is greater than <paramref name="a"/>, negative if less, and zero if they are equal.</returns>
virtual int CompareTo(MPEXPR_NAME^ a) sealed;
@ -360,16 +365,16 @@ namespace MPIR
virtual bool Equals(Object^ a) override sealed;
/// <summary>Compares two numbers approximately, taking into account <paramref name="precision"/> most significant bits of the mantissa.
/// <para>If any argument is an expression, it is evaluated into a temporary variable before the comparison is performed.
/// <para>If any argument is an expression, it is evaluated into a temporary variable with the specified <paramref name="precision"/> before the comparison is performed.
/// </para>In the future values like 1000 and 0111 may be considered the same to 3 bits (on the basis that their difference is that small).
/// </summary>
/// <param name="a">Value to compare the source with</param>
/// <param name="precision">The number of most significant bits that must match for the two numbers to be considered equal</param>
/// <returns>true if the values of the source and <paramref name="a"/> are equal to <paramref name="precision"/>, false otherwise.</returns>
bool Equals(MPEXPR_NAME^ a, mp_bitcnt_t precision) { IN_CONTEXT(this, a); return MP(eq)(CTXT(0), CTXT(1), precision) != 0; }
bool Equals(MPEXPR_NAME^ a, mp_bitcnt_t precision) { IN_SPECIFIC_CONTEXT(precision, this, a); return MP(eq)(CTXT(0), CTXT(1), precision) != 0; }
/// <summary>Computes the hash code of the source value.
/// <para>If called on an expression, it is evaluated into a temporary variable before the comparison is performed.
/// <para>If called on an expression, it is evaluated into a temporary variable with the current default float precision before the calculation is performed.
/// </para>Multi-precision classes are mutable with value semantics. The hash code is based on the value, and will change if the value changes.
/// For this reason, the value of an object must not be modified while the object is contained in a hash table.</summary>
/// <returns>a signed integer hash code for the value.</returns>
@ -712,10 +717,10 @@ namespace MPIR
static bool operator == (double b, MPEXPR_NAME^ a) { return !IS_NULL(a) && a->CompareTo(b) == 0; }
/// <summary>Calculates the sign (+1, 0, or -1) of the source value.
/// <para>If the source is an expression, it is evaluated into a temporary variable before the sign is computed.
/// <para>If the source is an expression, it is evaluated into a temporary variable with the current default float precision before the sign is computed.
/// </para></summary>
/// <returns>+1 if the source is positive, -1 if negative, and 0 if zero.</returns>
int Sign() { IN_CONTEXT(this); return MP(sgn)(CTXT(0)); }
int Sign();
/// <summary>Compares two numbers.
/// <para>If any argument is an expression, it is evaluated into a temporary variable before the comparison is performed.

View File

@ -35,6 +35,7 @@ using namespace System::Runtime::InteropServices;
#undef CTXT
#undef ASSIGN_TO
#undef Mpt
#undef SET_CONTEXT_PRECISION
#endif
#define SPECIALIZE_EXPRESSIONS
#define Mpt Int
@ -47,6 +48,7 @@ using namespace System::Runtime::InteropServices;
#define MPEXPR(x) LIT(MPTYPE_NAME)##x##Expression
#define CTXT(x) context.IntArgs[x]
#define ASSIGN_TO CONCAT(AssignTo, LIT(MPTYPE_NAME))
#define SET_CONTEXT_PRECISION
#include "ExpressionMacros.h"
extern __mpz_struct HugeIntConst1;
@ -1665,7 +1667,7 @@ namespace MPIR
/// <summary>
/// Sets the value of the integer object. Any fractional portion is truncated.
/// <para>Do not change the value of an object while it is contained in a hash table, because that changes its hash code.
/// </para></summary>
/// </para>If the argument is an expression, it is evaluated with the current default float precision.</summary>
/// <param name="value">new value for the object</param>
void SetTo(FloatExpression^ value);

View File

@ -36,6 +36,7 @@ using namespace System::Runtime::InteropServices;
#undef CTXTI
#undef ASSIGN_TO
#undef Mpt
#undef SET_CONTEXT_PRECISION
#endif
#define SPECIALIZE_EXPRESSIONS
#define Mpt Rat
@ -49,6 +50,7 @@ using namespace System::Runtime::InteropServices;
#define CTXT(x) context.RationalArgs[x]
#define CTXTI(x) context.IntArgs[x]
#define ASSIGN_TO CONCAT(AssignTo, LIT(MPTYPE_NAME))
#define SET_CONTEXT_PRECISION
#include "ExpressionMacros.h"
namespace MPIR
@ -1177,7 +1179,7 @@ namespace MPIR
/// <summary>
/// Sets the value of the rational object. There is no rounding, this conversion is exact.
/// <para>Do not change the value of an object while it is contained in a hash table, because that changes its hash code.
/// </para></summary>
/// </para>If the argument is an expression, it is evaluated with the current default float precision.</summary>
/// <param name="value">new value for the object</param>
void SetTo(FloatExpression^ value);

View File

@ -45,7 +45,6 @@ namespace MPIR
MP(rrandomb)(destination, Left->_value, BITS_TO_LIMBS(MP(get_prec)(destination)), Right);
}
//TODO implement "precision of destination" for context ops
//TODO investigate implementing raw IO for floats
DEFINE_ASSIGNMENT_PROLOG(RandomLimbsChunky)