6663b6c9
 adorian
projet complet av...
|
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
#include <poincare/derivative.h>
#include <poincare/symbol.h>
#include <poincare/complex.h>
#include <poincare/simplification_engine.h>
#include <poincare/undefined.h>
#include <cmath>
extern "C" {
#include <assert.h>
#include <float.h>
}
namespace Poincare {
Expression::Type Derivative::type() const {
return Type::Derivative;
}
Expression * Derivative::clone() const {
Derivative * a = new Derivative(m_operands, true);
return a;
}
Expression * Derivative::shallowReduce(Context& context, AngleUnit angleUnit) {
Expression * e = Expression::shallowReduce(context, angleUnit);
if (e != this) {
return e;
}
#if MATRIX_EXACT_REDUCING
if (operand(0)->type() == Type::Matrix || operand(1)->type() == Type::Matrix) {
return replaceWith(new Undefined(), true);
}
#endif
// TODO: to be implemented diff(+) -> +diff() etc
return this;
}
template<typename T>
Expression * Derivative::templatedApproximate(Context& context, AngleUnit angleUnit) const {
static T min = sizeof(T) == sizeof(double) ? DBL_MIN : FLT_MIN;
static T epsilon = sizeof(T) == sizeof(double) ? DBL_EPSILON : FLT_EPSILON;
VariableContext<T> xContext = VariableContext<T>('x', &context);
Symbol xSymbol('x');
Expression * xInput = operand(1)->approximate<T>(context, angleUnit);
T x = xInput->type() == Type::Complex ? static_cast<Complex<T> *>(xInput)->toScalar() : NAN;
delete xInput;
Complex<T> e = Complex<T>::Float(x);
xContext.setExpressionForSymbolName(&e, &xSymbol, xContext);
Expression * fInput = operand(0)->approximate<T>(xContext, angleUnit);
T functionValue = fInput->type() == Type::Complex ? static_cast<Complex<T> *>(fInput)->toScalar() : NAN;
delete fInput;
// No complex/matrix version of Derivative
if (std::isnan(x) || std::isnan(functionValue)) {
return new Complex<T>(Complex<T>::Float(NAN));
}
T error, result;
T h = k_minInitialRate;
do {
result = riddersApproximation(xContext, angleUnit, x, h, &error);
h /= 10.0;
} while ((std::fabs(error/result) > k_maxErrorRateOnApproximation || std::isnan(error)) && h >= epsilon);
/* if the error is too big regarding the value, do not return the answer */
if (std::fabs(error/result) > k_maxErrorRateOnApproximation || std::isnan(error)) {
return new Complex<T>(Complex<T>::Float(NAN));
}
if (std::fabs(error) < min) {
return new Complex<T>(Complex<T>::Float(result));
}
error = std::pow((T)10, std::floor(std::log10(std::fabs(error)))+2);
return new Complex<T>(Complex<T>::Float(std::round(result/error)*error));
}
template<typename T>
T Derivative::growthRateAroundAbscissa(T x, T h, VariableContext<T> xContext, AngleUnit angleUnit) const {
Symbol xSymbol('x');
Complex<T> e = Complex<T>::Float(x + h);
xContext.setExpressionForSymbolName(&e, &xSymbol, xContext);
Expression * fInput = operand(0)->approximate<T>(xContext, angleUnit);
T expressionPlus = fInput->type() == Type::Complex ? static_cast<Complex<T> *>(fInput)->toScalar() : NAN;
delete fInput;
e = Complex<T>::Float(x-h);
xContext.setExpressionForSymbolName(&e, &xSymbol, xContext);
fInput = operand(0)->approximate<T>(xContext, angleUnit);
T expressionMinus = fInput->type() == Type::Complex ? static_cast<Complex<T> *>(fInput)->toScalar() : NAN;
delete fInput;
return (expressionPlus - expressionMinus)/(2*h);
}
template<typename T>
T Derivative::riddersApproximation(VariableContext<T> xContext, AngleUnit angleUnit, T x, T h, T * error) const {
/* Ridders' Algorithm
* Blibliography:
* - Ridders, C.J.F. 1982, Advances in Engineering Software, vol. 4, no. 2,
* pp. 75–76. */
*error = sizeof(T) == sizeof(float) ? FLT_MAX : DBL_MAX;
// Initialize hh
assert(h != 0.0);
/* Make hh an exactly representable number */
volatile T temp = x+h;
T hh = temp - x;
/* a is matrix storing the function extrapolations for different stepsizes at
* different order */
T a[10][10];
for (int i = 0; i < 10; i++) {
for (int j = 0; j < 10; j++) {
a[i][j] = 1;
}
}
a[0][0] = growthRateAroundAbscissa(x, hh, xContext, angleUnit);
T ans = 0;
T errt = 0;
/* Loop on i: change the step size */
for (int i = 1; i < 10; i++) {
hh /= k_rateStepSize;
/* Make hh an exactly representable number */
volatile T temp = x+hh;
hh = temp - x;
a[0][i] = growthRateAroundAbscissa(x, hh, xContext, angleUnit);
T fac = k_rateStepSize*k_rateStepSize;
/* Loop on j: compute extrapolation for several orders */
for (int j = 1; j < 10; j++) {
a[j][i] = (a[j-1][i]*fac-a[j-1][i-1])/(fac-1);
fac = k_rateStepSize*k_rateStepSize*fac;
errt = std::fabs(a[j][i]-a[j-1][i]) > std::fabs(a[j][i]-a[j-1][i-1]) ? std::fabs(a[j][i]-a[j-1][i]) : std::fabs(a[j][i]-a[j-1][i-1]);
/* Update error and answer if error decreases */
if (errt < *error) {
*error = errt;
ans = a[j][i];
}
}
/* If higher extrapolation order significantly increases the error, return
* early */
if (std::fabs(a[i][i]-a[i-1][i-1]) > 2*(*error)) {
break;
}
}
return ans;
}
}
|