fixed signed overflow in ldexpl and implemented scalbln(f/l) in asm

Author: ZERICO2005

src/libc/ldexpl.c

@@ -42,6 +42,14 @@ static long double generate_ldexpl_mult(int expon) {
  * subnormal values.
  */
 static long double _ldexpl_c_positive(long double x, int expon) {
+    /* clamps the exponent to avoid signed overflow bugs */
+    /* hopefully these can also remove call pe, __setflag from the assembly output */
+    if (expon > 4095) {
+        expon = 4095;
+    }
+    if (expon < -4096) {
+        expon = -4096;
+    }
     F64_pun val;
     val.flt = x;
     /* expon == 0 || iszero(x) */

src/libc/scalbln.c

@@ -0,0 +1,32 @@
+	assume	adl=1
+
+	section	.text
+
+	public	_scalblnf, _scalbln
+
+; float _scalblnf(float, long)
+_scalbln:
+_scalblnf:
+	; We need to cast the int32_t scale to int24_t via saturation.
+	; scalblnf will always overflow/underflow if the scale is larger than +-280.
+	ld	hl, 12
+	add	hl, sp
+	ld	a, (hl)
+	dec	hl
+	; If bits [16, 23] match bits [24, 31] then we don't need to do anything.
+	; It implies that the scale is either already [-65536, +65535] and can
+	; be safely truncated, or that the original and truncated scale values
+	; are both larger than +-65535 and have the same sign.
+.overflow:
+	cp	a, (hl)
+	jp	z, _scalbnf
+	; We need modify the scale value to ensure that overflow/underflow still occurs.
+	; Ensure that bits[16, 23] are not all zeros/ones to set the scale to a large magnitude.
+	or	a, $03
+	dec	a
+	; bit 23 = bit 31, bit 16 is cleared, and bit 17 is set
+	ld	(hl), a	; store the new scale
+	; cp a, a	; set Z flag
+	jr	.overflow
+
+	extern	_scalbnf

src/libc/scalblnf.src

@@ -0,0 +1,31 @@
+	assume	adl=1
+
+	section	.text
+
+	public	_scalblnl
+
+; long double _scalblnl(long double, long)
+_scalblnl:
+	; We need to cast the int32_t scale to int24_t via saturation.
+	; scalblnl will always overflow/underflow if the scale is larger than +-2100.
+	ld	hl, 15
+	add	hl, sp
+	ld	a, (hl)
+	dec	hl
+	; If bits [16, 23] match bits [24, 31] then we don't need to do anything.
+	; It implies that the scale is either already [-65536, +65535] and can
+	; be safely truncated, or that the original and truncated scale values
+	; are both larger than +-65535 and have the same sign.
+.overflow:
+	cp	a, (hl)
+	jp	z, _scalbnl
+	; We need modify the scale value to ensure that overflow/underflow still occurs.
+	; Ensure that bits[16, 23] are not all zeros/ones to set the scale to a large magnitude.
+	or	a, $03
+	dec	a
+	; bit 23 = bit 31, bit 16 is cleared, and bit 17 is set
+	ld	(hl), a	; store the new scale
+	; cp a, a
+	jr	.overflow
+
+	extern	_scalbnl

src/libc/scalblnl.src

@@ -2,6 +2,8 @@
 #include <stddef.h>
 #include <stdint.h>
 #include <stdio.h>
+#include <limits.h>
+#include <string.h>
 #include <math.h>
 #include <fenv.h>
 #include <errno.h>
@@ -92,8 +94,58 @@ size_t run_test(void) {
     return SIZE_MAX;
 }
 
+static int32_t rand_i32_expon() {
+    union {
+        struct {
+            uint16_t part[2];
+        };
+        int32_t full;
+    } ret;
+    ret.part[0] = rand() & 0x03FF;
+    ret.part[1] = rand() & 0xC0C0;
+    return ret.full;
+}
+
+static float rand_f32() {
+    union {
+        struct {
+            uint16_t part[2];
+        };
+        float full;
+    } ret;
+    ret.part[0] = rand();
+    ret.part[1] = rand();
+    return ret.full;
+}
+
+static int clamp_exponent(long expon) {
+    if (expon > INT_MAX) {
+        return INT_MAX;
+    }
+    if (expon < INT_MIN) {
+        return INT_MIN;
+    }
+    return (int)expon;
+}
+
+void run_edge_case(void) {
+    srand(0x7184CE);
+    float input, output_i24, output_i32;
+    for (size_t i = 0; i < 512; i++) {
+        input = rand_f32();
+        int32_t expon = rand_i32_expon();
+        output_i24 = scalbnf(input, clamp_exponent(expon));
+        output_i32 = scalblnf(input, expon);
+        if (memcmp(&output_i24, &output_i32, sizeof(float)) != 0) {
+            fputs("Failed edge case\n", stdout);
+            return;
+        }
+    }
+}
+
 int main(void) {
     os_ClrHome();
+    run_edge_case();
     size_t fail_index = run_test();
     if (fail_index == SIZE_MAX) {
         fputs("All tests passed", stdout);

test/floating_point/float32_ldexp/src/main.c

@@ -2,6 +2,8 @@
 #include <stddef.h>
 #include <stdint.h>
 #include <stdio.h>
+#include <limits.h>
+#include <string.h>
 #include <math.h>
 #include <assert.h>
 #include <ti/screen.h>
@@ -40,8 +42,60 @@ size_t run_test(void) {
     return SIZE_MAX;
 }
 
+static int32_t rand_i32_expon() {
+    union {
+        struct {
+            uint16_t part[2];
+        };
+        int32_t full;
+    } ret;
+    ret.part[0] = rand() & 0x1FFF;
+    ret.part[1] = rand() & 0xC0C0;
+    return ret.full;
+}
+
+static long double rand_f64() {
+    union {
+        struct {
+            uint16_t part[4];
+        };
+        long double full;
+    } ret;
+    ret.part[0] = rand();
+    ret.part[1] = rand();
+    ret.part[2] = rand();
+    ret.part[3] = rand();
+    return ret.full;
+}
+
+static int clamp_exponent(long expon) {
+    if (expon > INT_MAX) {
+        return INT_MAX;
+    }
+    if (expon < INT_MIN) {
+        return INT_MIN;
+    }
+    return (int)expon;
+}
+
+void run_edge_case(void) {
+    srand(0x7184CE);
+    long double input, output_i24, output_i32;
+    for (size_t i = 0; i < 512; i++) {
+        input = rand_f64();
+        int32_t expon = rand_i32_expon();
+        output_i24 = scalbnl(input, clamp_exponent(expon));
+        output_i32 = scalblnl(input, expon);
+        if (memcmp(&output_i24, &output_i32, sizeof(long double)) != 0) {
+            fputs("Failed edge case\n", stdout);
+            return;
+        }
+    }
+}
+
 int main(void) {
     os_ClrHome();
+    run_edge_case();
     size_t fail_index = run_test();
     if (fail_index == SIZE_MAX) {
         fputs("All tests passed", stdout);