[poincare] BinomialCoefficient

This commit is contained in:
Léa Saviot
2018-08-27 10:21:48 +02:00
parent 48a6d92676
commit b3fdf28ddc
6 changed files with 133 additions and 91 deletions

View File

@@ -40,6 +40,7 @@ objs += $(addprefix poincare/src/,\
arc_sine.o\
arc_tangent.o\
arithmetic.o\
binomial_coefficient.o\
complex.o\
cosine.o\
decimal.o\
@@ -214,6 +215,7 @@ tests += $(addprefix poincare/test/,\
tree/tree_by_reference.cpp\
expression.cpp\
addition.cpp\
function.cpp\
helper.cpp\
integer.cpp\
logarithm.cpp\

View File

@@ -110,6 +110,7 @@
#include <poincare/arc_cosine.h>
#include <poincare/arc_sine.h>
#include <poincare/arc_tangent.h>
#include <poincare/binomial_coefficient.h>
#include <poincare/cosine.h>
#include <poincare/sine.h>
#include <poincare/tangent.h>

View File

@@ -1,30 +1,56 @@
#ifndef POINCARE_BINOMIAL_COEFFICIENT_H
#define POINCARE_BINOMIAL_COEFFICIENT_H
#include <poincare/evaluation.h>
#include <poincare/approximation_helper.h>
#include <poincare/expression.h>
#include <poincare/layout_helper.h>
#include <poincare/static_hierarchy.h>
namespace Poincare {
class BinomialCoefficient : public StaticHierarchy<2> {
using StaticHierarchy<2>::StaticHierarchy;
class BinomialCoefficientNode : public ExpressionNode {
public:
Type type() const override;
static BinomialCoefficientNode * FailedAllocationStaticNode();
BinomialCoefficientNode * failedAllocationStaticNode() override { return FailedAllocationStaticNode(); }
// TreeNode
size_t size() const override { return sizeof(BinomialCoefficientNode); }
int numberOfChildren() const override { return 2; }
#if POINCARE_TREE_LOG
virtual void logNodeName(std::ostream & stream) const override {
stream << "BinomialCoefficient";
}
#endif
// ExpressionNode
// Properties
Type type() const override{ return Type::BinomialCoefficient; }
template<typename T> static T compute(T k, T n);
private:
constexpr static int k_maxNValue = 300;
/* Layout */
LayoutRef createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override;
int serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override {
return SerializationHelper::Prefix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, "binomial");
}
/* Simplification */
// Layout
LayoutReference createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override;
int serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const override;
// Simplification
Expression shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const override;
/* Evaluation */
// Evaluation
Evaluation<float> approximate(SinglePrecision p, Context& context, Preferences::AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
Evaluation<double> approximate(DoublePrecision p, Context& context, Preferences::AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
template<typename T> Evaluation<T> templatedApproximate(Context& context, Preferences::AngleUnit angleUnit) const;
template<typename T> Complex<T> templatedApproximate(Context& context, Preferences::AngleUnit angleUnit) const;
};
class BinomialCoefficient : public Expression {
public:
BinomialCoefficient() : Expression(TreePool::sharedPool()->createTreeNode<BinomialCoefficientNode>()) {}
BinomialCoefficient(const BinomialCoefficientNode * n) : Expression(n) {}
BinomialCoefficient(Expression child1, Expression child2) : BinomialCoefficient() {
replaceChildAtIndexInPlace(0, child1);
replaceChildAtIndexInPlace(1, child2);
}
// Expression
Expression shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const;
private:
constexpr static int k_maxNValue = 300;
};
}

View File

@@ -1,96 +1,45 @@
#include <poincare/binomial_coefficient.h>
#include <poincare/undefined.h>
#include <poincare/rational.h>
#include <poincare/binomial_coefficient_layout_node.h>
extern "C" {
#include <poincare/rational.h>
#include <poincare/serialization_helper.h>
#include <poincare/undefined.h>
#include <stdlib.h>
#include <assert.h>
}
#include <cmath>
namespace Poincare {
ExpressionNode::Type BinomialCoefficient::type() const {
return Type::BinomialCoefficient;
BinomialCoefficientNode * BinomialCoefficientNode::FailedAllocationStaticNode() {
static AllocationFailureExpressionNode<BinomialCoefficientNode> failure;
TreePool::sharedPool()->registerStaticNodeIfRequired(&failure);
return &failure;
}
Expression * BinomialCoefficient::clone() const {
BinomialCoefficient * b = new BinomialCoefficient(m_operands, true);
return b;
Expression BinomialCoefficientNode::shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const {
return BinomialCoefficient(this).shallowReduce(context, angleUnit);
}
Expression BinomialCoefficient::shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const {
Expression e = Expression::defaultShallowReduce(context, angleUnit);
if (e.isUndefinedOrAllocationFailure()) {
return e;
}
Expression * op0 = childAtIndex(0);
Expression * op1 = childAtIndex(1);
#if MATRIX_EXACT_REDUCING
if (op0->type() == Type::Matrix || op1->type() == Type::Matrix) {
return replaceWith(new Undefined(), true);
}
#endif
if (op0->type() == Type::Rational) {
Rational * r0 = static_cast<Rational *>(op0);
if (!r0->denominator().isOne() || r0->numerator().isNegative()) {
return replaceWith(new Undefined(), true);
}
}
if (op1->type() == Type::Rational) {
Rational * r1 = static_cast<Rational *>(op1);
if (!r1->denominator().isOne() || r1->numerator().isNegative()) {
return replaceWith(new Undefined(), true);
}
}
if (op0->type() != Type::Rational || op1->type() != Type::Rational) {
return this;
}
Rational * r0 = static_cast<Rational *>(op0);
Rational * r1 = static_cast<Rational *>(op1);
Integer n = r0->numerator();
Integer k = r1->numerator();
if (n.isLowerThan(k)) {
return replaceWith(new Undefined(), true);
}
/* if n is too big, we do not reduce to avoid too long computation.
* The binomial coefficient will be evaluate approximatively later */
if (Integer(k_maxNValue).isLowerThan(n)) {
return this;
}
Rational result(1);
Integer kBis = Integer::Subtraction(n, k);
k = kBis.isLowerThan(k) ? kBis : k;
int clippedK = k.extractedInt(); // Authorized because k < n < k_maxNValue
for (int i = 0; i < clippedK; i++) {
Rational factor = Rational(Integer::Subtraction(n, Integer(i)), Integer::Subtraction(k, Integer(i)));
result = Rational::Multiplication(result, factor);
}
return replaceWith(new Rational(result), true);
}
LayoutRef BinomialCoefficient::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
LayoutReference BinomialCoefficientNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
return BinomialCoefficientLayoutRef(
childAtIndex(0)->createLayout(floatDisplayMode, numberOfSignificantDigits),
childAtIndex(1)->createLayout(floatDisplayMode, numberOfSignificantDigits));
}
template<typename T>
Complex<T> * BinomialCoefficient::templatedApproximate(Context& context, Preferences::AngleUnit angleUnit) const {
Evaluation<T> * nInput = childAtIndex(0)->privateApproximate(T(), context, angleUnit);
Evaluation<T> * kInput = childAtIndex(1)->privateApproximate(T(), context, angleUnit);
T n = nInput->toScalar();
T k = kInput->toScalar();
delete nInput;
delete kInput;
return new Complex<T>(compute(k, n));
int BinomialCoefficientNode::serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
return SerializationHelper::Prefix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, "binomial");
}
template<typename T>
Complex<T> BinomialCoefficientNode::templatedApproximate(Context& context, Preferences::AngleUnit angleUnit) const {
Evaluation<T> nInput = childAtIndex(0)->approximate(T(), context, angleUnit);
Evaluation<T> kInput = childAtIndex(1)->approximate(T(), context, angleUnit);
T n = nInput.toScalar();
T k = kInput.toScalar();
return Complex<T>(compute(k, n));
}
template<typename T>
T BinomialCoefficient::compute(T k, T n) {
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) {
return NAN;
@@ -105,7 +54,58 @@ T BinomialCoefficient::compute(T k, T n) {
return std::round(result);
}
template double BinomialCoefficient::compute(double k, double n);
template float BinomialCoefficient::compute(float k, float n);
Expression BinomialCoefficient::shallowReduce(Context& context, Preferences::AngleUnit angleUnit) const {
Expression e = Expression::defaultShallowReduce(context, angleUnit);
if (e.isUndefinedOrAllocationFailure()) {
return e;
}
Expression op0 = childAtIndex(0);
Expression op1 = childAtIndex(1);
#if MATRIX_EXACT_REDUCING
if (op0.type() == ExpressionNode::Type::Matrix || op1.type() == ExpressionNode::Type::Matrix) {
return Undefined();
}
#endif
if (op0.type() == ExpressionNode::Type::Rational) {
Rational r0 = static_cast<Rational>(op0);
if (!r0.integerDenominator().isOne() || r0.integerDenominator().isNegative()) {
return Undefined();
}
}
if (op1.type() == ExpressionNode::Type::Rational) {
Rational r1 = static_cast<Rational>(op1);
if (!r1.integerDenominator().isOne() || r1.integerDenominator().isNegative()) {
return Undefined();
}
}
if (op0.type() != ExpressionNode::Type::Rational || op1.type() != ExpressionNode::Type::Rational) {
return *this;
}
Rational r0 = static_cast<Rational>(op0);
Rational r1 = static_cast<Rational>(op1);
Integer n = r0.signedIntegerNumerator();
Integer k = r1.signedIntegerNumerator();
if (n.isLowerThan(k)) {
return Undefined();
}
/* 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)) {
return *this;
}
Rational result(1);
Integer kBis = Integer::Subtraction(n, k);
k = kBis.isLowerThan(k) ? kBis : k;
int clippedK = k.extractedInt(); // Authorized because k < n < k_maxNValue
for (int i = 0; i < clippedK; i++) {
Rational factor = Rational(Integer::Subtraction(n, Integer(i)), Integer::Subtraction(k, Integer(i)));
result = Rational::Multiplication(result, factor);
}
return Rational(result);
}
template double BinomialCoefficientNode::compute(double k, double n);
template float BinomialCoefficientNode::compute(float k, float n);
}

View File

@@ -123,7 +123,8 @@ asin { poincare_expression_yylval.expression = ArcSine(); return FUNCTION; }
asinh { poincare_expression_yylval.expression = HyperbolicArcSine(); return FUNCTION; }
atan { poincare_expression_yylval.expression = ArcTangent(); return FUNCTION; }
atanh { poincare_expression_yylval.expression = HyperbolicArcTangent(); return FUNCTION; }
/*binomial { poincare_expression_yylval.expression = new BinomialCoefficient(); return FUNCTION; }
binomial { poincare_expression_yylval.expression = BinomialCoefficient(); return FUNCTION; }
/*
ceil { poincare_expression_yylval.expression = new Ceiling(); return FUNCTION; }
confidence { poincare_expression_yylval.expression = new ConfidenceInterval(); return FUNCTION; }
diff { poincare_expression_yylval.expression = new Derivative(); return FUNCTION; }

View File

@@ -8,17 +8,20 @@
using namespace Poincare;
template<typename T>
void assert_exp_is_bounded(Expression * exp, T lowBound, T upBound, bool upBoundIncluded = false) {
void assert_exp_is_bounded(Expression exp, T lowBound, T upBound, bool upBoundIncluded = false) {
GlobalContext globalContext;
T result = exp->approximateToScalar<T>(globalContext, Radian);
T result = exp.approximateToScalar<T>(globalContext, Radian);
assert(result >= lowBound);
assert(result < upBound || (result == upBound && upBoundIncluded));
}
QUIZ_CASE(poincare_parse_function) {
#if 0
assert_parsed_expression_type("abs(-1)", ExpressionNode::Type::AbsoluteValue);
assert_parsed_expression_type("arg(2+I)", ExpressionNode::Type::ComplexArgument);
#endif
assert_parsed_expression_type("binomial(10, 4)", ExpressionNode::Type::BinomialCoefficient);
#if 0
assert_parsed_expression_type("ceil(0.2)", ExpressionNode::Type::Ceiling);
assert_parsed_expression_type("diff(2*x, 2)", ExpressionNode::Type::Derivative);
#if MATRICES_ARE_DEFINED
@@ -57,10 +60,12 @@ QUIZ_CASE(poincare_parse_function) {
assert_parsed_expression_type("transpose([[1,2,3][4,5,6][7,8,9]])", ExpressionNode::Type::MatrixTranspose);
#endif
assert_parsed_expression_type("6!", ExpressionNode::Type::Factorial);
#endif
}
QUIZ_CASE(poincare_function_evaluate) {
#if 0
assert_parsed_expression_evaluates_to<float>("abs(-1)", "1");
assert_parsed_expression_evaluates_to<double>("abs(-1)", "1");
@@ -73,8 +78,10 @@ QUIZ_CASE(poincare_function_evaluate) {
assert_parsed_expression_evaluates_to<float>("abs([[3+2I,3+4I][5+2I,3+2I]])", "[[3.605551,5][5.385165,3.605551]]");
assert_parsed_expression_evaluates_to<double>("abs([[3+2I,3+4I][5+2I,3+2I]])", "[[3.605551275464,5][5.3851648071345,3.605551275464]]");
#endif
assert_parsed_expression_evaluates_to<float>("binomial(10, 4)", "210");
assert_parsed_expression_evaluates_to<double>("binomial(10, 4)", "210");
#if 0
assert_parsed_expression_evaluates_to<float>("ceil(0.2)", "1");
assert_parsed_expression_evaluates_to<double>("ceil(0.2)", "1");
@@ -223,13 +230,17 @@ QUIZ_CASE(poincare_function_evaluate) {
assert_exp_is_bounded(exp, 4.0f, 45.0f, true);
assert_exp_is_bounded(exp, 4.0, 45.0, true);
delete exp;
#endif
}
QUIZ_CASE(poincare_function_simplify) {
#if 0
assert_parsed_expression_simplify_to("abs(P)", "P");
assert_parsed_expression_simplify_to("abs(-P)", "P");
#endif
assert_parsed_expression_simplify_to("binomial(20,3)", "1140");
assert_parsed_expression_simplify_to("binomial(20,10)", "184756");
#if 0
assert_parsed_expression_simplify_to("ceil(-1.3)", "-1");
assert_parsed_expression_simplify_to("conj(1/2)", "1/2");
assert_parsed_expression_simplify_to("quo(19,3)", "6");
@@ -261,4 +272,5 @@ QUIZ_CASE(poincare_function_simplify) {
assert_parsed_expression_simplify_to("permute(99,4)", "90345024");
assert_parsed_expression_simplify_to("permute(20,-10)", "undef");
assert_parsed_expression_simplify_to("re(1/2)", "1/2");
#endif
}