shishantbiswas/c3c
1
1
Fork 0
mirror of https://github.com/c3lang/c3c.git synced 2026-02-27 12:01:16 +00:00
Code Issues Packages Projects Releases Wiki Activity

Adding sincos / libc tan/tanf.

Browse Source
This commit is contained in:
Christoffer Lerno
2023-02-02 19:29:29 +01:00
parent 0f4d20f168
commit 3b3773663a
5 changed files with 320 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

83
lib/std/math/math_nolibc/__tan.c3 Normal file
View File

@@ -0,0 +1,83 @@
module std::math::nolibc;
$if (!env::COMPILER_LIBC_AVAILABLE):
/* origin: FreeBSD /usr/src/lib/msun/src/k_tan.c */
/*
* ====================================================
* Copyright 2004 Sun Microsystems, Inc. All Rights Reserved.
*
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
const double[*] TAN_T = {
3.33333333333334091986e-01, /* 3FD55555, 55555563 */
1.33333333333201242699e-01, /* 3FC11111, 1110FE7A */
5.39682539762260521377e-02, /* 3FABA1BA, 1BB341FE */
2.18694882948595424599e-02, /* 3F9664F4, 8406D637 */
8.86323982359930005737e-03, /* 3F8226E3, E96E8493 */
3.59207910759131235356e-03, /* 3F6D6D22, C9560328 */
1.45620945432529025516e-03, /* 3F57DBC8, FEE08315 */
5.88041240820264096874e-04, /* 3F4344D8, F2F26501 */
2.46463134818469906812e-04, /* 3F3026F7, 1A8D1068 */
7.81794442939557092300e-05, /* 3F147E88, A03792A6 */
7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */
-1.85586374855275456654e-05, /* BEF375CB, DB605373 */
2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */
};
fn double __tan(double x, double y, int odd) @extern("__tan") @weak
{
const double PIO4 = 7.85398163397448278999e-01; /* 3FE921FB, 54442D18 */
const double PIO4LO = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */
uint hx = (uint)(bitcast(x, ulong) >> 32);
bool big = (hx &0x7fffffff) >= 0x3FE59428; // |x| >= 0.6744
int sign = void;
if (big)
{
sign = hx >> 31;
if (sign)
{
x = -x;
y = -y;
}
x = (PIO4 - x) + (PIO4LO - y);
y = 0.0;
}
double z = x * x;
double w = z * z;
/*
* Break x^5*(T[1]+x^2*T[2]+...) into
* x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
* x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
*/
double r = TAN_T[1] + w * (TAN_T[3] + w * (TAN_T[5] + w * (TAN_T[7] + w * (TAN_T[9] + w * TAN_T[11]))));
double v = z * (TAN_T[2] + w * (TAN_T[4] + w * (TAN_T[6] + w * (TAN_T[8] + w * (TAN_T[10] + w * TAN_T[12])))));
double s = z * x;
r = y + z * (s * (r + v) + y) + s * TAN_T[0];
w = x + r;
if (big)
{
s = 1 - 2 * odd;
v = s - 2.0 * (x + (r - w*w/(w + s)));
return sign ? -v : v;
}
if (!odd) return w;
// -1.0/(x+r) has up to 2ulp error, so compute it accurately
// Clear low word
ulong d = bitcast(w, ulong) & 0xFFFF_FFFF_0000_0000;
double w0 = bitcast(d, double);
v = r - (w0 - x); // w0+v = r+x
double a = -1.0 / w;
d = bitcast(a, ulong) & 0xFFFF_FFFF_0000_0000;
double a0 = bitcast(d, double);
return a0 + a * (1.0 + a0 * w0 + a0 * v);
}
$endif;

56
lib/std/math/math_nolibc/__tandf.c3 Normal file
View File

@@ -0,0 +1,56 @@
module std::math::nolibc;
$if (!env::COMPILER_LIBC_AVAILABLE):
/* origin: FreeBSD /usr/src/lib/msun/src/k_tanf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
* Optimized by Bruce D. Evans.
*/
/*
* ====================================================
* Copyright 2004 Sun Microsystems, Inc. All Rights Reserved.
*
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
// |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]).
const double[*] TANDF = {
0x15554d3418c99f.0p-54, /* 0.333331395030791399758 */
0x1112fd38999f72.0p-55, /* 0.133392002712976742718 */
0x1b54c91d865afe.0p-57, /* 0.0533812378445670393523 */
0x191df3908c33ce.0p-58, /* 0.0245283181166547278873 */
0x185dadfcecf44e.0p-61, /* 0.00297435743359967304927 */
0x1362b9bf971bcd.0p-59, /* 0.00946564784943673166728 */
};
fn float __tandf(double x, int odd) @extern("__tandf") @weak
{
double z = x * x;
/*
* Split up the polynomial into small independent terms to give
* opportunities for parallel evaluation. The chosen splitting is
* micro-optimized for Athlons (XP, X64). It costs 2 multiplications
* relative to Horner's method on sequential machines.
*
* We add the small terms from lowest degree up for efficiency on
* non-sequential machines (the lowest degree terms tend to be ready
* earlier). Apart from this, we don't care about order of
* operations, and don't need to to care since we have precision to
* spare. However, the chosen splitting is good for accuracy too,
* and would give results as accurate as Horner's method if the
* small terms were added from highest degree down.
*/
double r = TANDF[4] + z * TANDF[5];
double t = TANDF[2] + z * TANDF[3];
double w = z * z;
double s = z * x;
double u = TANDF[0] + z * TANDF[1];
r = (x + s * u) + (s * w) * (t + w * r);
return (float)(odd ? -1.0 / r : r);
}
$endif;

