diff --git a/src/libc/frexpf.src b/src/libc/frexpf.src index 9b5e2682a..d7d723893 100644 --- a/src/libc/frexpf.src +++ b/src/libc/frexpf.src @@ -15,6 +15,7 @@ else _frexp: _frexpf: ld iy, 0 + lea de, iy - 127 ; bias ; $FFFF81 add iy, sp ld hl, (iy + 3) add hl, hl @@ -29,40 +30,39 @@ _frexpf: rr b sbc hl, hl ld l, a - ld de, -127 ; bias +; ld de, -127 ; bias add hl, de +.ret_zero: + ex de, hl .ret_subnormal: res 7, (iy + 5) -.ret_zero: .ret_self: - ld de, (iy + 9) ; int *expon - ex de, hl + ld hl, (iy + 9) ; int *expon ld (hl), de ld hl, (iy + 3) ; mantissa ld e, b ; exponent ret .inf_nan: - ld hl, $7FFFFF ; INT_MAX + ld de, $7FFFFF ; INT_MAX jr .ret_self .maybe_subnormal: - add hl, bc - or a, a - sbc hl, bc - jr z, .ret_zero + ld e, d ; ld de, -1 + add hl, de + inc hl ; restore HL + jr nc, .ret_zero ; input: HL output: A call __ictlz ld c, a call __ishl ld (iy + 3), hl - scf - sbc hl, hl + sub a, 131 ; (127 + 3) + 1? idk where this magic number comes from cpl - add a, 131 ; 127 + 3 + 1? idk where this magic number comes from - ld l, a - ld a, b ; exponent - xor a, $3F + ld e, a ; DE was -1 from before + ; copy exponent of 0.5f + ld a, b + xor a, $3F ld b, a jr .ret_subnormal diff --git a/src/libc/ldexpf.c b/src/libc/ldexpf.c deleted file mode 100644 index 81e39be6a..000000000 --- a/src/libc/ldexpf.c +++ /dev/null @@ -1,211 +0,0 @@ -/** - * @remarks Zilog's implementation of ldexpf is probably smaller since it does - * not handle subnormals. Once subnormal support in _fmul is added, switch to - * the other routine. - */ -#if 1 - -/************************************************************************/ -/* */ -/* Copyright (C)1987-2008 by */ -/* Zilog, Inc. */ -/* */ -/* San Jose, California */ -/* */ -/************************************************************************/ -#include -#include -#include - -typedef union -{ - float value; - unsigned long bits; -} Ieee754; - -enum -{ // The IEEE 754 format is: - // SEEEEEEE EMMMMMMM MMMMMMMM MMMMMMMM - // (with an implicit mantissa high-order 1-bit.) - mastissa_shift = 0, - mastissa_bits = 23, - exponent_bits = 8, - exponent_shift = mastissa_shift + mastissa_bits, - sign_shift = exponent_shift + exponent_bits, - exponent_mask = (1 << exponent_bits) -1, // shifted = 0x7F800000 - exponent_max = (1 << exponent_bits) -1, - exponent_base = 127, - exponent_min = 0 -}; - -// ldexp - Standard C library routine -// ldexp returns the argument multiplied by an integral (positive or -// negative) power of two. -// -// Arguments: -// value - the floating point argument -// power - the power of two to be used -// -// Returns: -// - the argument multiplied by an integral power of two -// -float _ldexpf_c(float value, int power) -{ - Ieee754 floating; - int exponent; - int powerplusexponent; - - if ( value == 0.0 || !isfinite(value) ) return value; - - floating.value = value; - exponent = (floating.bits >> exponent_shift) & exponent_mask; - powerplusexponent = power + exponent; - - if ( powerplusexponent > exponent_max ) - { - errno = ERANGE; - return (floating.bits & (1L << sign_shift)) == 0 ? HUGE_VALF : - HUGE_VALF; - } - if ( powerplusexponent <= exponent_min ) // CR 3964 - { - errno = ERANGE; - return 0.0; - } - floating.bits += (long) power << exponent_shift; // adjust exponent - return floating.value; -} - -#else - -#include -#include -#include - -typedef union F32_pun { - float flt; - uint32_t bin; - struct { - uint24_t HL; - uint8_t E; - } reg; -} F32_pun; - -#define Float32_exp_bias 127 -#define Float32_ilogb_subnorm_max -127 -#define Float32_max_exp 127 -#define Float32_norm_min_exp -126 -#define Float32_biased_inf_nan_exp 255 - -#define Float32_mantissa_bits 23 -#define Float32_exponent_bits 8 -#define Float32_sign_bits 1 -#define Float32_implicit_mantissa_bits 1 - -/** Generates a normalized constant that is 2^expon */ -static float generate_ldexpf_mult(int expon) { - F32_pun ret; - ret.bin = (uint32_t)(expon + Float32_exp_bias) << Float32_mantissa_bits; - return ret.flt; -} - -/** - * @note x is assumed to be positive - * @remarks Cases where the result/output is subnormal rely on _fmul handling - * subnormals correctly. - */ -float _ldexpf_c_positive(float x, int expon) { - F32_pun val; - val.flt = x; - int x_class = fpclassify(val.flt); - switch(x_class) { - case FP_INFINITE: - case FP_NAN: - case FP_ZERO: { - // return unmodified - return val.flt; - } - case FP_NORMAL: { - val.bin <<= 1; - int x_exp = val.reg.E; - x_exp -= Float32_exp_bias; - x_exp += expon; - // overflow - if (x_exp > Float32_max_exp) { - errno = ERANGE; - return HUGE_VALF; - } - // normalized - if (x_exp >= Float32_norm_min_exp) { - // Copy the new exponent - x_exp += Float32_exp_bias; - val.reg.E = (uint8_t)x_exp; - val.bin >>= 1; - return val.flt; - } - // make subnormal (with correct rounding) - int subnorm_exp = x_exp - Float32_norm_min_exp; - if (subnorm_exp < -(Float32_mantissa_bits + Float32_implicit_mantissa_bits)) { - // rounds to zero (round to nearest ties to even) - errno = ERANGE; - return 0.0f; - } - // Copy the exponent of the minimum normalized value - val.reg.E = 0x01; - val.bin >>= 1; - // precision may be lost - val.flt *= generate_ldexpf_mult(subnorm_exp); - return val.flt; - } - case FP_SUBNORMAL: { - // make normalized (with correct rounding) - if (expon < -Float32_mantissa_bits) { - // rounds to zero (round to nearest ties to even) - errno = ERANGE; - return 0.0f; - } - - int clz_result = __builtin_clz(val.reg.HL); - int shift_amount = clz_result - Float32_implicit_mantissa_bits + 1; - - // fallsback to multiplication if the value will stay subnormal - if (expon - shift_amount < 0) { - // precision may be lost - val.flt *= generate_ldexpf_mult(expon); - return val.flt; - } - - expon -= shift_amount; - if (expon >= Float32_biased_inf_nan_exp - 1) { - // overflow - errno = ERANGE; - return HUGE_VALF; - } - // Shift everything such that the MSB of the mantissa is in the LSB of the exponent - val.reg.HL <<= shift_amount; - // Add the exponent - #if 1 - if ((expon & 1)) { - val.reg.HL += 0x800000; - expon++; - } - val.reg.E += ((uint8_t)expon >> 1); - return val.flt; - #else - val.bin <<= 1; - val.reg.E += (uint8_t)expon; - val.bin >>= 1; - return val.flt; - #endif - } - default: - __builtin_unreachable(); - } -} - -float _ldexpf_c(float x, int expon) { - return copysignf(_ldexpf_c_positive(fabsf(x), expon), x); -} - -#endif - -double _ldexp_c(double, int) __attribute__((alias("_ldexpf_c"))); diff --git a/src/libc/ldexpf.src b/src/libc/ldexpf.src index edb1601c6..74bd1c43c 100644 --- a/src/libc/ldexpf.src +++ b/src/libc/ldexpf.src @@ -16,12 +16,198 @@ _scalbn := _ldexpf else -_ldexpf := __ldexpf_c -_scalbnf := __ldexpf_c -_ldexp := __ldexp_c -_scalbn := __ldexp_c +if 1 - extern __ldexpf_c - extern __ldexp_c +; NOTE: since Ti floats are used, negative zero will not be returned unless the +; input was negative zero. +; +; normal inputs are handled correctly, unless the output is subnormal +; subnormal inputs/outputs return zero and set ERANGE and FE_INEXACT +; zero/infinite/NaN inputs are handled correctly +_ldexpf: +_ldexp: +_scalbnf: +_scalbn: + ld iy, 0 + lea bc, iy + 0 + add iy, sp + ld hl, (iy + 3) ; mant + add hl, hl + ld a, (iy + 6) ; expon + adc a, a + jr z, .maybe_subnormal + inc a + jr z, .ret_self ; inf NaN + dec a + ex de, hl + ld hl, (iy + 9) ; scale + ld c, a + add hl, bc ; add expon + ld a, l + bit 7, (iy + 11) ; scale sign + jr z, .scale_up +.scale_down: + ; test signbit + add hl, hl + jr nc, .finish +.underflow_to_zero: + ld hl, ___fe_cur_env + set 5, (hl) ; FE_INEXACT +if 0 + ld a, (iy + 6) + and a, $80 ; copysign +else + xor a, a ; avoid returning negative zero with Ti's floats +end if + sbc hl, hl +.common_erange: + ld bc, 5 ; ERANGE + ld (_errno), bc + ld e, a + ret +.overflow_to_inf: + ld hl, $800000 + ld a, (iy + 6) + or a, $7F ; copysign + jr .common_erange + +.scale_up: + ld bc, -255 + add hl, bc + jr c, .overflow_to_inf +.finish: + ld l, a + ex de, hl + ; signbit(A) E:UHL >>= 1 + ld a, (iy + 6) ; expon + push hl + rla + rr e + rr (iy - 1) + pop hl + rr h + rr l + ret + +.maybe_subnormal: + dec bc ; BC is now -1 + add hl, bc + jr c, .underflow_to_zero + ; return zero +.ret_self: + ld hl, (iy + 3) + ld e, (iy + 6) + ret + + +else + +; normal inputs are handled correctly, unless the output is subnormal +; subnormal inputs are handled correctly for positive scaling values +; subnormal outputs return zero and set ERANGE and FE_INEXACT for negative scaling values +; zero/infinite/NaN inputs are handled correctly +_ldexpf: +_ldexp: +_scalbnf: +_scalbn: + ld iy, 0 + lea bc, iy + 0 + add iy, sp + ld hl, (iy + 3) ; mant + add hl, hl + ld e, (iy + 6) ; expon + ld a, e + rl e + jr z, .maybe_subnormal + inc e + jr z, .ret_self ; inf NaN + dec e + ld c, e + ex de, hl + ld hl, (iy + 9) ; scale + add hl, bc ; add expon + ld a, l + bit 7, (iy + 11) ; scale sign + jr z, .scale_up +.scale_down: + ; test signbit + push hl + add hl, hl + pop hl + jr nc, .finish +; jr .underflow_to_zero +.underflow_to_zero: + ld hl, ___fe_cur_env + set 5, (hl) ; FE_INEXACT + ld a, (iy + 6) + and a, $80 ; copysign + sbc hl, hl +.common_erange: + ld bc, 5 ; ERANGE + ld (_errno), bc + ld e, a + ret +.overflow_to_inf: + ld hl, $800000 + ld a, (iy + 6) + or a, $7F ; copysign + jr .common_erange + +.scale_up: + ld bc, -255 + add hl, bc + jr c, .overflow_to_inf +.finish: + ld l, a + ex de, hl +.finish_subnormal: + ; signbit(A) E:UHL >>= 1 + ld a, (iy + 6) ; expon + push hl + rla + rr e + rr (iy - 1) + pop hl + rr h + rr l + ret + +.maybe_subnormal: + dec bc ; BC is now -1 + add hl, bc +; jr c, .underflow_to_zero + jr c, .subnormal + ; return zero +.ret_self: + ld hl, (iy + 3) + ld e, (iy + 6) + ret + +.subnormal: + ; BC is -1 here + bit 7, (iy + 11) ; scale sign + jr nz, .underflow_to_zero + ld de, (iy + 9) ; scale + ld hl, (iy + 3) ; mant +.norm_loop: + add hl, hl + jr c, .normalized + ex de, hl + add hl, bc ; --scale + ex de, hl + jr c, .norm_loop +; .still_subnormal: + ld e, 0 + jr .finish_subnormal +.normalized: + inc de + ex de, hl + ld a, l + jr .scale_up + +end if + + extern ___fe_cur_env + extern _errno end if diff --git a/test/floating_point/float32_frexp/src/main.c b/test/floating_point/float32_frexp/src/main.c index 62d4e2f51..e8134b33c 100644 --- a/test/floating_point/float32_frexp/src/main.c +++ b/test/floating_point/float32_frexp/src/main.c @@ -31,8 +31,10 @@ size_t run_test(void) { result.flt = frexpf(input[i], &expon); if (result.bin != output[i].frac.bin || expon != output[i].expon) { if (!(isnan(result.flt) && isnan(output[i].frac.flt))) { - // printf("G: %08lX %d\n", result.bin, expon); - // printf("T: %08lX %d\n", output[i].frac.bin, output[i].expon); + #if 0 + printf("G: %08lX %d\n", result.bin, expon); + printf("T: %08lX %d\n", output[i].frac.bin, output[i].expon); + #endif return i; } } diff --git a/test/floating_point/float32_ldexp/src/main.c b/test/floating_point/float32_ldexp/src/main.c index 4a28f5332..1256e0986 100644 --- a/test/floating_point/float32_ldexp/src/main.c +++ b/test/floating_point/float32_ldexp/src/main.c @@ -8,9 +8,6 @@ #include #include -/* enable if the toolchain is configured to use the subnormal compliant ldexpf */ -#if 0 - #include "f32_ldexp_LUT.h" #define ARRAY_LENGTH(x) (sizeof(x) / sizeof(x[0])) @@ -30,10 +27,26 @@ size_t run_test(void) { for (size_t i = 0; i < length; i++) { F32_pun result; result.flt = ldexpf(input[i].value, input[i].expon); + // ignoring subnormal inputs for now + if (issubnormal(input[i].value) || issubnormal(output[i].flt)) { + continue; + } if (result.bin != output[i].bin) { - if (!(isnan(result.flt) && isnan(output[i].flt))) { + // ignore NaN's with differing payloads + // treat signed zeros as equal for now + if ( + (!(isnan(result.flt) && isnan(output[i].flt))) && + (!(iszero(result.flt) && iszero(output[i].flt))) + ) { /* Float multiplication does not handle subnormals yet */ if (!(iszero(result.flt) && issubnormal(output[i].flt))) { + #if 0 + printf( + "%zu:\nI: %08lX %+d\nG: %08lX\nT: %08lX\n", + i, *(uint32_t*)(void*)&(input[i].value), input[i].expon, + result.bin, output[i].bin + ); + #endif return i; } } @@ -57,14 +70,3 @@ int main(void) { return 0; } - -#else - -int main(void) { - os_ClrHome(); - printf("All tests passed"); - while (!os_GetCSC()); - return 0; -} - -#endif