UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-03-18 21:30:38 +01:00
Code Issues Packages Projects Releases Wiki Activity

[poincare] Allow more BinomialCoefficient exact results

Browse Source 
Change-Id: I9dbbf9471ee6d9e12fe4861a5b11990858382562
This commit is contained in:
Hugo Saint-Vignes
2020-09-25 11:52:23 +02:00
committed by Émilie Feral
parent 857ca4eecd
commit a89878a24d
3 changed files with 33 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
poincare/src/binomial_coefficient.cpp
View File

@@ -87,32 +87,40 @@ Expression BinomialCoefficient::shallowReduce(Context * context) {
}
if (!r0.isInteger()) {
// Generalized binomial coefficient
// Generalized binomial coefficient (n is not an integer)
return *this;
}
Integer n = r0.signedIntegerNumerator();
Integer k = r1.signedIntegerNumerator();
/* Check for situations where there should be no reduction in order to avoid
* too long computation and a huge result. The binomial coefficient will be
* approximatively evaluated later. */
if (n.isLowerThan(k)) {
// Generalized binomial coefficient
return *this;
}
/* If n is too big, we do not reduce in order to avoid too long computation.
* The binomial coefficient will be approximatively evaluated later. */
if (Integer(k_maxNValue).isLowerThan(n)) {
// Generalized binomial coefficient (n < k)
if (!n.isNegative()) {
// When n is an integer and 0 <= n < k, binomial(n,k) is 0.
return Rational::Builder(0);
}
if (Integer(k_maxNValue).isLowerThan(Integer::Subtraction(k, n))) {
return *this;
}
} else if (Integer(k_maxNValue).isLowerThan(n)) {
return *this;
}
Rational result = Rational::Builder(1);
Integer kBis = Integer::Subtraction(n, k);
k = kBis.isLowerThan(k) ? kBis : k;
int clippedK = k.extractedInt(); // Authorized because k < n < k_maxNValue
// Take advantage of symmetry if n >= k
k = !n.isLowerThan(k) && kBis.isLowerThan(k) ? kBis : k;
int clippedK = k.extractedInt(); // Authorized because k < k_maxNValue
for (int i = 0; i < clippedK; i++) {
Integer nMinusI = Integer::Subtraction(n, Integer(i));
Integer kMinusI = Integer::Subtraction(k, Integer(i));
Rational factor = Rational::Builder(nMinusI, kMinusI);
result = Rational::Multiplication(result, factor);
}
// As we cap the n < k_maxNValue = 300, result < binomial(300, 150) ~2^89
// As we cap the n < k_maxNValue = 300, result < binomial(300, 150) ~10^89
// If n was negative, k - n < k_maxNValue, result < binomial(-150,150) ~10^88
assert(!result.numeratorOrDenominatorIsInfinity());
replaceWithInPlace(result);
return std::move(result);

19
poincare/test/approximation.cpp
View File

@@ -249,6 +249,16 @@ QUIZ_CASE(poincare_approximation_function) {
assert_expression_approximates_to<float>("binomial(10, 4)", "210");
assert_expression_approximates_to<double>("binomial(10, 4)", "210");
assert_expression_approximates_to<float>("binomial(12, 3)", "220");
assert_expression_approximates_to<double>("binomial(12, 3)", "220");
assert_expression_approximates_to<float>("binomial(-4.6, 3)", "-28.336");
assert_expression_approximates_to<double>("binomial(-4.6, 3)", "-28.336");
assert_expression_approximates_to<float>("binomial(π, 3)", "1.280108");
assert_expression_approximates_to<double>("binomial(π, 3)", "1.2801081307019");
assert_expression_approximates_to<float>("binomial(7, 9)", "0");
assert_expression_approximates_to<double>("binomial(7, 9)", "0");
assert_expression_approximates_to<float>("binomial(-7, 9)", "-5005");
assert_expression_approximates_to<double>("binomial(-7, 9)", "-5005");
assert_expression_approximates_to<float>("binompdf(4.4, 9, 0.7)", "0.0735138", Degree, Metric, Cartesian, 6); // FIXME: precision problem
assert_expression_approximates_to<double>("binompdf(4.4, 9, 0.7)", "0.073513818");
@@ -461,15 +471,6 @@ QUIZ_CASE(poincare_approximation_function) {
assert_expression_approximation_is_bounded("randint(4,45)", 4.0f, 45.0f, true);
assert_expression_approximation_is_bounded("randint(4,45)", 4.0, 45.0, true);
assert_expression_approximates_to<float>("binomial(12, 3)", "220");
assert_expression_approximates_to<double>("binomial(12, 3)", "220");
assert_expression_approximates_to<float>("binomial(-4.6, 3)", "-28.336");
assert_expression_approximates_to<double>("binomial(-4.6, 3)", "-28.336");
assert_expression_approximates_to<float>("binomial(π, 3)", "1.280108");
assert_expression_approximates_to<double>("binomial(π, 3)", "1.2801081307019");
}
QUIZ_CASE(poincare_approximation_trigonometry_functions) {

5
poincare/test/simplification.cpp
View File

@@ -687,6 +687,11 @@ QUIZ_CASE(poincare_simplification_function) {
assert_parsed_expression_simplify_to("arg(1+𝐢)", "π/4");
assert_parsed_expression_simplify_to("binomial(20,3)", "1140");
assert_parsed_expression_simplify_to("binomial(20,10)", "184756");
assert_parsed_expression_simplify_to("binomial(10,20)", "0");
assert_parsed_expression_simplify_to("binomial(-10,10)", "92378");
assert_parsed_expression_simplify_to("binomial(2.5,3)", "binomial(5/2,3)");
assert_parsed_expression_simplify_to("binomial(-200,120)", "binomial(-200,120)");
assert_parsed_expression_simplify_to("binomial(400,1)", "binomial(400,1)");
assert_parsed_expression_simplify_to("ceil(-1.3)", "-1");
assert_parsed_expression_simplify_to("ceil(2π)", "7");
assert_parsed_expression_simplify_to("ceil(123456789012345678901234567892/3)", "41152263004115226300411522631");