Skip to content

add function sncs1cs to get accurate 1-cosine #44

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
nsajko wants to merge 1 commit into
base: master
Choose a base branch
from
Open

add function sncs1cs to get accurate 1-cosine #44

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
152 changes: 152 additions & 0 deletions mathc.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,140 @@ the following restrictions:

#include "mathc.h"

#if !defined(MATHC_FLOATING_POINT_TYPE) && defined(MATHC_USE_FLOATING_POINT)
/* The following code (sincos1cos and float and double versions of
* sncs1cs) is Copyright © 1985, 1995, 2000 Stephen L. Moshier and
* Copyright © 2020 Neven Sajko. The intention is to get accurate
* 1-cosine, while also getting the sine and cosine as a bonus. The
* implementation is derived from the Cephes Math Library's sin.c and
* sinf.c. To be more specific, I took Stephen Moshier's sin, cos, sinf
* and cosf (without changing the polynomials) and adapted them to give
* all three required function values (in double and float versions),
* without unnecessary accuracy losses.
*
* sncs1cs is not correct for values of x of huge magnitude. That can
* be fixed by more elaborate range reduction.
*/

typedef struct {
/* Sine, cosine, 1-cosine */
mfloat_t sin, cos, omc;
} sincos1cos;

static
sincos1cos sncs1cs(mfloat_t x)
{
const mfloat_t fourOverPi = 1.27323954473516268615;

mfloat_t y, z, zz;
mint_t j, sign = 1, csign = 1;
sincos1cos r;

/* Handle +-0. */
if (x == (mfloat_t)0) {
r.sin = x;
r.cos = 1;
r.omc = 0;
return r;
}
if (isnan(x)) {
r.sin = r.cos = r.omc = x;
return r;
}
if (isinf(x)) {
r.sin = r.cos = r.omc = x - x;
return r;
}
if (x < 0) {
sign = -1;
x = -x;
}
j = (mint_t)(x * fourOverPi);
y = (mfloat_t)j;
/* map zeros to origin */
if ((j & 1)) {
j += 1;
y += 1;
}
j = j & 7; /* octant modulo one turn */
/* reflect in x axis */
if (j > 3) {
sign = -sign;
csign = -csign;
j -= 4;
}
if (j > 1) {
csign = -csign;
}

# if defined(MATHC_USE_SINGLE_FLOATING_POINT)
const float sc[] = {-1.9515295891E-4, 8.3321608736E-3, -1.6666654611E-1};

const float cc[] = {2.443315711809948E-5, -1.388731625493765E-3, 4.166664568298827E-2};

/* These are for a 24-bit significand: */
const float DP1 = 0.78515625;
const float DP2 = 2.4187564849853515625e-4;
const float DP3 = 3.77489497744594108e-8;

/* Extended precision modular arithmetic */
z = ((x - y * DP1) - y * DP2) - y * DP3;
zz = z * z;
r.sin = z + zz*z*((sc[0]*zz + sc[1])*zz + sc[2]);
r.omc = (mfloat_t)0.5*zz - zz*zz*((cc[0]*zz + cc[1])*zz + cc[2]);
# elif defined(MATHC_USE_DOUBLE_FLOATING_POINT)
const double sc[] = {
1.58962301576546568060E-10,
-2.50507477628578072866E-8,
2.75573136213857245213E-6,
-1.98412698295895385996E-4,
8.33333333332211858878E-3,
-1.66666666666666307295E-1,
};

const double cc[] = {
-1.13585365213876817300E-11,
2.08757008419747316778E-9,
-2.75573141792967388112E-7,
2.48015872888517045348E-5,
-1.38888888888730564116E-3,
4.16666666666665929218E-2,
};

const double DP1 = 7.85398125648498535156E-1;
const double DP2 = 3.77489470793079817668E-8;
const double DP3 = 2.69515142907905952645E-15;

/* Extended precision modular arithmetic */
z = ((x - y * DP1) - y * DP2) - y * DP3;
zz = z * z;
r.sin = z + zz*z*(((((sc[0]*zz + sc[1])*zz + sc[2])*zz + sc[3])*zz + sc[4])*zz + sc[5]);
r.omc = (mfloat_t)0.5*zz - zz*zz*(((((cc[0]*zz + cc[1])*zz + cc[2])*zz + cc[3])*zz + cc[4])*zz + cc[5]);
# else
# error "Unknown preprocessing configuration."
# endif
if (j == 1 || j == 2) {
if (csign < 0) {
r.sin = -r.sin;
}
r.cos = r.sin;
r.sin = 1 - r.omc;
r.omc = 1 - r.cos;
} else {
if (csign < 0) {
r.cos = r.omc - 1;
r.omc = 1 - r.cos;
} else {
r.cos = 1 - r.omc;
}
}
if (sign < 0) {
r.sin = -r.sin;
}
return r;
}
#endif

#if defined(MATHC_USE_INT)
mint_t clampi(mint_t value, mint_t min, mint_t max)
{
Expand Down Expand Up @@ -2520,9 +2654,15 @@ mfloat_t *mat3_rotation_z(mfloat_t *result, mfloat_t f)

mfloat_t *mat3_rotation_axis(mfloat_t *result, mfloat_t *v0, mfloat_t f)
{
#if !defined(MATHC_FLOATING_POINT_TYPE)
// Get higher accuracy 1-cosine when mfloat_t is known.
sincos1cos sc1c = sncs1cs(f);
mfloat_t c = sc1c.cos, s = sc1c.sin, one_c = sc1c.omc;
#else
mfloat_t c = MCOS(f);
mfloat_t s = MSIN(f);
mfloat_t one_c = MFLOAT_C(1.0) - c;
#endif
mfloat_t x = v0[0];
mfloat_t y = v0[4];
mfloat_t z = v0[8];
Expand Down Expand Up @@ -2977,9 +3117,15 @@ mfloat_t *mat4_rotation_z(mfloat_t *result, mfloat_t f)

mfloat_t *mat4_rotation_axis(mfloat_t *result, mfloat_t *v0, mfloat_t f)
{
#if !defined(MATHC_FLOATING_POINT_TYPE)
// Get higher accuracy 1-cosine when mfloat_t is known.
sincos1cos sc1c = sncs1cs(f);
mfloat_t c = sc1c.cos, s = sc1c.sin, one_c = sc1c.omc;
#else
mfloat_t c = MCOS(f);
mfloat_t s = MSIN(f);
mfloat_t one_c = MFLOAT_C(1.0) - c;
#endif
mfloat_t x = v0[0];
mfloat_t y = v0[1];
mfloat_t z = v0[2];
Expand Down Expand Up @@ -6753,7 +6899,13 @@ mfloat_t sine_ease_in(mfloat_t f)

mfloat_t sine_ease_in_out(mfloat_t f)
{
#if !defined(MATHC_FLOATING_POINT_TYPE)
// Get higher accuracy 1-cosine when mfloat_t is known.
sincos1cos sc1c = sncs1cs(f * MPI);
return MFLOAT_C(0.5) * sc1c.omc;
#else
return MFLOAT_C(0.5) * (MFLOAT_C(1.0) - MCOS(f * MPI));
#endif
}

mfloat_t circular_ease_out(mfloat_t f)
Expand Down