summaryrefslogtreecommitdiffstats
path: root/mathfuncs_sin.h
diff options
context:
space:
mode:
Diffstat (limited to 'mathfuncs_sin.h')
-rw-r--r--mathfuncs_sin.h233
1 files changed, 38 insertions, 195 deletions
diff --git a/mathfuncs_sin.h b/mathfuncs_sin.h
index 107e6eb..9e1b900 100644
--- a/mathfuncs_sin.h
+++ b/mathfuncs_sin.h
@@ -15,40 +15,52 @@ namespace vecmathlib {
template<typename realvec_t>
realvec_t mathfuncs<realvec_t>::vml_sin(realvec_t x)
{
- typename realvec_t::scalar_t eps __attribute__((__unused__)) = 0.000001;
+ // Rescale input
+ x *= RV(1.0/(2.0*M_PI));
// Reduce range: sin(x) = sin(x + 2pi)
- x = remainder(x, RV(2.0*M_PI));
- assert(all(x >= RV(-(1.0+eps)*M_PI) && x <= RV((1.0+eps)*M_PI)));
+ x -= round(x);
+ VML_ASSERT(all(x >= RV(-0.5) && x <= RV(+0.5)));
// Reduce range: sin(x) = -sin(-x)
realvec_t sign = x;
x = fabs(x);
- assert(all(x >= RV(0.0) && x <= RV((1.0+eps)*M_PI)));
+ VML_ASSERT(all(x >= RV(0.0) && x <= RV(0.5)));
- // Reduce range: sin(x) = sin(pi-x)
- x = ifthen(x > RV(M_PI_2), RV(M_PI) - x, x);
- assert(all(x >= RV(0.0-eps) && x <= RV(M_PI_2)));
+ // Reduce range: sin(x) = sin(pi - x)
+ x = fmin(x, RV(0.5)-x);
+ VML_ASSERT(all(x >= RV(0.0) && x <= RV(0.25)));
- // Reduce range: cos(x) = sin(pi/2 - x)
- boolvec_t use_cos = x > RV(0.5*M_PI_2);
- x = ifthen(use_cos, RV(M_PI_2) - x, x);
-
- // Taylor expansion
- // sin(x) = Sum[n=0,nmax,n%2!=0] (-1)^((n-1)/2) x^n / n!
- // cos(x) = Sum[n=0,nmax,n%2==0] (-1)^((n-1)/2) x^n / n!
- realvec_t noffset = ifthen(use_cos, RV(-1.0), RV(0.0));
- int const nmax = 15;
+ // Polynomial expansion
realvec_t x2 = x*x;
- realvec_t y = ifthen(use_cos, RV(1.0), x);
- realvec_t r = y;
- for (int n=3; n<nmax; n+=2) {
- // sin: (n-1)*n
- // cos: n*(n+1)
- realvec_t rn = RV(FP::convert_float(n)) + noffset;
- y *= - x2 * rcp((rn-RV(1.0))*rn);
- r += y;
+ realvec_t r;
+ switch (sizeof(real_t)) {
+ case 4:
+ // float, error=9.63199952895688527467513447362e-9
+ r = RV(+39.6507114184650539767708004425);
+ r = fma(r, x2, RV(-76.564420291497946510082570076));
+ r = fma(r, x2, RV(+81.601622084783653475948476293));
+ r = fma(r, x2, RV(-41.3416645219970165527584702233));
+ r = fma(r, x2, RV(+6.2831851984651420764055506572));
+ break;
+ case 8:
+ // double, error=5.7727879332805926019333251099475986171199401686760527068393e-19
+ x*=RV(4.0);
+ x2*=RV(16.0);
+ r = RV(+5.86762142331525639276611151945042252894600631928465887664953969158e-12);
+ r = fma(r, x2, RV(-6.6841794237726945077103730471119313610778831290978967922327437499159e-10));
+ r = fma(r, x2, RV(+5.69213209691935763964803333647156625719770522961368e-8));
+ r = fma(r, x2, RV(-3.598842978568308467678745880795603458651065024316994e-6));
+ r = fma(r, x2, RV(+0.000160441184690020871131343239138998965356942754919962001));
+ r = fma(r, x2, RV(-0.0046817541352970526059161105189653385731022565875740295038));
+ r = fma(r, x2, RV(+0.0796926262461644478160657529010052899007494340748685273131));
+ r = fma(r, x2, RV(-0.6459640975062461124233071963502627016003287601029162984167));
+ r = fma(r, x2, RV(+1.5707963267948966169881546168049607518586729643258729850326));
+ break;
+ default:
+ __builtin_unreachable();
}
+ r *= x;
// Undo range reduction
r = copysign(r, sign);
@@ -65,178 +77,9 @@ namespace vecmathlib {
template<typename realvec_t>
realvec_t mathfuncs<realvec_t>::vml_tan(realvec_t x)
{
- // Reduce range: tan(x) = tan(x + pi)
- x = remainder(x, RV(M_PI));
- assert(all(x >= RV(-M_PI_2) && x <= RV(M_PI_2)));
-
- // Reduce range: tan(x) = -tan(-x)
- realvec_t sign = x;
- x = fabs(x);
- assert(all(x >= RV(0.0) && x <= RV(M_PI_2)));
-
- // Reduce range: cot(x) = -tan(x + pi/2)
- boolvec_t use_cot = x > RV(0.5 * M_PI_2);
- x = ifthen(use_cot, RV(M_PI_2) -x, x);
- assert(all(x >= RV(0.0) && x <= RV(0.5 * M_PI_2)));
-
- // Calculate tan
- realvec_t r = sin(x) / cos(x);
-
- // Undo range reduction
- r = ifthen(use_cot, -rcp(r), r);
-
- // Undo range reducion
- r = copysign(r, sign);
-
- return r;
- }
-
- template<typename realvec_t>
- realvec_t mathfuncs<realvec_t>::vml_sin_chebyshev_single(realvec_t x)
- {
- typename realvec_t::scalar_t const eps __attribute__((__unused__)) = 0.000001;
- realvec_t r, s, t;
-
- // Reduce range: sin(x) = sin(x + 2PI)
- t = x + RV(M_PI);
- t = t * RV(1.0 / (2.0 * M_PI));
- t = floor(t);
- t = x - (t * RV(2.0 * M_PI));
- assert(all(t >= RV(-(1.0 + eps) * M_PI) && t <= RV((1.0 + eps) * M_PI)));
-
- // Reduce range: -sin(x) = -sin(-x)
- s = copysign(RV(1.0), t);
- t = fabs(t);
-
- // Reduce range: sin(x) = sin(pi-x)
- x = ifthen(t > RV(M_PI_2), RV(M_PI) - t, t);
- assert(all(x >= RV(0.0 - eps) && x <= RV(M_PI_2)));
-
- // Evaluate Chebyshev polynomial expansion r = s*T(x)
- // where T(x) = x + c1*x^3 + c2*x^5 + c3*x^7 + ...
- t = x * x;
- x = x * s;
- s = t * x;
- r = fma(RV(-2.2636462059494963e-09), t, RV(1.9761453001926614e-08));
- r = fma(r, t, RV(-9.1621870444937360e-08));
- r = fma(r, t, RV(2.8673509297532770e-06));
- r = fma(r, t, RV(-1.9850777498822076e-04));
- r = fma(r, t, RV(8.3333708144611930e-03));
- r = fma(r, t, RV(-1.6666667163372040e-01));
- r = fma(r, s, x);
- return r;
- }
-
- template<typename realvec_t>
- realvec_t mathfuncs<realvec_t>::vml_sin_chebyshev_double(realvec_t x)
- {
- typename realvec_t::scalar_t const eps __attribute__((__unused__)) = 0.000001;
- realvec_t r, s, t;
-
- // Reduce range: sin(x) = sin(x + 2PI)
- t = x + RV(M_PI);
- t = t * RV(1.0 / (2.0 * M_PI));
- t = floor(t);
- t = x - (t * RV(2.0 * M_PI));
- assert(all(t >= RV(-(1.0 + eps) * M_PI) && t <= RV((1.0 + eps) * M_PI)));
-
- // Reduce range: -sin(x) = -sin(-x)
- s = copysign(RV(1.0), t);
- t = fabs(t);
-
- // Reduce range: sin(x) = sin(pi-x)
- x = ifthen(t > RV(M_PI_2), RV(M_PI) - t, t);
- assert(all(x >= RV(0.0 - eps) && x <= RV(M_PI_2)));
-
- // Evaluate Chebyshev polynomial expansion r = s*T(x)
- // where T(x) = x + c1*x^3 + c2*x^5 + c3*x^7 + ...
- t = x * x;
- x = x * s;
- s = t * x;
- r = fma(RV(1.1888871779171205e-23), t, RV(-1.6213346583948200e-21));
- r = fma(r, t, RV(9.4674830124704450e-19));
- r = fma(r, t, RV(-3.6586864533554100e-16));
- r = fma(r, t, RV(-7.5815669263036780e-13));
- r = fma(r, t, RV(1.6058175109947732e-10));
- r = fma(r, t, RV(-2.5052101017560582e-08));
- r = fma(r, t, RV(2.7557319185249400e-06));
- r = fma(r, t, RV(-1.9841269841152493e-04));
- r = fma(r, t, RV(8.3333333333331560e-03));
- r = fma(r, t, RV(-1.6666666666666666e-01));
- r = fma(r, s, x);
- return r;
- }
-
- template<typename realvec_t>
- realvec_t mathfuncs<realvec_t>::vml_cos_chebyshev_single(realvec_t x)
- {
- typename realvec_t::boolvec_t m;
- realvec_t r, s, t;
-
- // Reduce range: cos(x) = cos(x + 2PI)
- t = x + RV(M_PI);
- t = t * RV(1.0 / (2.0 * M_PI));
- t = floor(t);
- t = x - (t * RV(2.0 * M_PI));
-
- // Reduce range: cos(x) = sin(PI/2 - x)
- t = fabs(t);
- m = t > RV(M_PI_2);
- x = RV(M_PI) - t;
- s = ifthen(m, RV(-1.0), RV(1.0));
- x = ifthen(m, x, t);
-
- // Evaluate Chebyshev polynomial s*T(x) where
- // T(x) = 1 + c1*x^2 + c2*x^4 + c3*x^6 + ...
- t = x * x;
- x = t * s;
- r = fma(RV(-1.0986536764839979e-11), t, RV(2.0856702467901100e-09));
- r = fma(r, t, RV(-2.7556891974788950e-07));
- r = fma(r, t, RV(2.4801582421938170e-05));
- r = fma(r, t, RV(-1.3888888860984269e-03));
- r = fma(r, t, RV(4.1666666666056330e-02));
- r = fma(r, t, RV(-5.0000000000000000e-01));
- r = fma(r, x, s);
- return r;
- }
-
- template<typename realvec_t>
- realvec_t mathfuncs<realvec_t>::vml_cos_chebyshev_double(realvec_t x)
- {
- typename realvec_t::boolvec_t m;
- realvec_t r, s, t;
-
- // Reduce range: cos(x) = cos(x + 2PI)
- t = x + RV(M_PI);
- t = t * RV(1.0 / (2.0 * M_PI));
- t = floor(t);
- t = x - (t * RV(2.0 * M_PI));
-
- // Reduce range: cos(x) = sin(PI/2 - x)
- t = fabs(t);
- m = t > RV(M_PI_2);
- x = RV(M_PI) - t;
- s = ifthen(m, RV(-1.0), RV(1.0));
- x = ifthen(m, x, t);
-
- // Evaluate Chebyshev polynomial s*T(x) where
- // T(x) = 1 + c1*x^2 + c2*x^4 + c3*x^6 + ...
- t = x * x;
- x = t * s;
- r = fma(RV(-8.6512994843471700e-22), t, RV(4.1086675770914360e-19));
- r = fma(r, t, RV(-1.5619143199049570e-16));
- r = fma(r, t, RV(4.7794771764282040e-14));
- r = fma(r, t, RV(-1.1470745595224050e-11));
- r = fma(r, t, RV(2.0876756987841530e-09));
- r = fma(r, t, RV(-2.7557319223985710e-07));
- r = fma(r, t, RV(2.4801587301587300e-05));
- r = fma(r, t, RV(-1.3888888888888890e-03));
- r = fma(r, t, RV(4.1666666666666664e-02));
- r = fma(r, t, RV(-5.0000000000000000e-01));
- r = fma(r, x, s);
- return r;
+ return sin(x) / cos(x);
}
-
+
}; // namespace vecmathlib
#endif // #ifndef MATHFUNCS_SIN_H
OpenPOWER on IntegriCloud