UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-01-19 00:37:25 +01:00
Code Issues Packages Projects Releases Wiki Activity

Performance fixes relating to floating point constants

Browse Source 
This replaces unnecessary double-precision soft-float operations with
single-precision floating-point operations, mainly by casting.

In a couple places I also replace a function call with a constant.
This commit is contained in:
Neven Sajko
2020-02-21 19:02:49 +00:00
committed by LeaNumworks
parent 068325d151
commit c92b770112
26 changed files with 46 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
apps/calculation/additional_outputs/complex_graph_cell.cpp
View File

@@ -39,7 +39,7 @@ void ComplexGraphView::drawRect(KDContext * ctx, KDRect rect) const {
* and the line of equation (real*t,imag*t).
* (a*cos(t), b*sin(t)) = (real*t,imag*t) --> tan(t) = sign(a)*sign(b) (± π)
* --> t = π/4 [π/2] according to sign(a) and sign(b). */
float th = real < 0.0f ? 3.0f*M_PI/4.0f : M_PI/4.0f;
float th = real < 0.0f ? 3.0f*(float)M_PI/4.0f : (float)M_PI/4.0f;
th = imag < 0.0f ? -th : th;
// Compute ellipsis parameters a and b
float factor = 5.0f;
@@ -48,7 +48,7 @@ void ComplexGraphView::drawRect(KDContext * ctx, KDRect rect) const {
// Avoid flat ellipsis for edge cases (for real = 0, the case imag = 0 is excluded)
if (real == 0.0f) {
a = 1.0f/factor;
th = imag < 0.0f ? -M_PI/2.0f : M_PI/2.0f;
th = imag < 0.0f ? -(float)M_PI/2.0f : (float)M_PI/2.0f;
}
std::complex<float> parameters(a,b);
drawCurve(ctx, rect, 0.0f, 1.0f, 0.01f,

2
apps/calculation/additional_outputs/trigonometry_model.h
View File

@@ -18,7 +18,7 @@ public:
float yMax() const override { return yCenter() + yHalfRange(); }
void setAngle(float f) { m_angle = f; }
float angle() const { return m_angle*M_PI/Poincare::Trigonometry::PiInAngleUnit(Poincare::Preferences::sharedPreferences()->angleUnit()); }
float angle() const { return m_angle*(float)M_PI/(float)Poincare::Trigonometry::PiInAngleUnit(Poincare::Preferences::sharedPreferences()->angleUnit()); }
private:
constexpr static float k_xHalfRange = 2.1f;
// We center the yRange around the semi-circle where the angle is

2
apps/probability/distribution/binomial_distribution.cpp
View File

@@ -66,7 +66,7 @@ double BinomialDistribution::cumulativeDistributiveInverseForProbability(double
}
double BinomialDistribution::rightIntegralInverseForProbability(double * probability) {
if (m_parameter1 == 0.0 && (m_parameter2 == 0.0 || m_parameter2 == 1.0)) {
if (m_parameter1 == 0.0f && (m_parameter2 == 0.0f || m_parameter2 == 1.0f)) {
return NAN;
}
if (*probability <= 0.0) {

2
apps/probability/distribution/chi_squared_distribution.cpp
View File

@@ -44,7 +44,7 @@ double ChiSquaredDistribution::cumulativeDistributiveFunctionAtAbscissa(double x
return 0.0;
}
double result = 0.0;
if (regularizedGamma(m_parameter1/2.0, x/2.0, k_regularizedGammaPrecision, k_maxRegularizedGammaIterations, &result)) {
if (regularizedGamma(m_parameter1/2.0f, x/2.0, k_regularizedGammaPrecision, k_maxRegularizedGammaIterations, &result)) {
return result;
}
return NAN;

2
apps/probability/distribution/student_distribution.cpp
View File

@@ -23,7 +23,7 @@ float StudentDistribution::evaluateAtAbscissa(float x) const {
}
bool StudentDistribution::authorizedValueAtIndex(float x, int index) const {
return x >= FLT_EPSILON && x <= 200.0; // We cannot draw the curve for x > 200 (coefficient() is too small)
return x >= FLT_EPSILON && x <= 200.0f; // We cannot draw the curve for x > 200 (coefficient() is too small)
}
double StudentDistribution::cumulativeDistributiveFunctionAtAbscissa(double x) const {

4
apps/shared/function_graph_view.cpp
View File

@@ -69,8 +69,8 @@ void FunctionGraphView::reloadBetweenBounds(float start, float end) {
if (start == end) {
return;
}
float pixelLowerBound = floatToPixel(Axis::Horizontal, start)-2.0;
float pixelUpperBound = floatToPixel(Axis::Horizontal, end)+4.0;
float pixelLowerBound = floatToPixel(Axis::Horizontal, start) - 2.0f;
float pixelUpperBound = floatToPixel(Axis::Horizontal, end) + 4.0f;
/* We exclude the banner frame from the dirty zone to avoid unnecessary
* redrawing */
KDRect dirtyZone(KDRect(pixelLowerBound, 0, pixelUpperBound-pixelLowerBound,

4
poincare/include/poincare/ieee754.h
View File

@@ -53,7 +53,7 @@ public:
}
static int exponentBase10(T f) {
static T k_log10base2 = 3.321928094887362347870319429489390175864831393024580612054;
if (f == 0.0) {
if (f == (T)0.0) {
return 0;
}
T exponentBase2 = exponent(f);
@@ -68,7 +68,7 @@ public:
* in -0.31 < x < 1, we get:
* e2 = [e1/log(10,2)] or e2 = [e1/log(10,2)]-1 depending on m1. */
int exponentBase10 = std::round(exponentBase2/k_log10base2);
if (std::pow(10.0, exponentBase10) > std::fabs(f)) {
if (std::pow((T)10.0, exponentBase10) > std::fabs(f)) {
exponentBase10--;
}
return exponentBase10;

4
poincare/src/approximation_helper.cpp
View File

@@ -19,7 +19,7 @@ template < typename T> T minimalNonNullMagnitudeOfParts(std::complex<T> c) {
T absRealInput = std::fabs(c.real());
T absImagInput = std::fabs(c.imag());
// If the magnitude of one part is null, ignore it
if (absRealInput == 0.0 || (absImagInput > 0.0 && absImagInput < absRealInput)) {
if (absRealInput == (T)0.0 || (absImagInput > (T)0.0 && absImagInput < absRealInput)) {
return absImagInput;
}
return absRealInput;
@@ -52,7 +52,7 @@ template <typename T> std::complex<T> ApproximationHelper::NeglectRealOrImaginar
}
T magnitude1 = minimalNonNullMagnitudeOfParts(input1);
T magnitude2 = minimalNonNullMagnitudeOfParts(input2);
T precision = 10.0*Expression::Epsilon<T>();
T precision = ((T)10.0)*Expression::Epsilon<T>();
if (isNegligeable(result.imag(), precision, magnitude1, magnitude2)) {
result.imag(0);
}

2
poincare/src/arc_cosine.cpp
View File

@@ -26,7 +26,7 @@ Expression ArcCosineNode::shallowReduce(ReductionContext reductionContext) {
template<typename T>
Complex<T> ArcCosineNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
std::complex<T> result;
if (c.imag() == 0 && std::fabs(c.real()) <= 1.0) {
if (c.imag() == 0 && std::fabs(c.real()) <= (T)1.0) {
/* acos: [-1;1] -> R
* In these cases we rather use std::acos(double) because acos on complexes
* is not as precise as pow on double in std library. For instance,

2
poincare/src/arc_sine.cpp
View File

@@ -26,7 +26,7 @@ Expression ArcSineNode::shallowReduce(ReductionContext reductionContext) {
template<typename T>
Complex<T> ArcSineNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
std::complex<T> result;
if (c.imag() == 0 && std::fabs(c.real()) <= 1.0) {
if (c.imag() == 0 && std::fabs(c.real()) <= (T)1.0) {
/* asin: [-1;1] -> R
* In these cases we rather use std::asin(double) because asin on complexes
* is not as precise as asin on double in std library. For instance,

2
poincare/src/arc_tangent.cpp
View File

@@ -22,7 +22,7 @@ int ArcTangentNode::serialize(char * buffer, int bufferSize, Preferences::PrintF
template<typename T>
Complex<T> ArcTangentNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
std::complex<T> result;
if (c.imag() == 0 && std::fabs(c.real()) <= 1.0) {
if (c.imag() == 0 && std::fabs(c.real()) <= (T)1.0) {
/* atan: R -> R
* In these cases we rather use std::atan(double) because atan on complexes
* is not as precise as atan on double in std library. For instance,

6
poincare/src/complex.cpp
View File

@@ -26,7 +26,7 @@ ComplexNode<T>::ComplexNode(std::complex<T> c) :
EvaluationNode<T>(),
std::complex<T>(c.real(), c.imag())
{
if (!std::isnan(c.imag()) && c.imag() != 0.0) {
if (!std::isnan(c.imag()) && c.imag() != (T)0.0) {
Expression::SetEncounteredComplex(true);
}
if (this->real() == -0) {
@@ -39,7 +39,7 @@ ComplexNode<T>::ComplexNode(std::complex<T> c) :
template<typename T>
T ComplexNode<T>::toScalar() const {
if (this->imag() == 0.0) {
if (this->imag() == (T)0.0) {
return this->real();
}
return NAN;
@@ -63,7 +63,7 @@ Expression ComplexNode<T>::complexToExpression(Preferences::ComplexFormat comple
Number::DecimalNumber<T>(std::fabs(tb)),
complexFormat,
(std::isnan(this->real()) || std::isnan(this->imag())),
ra == 0.0, std::fabs(ra) == 1.0, tb == 0.0, std::fabs(tb) == 1.0, ra < 0.0, tb < 0.0
ra == (T)0.0, std::fabs(ra) == (T)1.0, tb == (T)0.0, std::fabs(tb) == (T)1.0, ra < (T)0.0, tb < (T)0.0
);
}

2
poincare/src/complex_cartesian.cpp
View File

@@ -35,7 +35,7 @@ Complex<T> ComplexCartesianNode::templatedApproximate(Context * context, Prefere
assert(imagEvalution.type() == EvaluationNode<T>::Type::Complex);
std::complex<T> a = static_cast<Complex<T> &>(realEvaluation).stdComplex();
std::complex<T> b = static_cast<Complex<T> &>(imagEvalution).stdComplex();
if ((a.imag() != 0.0 && !std::isnan(a.imag())) || (b.imag() != 0.0 && !std::isnan(b.imag()))) {
if ((a.imag() != (T)0.0 && !std::isnan(a.imag())) || (b.imag() != (T)0.0 && !std::isnan(b.imag()))) {
/* a and b are supposed to be real (if they are not undefined). However,
* due to double precision limit, the approximation of the real part or the
* imaginary part can lead to complex values.

2
poincare/src/derivative.cpp
View File

@@ -56,7 +56,7 @@ Evaluation<T> DerivativeNode::templatedApproximate(Context * context, Preference
do {
T currentError;
T currentResult = riddersApproximation(context, complexFormat, angleUnit, evaluationArgument, h, &currentError);
h /= 10.0;
h /= (T)10.0;
if (std::isnan(currentError) || currentError > error) {
continue;
}

2
poincare/src/division.cpp
View File

@@ -47,7 +47,7 @@ Expression DivisionNode::shallowReduce(ReductionContext reductionContext) {
}
template<typename T> Complex<T> DivisionNode::compute(const std::complex<T> c, const std::complex<T> d, Preferences::ComplexFormat complexFormat) {
if (d.real() == 0.0 && d.imag() == 0.0) {
if (d.real() == (T)0.0 && d.imag() == (T)0.0) {
return Complex<T>::Undefined();
}
return Complex<T>::Builder(c/d);

8
poincare/src/integral.cpp
View File

@@ -128,8 +128,8 @@ IntegralNode::DetailedResult<T> IntegralNode::kronrodGaussQuadrature(T a, T b, C
T fv1[10];
T fv2[10];
T center = 0.5 * (a+b);
T halfLength = 0.5 * (b-a);
T center = (T)0.5 * (a+b);
T halfLength = (T)0.5 * (b-a);
T absHalfLength = std::fabs(halfLength);
DetailedResult<T> errorResult;
@@ -163,7 +163,7 @@ IntegralNode::DetailedResult<T> IntegralNode::kronrodGaussQuadrature(T a, T b, C
absKronrodIntegral += wKronrod[j] * (std::fabs(fval1) + std::fabs(fval2));
}
T halfKronrodIntegral = 0.5 * kronrodIntegral;
T halfKronrodIntegral = (T)0.5 * kronrodIntegral;
T kronrodIntegralDifference = wKronrod[10] * std::fabs(fCenter - halfKronrodIntegral);
for (int j = 0; j < 10; j++) {
kronrodIntegralDifference += wKronrod[j] * (std::fabs(fv1[j] - halfKronrodIntegral) + std::fabs(fv2[j] - halfKronrodIntegral));
@@ -176,7 +176,7 @@ IntegralNode::DetailedResult<T> IntegralNode::kronrodGaussQuadrature(T a, T b, C
T errorCoefficient = std::pow((T)(200*absError/kronrodIntegralDifference), (T)1.5);
absError = 1 > errorCoefficient ? kronrodIntegralDifference * errorCoefficient : kronrodIntegralDifference;
}
if (absKronrodIntegral > max/(50.0 * epsilon)) {
if (absKronrodIntegral > max/((T)50.0 * epsilon)) {
T minError = epsilon * 50 * absKronrodIntegral;
absError = absError > minError ? absError : minError;
}

2
poincare/src/logarithm.cpp
View File

@@ -279,7 +279,7 @@ Expression Logarithm::simpleShallowReduce(Context * context, Preferences::Comple
if (logDenominator.imag() != 0.0f || logDenominator.real() == 0.0f) {
result = Undefined::Builder();
}
isNegative = logDenominator.real() > 0.0;
isNegative = logDenominator.real() > 0.0f;
result = result.isUninitialized() ? Infinity::Builder(isNegative) : result;
replaceWithInPlace(result);
return result;

4
poincare/src/matrix.cpp
View File

@@ -270,7 +270,7 @@ void Matrix::ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns,
// No non-null coefficient in this column, skip
k++;
// Update determinant: det *= 0
if (determinant) { *determinant *= 0.0; }
if (determinant) { *determinant *= (T)0.0; }
} else {
// Swap row h and iPivot
if (iPivot != h) {
@@ -281,7 +281,7 @@ void Matrix::ArrayRowCanonize(T * array, int numberOfRows, int numberOfColumns,
array[h*numberOfColumns+col] = temp;
}
// Update determinant: det *= -1
if (determinant) { *determinant *= -1.0; }
if (determinant) { *determinant *= (T)-1.0; }
}
// Set to 1 array[h][k] by linear combination
T divisor = array[h*numberOfColumns+k];

8
poincare/src/normal_distribution.cpp
View File

@@ -5,6 +5,8 @@
#include <float.h>
#include <assert.h>
#define M_SQRT_2PI 2.506628274631000502415765284811
namespace Poincare {
template<typename T>
@@ -13,7 +15,7 @@ T NormalDistribution::EvaluateAtAbscissa(T x, T mu, T sigma) {
return NAN;
}
const float xMinusMuOverVar = (x - mu)/sigma;
return ((T)1.0)/(std::fabs(sigma) * std::sqrt(((T)2.0) * M_PI)) * std::exp(-((T)0.5) * xMinusMuOverVar * xMinusMuOverVar);
return ((T)1.0)/(std::fabs(sigma) * (T)M_SQRT_2PI) * std::exp(-((T)0.5) * xMinusMuOverVar * xMinusMuOverVar);
}
template<typename T>
@@ -95,7 +97,7 @@ T NormalDistribution::StandardNormalCumulativeDistributiveFunctionAtAbscissa(T a
if (abscissa < (T)0.0) {
return ((T)1.0) - StandardNormalCumulativeDistributiveFunctionAtAbscissa(-abscissa);
}
return ((T)0.5) + ((T)0.5) * std::erf(abscissa/std::sqrt(((T)2.0)));
return ((T)0.5) + ((T)0.5) * std::erf(abscissa/(T)M_SQRT2);
}
template<typename T>
@@ -113,7 +115,7 @@ T NormalDistribution::StandardNormalCumulativeDistributiveInverseForProbability(
if (probability < (T)0.5) {
return -StandardNormalCumulativeDistributiveInverseForProbability(((T)1.0)-probability);
}
return std::sqrt((T)2.0) * erfInv(((T)2.0) * probability - (T)1.0);
return (T)M_SQRT2 * erfInv(((T)2.0) * probability - (T)1.0);
}
template float NormalDistribution::EvaluateAtAbscissa<float>(float, float, float);

2
poincare/src/number.cpp
View File

@@ -72,7 +72,7 @@ Number Number::DecimalNumber(T f) {
return Undefined::Builder();
}
if (std::isinf(f)) {
return Infinity::Builder(f < 0.0);
return Infinity::Builder(f < (T)0.0);
}
return Decimal::Builder(f);
}

6
poincare/src/power.cpp
View File

@@ -136,7 +136,7 @@ Complex<T> PowerNode::computeNotPrincipalRealRootOfRationalPow(const std::comple
* where the principal root is real as these cases are handled generically
* later (for instance 1232^(1/8) which has a real principal root is not
* handled here). */
if (c.imag() == 0 && std::pow((T)-1.0, q) < 0.0) {
if (c.imag() == (T)0.0 && std::pow((T)-1.0, q) < (T)0.0) {
/* If c real and q odd integer (q odd if (-1)^q = -1), a real root does
* exist (which is not necessarily the principal root)!
* For q even integer, a real root does not necessarily exist (example:
@@ -148,7 +148,7 @@ Complex<T> PowerNode::computeNotPrincipalRealRootOfRationalPow(const std::comple
/* As q is odd, c^(p/q) = (sign(c)^(1/q))^p * |c|^(p/q)
* = sign(c)^p * |c|^(p/q)
* = -|c|^(p/q) iff c < 0 and p odd */
return c.real() < 0 && std::pow((T)-1.0, p) < 0.0 ? Complex<T>::Builder(-absCPowD.stdComplex()) : absCPowD;
return c.real() < (T)0.0 && std::pow((T)-1.0, p) < (T)0.0 ? Complex<T>::Builder(-absCPowD.stdComplex()) : absCPowD;
}
return Complex<T>::Undefined();
}
@@ -156,7 +156,7 @@ Complex<T> PowerNode::computeNotPrincipalRealRootOfRationalPow(const std::comple
template<typename T>
Complex<T> PowerNode::compute(const std::complex<T> c, const std::complex<T> d, Preferences::ComplexFormat complexFormat) {
std::complex<T> result;
if (c.imag() == 0.0 && d.imag() == 0.0 && c.real() != 0.0 && (c.real() > 0.0 || std::round(d.real()) == d.real())) {
if (c.imag() == (T)0.0 && d.imag() == (T)0.0 && c.real() != (T)0.0 && (c.real() > (T)0.0 || std::round(d.real()) == d.real())) {
/* pow: (R+, R) -> R+ (2^1.3 ~ 2.46)
* pow: (R-, N) -> R+ ((-2)^3 = -8)
* In these cases we rather use std::pow(double, double) because:

4
poincare/src/prediction_interval.cpp
View File

@@ -40,8 +40,8 @@ Evaluation<T> PredictionIntervalNode::templatedApproximate(Context * context, Pr
return Complex<T>::RealUndefined();
}
std::complex<T> operands[2];
operands[0] = std::complex<T>(p - 1.96*std::sqrt(p*(1.0-p))/std::sqrt(n));
operands[1] = std::complex<T>(p + 1.96*std::sqrt(p*(1.0-p))/std::sqrt(n));
operands[0] = std::complex<T>(p - (T)1.96*std::sqrt(p*((T)1.0-p))/std::sqrt(n));
operands[1] = std::complex<T>(p + (T)1.96*std::sqrt(p*((T)1.0-p))/std::sqrt(n));
return MatrixComplex<T>::Builder(operands, 1, 2);
}

2
poincare/src/print_float.cpp
View File

@@ -238,7 +238,7 @@ PrintFloat::TextLengths PrintFloat::ConvertFloatToTextPrivate(T f, char * buffer
// Correct the number of digits in mantissa after rounding
if (IEEE754<T>::exponentBase10(mantissa) >= numberOfSignificantDigits) {
mantissa = mantissa/10.0;
mantissa = mantissa / (T)10.0;
}
// Number of chars for the mantissa

4
poincare/src/randint.cpp
View File

@@ -49,10 +49,10 @@ template <typename T> Evaluation<T> RandintNode::templateApproximate(Context * c
if (std::isnan(a) || std::isnan(b) || std::isinf(a) || std::isinf(b)
|| a > b
|| a != (int)a || b != (int)b
|| (Expression::Epsilon<T>()*(b+1.0-a) > 1.0)) {
|| (Expression::Epsilon<T>()*(b+(T)1.0-a) > (T)1.0)) {
return Complex<T>::RealUndefined();
}
T result = std::floor(Random::random<T>()*(b+1.0-a)+a);
T result = std::floor(Random::random<T>()*(b+(T)1.0-a)+a);
return Complex<T>::Builder(result);
}

4
poincare/src/solver.cpp
View File

@@ -230,11 +230,11 @@ T Solver::CumulativeDistributiveInverseForNDefinedFunction(T * probability, Valu
do {
delta = std::fabs(*probability-p);
p += evaluation(k++, context, complexFormat, angleUnit, context1, context2, context3);
if (p >= k_maxProbability && std::fabs(*probability-1.0) <= delta) {
if (p >= k_maxProbability && std::fabs(*probability-(T)1.0) <= delta) {
*probability = (T)1.0;
return (T)(k-1);
}
} while (std::fabs(*probability-p) <= delta && k < k_maxNumberOfOperations && p < 1.0);
} while (std::fabs(*probability-p) <= delta && k < k_maxNumberOfOperations && p < (T)1.0);
p -= evaluation(--k, context, complexFormat, angleUnit, context1, context2, context3);
if (k == k_maxNumberOfOperations) {
*probability = (T)1.0;

4
poincare/src/trigonometry.cpp
View File

@@ -399,7 +399,7 @@ Expression Trigonometry::shallowReduceInverseFunction(Expression & e, Expressio
template <typename T>
std::complex<T> Trigonometry::ConvertToRadian(const std::complex<T> c, Preferences::AngleUnit angleUnit) {
if (angleUnit != Preferences::AngleUnit::Radian) {
return c * std::complex<T>(M_PI/Trigonometry::PiInAngleUnit(angleUnit));
return c * std::complex<T>((T)M_PI/(T)Trigonometry::PiInAngleUnit(angleUnit));
}
return c;
}
@@ -407,7 +407,7 @@ std::complex<T> Trigonometry::ConvertToRadian(const std::complex<T> c, Preferenc
template <typename T>
std::complex<T> Trigonometry::ConvertRadianToAngleUnit(const std::complex<T> c, Preferences::AngleUnit angleUnit) {
if (angleUnit != Preferences::AngleUnit::Radian) {
return c * std::complex<T>(Trigonometry::PiInAngleUnit(angleUnit)/M_PI);
return c * std::complex<T>((T)Trigonometry::PiInAngleUnit(angleUnit)/(T)M_PI);
}
return c;
}