122 lines
3.8 KiB
C++
122 lines
3.8 KiB
C++
/*
|
|
* Licensed under the Apache License, Version 2.0 (the "License");
|
|
* you may not use this file except in compliance with the License.
|
|
* You may obtain a copy of the License at
|
|
*
|
|
* http://www.apache.org/licenses/LICENSE-2.0
|
|
*
|
|
* Unless required by applicable law or agreed to in writing, software
|
|
* distributed under the License is distributed on an "AS IS" BASIS,
|
|
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|
* See the License for the specific language governing permissions and
|
|
* limitations under the License.
|
|
*
|
|
* Copyright 2018 Danny Robson <danny@nerdcruft.net>
|
|
*/
|
|
|
|
#ifndef CRUFT_UTIL_MATHS_FAST_HPP
|
|
#define CRUFT_UTIL_MATHS_FAST_HPP
|
|
|
|
#include "../maths.hpp"
|
|
#include "../debug.hpp"
|
|
|
|
namespace util::maths::fast {
|
|
float
|
|
estrin (float t, float a, float b, float c, float d)
|
|
{
|
|
const auto t2 = t * t;
|
|
return std::fma (b, t, a) + std::fma (d, t, c) * t2;
|
|
}
|
|
|
|
// approximates sin over the range [0; pi/4] using a polynomial.
|
|
//
|
|
// the relative error should be on the order of 1e-6.
|
|
//
|
|
// sollya:
|
|
// display=hexadecimal;
|
|
// canonical=on!;
|
|
// P = fpminimax(sin(x), [|3,5,7|], [|single...|],[0x1p-126;pi/4], , relative, x);
|
|
// P;
|
|
// dirtyinfnorm(sin(x)-P, [0;pi/4]);
|
|
// plot(P(x)-sin(x),[0x1p-126;pi/4]);
|
|
//
|
|
// using minimax(sin(x)-x, [|3,5,7,9|], ...) has higher accuracy, but
|
|
// probably not enough to offset the extra multiplications.
|
|
float
|
|
sin (float x) noexcept
|
|
{
|
|
CHECK_LT (x, 2 * pi<float>);
|
|
CHECK_GT (x, -2 * pi<float>);
|
|
|
|
// flip negative angles due to having an even function
|
|
float neg = x < 0 ? x *= -1, -1 : 1;
|
|
|
|
// reduce the angle due to the period
|
|
int periods = static_cast<int> (x / 2 / pi<float>);
|
|
x -= periods * 2 * pi<float>;
|
|
|
|
// 1,3,5,7: 1, -0x1.555546p-3f, 0x1.11077ap-7f, -0x1.9954e8p-13f;
|
|
// 1,3,5,7,9: 1, -0x1.555556p-3f, 0x1.11115cp-7f, -0x1.a0403ap-13f, 0x1.75dc26p-19f;
|
|
|
|
// sin(x) = x + Ax3 + Bx5 + Cx7
|
|
// sin(x) = x (1 + Ax2 + Bx4 + Cx6)
|
|
// sin(x) = x (1 + Ay1 + By2 + Cy3), y = x^2
|
|
|
|
return neg * x * estrin (x * x, 1, -0x1.555546p-3f, 0x1.11077ap-7f, -0x1.9954e8p-13f);
|
|
}
|
|
|
|
// calculates an approximation of e^x
|
|
//
|
|
// we split the supplied value into integer and fractional components as
|
|
// use the identity: e^(a+b) == e^a * e^b;
|
|
//
|
|
// integer powers can be computed efficiently using:
|
|
// e^x
|
|
// = 2^(x/ln2)
|
|
// = 2^kx, where k = 1/ln2
|
|
// we set the floating point exponent directly with kx
|
|
//
|
|
// the fractional component is approximated by a polynomial across [0;1]
|
|
//
|
|
// we force the first two terms as 1 and x because they are fantastically
|
|
// efficient to load in this context.
|
|
//
|
|
// using an order 5 polynomial is overkill, but we use something with 4
|
|
// coefficients because it should be pretty efficient to evaluate using
|
|
// SIMD.
|
|
//
|
|
// sollya:
|
|
// display=hexadecimal;
|
|
// canonical=on!;
|
|
// P = fpminimax(exp(x), [|2,3,4,5|], [|single...|], [0x1p-126;1], 1+x);
|
|
// P;
|
|
// dirtyinfnorm(P(x)-exp(x), [0;1]);
|
|
// plot(P(x)-exp(x),[0; 1]);
|
|
float
|
|
exp (float x)
|
|
{
|
|
union {
|
|
int32_t i;
|
|
float f;
|
|
} whole;
|
|
|
|
whole.i = static_cast<int32_t> (x);
|
|
float frac = x - whole.i;
|
|
|
|
whole.i *= 0b00000000101110001010101000111011;
|
|
//whole.i *= static_cast<int32_t> ((1<<(23)) / std::log(2));
|
|
whole.i += 0b00111111100000000000000000000000; //127 << 23;
|
|
|
|
frac = 1 + frac + frac * frac * estrin (frac,
|
|
0.4997530281543731689453125f,
|
|
0.16853784024715423583984375,
|
|
3.6950431764125823974609375e-2,
|
|
1.303750835359096527099609375e-2
|
|
);
|
|
|
|
return frac * whole.f;
|
|
}
|
|
}
|
|
|
|
#endif
|