51
lib/std/math/math_nolibc/sincosf.c3 → lib/std/math/math_nolibc/sincos.c3
View File

@@ -78,11 +78,9 @@ fn void sincosf(float x, float *sin, float *cos) @extern("sincosf") @weak
}
*sin = __sindf(sign ? x + S4PI2 : x - S4PI2);
*cos = __cosdf(sign ? x + S4PI2 : x - S4PI2);
return;
case ix >= 0x7f800000:
// sin(Inf or NaN) is NaN
*sin = *cos = x - x;
return;
default:
// general argument reduction needed
double y = void;
@@ -109,4 +107,53 @@ fn void sincosf(float x, float *sin, float *cos) @extern("sincosf") @weak
}
fn void sincos(double x, double *sin, double *cos) @extern("sincos") @weak
{
// High word of x.
uint ix = (uint)(bitcast(x, ulong) >> 32);
ix &= 0x7fffffff;
// |x| ~< pi/4
switch
{
// |x| ~< pi/4
case ix <= 0x3fe921fb:
if (ix < 0x3e46a09e)
{
// raise inexact if x!=0 and underflow if subnormal
force_eval_add(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f, 0);
*sin = x;
*cos = 1.0;
return;
}
// if |x| < 2**-27 * sqrt(2)
*sin = __sin(x, 0.0, 0);
*cos = __cos(x, 0.0);
// sincos(Inf or NaN) is NaN
case ix >= 0x7ff00000:
*sin = *cos = x - x;
default:
// argument reduction needed
double[2] y;
int n = __rem_pio2(x, &y);
double s = __sin(y[0], y[1], 1);
double c = __cos(y[0], y[1]);
switch (n & 3)
{
case 0:
*sin = s;
*cos = c;
case 1:
*sin = c;
*cos = -s;
case 2:
*sin = -s;
*cos = -c;
case 3:
default:
*sin = -c;
*cos = s;
}
}
}
$endif;

103
lib/std/math/math_nolibc/tan.c3 Normal file
View File

@@ -0,0 +1,103 @@
module std::math::nolibc;
$if (!env::COMPILER_LIBC_AVAILABLE):
/* origin: FreeBSD /usr/src/lib/msun/src/s_tan.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
fn double tan(double x) @extern("tan") @weak
{
uint ix = (uint)(bitcast(x, ulong) >> 32);
ix &= 0x7fffffff;
// |x| ~< pi/4
switch
{
case ix <= 0x3fe921fb:
if (ix < 0x3e400000)
{
// |x| < 2**-27 */
/* raise inexact if x!=0 and underflow if subnormal */
force_eval_add(ix < 0x00100000 ? x / 0x1p120f : x + 0x1p120f, 0);
return x;
}
return __tan(x, 0.0, 0);
// tan(Inf or NaN) is NaN
case ix >= 0x7ff00000:
return x - x;
default:
// argument reduction
double[2] y;
uint n = __rem_pio2(x, &y);
return __tan(y[0], y[1], n & 1);
}
}
/* origin: FreeBSD /usr/src/lib/msun/src/s_tanf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
* Optimized by Bruce D. Evans.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
fn float tanf(float x) @extern("tanf") @weak
{
uint ix = bitcast(x, uint);
uint sign = ix >> 31;
ix &= 0x7fffffff;
switch
{
// |x| ~<= pi/4
case ix <= 0x3f490fda:
// |x| < 2**-12
if (ix < 0x39800000)
{
/* raise inexact if x!=0 and underflow if subnormal */
force_eval_add(ix < 0x00800000 ? x / 0x1p120f : x + 0x1p120f, 0);
return x;
}
return __tandf(x, 0);
// |x| ~<= 5*pi/4
case ix <= 0x407b53d1:
// |x| ~<= 3pi/4
return ix <= 0x4016cbe3
? __tandf(sign ? x + S1PI2 : x - S1PI2, 1)
: __tandf(sign ? x + S2PI2 : x - S2PI2, 0);
// |x| ~<= 9*pi/4
case ix <= 0x40e231d5:
// |x| ~<= 7*pi/4 */
return ix <= 0x40afeddf
? __tandf(sign ? x + S3PI2 : x - S3PI2, 1)
: __tandf(sign ? x + S4PI2 : x - S4PI2, 0);
// tan(Inf or NaN) is NaN
case ix >= 0x7f800000:
return x - x;
default:
// argument reduction
double y;
uint n = __rem_pio2f(x, &y);
return __tandf(y, n & 1);
}
}
$endif;