binomial_coefficient.cpp
3.54 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#include <poincare/binomial_coefficient.h>
#include <poincare/undefined.h>
#include <poincare/rational.h>
#include "layout/binomial_coefficient_layout.h"
extern "C" {
#include <stdlib.h>
#include <assert.h>
}
#include <cmath>
namespace Poincare {
Expression::Type BinomialCoefficient::type() const {
return Type::BinomialCoefficient;
}
Expression * BinomialCoefficient::clone() const {
BinomialCoefficient * b = new BinomialCoefficient(m_operands, true);
return b;
}
Expression * BinomialCoefficient::shallowReduce(Context& context, AngleUnit angleUnit) {
Expression * e = Expression::shallowReduce(context, angleUnit);
if (e != this) {
return e;
}
Expression * op0 = editableOperand(0);
Expression * op1 = editableOperand(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);
}
ExpressionLayout * BinomialCoefficient::createLayout(PrintFloat::Mode floatDisplayMode, int numberOfSignificantDigits) const {
return new BinomialCoefficientLayout(
operand(0)->createLayout(floatDisplayMode, numberOfSignificantDigits),
operand(1)->createLayout(floatDisplayMode, numberOfSignificantDigits),
false);
}
template<typename T>
Complex<T> * BinomialCoefficient::templatedApproximate(Context& context, AngleUnit angleUnit) const {
Evaluation<T> * nInput = operand(0)->privateApproximate(T(), context, angleUnit);
Evaluation<T> * kInput = operand(1)->privateApproximate(T(), context, angleUnit);
T n = nInput->toScalar();
T k = kInput->toScalar();
delete nInput;
delete kInput;
return new Complex<T>(compute(k, n));
}
template<typename T>
T BinomialCoefficient::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;
}
T result = 1;
for (int i = 0; i < k; i++) {
result *= (n-(T)i)/(k-(T)i);
if (std::isinf(result) || std::isnan(result)) {
return result;
}
}
return std::round(result);
}
template double BinomialCoefficient::compute(double k, double n);
template float BinomialCoefficient::compute(float k, float n);
}