From 94c211f5bed1c18fe73a9ab3a4f14a76042f2fa3 Mon Sep 17 00:00:00 2001 From: ZERICO2005 <71151164+ZERICO2005@users.noreply.github.com> Date: Sat, 1 Mar 2025 17:16:53 -0700 Subject: [PATCH 1/2] Renamed math.h functions 'funcf.*' is float, 'funcl.*' is long double, 'func.*' is float and long double --- src/libc/{acos.src => acosf.src} | 0 src/libc/{acosh.c => acoshf.c} | 5 -- src/libc/acoshl.c | 6 +++ src/libc/{asin.c => asinf.c} | 0 src/libc/{asin.src => asinf.src} | 0 src/libc/{asinh.c => asinhf.c} | 5 -- src/libc/asinhl.c | 6 +++ src/libc/{atan2.c => atan2f.c} | 0 src/libc/{atan2.src => atan2f.src} | 0 src/libc/{atan.c => atanf.c} | 0 src/libc/{atan.src => atanf.src} | 0 src/libc/{atanh.c => atanhf.c} | 5 -- src/libc/atanhl.c | 6 +++ src/libc/{cbrt.c => cbrtf.c} | 0 src/libc/{cos.src => cosf.src} | 0 src/libc/{cosh.c => coshf.c} | 0 src/libc/{cosh.src => coshf.src} | 0 src/libc/erfcf.c | 10 ++++ src/libc/{erfc.c => erfcl.c} | 9 ---- src/libc/{erf.c => erff.c} | 5 -- src/libc/erfl.c | 6 +++ src/libc/{exp2.c => exp2f.c} | 0 src/libc/{exp.c => expf.c} | 0 src/libc/{exp.src => expf.src} | 0 src/libc/{expm1.c => expm1f.c} | 0 src/libc/{lgamma.c => lgammaf.c} | 0 src/libc/{llrint.c => llrintf.c} | 0 src/libc/{llround.c => llroundf.c} | 0 src/libc/{log10.c => log10f.c} | 0 src/libc/{log10.src => log10f.src} | 0 src/libc/{log1p.c => log1pf.c} | 0 src/libc/{log2.c => log2f.c} | 0 src/libc/{log.c => logf.c} | 0 src/libc/{log.src => logf.src} | 0 src/libc/{lrint.c => lrintf.c} | 0 src/libc/{lround.c => lroundf.c} | 0 src/libc/nextafterf.c | 58 ++++++++++++++++++++++ src/libc/{nextafter.c => nextafterl.c} | 57 --------------------- src/libc/{nextupdown.c => nextupdownf.c} | 62 ----------------------- src/libc/nextupdownl.c | 63 ++++++++++++++++++++++++ src/libc/{pow.c => powf.c} | 0 src/libc/{pow.src => powf.src} | 0 src/libc/{remainder.c => remainderf.c} | 0 src/libc/{remquo.c => remquof.c} | 0 src/libc/{rint.c => rintf.c} | 0 src/libc/{sin.c => sinf.c} | 0 src/libc/{sin.src => sinf.src} | 0 src/libc/{sinh.c => sinhf.c} | 0 src/libc/{sinh.src => sinhf.src} | 0 src/libc/{tan.c => tanf.c} | 0 src/libc/{tan.src => tanf.src} | 0 src/libc/{tanh.c => tanhf.c} | 0 src/libc/{tanh.src => tanhf.src} | 0 src/libc/{tgamma.c => tgammaf.c} | 0 54 files changed, 155 insertions(+), 148 deletions(-) rename src/libc/{acos.src => acosf.src} (100%) rename src/libc/{acosh.c => acoshf.c} (64%) create mode 100644 src/libc/acoshl.c rename src/libc/{asin.c => asinf.c} (100%) rename src/libc/{asin.src => asinf.src} (100%) rename src/libc/{asinh.c => asinhf.c} (64%) create mode 100644 src/libc/asinhl.c rename src/libc/{atan2.c => atan2f.c} (100%) rename src/libc/{atan2.src => atan2f.src} (100%) rename src/libc/{atan.c => atanf.c} (100%) rename src/libc/{atan.src => atanf.src} (100%) rename src/libc/{atanh.c => atanhf.c} (63%) create mode 100644 src/libc/atanhl.c rename src/libc/{cbrt.c => cbrtf.c} (100%) rename src/libc/{cos.src => cosf.src} (100%) rename src/libc/{cosh.c => coshf.c} (100%) rename src/libc/{cosh.src => coshf.src} (100%) create mode 100644 src/libc/erfcf.c rename src/libc/{erfc.c => erfcl.c} (54%) rename src/libc/{erf.c => erff.c} (65%) create mode 100644 src/libc/erfl.c rename src/libc/{exp2.c => exp2f.c} (100%) rename src/libc/{exp.c => expf.c} (100%) rename src/libc/{exp.src => expf.src} (100%) rename src/libc/{expm1.c => expm1f.c} (100%) rename src/libc/{lgamma.c => lgammaf.c} (100%) rename src/libc/{llrint.c => llrintf.c} (100%) rename src/libc/{llround.c => llroundf.c} (100%) rename src/libc/{log10.c => log10f.c} (100%) rename src/libc/{log10.src => log10f.src} (100%) rename src/libc/{log1p.c => log1pf.c} (100%) rename src/libc/{log2.c => log2f.c} (100%) rename src/libc/{log.c => logf.c} (100%) rename src/libc/{log.src => logf.src} (100%) rename src/libc/{lrint.c => lrintf.c} (100%) rename src/libc/{lround.c => lroundf.c} (100%) create mode 100644 src/libc/nextafterf.c rename src/libc/{nextafter.c => nextafterl.c} (52%) rename src/libc/{nextupdown.c => nextupdownf.c} (52%) create mode 100644 src/libc/nextupdownl.c rename src/libc/{pow.c => powf.c} (100%) rename src/libc/{pow.src => powf.src} (100%) rename src/libc/{remainder.c => remainderf.c} (100%) rename src/libc/{remquo.c => remquof.c} (100%) rename src/libc/{rint.c => rintf.c} (100%) rename src/libc/{sin.c => sinf.c} (100%) rename src/libc/{sin.src => sinf.src} (100%) rename src/libc/{sinh.c => sinhf.c} (100%) rename src/libc/{sinh.src => sinhf.src} (100%) rename src/libc/{tan.c => tanf.c} (100%) rename src/libc/{tan.src => tanf.src} (100%) rename src/libc/{tanh.c => tanhf.c} (100%) rename src/libc/{tanh.src => tanhf.src} (100%) rename src/libc/{tgamma.c => tgammaf.c} (100%) diff --git a/src/libc/acos.src b/src/libc/acosf.src similarity index 100% rename from src/libc/acos.src rename to src/libc/acosf.src diff --git a/src/libc/acosh.c b/src/libc/acoshf.c similarity index 64% rename from src/libc/acosh.c rename to src/libc/acoshf.c index 34ef3a415..b42f45646 100644 --- a/src/libc/acosh.c +++ b/src/libc/acoshf.c @@ -6,8 +6,3 @@ float acoshf(float x) } double acosh(double) __attribute__((alias("acoshf"))); - -long double acoshl(long double x) -{ - return logl(x + sqrtl(x * x - 1)); -} diff --git a/src/libc/acoshl.c b/src/libc/acoshl.c new file mode 100644 index 000000000..d30cdad42 --- /dev/null +++ b/src/libc/acoshl.c @@ -0,0 +1,6 @@ +#include + +long double acoshl(long double x) +{ + return logl(x + sqrtl(x * x - 1)); +} diff --git a/src/libc/asin.c b/src/libc/asinf.c similarity index 100% rename from src/libc/asin.c rename to src/libc/asinf.c diff --git a/src/libc/asin.src b/src/libc/asinf.src similarity index 100% rename from src/libc/asin.src rename to src/libc/asinf.src diff --git a/src/libc/asinh.c b/src/libc/asinhf.c similarity index 64% rename from src/libc/asinh.c rename to src/libc/asinhf.c index deb268d2b..d3c41bde7 100644 --- a/src/libc/asinh.c +++ b/src/libc/asinhf.c @@ -6,8 +6,3 @@ float asinhf(float x) } double asinh(double) __attribute__((alias("asinhf"))); - -long double asinhl(long double x) -{ - return logl(x + sqrtl(x * x + 1)); -} diff --git a/src/libc/asinhl.c b/src/libc/asinhl.c new file mode 100644 index 000000000..603ad0889 --- /dev/null +++ b/src/libc/asinhl.c @@ -0,0 +1,6 @@ +#include + +long double asinhl(long double x) +{ + return logl(x + sqrtl(x * x + 1)); +} diff --git a/src/libc/atan2.c b/src/libc/atan2f.c similarity index 100% rename from src/libc/atan2.c rename to src/libc/atan2f.c diff --git a/src/libc/atan2.src b/src/libc/atan2f.src similarity index 100% rename from src/libc/atan2.src rename to src/libc/atan2f.src diff --git a/src/libc/atan.c b/src/libc/atanf.c similarity index 100% rename from src/libc/atan.c rename to src/libc/atanf.c diff --git a/src/libc/atan.src b/src/libc/atanf.src similarity index 100% rename from src/libc/atan.src rename to src/libc/atanf.src diff --git a/src/libc/atanh.c b/src/libc/atanhf.c similarity index 63% rename from src/libc/atanh.c rename to src/libc/atanhf.c index 8ae497138..c62c075dc 100644 --- a/src/libc/atanh.c +++ b/src/libc/atanhf.c @@ -6,8 +6,3 @@ float atanhf(float x) } double atanh(double) __attribute__((alias("atanhf"))); - -long double atanhl(long double x) -{ - return .5l * logl((1 + x) / (1 - x)); -} diff --git a/src/libc/atanhl.c b/src/libc/atanhl.c new file mode 100644 index 000000000..393201d4b --- /dev/null +++ b/src/libc/atanhl.c @@ -0,0 +1,6 @@ +#include + +long double atanhl(long double x) +{ + return .5l * logl((1 + x) / (1 - x)); +} diff --git a/src/libc/cbrt.c b/src/libc/cbrtf.c similarity index 100% rename from src/libc/cbrt.c rename to src/libc/cbrtf.c diff --git a/src/libc/cos.src b/src/libc/cosf.src similarity index 100% rename from src/libc/cos.src rename to src/libc/cosf.src diff --git a/src/libc/cosh.c b/src/libc/coshf.c similarity index 100% rename from src/libc/cosh.c rename to src/libc/coshf.c diff --git a/src/libc/cosh.src b/src/libc/coshf.src similarity index 100% rename from src/libc/cosh.src rename to src/libc/coshf.src diff --git a/src/libc/erfcf.c b/src/libc/erfcf.c new file mode 100644 index 000000000..fdb427aa1 --- /dev/null +++ b/src/libc/erfcf.c @@ -0,0 +1,10 @@ +#include + +float erfcf(float x) +{ + static const float p = 0.47047f, a1 = 0.3480242f, a2 = -0.0958798f, a3 = 0.7478556f; + const float t = 1.0f / (1.0f + p * fabsf(x)); + return copysignf(t * (a1 + t * (a2 + t * a3)) * expf(-x * x), x); +} + +double erfc(double) __attribute__((alias("erfcf"))); diff --git a/src/libc/erfc.c b/src/libc/erfcl.c similarity index 54% rename from src/libc/erfc.c rename to src/libc/erfcl.c index a6ca51406..bc9d0ea93 100644 --- a/src/libc/erfc.c +++ b/src/libc/erfcl.c @@ -1,14 +1,5 @@ #include -float erfcf(float x) -{ - static const float p = 0.47047f, a1 = 0.3480242f, a2 = -0.0958798f, a3 = 0.7478556f; - const float t = 1.0f / (1.0f + p * fabsf(x)); - return copysignf(t * (a1 + t * (a2 + t * a3)) * expf(-x * x), x); -} - -double erfc(double) __attribute__((alias("erfcf"))); - long double erfcl(long double x) { static const long double p = 0.3275911l, a1 = 0.254829592l, a2 = -0.284496736l, a3 = 1.421413741l, a4 = -1.453152027l, a5 = 1.061405429l; diff --git a/src/libc/erf.c b/src/libc/erff.c similarity index 65% rename from src/libc/erf.c rename to src/libc/erff.c index 74bb7ae33..421688f2b 100644 --- a/src/libc/erf.c +++ b/src/libc/erff.c @@ -6,8 +6,3 @@ float erff(float x) } double erf(double) __attribute__((alias("erff"))); - -long double erfl(long double x) -{ - return 1 - erfcl(x); -} diff --git a/src/libc/erfl.c b/src/libc/erfl.c new file mode 100644 index 000000000..a7027954f --- /dev/null +++ b/src/libc/erfl.c @@ -0,0 +1,6 @@ +#include + +long double erfl(long double x) +{ + return 1 - erfcl(x); +} diff --git a/src/libc/exp2.c b/src/libc/exp2f.c similarity index 100% rename from src/libc/exp2.c rename to src/libc/exp2f.c diff --git a/src/libc/exp.c b/src/libc/expf.c similarity index 100% rename from src/libc/exp.c rename to src/libc/expf.c diff --git a/src/libc/exp.src b/src/libc/expf.src similarity index 100% rename from src/libc/exp.src rename to src/libc/expf.src diff --git a/src/libc/expm1.c b/src/libc/expm1f.c similarity index 100% rename from src/libc/expm1.c rename to src/libc/expm1f.c diff --git a/src/libc/lgamma.c b/src/libc/lgammaf.c similarity index 100% rename from src/libc/lgamma.c rename to src/libc/lgammaf.c diff --git a/src/libc/llrint.c b/src/libc/llrintf.c similarity index 100% rename from src/libc/llrint.c rename to src/libc/llrintf.c diff --git a/src/libc/llround.c b/src/libc/llroundf.c similarity index 100% rename from src/libc/llround.c rename to src/libc/llroundf.c diff --git a/src/libc/log10.c b/src/libc/log10f.c similarity index 100% rename from src/libc/log10.c rename to src/libc/log10f.c diff --git a/src/libc/log10.src b/src/libc/log10f.src similarity index 100% rename from src/libc/log10.src rename to src/libc/log10f.src diff --git a/src/libc/log1p.c b/src/libc/log1pf.c similarity index 100% rename from src/libc/log1p.c rename to src/libc/log1pf.c diff --git a/src/libc/log2.c b/src/libc/log2f.c similarity index 100% rename from src/libc/log2.c rename to src/libc/log2f.c diff --git a/src/libc/log.c b/src/libc/logf.c similarity index 100% rename from src/libc/log.c rename to src/libc/logf.c diff --git a/src/libc/log.src b/src/libc/logf.src similarity index 100% rename from src/libc/log.src rename to src/libc/logf.src diff --git a/src/libc/lrint.c b/src/libc/lrintf.c similarity index 100% rename from src/libc/lrint.c rename to src/libc/lrintf.c diff --git a/src/libc/lround.c b/src/libc/lroundf.c similarity index 100% rename from src/libc/lround.c rename to src/libc/lroundf.c diff --git a/src/libc/nextafterf.c b/src/libc/nextafterf.c new file mode 100644 index 000000000..65f4e5ee8 --- /dev/null +++ b/src/libc/nextafterf.c @@ -0,0 +1,58 @@ +#include +#include +#include + +typedef union F32_pun { + float flt; + uint32_t bin; +} F32_pun; + +#define Float32_pos_subnorm_min UINT32_C(0x00000001) +#define Float32_neg_subnorm_min UINT32_C(0x80000001) + +float nextafterf(float x, float y) { + F32_pun arg_x, arg_y; + arg_x.flt = x; + arg_y.flt = y; + + if (isnan(y)) { + return y; + } + if (isnan(x)) { + return x; + } + if (arg_x.bin == arg_y.bin) { + return y; + } + + if (iszero(x)) { + if (iszero(y)) { + // special case where `+0.0 --> -0.0` and `-0.0 --> +0.0` + return y; + } + feraiseexcept(FE_INEXACT | FE_UNDERFLOW); + F32_pun ret; + ret.bin = signbit(y) ? Float32_neg_subnorm_min : Float32_pos_subnorm_min; + return ret.flt; + } + if (isless(x, y) != signbit(x)) { + // Towards positive/negative infinity + arg_x.bin++; + } else { + // Towards negative/positive zero + arg_x.bin--; + } + if (isnormal(arg_x.flt)) { + return arg_x.flt; + } + if (isinf(arg_x.flt)) { + // overflow to infinity + feraiseexcept(FE_INEXACT | FE_OVERFLOW); + } else { + // result is subnormal or zero + feraiseexcept(FE_INEXACT | FE_UNDERFLOW); + } + return arg_x.flt; +} + +double nextafter(double, double) __attribute__((alias("nextafterf"))); diff --git a/src/libc/nextafter.c b/src/libc/nextafterl.c similarity index 52% rename from src/libc/nextafter.c rename to src/libc/nextafterl.c index 056755aa0..b478c3c3a 100644 --- a/src/libc/nextafter.c +++ b/src/libc/nextafterl.c @@ -2,63 +2,6 @@ #include #include -typedef union F32_pun { - float flt; - uint32_t bin; -} F32_pun; - -#define Float32_pos_subnorm_min UINT32_C(0x00000001) -#define Float32_neg_subnorm_min UINT32_C(0x80000001) - -float nextafterf(float x, float y) { - F32_pun arg_x, arg_y; - arg_x.flt = x; - arg_y.flt = y; - - if (isnan(y)) { - return y; - } - if (isnan(x)) { - return x; - } - if (arg_x.bin == arg_y.bin) { - return y; - } - - if (iszero(x)) { - if (iszero(y)) { - // special case where `+0.0 --> -0.0` and `-0.0 --> +0.0` - return y; - } - feraiseexcept(FE_INEXACT | FE_UNDERFLOW); - F32_pun ret; - ret.bin = signbit(y) ? Float32_neg_subnorm_min : Float32_pos_subnorm_min; - return ret.flt; - } - if (isless(x, y) != signbit(x)) { - // Towards positive/negative infinity - arg_x.bin++; - } else { - // Towards negative/positive zero - arg_x.bin--; - } - if (isnormal(arg_x.flt)) { - return arg_x.flt; - } - if (isinf(arg_x.flt)) { - // overflow to infinity - feraiseexcept(FE_INEXACT | FE_OVERFLOW); - } else { - // result is subnormal or zero - feraiseexcept(FE_INEXACT | FE_UNDERFLOW); - } - return arg_x.flt; -} - -double nextafter(double, double) __attribute__((alias("nextafterf"))); - - - typedef union F64_pun { long double flt; uint64_t bin; diff --git a/src/libc/nextupdown.c b/src/libc/nextupdownf.c similarity index 52% rename from src/libc/nextupdown.c rename to src/libc/nextupdownf.c index f171c1407..dc5d33052 100644 --- a/src/libc/nextupdown.c +++ b/src/libc/nextupdownf.c @@ -64,65 +64,3 @@ float nextdownf(float x) { double nextup(double) __attribute__((alias("nextupf"))); double nextdown(double) __attribute__((alias("nextdownf"))); - - - -typedef union F64_pun { - long double flt; - uint64_t bin; -} F64_pun; - -#define Float64_pos_subnorm_min UINT64_C(0x0000000000000001) -#define Float64_neg_subnorm_min UINT64_C(0x8000000000000001) -#define Float64_pos_inf UINT64_C(0x7FF0000000000000) -#define Float64_neg_inf UINT64_C(0xFFF0000000000000) - -long double nextupl(long double x) { - F64_pun arg_x; - arg_x.flt = x; - - if (isnan(x) || arg_x.bin == Float64_pos_inf) { - // return unmodified - return x; - } - - if (iszero(x)) { - F64_pun ret; - ret.bin = Float64_pos_subnorm_min; - return ret.flt; - } - - if (signbit(x)) { - // Towards negative zero - arg_x.bin--; - } else { - // Towards positive infinity - arg_x.bin++; - } - return arg_x.flt; -} - -long double nextdownl(long double x) { - F64_pun arg_x; - arg_x.flt = x; - - if (isnan(x) || arg_x.bin == Float64_neg_inf) { - // return unmodified - return x; - } - - if (iszero(x)) { - F64_pun ret; - ret.bin = Float64_neg_subnorm_min; - return ret.flt; - } - - if (signbit(x)) { - // Towards negative infinity - arg_x.bin++; - } else { - // Towards positive zero - arg_x.bin--; - } - return arg_x.flt; -} diff --git a/src/libc/nextupdownl.c b/src/libc/nextupdownl.c new file mode 100644 index 000000000..f2b98cce3 --- /dev/null +++ b/src/libc/nextupdownl.c @@ -0,0 +1,63 @@ +#include +#include +#include + +typedef union F64_pun { + long double flt; + uint64_t bin; +} F64_pun; + +#define Float64_pos_subnorm_min UINT64_C(0x0000000000000001) +#define Float64_neg_subnorm_min UINT64_C(0x8000000000000001) +#define Float64_pos_inf UINT64_C(0x7FF0000000000000) +#define Float64_neg_inf UINT64_C(0xFFF0000000000000) + +long double nextupl(long double x) { + F64_pun arg_x; + arg_x.flt = x; + + if (isnan(x) || arg_x.bin == Float64_pos_inf) { + // return unmodified + return x; + } + + if (iszero(x)) { + F64_pun ret; + ret.bin = Float64_pos_subnorm_min; + return ret.flt; + } + + if (signbit(x)) { + // Towards negative zero + arg_x.bin--; + } else { + // Towards positive infinity + arg_x.bin++; + } + return arg_x.flt; +} + +long double nextdownl(long double x) { + F64_pun arg_x; + arg_x.flt = x; + + if (isnan(x) || arg_x.bin == Float64_neg_inf) { + // return unmodified + return x; + } + + if (iszero(x)) { + F64_pun ret; + ret.bin = Float64_neg_subnorm_min; + return ret.flt; + } + + if (signbit(x)) { + // Towards negative infinity + arg_x.bin++; + } else { + // Towards positive zero + arg_x.bin--; + } + return arg_x.flt; +} diff --git a/src/libc/pow.c b/src/libc/powf.c similarity index 100% rename from src/libc/pow.c rename to src/libc/powf.c diff --git a/src/libc/pow.src b/src/libc/powf.src similarity index 100% rename from src/libc/pow.src rename to src/libc/powf.src diff --git a/src/libc/remainder.c b/src/libc/remainderf.c similarity index 100% rename from src/libc/remainder.c rename to src/libc/remainderf.c diff --git a/src/libc/remquo.c b/src/libc/remquof.c similarity index 100% rename from src/libc/remquo.c rename to src/libc/remquof.c diff --git a/src/libc/rint.c b/src/libc/rintf.c similarity index 100% rename from src/libc/rint.c rename to src/libc/rintf.c diff --git a/src/libc/sin.c b/src/libc/sinf.c similarity index 100% rename from src/libc/sin.c rename to src/libc/sinf.c diff --git a/src/libc/sin.src b/src/libc/sinf.src similarity index 100% rename from src/libc/sin.src rename to src/libc/sinf.src diff --git a/src/libc/sinh.c b/src/libc/sinhf.c similarity index 100% rename from src/libc/sinh.c rename to src/libc/sinhf.c diff --git a/src/libc/sinh.src b/src/libc/sinhf.src similarity index 100% rename from src/libc/sinh.src rename to src/libc/sinhf.src diff --git a/src/libc/tan.c b/src/libc/tanf.c similarity index 100% rename from src/libc/tan.c rename to src/libc/tanf.c diff --git a/src/libc/tan.src b/src/libc/tanf.src similarity index 100% rename from src/libc/tan.src rename to src/libc/tanf.src diff --git a/src/libc/tanh.c b/src/libc/tanhf.c similarity index 100% rename from src/libc/tanh.c rename to src/libc/tanhf.c diff --git a/src/libc/tanh.src b/src/libc/tanhf.src similarity index 100% rename from src/libc/tanh.src rename to src/libc/tanhf.src diff --git a/src/libc/tgamma.c b/src/libc/tgammaf.c similarity index 100% rename from src/libc/tgamma.c rename to src/libc/tgammaf.c From d1ef40a07ede63c82a78de5fe176f1ba9c95ef29 Mon Sep 17 00:00:00 2001 From: ZERICO2005 <71151164+ZERICO2005@users.noreply.github.com> Date: Sun, 2 Mar 2025 17:43:02 -0700 Subject: [PATCH 2/2] Fixed formatting and indentation in math.h functions, corrected softfloat linker.mk, fixed lgammaf softlock bug, fixed erfcf/erfcl for negative inputs, and fixed the zilog copyright headers so I can sleep peacefully --- src/libc/asinf.c | 46 ++++---- src/libc/atan2f.c | 38 +++--- src/libc/atanf.c | 66 +++++------ src/libc/atanhl.c | 2 +- src/libc/bsearch.c | 108 ++++++++--------- src/libc/coshf.c | 50 ++++---- src/libc/difftime.c | 2 +- src/libc/erfcf.c | 16 ++- src/libc/erfcl.c | 20 +++- src/libc/expf.c | 70 +++++------ src/libc/floorf.c | 33 +++--- src/libc/frexpf.c | 4 +- src/libc/lgammaf.c | 30 +++-- src/libc/log10f.c | 2 +- src/libc/log1pf.c | 2 +- src/libc/log2f.c | 2 +- src/libc/logf.c | 72 ++++++------ src/libc/nextupdownf.c | 3 +- src/libc/powf.c | 50 ++++---- src/libc/qsort.c | 156 +++++++++++++------------ src/libc/roundeven.c | 2 +- src/libc/signbitf.src | 14 +++ src/libc/{signbit.src => signbitl.src} | 11 +- src/libc/sinf.c | 93 ++++++++------- src/libc/sinhf.c | 78 ++++++------- src/libc/strtof.c | 25 ++-- src/libc/strtold.c | 22 ++-- src/libc/tanf.c | 115 +++++++++--------- src/libc/tanhf.c | 38 +++--- src/libc/tgammaf.c | 8 +- src/linker.mk | 2 +- 31 files changed, 614 insertions(+), 566 deletions(-) create mode 100644 src/libc/signbitf.src rename src/libc/{signbit.src => signbitl.src} (51%) diff --git a/src/libc/asinf.c b/src/libc/asinf.c index 03f647fdb..8a24fe5a8 100644 --- a/src/libc/asinf.c +++ b/src/libc/asinf.c @@ -5,33 +5,33 @@ * arctan is called after appropriate range reduction. */ -#include -#include +#include +#include #define pio2 1.57079632679490 float _asinf_c(float arg) { - float sign, temp; - - sign = 1.; - if(arg < 0) { - arg = -arg; - sign = -1.; - } - - if(arg > 1.) { - errno = EDOM; - return(0.); - } - - temp = sqrt(1. - arg*arg); - if(arg > 0.7) { - temp = pio2 - atan(temp/arg); - } else { - temp = atan(arg/temp); - } - - return(sign*temp); + float sign, temp; + + sign = 1.; + if(arg < 0) { + arg = -arg; + sign = -1.; + } + + if(arg > 1.) { + errno = EDOM; + return(0.); + } + + temp = sqrtf(1. - arg*arg); + if(arg > 0.7) { + temp = pio2 - atanf(temp/arg); + } else { + temp = atanf(arg/temp); + } + + return(sign*temp); } double _asin_c(double) __attribute__((alias("_asinf_c"))); diff --git a/src/libc/atan2f.c b/src/libc/atan2f.c index 046e4eb72..48957946c 100644 --- a/src/libc/atan2f.c +++ b/src/libc/atan2f.c @@ -1,28 +1,28 @@ /* Copyright (c) 2000-2008 Zilog, Inc. */ #include -#define pio2 1.57079632679489e0 +#define pio2 1.57079632679489e0 float _atan2f_c(float arg1,float arg2) { - float satan(float); + float satan(float); - if((arg1+arg2)==arg1) { - if(arg1 >= 0.) { - return(pio2); - } else { - return(-pio2); - } - } else if(arg2 < 0.) { - if(arg1 >= 0.) { - return(pio2+pio2 - satan(-arg1/arg2)); - } else { - return(-pio2-pio2 + satan(arg1/arg2)); - } - } else if(arg1 > 0) { - return(satan(arg1/arg2)); - } else { - return(-satan(-arg1/arg2)); - } + if((arg1+arg2)==arg1) { + if(arg1 >= 0.) { + return(pio2); + } else { + return(-pio2); + } + } else if(arg2 < 0.) { + if(arg1 >= 0.) { + return(pio2+pio2 - satan(-arg1/arg2)); + } else { + return(-pio2-pio2 + satan(arg1/arg2)); + } + } else if(arg1 > 0) { + return(satan(arg1/arg2)); + } else { + return(-satan(-arg1/arg2)); + } } double _atan2_c(double, double) __attribute__((alias("_atan2f_c"))); diff --git a/src/libc/atanf.c b/src/libc/atanf.c index 31eb395ab..152d21b03 100644 --- a/src/libc/atanf.c +++ b/src/libc/atanf.c @@ -15,20 +15,20 @@ */ #include -#define sq2p1 2.41421356237309e0 -#define sq2m1 0.414213562373095e0 -#define pio2 1.57079632679489e0 -#define pio4 0.785398163397448e0 -#define p4 0.161536412982230e2 -#define p3 0.268425481955040e3 -#define p2 0.115302935154049e4 -#define p1 0.178040631643320e4 -#define p0 0.896785974036639e3 -#define q4 0.589569705084446e2 -#define q3 0.536265374031215e3 -#define q2 0.166678381488163e4 -#define q1 0.207933497444541e4 -#define q0 0.896785974036639e3 +#define sq2p1 2.41421356237309e0 +#define sq2m1 0.414213562373095e0 +#define pio2 1.57079632679489e0 +#define pio4 0.785398163397448e0 +#define p4 0.161536412982230e2 +#define p3 0.268425481955040e3 +#define p2 0.115302935154049e4 +#define p1 0.178040631643320e4 +#define p0 0.896785974036639e3 +#define q4 0.589569705084446e2 +#define q3 0.536265374031215e3 +#define q2 0.166678381488163e4 +#define q1 0.207933497444541e4 +#define q0 0.896785974036639e3 /** * atan makes its argument positive and @@ -36,13 +36,13 @@ */ float _atanf_c(float arg) { - float satan(float); + float satan(float); - if(arg>0) { - return(satan(arg)); - } else { - return(-satan(-arg)); - } + if(arg>0) { + return(satan(arg)); + } else { + return(-satan(-arg)); + } } double _atan_c(double) __attribute__((alias("_atanf_c"))); @@ -59,13 +59,13 @@ double _atan_c(double) __attribute__((alias("_atanf_c"))); */ static float xatan(float arg) { - float argsq; - float value; + float argsq; + float value; - argsq = arg*arg; - value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0); - value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0); - return(value*arg); + argsq = arg*arg; + value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0); + value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0); + return(value*arg); } /** @@ -74,11 +74,11 @@ static float xatan(float arg) { */ float satan(float arg) { - if(arg < sq2m1) { - return(xatan(arg)); - } else if(arg > sq2p1) { - return(pio2 - xatan(1.0/arg)); - } else { - return(pio4 + xatan((arg-1.0)/(arg+1.0))); - } + if(arg < sq2m1) { + return(xatan(arg)); + } else if(arg > sq2p1) { + return(pio2 - xatan(1.0/arg)); + } else { + return(pio4 + xatan((arg-1.0)/(arg+1.0))); + } } diff --git a/src/libc/atanhl.c b/src/libc/atanhl.c index 393201d4b..e3c0a756c 100644 --- a/src/libc/atanhl.c +++ b/src/libc/atanhl.c @@ -2,5 +2,5 @@ long double atanhl(long double x) { - return .5l * logl((1 + x) / (1 - x)); + return .5L * logl((1 + x) / (1 - x)); } diff --git a/src/libc/bsearch.c b/src/libc/bsearch.c index e43eb543a..77cbb8c7d 100644 --- a/src/libc/bsearch.c +++ b/src/libc/bsearch.c @@ -1,11 +1,12 @@ /************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C)1987-2008 by */ +/* Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ + #include #include @@ -14,62 +15,63 @@ * bsearch - binary search * * Inputs: -* key - key to search for -* base - base of array to be sorted -* num - number of elements -* width - size of each element -* comp - pointer to function for comparison +* key - key to search for +* base - base of array to be sorted +* num - number of elements +* width - size of each element +* comp - pointer to function for comparison * * Returns: -* nothing +* nothing * *************************************************/ -void *bsearch(void *keyp, void *ptr, size_t num, size_t width, - int (*comp)(const void *, const void *)) -{ - char *key = keyp; - char *base = ptr; - unsigned int mid; - unsigned int low; - unsigned int high; - int d; - unsigned int pmid; - char *addr; +void *bsearch( + void *keyp, void *ptr, size_t num, size_t width, + int (*comp)(const void *, const void *) +) { + char *key = keyp; + char *base = ptr; + unsigned int mid; + unsigned int low; + unsigned int high; + int d; + unsigned int pmid; + char *addr; - /* make high the nearest power of two for */ - /* efficiency and to ensure we always */ - /* terminate properly. */ - /* Note that high is always higher */ - /* than it should be so that we will */ - /* not fail to find the last entry in the */ - /* table. */ + /* make high the nearest power of two for */ + /* efficiency and to ensure we always */ + /* terminate properly. */ + /* Note that high is always higher */ + /* than it should be so that we will */ + /* not fail to find the last entry in the */ + /* table. */ - high = 0x0001; - while (high <= num) - high <<= 1; - low = 0; - mid = 0; + high = 0x0001; + while (high <= num) + high <<= 1; + low = 0; + mid = 0; - /* begin the search */ + /* begin the search */ - for(;;) { - pmid = mid; - mid = ((high - low) >> 1) + low; + for(;;) { + pmid = mid; + mid = ((high - low) >> 1) + low; - if (pmid == mid) /* terminate because we're */ - return(NULL); /* oscilating. */ + if (pmid == mid) /* terminate because we're */ + return(NULL); /* oscilating. */ - if (mid >= num) { /* we're above the array. */ - high = mid; /* pretend element is larger*/ - continue; /* than the key. */ - } + if (mid >= num) { /* we're above the array. */ + high = mid; /* pretend element is larger*/ + continue; /* than the key. */ + } - d = (*comp)(key,addr = base + mid * width); - if (d == 0) /* we found it */ - return(addr); - if (d < 0) /* key is less than mid, */ - high = mid; /* set high to mid. */ - else /* key is greater than mid, */ - low = mid; /* set low to mid. */ - } + d = (*comp)(key,addr = base + mid * width); + if (d == 0) /* we found it */ + return(addr); + if (d < 0) /* key is less than mid, */ + high = mid; /* set high to mid. */ + else /* key is greater than mid, */ + low = mid; /* set low to mid. */ + } } diff --git a/src/libc/coshf.c b/src/libc/coshf.c index fe39354b0..d6516299e 100644 --- a/src/libc/coshf.c +++ b/src/libc/coshf.c @@ -1,23 +1,23 @@ /************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C)1987-2008 by */ +/* Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ /* - sinh(arg) returns the hyperbolic sine of its floating- - point argument. + sinh(arg) returns the hyperbolic sine of its floating- + point argument. - The exponential function is called for arguments - greater in magnitude than 0.5. + The exponential function is called for arguments + greater in magnitude than 0.5. - A series is used for arguments smaller in magnitude than 0.5. - The coefficients are #2029 from Hart & Cheney. (20.36D) + A series is used for arguments smaller in magnitude than 0.5. + The coefficients are #2029 from Hart & Cheney. (20.36D) - cosh(arg) is computed from the exponential function for - all arguments. + cosh(arg) is computed from the exponential function for + all arguments. */ #include @@ -31,21 +31,21 @@ #define q2 -0.173678953558234e+3 float _coshf_c(float arg) { - float val; + float val; - if(arg < 0) { - arg = -arg; - } + if(arg < 0) { + arg = -arg; + } - val = expf(arg); + val = expf(arg); - if(arg > 21.) { - return val/2; - } + if(arg > 21.) { + return val/2; + } - val += expf(-arg); - val /= 2; - return val; + val += expf(-arg); + val /= 2; + return val; } double _cosh_c(double) __attribute__((alias("_coshf_c"))); diff --git a/src/libc/difftime.c b/src/libc/difftime.c index c98ec889a..f25b312e8 100644 --- a/src/libc/difftime.c +++ b/src/libc/difftime.c @@ -2,6 +2,6 @@ double difftime(time_t end, time_t beginning) { - /* assuming typedef unsigned long time_t */ + /* assuming typedef unsigned long time_t */ return (double)((signed long)(end - beginning)); } diff --git a/src/libc/erfcf.c b/src/libc/erfcf.c index fdb427aa1..42a38bc15 100644 --- a/src/libc/erfcf.c +++ b/src/libc/erfcf.c @@ -1,10 +1,22 @@ #include +/** + * Algorithm from: + * https://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions + */ float erfcf(float x) { - static const float p = 0.47047f, a1 = 0.3480242f, a2 = -0.0958798f, a3 = 0.7478556f; + static const float + p = 0.47047f, + a1 = 0.3480242f, + a2 = -0.0958798f, + a3 = 0.7478556f; const float t = 1.0f / (1.0f + p * fabsf(x)); - return copysignf(t * (a1 + t * (a2 + t * a3)) * expf(-x * x), x); + float ret = t * (a1 + t * (a2 + t * a3)) * expf(-x * x); + if (signbit(x)) { + ret = 2.0f - ret; + } + return ret; } double erfc(double) __attribute__((alias("erfcf"))); diff --git a/src/libc/erfcl.c b/src/libc/erfcl.c index bc9d0ea93..d2ce8eac5 100644 --- a/src/libc/erfcl.c +++ b/src/libc/erfcl.c @@ -1,8 +1,22 @@ #include +/** + * Algorithm from: + * https://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions + */ long double erfcl(long double x) { - static const long double p = 0.3275911l, a1 = 0.254829592l, a2 = -0.284496736l, a3 = 1.421413741l, a4 = -1.453152027l, a5 = 1.061405429l; - const long double t = 1.0l / (1.0l + p * fabsl(x)); - return copysignl(t * (a1 + t * (a2 + t * (a3 + t * (a4 + t * a5)))) * expl(-x * x), x); + static const long double + p = 0.3275911L, + a1 = 0.254829592L, + a2 = -0.284496736L, + a3 = 1.421413741L, + a4 = -1.453152027L, + a5 = 1.061405429L; + const long double t = 1.0L / (1.0L + p * fabsl(x)); + long double ret = t * (a1 + t * (a2 + t * (a3 + t * (a4 + t * a5)))) * expl(-x * x); + if (signbit(x)) { + ret = 2.0L - ret; + } + return ret; } diff --git a/src/libc/expf.c b/src/libc/expf.c index 7616cd492..d4dc3da00 100644 --- a/src/libc/expf.c +++ b/src/libc/expf.c @@ -1,50 +1,50 @@ /************************************************************************/ -/* */ -/* Copyright (C) 2000-2008 Zilog, Inc. */ -/* */ +/* */ +/* Copyright (C) 2000-2008 Zilog, Inc. */ +/* */ /************************************************************************/ /* - exp returns the exponential function of its - floating-point argument. + exp returns the exponential function of its + floating-point argument. - The coefficients are #1069 from Hart and Cheney. (22.35D) + The coefficients are #1069 from Hart and Cheney. (22.35D) */ #include #include -#define p0 0.20803843466947e7 -#define p1 0.30286971697440e5 -#define p2 0.60614853300611e2 -#define q0 0.60027203602388e7 -#define q1 0.32772515180829e6 -#define q2 0.17492876890931e4 -#define log2e 1.44269504088896 -#define sqrt2 1.41421356237310 -#define maxf 10000 +#define p0 0.20803843466947e7 +#define p1 0.30286971697440e5 +#define p2 0.60614853300611e2 +#define q0 0.60027203602388e7 +#define q1 0.32772515180829e6 +#define q2 0.17492876890931e4 +#define log2e 1.44269504088896 +#define sqrt2 1.41421356237310 +#define maxf 10000 float _expf_c(float arg) { - float fraction; - float temp1, temp2, xsq; - int ent; + float fraction; + float temp1, temp2, xsq; + int ent; - if ( arg == 0.0 ){ - return 1.0; - } - if ( arg < -maxf ){ - return 0.0; - } - if ( arg > maxf ){ - errno = ERANGE; - return HUGE_VALF; - } - arg *= log2e; - ent = floor( arg ); - fraction = arg - ent - 0.5; - xsq = fraction * fraction; - temp1 = ((p2 * xsq + p1) * xsq + p0) * fraction; - temp2 = ((1.0 * xsq + q2) * xsq + q1) * xsq + q0; - return ldexp( sqrt2 * (temp2+temp1) / (temp2-temp1), ent ); + if ( arg == 0.0 ){ + return 1.0; + } + if ( arg < -maxf ){ + return 0.0; + } + if ( arg > maxf ){ + errno = ERANGE; + return HUGE_VALF; + } + arg *= log2e; + ent = floor( arg ); + fraction = arg - ent - 0.5; + xsq = fraction * fraction; + temp1 = ((p2 * xsq + p1) * xsq + p0) * fraction; + temp2 = ((1.0 * xsq + q2) * xsq + q1) * xsq + q0; + return ldexp( sqrt2 * (temp2+temp1) / (temp2-temp1), ent ); } double _exp_c(double) __attribute__((alias("_expf_c"))); diff --git a/src/libc/floorf.c b/src/libc/floorf.c index 582058c6e..a5e5ea5fd 100644 --- a/src/libc/floorf.c +++ b/src/libc/floorf.c @@ -1,25 +1,26 @@ /************************************************************************/ -/* */ -/* Copyright (C) 2000-2008 Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C) 2000-2008 Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ #include float _floorf_c(float d) { - float fraction; + float fraction; - if (d<0.0) { - d = -d; - fraction = modff(d, &d); - if (fraction != 0.0) - d += 1; - d = -d; - } else { - fraction = modff(d, &d); - } - return(d); + if (d<0.0) { + d = -d; + fraction = modff(d, &d); + if (fraction != 0.0) { + d += 1; + } + d = -d; + } else { + fraction = modff(d, &d); + } + return(d); } double _floor_c(double) __attribute__((alias("_floorf_c"))); diff --git a/src/libc/frexpf.c b/src/libc/frexpf.c index 4a549a74a..4fe24b048 100644 --- a/src/libc/frexpf.c +++ b/src/libc/frexpf.c @@ -5,9 +5,9 @@ #if 1 /************************************************************************/ /* */ -/* Copyright (C) 2000-2008 Zilog, Inc. */ +/* Copyright (C) 2000-2008 Zilog, Inc. */ /* */ -/* San Jose, California */ +/* San Jose, California */ /* */ /************************************************************************/ diff --git a/src/libc/lgammaf.c b/src/libc/lgammaf.c index 3fb425833..07fb48c37 100644 --- a/src/libc/lgammaf.c +++ b/src/libc/lgammaf.c @@ -10,7 +10,7 @@ #define N 8 -#define B0 1 /* Bernoulli numbers */ +#define B0 1.0 /* Bernoulli numbers */ #define B1 (-1.0 / 2.0) #define B2 ( 1.0 / 6.0) #define B4 (-1.0 / 30.0) @@ -23,15 +23,27 @@ float lgammaf(float x) { /* the natural logarithm of the Gamma function. */ float v, w; + v = 1.0; - v = 1; - while (x < N) { v *= x; x++; } - w = 1 / (x * x); - return ((((((((B16 / (16 * 15)) * w + (B14 / (14 * 13))) * w - + (B12 / (12 * 11))) * w + (B10 / (10 * 9))) * w - + (B8 / ( 8 * 7))) * w + (B6 / ( 6 * 5))) * w - + (B4 / ( 4 * 3))) * w + (B2 / ( 2 * 1))) / x - + 0.5 * M_LOG_2M_PI - log(v) - x + (x - 0.5) * log(x); + /** + * This loop will take forever to terminate if `x < -100.0f`, so we have a + * maximum iteration count to ensure that the loop will terminate in a + * reasonable amount of time. `v` should overflow when `x < -33.0f` + */ + const int maximum_iter = 36 + N; + for (int iter = 0; iter < maximum_iter; iter++) { + if (x < (float)N) { + break; + } + v *= x; + x++; + } + w = 1.0 / (x * x); + return ((((((((B16 / (16.0 * 15.0)) * w + (B14 / (14.0 * 13.0))) * w + + (B12 / (12.0 * 11.0))) * w + (B10 / (10.0 * 9.0))) * w + + (B8 / ( 8.0 * 7.0))) * w + (B6 / ( 6.0 * 5.0))) * w + + (B4 / ( 4.0 * 3.0))) * w + (B2 / ( 2.0 * 1.0))) / x + + 0.5 * (float)M_LOG_2M_PI - logf(v) - x + (x - 0.5) * logf(x); } double lgamma(double) __attribute__((alias("lgammaf"))); diff --git a/src/libc/log10f.c b/src/libc/log10f.c index 2eef9e58b..36c815861 100644 --- a/src/libc/log10f.c +++ b/src/libc/log10f.c @@ -2,7 +2,7 @@ float _log10f_c(float x) { - return logf(x) / M_LN10; + return logf(x) * (float)(M_LOG10E); } double _log10_c(double) __attribute__((alias("_log10f_c"))); diff --git a/src/libc/log1pf.c b/src/libc/log1pf.c index 99f596665..f4123ead8 100644 --- a/src/libc/log1pf.c +++ b/src/libc/log1pf.c @@ -1,7 +1,7 @@ #include float log1pf(float x) { - if (fabs(x) <= 0.125) { + if (fabsf(x) <= 0.125) { // pade(series(ln(1+x),x=0,6,polynom),x,5,3) // (-57*x**2-90*x)/(x**3-21*x**2-102*x-90) // relative error less than 1e-7 diff --git a/src/libc/log2f.c b/src/libc/log2f.c index 6ffe86aa3..0e3d17103 100644 --- a/src/libc/log2f.c +++ b/src/libc/log2f.c @@ -2,7 +2,7 @@ float log2f(float x) { - return logf(x) * (float)(1 / M_LN2); + return logf(x) * (float)(M_LOG2E); } double log2(double) __attribute__((alias("log2f"))); diff --git a/src/libc/logf.c b/src/libc/logf.c index f6b735215..6bff1d1a2 100644 --- a/src/libc/logf.c +++ b/src/libc/logf.c @@ -1,53 +1,53 @@ /************************************************************************/ -/* */ -/* Copyright (C)1999-2008 by Zilog, Inc. */ +/* */ +/* Copyright (C) 1999-2008 by Zilog, Inc. */ /* */ /************************************************************************/ /* - log returns the natural logarithm of its floating - point argument. + log returns the natural logarithm of its floating + point argument. - The coefficients are #2705 from Hart & Cheney. (19.38D) + The coefficients are #2705 from Hart & Cheney. (19.38D) - It calls frexp. + It calls frexp. */ #include #include -#define log2 0.693147180559945e0 -#define ln10 2.30258509299404 -#define sqrto2 0.707106781186548e0 -#define p0 -0.240139179559211e2 -#define p1 0.309572928215377e2 -#define p2 -0.963769093368687e1 -#define p3 0.421087371217980e0 -#define q0 -0.120069589779605e2 -#define q1 0.194809660700890e2 -#define q2 -0.891110902798312e1 +#define ln2 0.693147180559945e0 +#define ln10 2.30258509299404 +#define sqrto2 0.707106781186548e0 +#define p0 -0.240139179559211e2 +#define p1 0.309572928215377e2 +#define p2 -0.963769093368687e1 +#define p3 0.421087371217980e0 +#define q0 -0.120069589779605e2 +#define q1 0.194809660700890e2 +#define q2 -0.891110902798312e1 float _logf_c(float arg) { - double x,z, zsq, temp; - int exp; - - if (arg <= 0.0) { - errno = EDOM; - return -HUGE_VALF; - } - x = frexpf(arg, & exp); - if ( x < sqrto2 ){ - x *= 2; - exp--; - } - - z = (x-1)/(x+1); - zsq = z*z; - - temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0; - temp = temp/(((1.0*zsq + q2)*zsq + q1)*zsq + q0); - temp = temp*z + exp*log2; - return temp; + float x, z, zsq, temp; + int expon; + + if (arg <= 0.0) { + errno = EDOM; + return -HUGE_VALF; + } + x = frexpf(arg, & expon); + if ( x < sqrto2 ){ + x *= 2; + expon--; + } + + z = (x-1)/(x+1); + zsq = z*z; + + temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0; + temp = temp/(((1.0*zsq + q2)*zsq + q1)*zsq + q0); + temp = temp*z + expon*ln2; + return temp; } double _log_c(double) __attribute__((alias("_logf_c"))); diff --git a/src/libc/nextupdownf.c b/src/libc/nextupdownf.c index dc5d33052..89d1358c5 100644 --- a/src/libc/nextupdownf.c +++ b/src/libc/nextupdownf.c @@ -37,6 +37,8 @@ float nextupf(float x) { return arg_x.flt; } +double nextup(double) __attribute__((alias("nextupf"))); + float nextdownf(float x) { F32_pun arg_x; arg_x.flt = x; @@ -62,5 +64,4 @@ float nextdownf(float x) { return arg_x.flt; } -double nextup(double) __attribute__((alias("nextupf"))); double nextdown(double) __attribute__((alias("nextdownf"))); diff --git a/src/libc/powf.c b/src/libc/powf.c index ebaf6bd3b..f3392fd08 100644 --- a/src/libc/powf.c +++ b/src/libc/powf.c @@ -1,34 +1,34 @@ /************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C)1987-2008 by */ +/* Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ -#include -#include +#include +#include float _powf_c(float arg1, float arg2) { - float result; - long temp; + float result; + long temp; - if ( arg1 > 0.0 ){ - return expf( arg2 * log( arg1 ) ); - } - if ( arg1 < 0.0 ){ - temp = arg2; - if ( temp == arg2 ){ - result = expf( arg2 * log( -arg1 ) ); - return temp & 1 ? -result : result; - } - errno = EDOM; - } - if ( arg2 <= 0.0 ){ - errno = EDOM; - } - return 0.0; + if ( arg1 > 0.0 ){ + return expf( arg2 * logf( arg1 ) ); + } + if ( arg1 < 0.0 ){ + temp = arg2; + if ( temp == arg2 ){ + result = expf( arg2 * logf( -arg1 ) ); + return temp & 1 ? -result : result; + } + errno = EDOM; + } + if ( arg2 <= 0.0 ){ + errno = EDOM; + } + return 0.0; } double _pow_c(double, double) __attribute__((alias("_powf_c"))); diff --git a/src/libc/qsort.c b/src/libc/qsort.c index 348fc6d5e..335a064ec 100644 --- a/src/libc/qsort.c +++ b/src/libc/qsort.c @@ -1,11 +1,12 @@ /************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C)1987-2008 by */ +/* Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ + #include #include @@ -14,89 +15,90 @@ * qsort - Quick sort * * Inputs: -* base - base of array to be sorted -* nel - number of elements -* size - size of each element -* compar - pointer to function for comparison +* base - base of array to be sorted +* nel - number of elements +* size - size of each element +* compar - pointer to function for comparison * * Returns: -* nothing +* nothing * *************************************************/ static void swapmem(char *a,char *b,size_t size) { - register char t; - register size_t i; + register char t; + register size_t i; - for (i=0;il = i; - sp->r = r; - ++sp; - } - r = j; /* continue sorting left partition */ - } - else { - if (base < j) { /* stack request for left partition */ - sp->l = base; - sp->r = j; - ++sp; - } - base = i; /* continue sorting right partition */ - } - } while (base < r); + sp = stack; + r = base + (nel-1)*size; + for (;;) { + do { + x = base + (r-base)/size/2 * size; + i = base; + j = r; + do { + while ((*compar)(i,x) < 0) + i += size; + while ((*compar)(x,j) < 0) + j -= size; + if (i < j) { + swapmem(i, j, size); + if (i == x) + x = j; + else if (j == x) + x = i; + } + if (i <= j) { + i += size; + j -= size; + } + } while (i <= j); + if (j-base < r-i) { + if (i < r) { /* stack request for right partition */ + sp->l = i; + sp->r = r; + ++sp; + } + r = j; /* continue sorting left partition */ + } + else { + if (base < j) { /* stack request for left partition */ + sp->l = base; + sp->r = j; + ++sp; + } + base = i; /* continue sorting right partition */ + } + } while (base < r); - if (sp <= stack) - break; - --sp; - base = sp->l; - r = sp->r; - } + if (sp <= stack) + break; + --sp; + base = sp->l; + r = sp->r; + } } diff --git a/src/libc/roundeven.c b/src/libc/roundeven.c index 957a401a9..b2c8e7318 100644 --- a/src/libc/roundeven.c +++ b/src/libc/roundeven.c @@ -5,7 +5,7 @@ float roundevenf(float x) float i, f = modff(x, &i); if (!isgreaterequal(f, .5f)) return i; if (f == .5f) { - float half = ldexp(i, -1); + float half = ldexpf(i, -1); if (truncf(half) == half) return i; } return signbit(i) ? i - 1 : i + 1; diff --git a/src/libc/signbitf.src b/src/libc/signbitf.src new file mode 100644 index 000000000..3be7e06f9 --- /dev/null +++ b/src/libc/signbitf.src @@ -0,0 +1,14 @@ + assume adl=1 + + section .text + + public __signbitf + +; bool _signbitf(float) +__signbitf: + ld hl, 6 + add hl, sp + ld a, (hl) + rla + sbc a, a + ret diff --git a/src/libc/signbit.src b/src/libc/signbitl.src similarity index 51% rename from src/libc/signbit.src rename to src/libc/signbitl.src index e48ae18fe..120a32844 100644 --- a/src/libc/signbit.src +++ b/src/libc/signbitl.src @@ -2,16 +2,7 @@ section .text - public __signbitf, __signbitl - -; bool _signbitf(float) -__signbitf: - ld hl, 6 - add hl, sp - ld a, (hl) - rla - sbc a, a - ret + public __signbitl ; bool _signbitl(long double) __signbitl: diff --git a/src/libc/sinf.c b/src/libc/sinf.c index 2395c99e1..b8544f5c6 100644 --- a/src/libc/sinf.c +++ b/src/libc/sinf.c @@ -1,65 +1,64 @@ /************************************************************************/ -/* */ +/* */ /* Copyright (C) 1999-2008 by Zilog, Inc. */ -/* */ +/* */ /************************************************************************/ /* - C program for floating point sin/cos. - Calls modf. - There are no error exits. - Coefficients are #3370 from Hart & Cheney (18.80D). + C program for floating point sin/cos. + Calls modf. + There are no error exits. + Coefficients are #3370 from Hart & Cheney (18.80D). */ #include - -#define twoopi 0.636619772367581 -#define p0 0.135788409787738e8 -#define p1 -0.494290810090284e7 -#define p2 0.440103053537527e6 -#define p3 -0.138472724998245e5 -#define p4 0.145968840666577e3 -#define q0 0.864455865292253e7 -#define q1 0.408179225234330e6 -#define q2 0.946309610153821e4 -#define q3 0.132653490878614e3 +#define two_over_pi 0.636619772367581 +#define p0 0.135788409787738e8 +#define p1 -0.494290810090284e7 +#define p2 0.440103053537527e6 +#define p3 -0.138472724998245e5 +#define p4 0.145968840666577e3 +#define q0 0.864455865292253e7 +#define q1 0.408179225234330e6 +#define q2 0.946309610153821e4 +#define q3 0.132653490878614e3 float sinus(float arg, int quad) { - float e, f; - int k; - float ysq; - float x,y; - float temp1, temp2; + float e, f; + int k; + float ysq; + float x,y; + float temp1, temp2; - x = arg; - if(x<0) { - x = -x; - quad = quad + 2; - } - x = x*twoopi; /*underflow?*/ - if(x>32764){ - y = modff(x,&e); - e = e + quad; - modff(0.25*e,&f); - quad = e - 4*f; - }else{ - k = x; - y = x - k; - quad = (quad + k) & 03; - } - if (quad & 01) - y = 1-y; - if(quad > 1) - y = -y; + x = arg; + if(x<0) { + x = -x; + quad = quad + 2; + } + x = x*two_over_pi; /*underflow?*/ + if(x>32764){ + y = modff(x,&e); + e = e + quad; + modff(0.25*e,&f); + quad = e - 4*f; + }else{ + k = x; + y = x - k; + quad = (quad + k) & 03; + } + if (quad & 01) + y = 1-y; + if(quad > 1) + y = -y; - ysq = y*y; - temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y; - temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0); - return(temp1/temp2); + ysq = y*y; + temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y; + temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0); + return(temp1/temp2); } float _sinf_c(float arg) { - return sinus(arg, 0); + return sinus(arg, 0); } double _sin_c(double) __attribute__((alias("_sinf_c"))); diff --git a/src/libc/sinhf.c b/src/libc/sinhf.c index 079c4c657..e94fada98 100644 --- a/src/libc/sinhf.c +++ b/src/libc/sinhf.c @@ -1,23 +1,23 @@ /************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C)1987-2008 by */ +/* Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ /* - sinh(arg) returns the hyperbolic sine of its floating- - point argument. + sinh(arg) returns the hyperbolic sine of its floating- + point argument. - The exponential function is called for arguments - greater in magnitude than 0.5. + The exponential function is called for arguments + greater in magnitude than 0.5. - A series is used for arguments smaller in magnitude than 0.5. - The coefficients are #2029 from Hart & Cheney. (20.36D) + A series is used for arguments smaller in magnitude than 0.5. + The coefficients are #2029 from Hart & Cheney. (20.36D) - cosh(arg) is computed from the exponential function for - all arguments. + cosh(arg) is computed from the exponential function for + all arguments. */ #include @@ -31,31 +31,31 @@ #define q2 -0.173678953558234e+3 float _sinhf_c(float arg) { - float temp, argsq; - register int sign; - - sign = 1; - if(arg < 0) { - arg = -arg; - sign = -1; - } - - if(arg > 21.) { - temp = expf(arg)/2; - if (sign>0) - return(temp); - else - return(-temp); - } - - if(arg > 0.5) { - return(sign*(expf(arg) - expf(-arg))/2); - } - - argsq = arg*arg; - temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg; - temp /= (((argsq+q2)*argsq+q1)*argsq+q0); - return(sign*temp); + float temp, argsq; + register int sign; + + sign = 1; + if(arg < 0) { + arg = -arg; + sign = -1; + } + + if(arg > 21.) { + temp = expf(arg)/2; + if (sign>0) + return(temp); + else + return(-temp); + } + + if(arg > 0.5) { + return(sign*(expf(arg) - expf(-arg))/2); + } + + argsq = arg*arg; + temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg; + temp /= (((argsq+q2)*argsq+q1)*argsq+q0); + return(sign*temp); } double _sinh_c(double, double *) __attribute__((alias("_sinhf_c"))); diff --git a/src/libc/strtof.c b/src/libc/strtof.c index 7ba41b1a4..912cfdd24 100644 --- a/src/libc/strtof.c +++ b/src/libc/strtof.c @@ -1,10 +1,10 @@ /************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C)1987-2008 by */ +/* Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ #include @@ -24,18 +24,17 @@ typedef union F32_pun { /************************************************* * -* strtod - string to double conversion +* strtof - string to float conversion * * Inputs: -* str - pointer to the character string -* endptr - pointer to pointer to char to -* put address of first char past -* the end of the string -- or NULL +* str - pointer to the character string +* endptr - pointer to pointer to char to +* put address of first char past +* the end of the string -- or NULL * Returns: -* the value of the number +* the value of the number * *************************************************/ - /** * @remarks `*str >= '0' && *str <= '9'` is smaller than calls to `isdigit(*str)` * @todo Add support for INF INFINITY NAN NAN(...) diff --git a/src/libc/strtold.c b/src/libc/strtold.c index 219cb1c61..ac4e86b7b 100644 --- a/src/libc/strtold.c +++ b/src/libc/strtold.c @@ -1,12 +1,12 @@ #if 1 /************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C)1987-2008 by */ +/* Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ #include @@ -29,12 +29,12 @@ typedef union F64_pun { * strtold - string to long double conversion * * Inputs: -* str - pointer to the character string -* endptr - pointer to pointer to char to -* put address of first char past -* the end of the string -- or NULL +* str - pointer to the character string +* endptr - pointer to pointer to char to +* put address of first char past +* the end of the string -- or NULL * Returns: -* the value of the number +* the value of the number * *************************************************/ diff --git a/src/libc/tanf.c b/src/libc/tanf.c index 768d3d96e..eff712c77 100644 --- a/src/libc/tanf.c +++ b/src/libc/tanf.c @@ -1,79 +1,80 @@ /************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C)1987-2008 by */ +/* Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ /* - floating point tangent + floating point tangent - A series is used after range reduction. - Coefficients are #4285 from Hart & Cheney. (19.74D) + A series is used after range reduction. + Coefficients are #4285 from Hart & Cheney. (19.74D) */ #include #include -#define invpi 1.27323954473516 -#define p0 -0.130682026475483e+5 -#define p1 0.105597090171495e+4 -#define p2 -0.155068565348327e+2 -#define p3 0.342255438724100e-1 -#define p4 0.338663864267717e-4 -#define q0 -0.166389523894712e+5 -#define q1 0.476575136291648e+4 -#define q2 -0.155503316403171e+3 +#define four_mul_invpi 1.27323954473516 + +#define p0 -0.130682026475483e+5 +#define p1 0.105597090171495e+4 +#define p2 -0.155068565348327e+2 +#define p3 0.342255438724100e-1 +#define p4 0.338663864267717e-4 +#define q0 -0.166389523894712e+5 +#define q1 0.476575136291648e+4 +#define q2 -0.155503316403171e+3 float _tanf_c(float arg) { - float sign, temp, e, x, xsq; - int flag, i; + float sign, temp, e, x, xsq; + int flag, i; - flag = 0; - sign = 1.; - if(arg < 0.){ - arg = -arg; - sign = -1.; - } - arg = arg*invpi; /*overflow?*/ + flag = 0; + sign = 1.; + if(arg < 0.){ + arg = -arg; + sign = -1.; + } + arg = arg*four_mul_invpi; /*overflow?*/ x = modff(arg, &e); - i = e; - switch(i%4) { - case 1: - x = 1. - x; - flag = 1; - break; + i = e; + switch(i%4) { + case 1: + x = 1. - x; + flag = 1; + break; - case 2: - sign = - sign; - flag = 1; - break; + case 2: + sign = - sign; + flag = 1; + break; - case 3: - x = 1. - x; - sign = - sign; - break; + case 3: + x = 1. - x; + sign = - sign; + break; - case 0: - break; - } + case 0: + break; + } - xsq = x*x; - temp = ((((p4*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0)*x; - temp = temp/(((1.0*xsq+q2)*xsq+q1)*xsq+q0); + xsq = x*x; + temp = ((((p4*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0)*x; + temp = temp/(((1.0*xsq+q2)*xsq+q1)*xsq+q0); - if(flag == 1) { - if(temp == 0.) { - errno = ERANGE; - if (sign>0) - return(HUGE_VALF); - return(-HUGE_VALF); - } - temp = 1./temp; - } - return(sign*temp); + if(flag == 1) { + if(temp == 0.) { + errno = ERANGE; + if (sign>0) + return(HUGE_VALF); + return(-HUGE_VALF); + } + temp = 1./temp; + } + return(sign*temp); } double _tan_c(double) __attribute__((alias("_tanf_c"))); diff --git a/src/libc/tanhf.c b/src/libc/tanhf.c index fd867cb83..bef020ccd 100644 --- a/src/libc/tanhf.c +++ b/src/libc/tanhf.c @@ -1,35 +1,35 @@ /************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ +/* */ +/* Copyright (C)1987-2008 by */ +/* Zilog, Inc. */ +/* */ +/* San Jose, California */ +/* */ /************************************************************************/ /* - tanh(arg) computes the hyperbolic tangent of its floating - point argument. + tanh(arg) computes the hyperbolic tangent of its floating + point argument. - sinh and cosh are called except for large arguments, which - would cause overflow improperly. + sinh and cosh are called except for large arguments, which + would cause overflow improperly. */ #include float _tanhf_c(float arg) { - float sign; + float sign; - sign = 1.; - if(arg < 0.){ - arg = -arg; - sign = -1.; - } + sign = 1.; + if(arg < 0.){ + arg = -arg; + sign = -1.; + } - if(arg > 21.) - return(sign); + if(arg > 21.) + return(sign); - return(sign*sinhf(arg)/coshf(arg)); + return(sign*sinhf(arg)/coshf(arg)); } double _tanh_c(double) __attribute__((alias("_tanhf_c"))); diff --git a/src/libc/tgammaf.c b/src/libc/tgammaf.c index 589efa91f..b05a006dc 100644 --- a/src/libc/tgammaf.c +++ b/src/libc/tgammaf.c @@ -11,19 +11,19 @@ float tgammaf(float x) { /* Gamma function */ if (x == 0.0) { /* Pole Error */ errno = ERANGE; - return 1/x < 0 ? -HUGE_VAL : HUGE_VAL; + return signbit(x) ? -HUGE_VALF : HUGE_VALF; } if (x < 0) { int sign; - static float zero = 0.0; + static float zero = 0.0; float i, f; f = modff(-x, &i); if (f == 0.0) { /* Domain Error */ errno = EDOM; - return zero/zero; + return zero/zero; /* probably better to return NAN here */ } sign = (fmodf(i, 2.0) != 0.0) ? 1 : -1; - return sign * M_PI / (sinf(M_PI * f) * expf(lgammaf(1 - x))); + return (sign * M_PI) / (sinf(M_PI * f) * expf(lgammaf(1 - x))); } return expf(lgammaf(x)); } diff --git a/src/linker.mk b/src/linker.mk index 2e0b0b7dc..d99b6108d 100644 --- a/src/linker.mk +++ b/src/linker.mk @@ -63,8 +63,8 @@ linker_script: $(STATIC_FILES) $(LINKED_FILES) $(SHARED_FILES) $(Q)$(call APPEND_FILES,source ,ce,$(sort $(CE_FILES))) $(Q)$(call APPEND,if HAS_LIBC) $(Q)$(call APPEND_FILES, source ,libc,$(sort $(LIBC_FILES))) + $(Q)$(call APPEND_FILES, source ,softfloat,$(sort $(SOFTFLOAT_FILES))) $(Q)$(call APPEND,end if) $(Q)$(call APPEND,if HAS_LIBCXX) $(Q)$(call APPEND_FILES, source ,libcxx,$(sort $(LIBCXX_FILES))) - $(Q)$(call APPEND_FILES, source ,softfloat,$(sort $(SOFTFLOAT_FILES))) $(Q)$(call APPEND,end if)