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

[poincare] Generalize Binomial Coefficient n to any real

Browse Source 
Change-Id: I3cda882e67db4becfe232c6428e14905e922662b
This commit is contained in:
Hugo Saint-Vignes
2020-05-29 15:42:35 +02:00
committed by Émilie Feral
parent ada369bf08
commit 735fc6e1bb
2 changed files with 30 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
poincare/src/binomial_coefficient.cpp
View File

@@ -40,18 +40,23 @@ Complex<T> BinomialCoefficientNode::templatedApproximate(Context * context, Pref
template<typename T>
T BinomialCoefficientNode::compute(T k, T n) {
k = k > (n-k) ? n-k : k;
if (std::isnan(n) || std::isnan(k) || n != std::round(n) || k != std::round(k) || k > n || k < 0 || n < 0) {
if (std::isnan(n) || std::isnan(k) || k != std::round(k) || k < 0) {
return NAN;
}
// Generalized definition allows any n value
bool generalized = (n != std::round(n) || n < k);
// Take advantage of symmetry
k = (!generalized && k > (n - k)) ? n - k : k;
T result = 1;
for (int i = 0; i < k; i++) {
result *= (n-(T)i)/(k-(T)i);
result *= (n - (T)i) / (k - (T)i);
if (std::isinf(result) || std::isnan(result)) {
return result;
}
}
return std::round(result);
// If not generalized, the output must be round
return generalized ? result : std::round(result);
}
@@ -70,28 +75,27 @@ Expression BinomialCoefficient::shallowReduce(Context * context) {
return replaceWithUndefinedInPlace();
}
if (c0.type() == ExpressionNode::Type::Rational) {
Rational r0 = static_cast<Rational&>(c0);
if (!r0.isInteger() || r0.isNegative()) {
return replaceWithUndefinedInPlace();
}
}
if (c1.type() == ExpressionNode::Type::Rational) {
Rational r1 = static_cast<Rational&>(c1);
if (!r1.isInteger() || r1.isNegative()) {
return replaceWithUndefinedInPlace();
}
}
if (c0.type() != ExpressionNode::Type::Rational || c1.type() != ExpressionNode::Type::Rational) {
return *this;
}
Rational r0 = static_cast<Rational&>(c0);
Rational r1 = static_cast<Rational&>(c1);
if (!r1.isInteger() || r1.isNegative()) {
return replaceWithUndefinedInPlace();
}
if (!r0.isInteger()) {
// Generalized binomial coefficient
return *this;
}
Integer n = r0.signedIntegerNumerator();
Integer k = r1.signedIntegerNumerator();
if (n.isLowerThan(k)) {
return replaceWithUndefinedInPlace();
// 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. */

9
poincare/test/approximation.cpp
View File

@@ -426,6 +426,15 @@ 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) {