diff --git a/lib/std/math/math.c3 b/lib/std/math/math.c3 index 173d7be13..4c3e1ab29 100644 --- a/lib/std/math/math.c3 +++ b/lib/std/math/math.c3 @@ -134,6 +134,20 @@ macro atan2(x, y) $endif; } +/** + * @require values::@is_int(x) || values::@is_float(x) "Expected an integer or floating point value" + * @checked (*y)[0] = x, y.len + * @require y.len == 2 + **/ +macro @sincos(x, y) +{ + $if ($typeof(y[0]).typeid == float.typeid): + return _sincosf(x, y); + $else: + return _sincos(x, y); + $endif; +} + /** * @require values::@is_int(x) || values::@is_float(x) "Expected an integer or floating point value" * @checked x + y @@ -355,7 +369,17 @@ macro sqrt(x) => $$sqrt(values::promote_int(x)); /** * @require values::@is_promotable_to_floatlike(x) `The input must be a number or a float vector` **/ -macro tan(x) => sin(x) / cos(x); +macro tan(x) +{ + var $Type = $typeof(x); + $if (types.is_vector($Type)): + return $$sin(x) / $$cos(x); + $elif ($Type.typeid == float.typeid): + return _tanf(x); + $else: + return tan(x); + $endif; +} /** * @require values::@is_promotable_to_floatlike(x) `The input must be a number or a float vector` @@ -615,3 +639,7 @@ extern fn double _atan(double x) @extname("atan"); extern fn float _atanf(float x) @extname("atanf"); extern fn double _atan2(double, double) @extname("atan2"); extern fn float _atan2f(float, float) @extname("atan2f"); +extern fn void _sincos(double, double*) @extname("sincos"); +extern fn void _sincosf(float, float*) @extname("sincosf"); +extern fn void _tan(double x, double y) @extern("tan"); +extern fn void _tanf(double x, double y) @extern("tanf"); \ No newline at end of file diff --git a/lib/std/math/math_nolibc/__tan.c3 b/lib/std/math/math_nolibc/__tan.c3 new file mode 100644 index 000000000..c6b764000 --- /dev/null +++ b/lib/std/math/math_nolibc/__tan.c3 @@ -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; \ No newline at end of file diff --git a/lib/std/math/math_nolibc/__tandf.c3 b/lib/std/math/math_nolibc/__tandf.c3 new file mode 100644 index 000000000..ea87625b5 --- /dev/null +++ b/lib/std/math/math_nolibc/__tandf.c3 @@ -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; \ No newline at end of file diff --git a/lib/std/math/math_nolibc/sincosf.c3 b/lib/std/math/math_nolibc/sincos.c3 similarity index 72% rename from lib/std/math/math_nolibc/sincosf.c3 rename to lib/std/math/math_nolibc/sincos.c3 index 11f7b839c..f819fd0ff 100644 --- a/lib/std/math/math_nolibc/sincosf.c3 +++ b/lib/std/math/math_nolibc/sincos.c3 @@ -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; \ No newline at end of file diff --git a/lib/std/math/math_nolibc/tan.c3 b/lib/std/math/math_nolibc/tan.c3 new file mode 100644 index 000000000..95fb73ce2 --- /dev/null +++ b/lib/std/math/math_nolibc/tan.c3 @@ -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; \ No newline at end of file