Compiler Explorer

Source code

#include <assert.h>
#include <memory.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

typedef uint8_t u8;
typedef uint32_t u32;
typedef int32_t i32;
typedef int64_t i64;
typedef uint64_t u64;
typedef i64 isize;
typedef float f32;
typedef double f64;

#define SWAP(T, a, b)                                                                              \
    do {                                                                                           \
        T swap = a;                                                                                \
        a = b;                                                                                     \
        b = swap;                                                                                  \
    } while (false)

#define ARRAY_SIZE(array_literal) sizeof(array_literal) / sizeof(array_literal[0])

typedef struct {
    char *data;
    isize size;
} StringView;

#define BIGINT_MAX_SIZE 128

void bigint_multiply(
    u32 const *term_a,
    isize term_a_size,
    u32 const *term_b,
    isize term_b_size,
    u32 *product,
    isize product_size
) {
    for (isize term_b_index = 0; term_b_index < term_b_size; term_b_index += 1) {
        u32 carry = 0;

for (isize term_a_index = 0; term_a_index < term_a_size; term_a_index += 1) {
            if (term_a_index + term_b_index >= product_size) {
                break;
            }

u64 digits_product = (u64) term_a[term_a_index] * term_b[term_b_index] +
                                 product[term_a_index + term_b_index] + carry;

product[term_a_index + term_b_index] = (u32) digits_product;
            carry = (u32) (digits_product >> 32);
        }

if (term_b_index + term_a_size < product_size) {
            product[term_b_index + term_a_size] = carry;
        }
    }
}

u32 bigint_divide_by_digit(u32 const *dividend, isize dividend_size, u32 divisor, u32 *quotient) {
    u32 remainder = 0;

for (isize dividend_index = dividend_size - 1; dividend_index >= 0; dividend_index -= 1) {
        u64 dividend_digits_pair = ((u64) remainder << 32) | dividend[dividend_index];

quotient[dividend_index] = (u32) (dividend_digits_pair / divisor);
        remainder = (u32) (dividend_digits_pair % divisor);
    }

return remainder;
}

void bigint_subtract_digit(
    u32 const *minuend, isize minuend_size, u32 subtrahend, u32 *difference
) {
    i64 borrow = (i64) subtrahend;

for (isize minuend_index = 0; minuend_index < minuend_size; minuend_index += 1) {
        i64 digits_difference = (i64) minuend[minuend_index] - borrow;

if (digits_difference >= 0) {
            difference[minuend_index] = (u32) digits_difference;
            break;
        } else {
            difference[minuend_index] = (u32) (digits_difference + ((i64) 1 << 32));
            borrow = 1;
        }
    }
}

void bigint_add_digit(u32 const *addend, isize addend_size, u32 addend_digit, u32 *sum) {
    u32 carry = addend_digit;

for (isize addend_index = 0; addend_index < addend_size; addend_index += 1) {
        u64 digits_sum = (u64) addend[addend_index] + carry;

sum[addend_index] = (u32) digits_sum;
        carry = (u32) (digits_sum >> 32);
    }
}

bool bigint_is_zero(u32 *bigint, isize bigint_size) {
    u32 *bigint_iter = bigint;
    while (bigint_iter < bigint + bigint_size) {
        if (*bigint_iter != 0) {
            return false;
        }

bigint_iter += 1;
    }

return true;
}

void bigint_u32_power(u32 base, u32 exponent, u32 *result, isize result_size) {
    assert(result_size <= BIGINT_MAX_SIZE);

u32 power_buffer[BIGINT_MAX_SIZE] = {0}, power_swap_buffer[BIGINT_MAX_SIZE] = {0};
    u32 *power = power_buffer, *power_swap = power_swap_buffer;
    power[0] = 1;

u32 base_buffer[BIGINT_MAX_SIZE] = {0}, base_swap_buffer[BIGINT_MAX_SIZE] = {0};
    u32 *base_bigint = base_buffer, *base_swap = base_swap_buffer;
    base_bigint[0] = base;

u32 current_exponent = exponent;

while (current_exponent > 0) {
        if (current_exponent % 2 == 1) {
            memset(power_swap, 0, result_size);
            bigint_multiply(power, result_size, base_bigint, result_size, power_swap, result_size);
            SWAP(u32 *, power, power_swap);
        }

memset(base_swap, 0, result_size);
        bigint_multiply(base_bigint, result_size, base_bigint, result_size, base_swap, result_size);
        SWAP(u32 *, base_bigint, base_swap);

current_exponent /= 2;
    }

memcpy(result, power, result_size * sizeof(u32));
}

void bigint_power_of_two(u32 exponent, u32 *result, isize result_size) {
    memset(result, 0, result_size * sizeof(u32));
    result[exponent / 32] = 1 << (exponent % 32);
}

StringView bigint_string_print(
    u32 const *bigint, isize bigint_size, char *output_buffer, isize output_buffer_size
) {
    assert(bigint_size <= BIGINT_MAX_SIZE);

u32 bigint_quotient[BIGINT_MAX_SIZE];
    memcpy(bigint_quotient, bigint, bigint_size * sizeof(u32));

char *output_buffer_iter = output_buffer + output_buffer_size - 1;

do {
        u32 remainder = bigint_divide_by_digit(bigint_quotient, bigint_size, 10, bigint_quotient);

*output_buffer_iter = '0' + remainder;
        output_buffer_iter -= 1;
    } while (!bigint_is_zero(bigint_quotient, bigint_size) && output_buffer_iter >= output_buffer);

return (StringView){
        .data = output_buffer_iter + 1,
        .size = output_buffer + output_buffer_size - (output_buffer_iter + 1),
    };
}

i32 bigint_compare(
    u32 const *left_bigint,
    isize left_bigint_size,
    u32 const *right_big_int,
    isize right_bigint_size
) {
    assert(left_bigint_size == right_bigint_size);

for (isize i = left_bigint_size - 1; i >= 0; i -= 1) {
        if (left_bigint[i] > right_big_int[i]) {
            return 1;
        }

if (left_bigint[i] < right_big_int[i]) {
            return -1;
        }
    }

return 0;
}

typedef struct {
    bool accept_smaller;
    bool accept_larger;
    bool break_tie_down;
} FloatPrintConfig;

#define F32_DECIMAL_MANTISSA_BIGINT_SIZE 12
#define F32_MANTISSA_LENGTH 23
#define F32_EXPONENT_BIAS 127
#define F32_SIGN_MASK 0x80000000
#define F32_EXPONENT_MASK 0xFF
#define F32_DECIMAL_MANTISSA_MAX_DIGITS 9

void f32_print(f32 f32_value, FloatPrintConfig config, char *result) {
    u32 f32_as_u32;
    ((u8 *) &f32_as_u32)[0] = ((u8 *) &f32_value)[0];
    ((u8 *) &f32_as_u32)[1] = ((u8 *) &f32_value)[1];
    ((u8 *) &f32_as_u32)[2] = ((u8 *) &f32_value)[2];
    ((u8 *) &f32_as_u32)[3] = ((u8 *) &f32_value)[3];

bool is_negative = (f32_as_u32 & F32_SIGN_MASK) != 0;
    u32 raw_mantissa = f32_as_u32 & ((1 << F32_MANTISSA_LENGTH) - 1);
    u32 raw_exponent = (f32_as_u32 >> F32_MANTISSA_LENGTH) & F32_EXPONENT_MASK;

if (raw_exponent == F32_EXPONENT_MASK) {
        if (raw_mantissa == 0) {
            sprintf(result, "%cinf", is_negative ? '-' : '+');
            return;
        }

sprintf(result, "NaN");
        return;
    }

if (raw_exponent == 0 && raw_mantissa == 0) {
        sprintf(result, "%c0.0", is_negative ? '-' : '+');
        return;
    }

u32 mantissa = (1 << F32_MANTISSA_LENGTH) + raw_mantissa;
    i32 exponent = (i32) raw_exponent - F32_EXPONENT_BIAS - F32_MANTISSA_LENGTH;
    if (raw_exponent == 0) {
        mantissa = raw_mantissa;
        exponent = 1 - F32_EXPONENT_BIAS - F32_MANTISSA_LENGTH;
    }

i32 binary_exponent = exponent - 2;

// Prev / current / next binary mantissa
    u32 u = 4 * mantissa - (raw_mantissa == 0 && raw_exponent > 1 ? 1 : 2);
    u32 v = 4 * mantissa;
    u32 w = 4 * mantissa + 2;

// Prev / current / next decimal mantissa
    u32 a_buffer[F32_DECIMAL_MANTISSA_BIGINT_SIZE] = {0};
    u32 a_swap_buffer[F32_DECIMAL_MANTISSA_BIGINT_SIZE] = {0};

u32 b_buffer[F32_DECIMAL_MANTISSA_BIGINT_SIZE] = {0};

u32 c_buffer[F32_DECIMAL_MANTISSA_BIGINT_SIZE] = {0};
    u32 c_swap_buffer[F32_DECIMAL_MANTISSA_BIGINT_SIZE] = {0};

u32 to_decimal_multiplier[F32_DECIMAL_MANTISSA_BIGINT_SIZE];
    if (binary_exponent >= 0) {
        bigint_power_of_two(
            binary_exponent, to_decimal_multiplier, ARRAY_SIZE(to_decimal_multiplier)
        );
    } else {
        bigint_u32_power(
            5, -binary_exponent, to_decimal_multiplier, ARRAY_SIZE(to_decimal_multiplier)
        );
    }

bigint_multiply(
        to_decimal_multiplier,
        ARRAY_SIZE(to_decimal_multiplier),
        &u,
        1,
        a_buffer,
        ARRAY_SIZE(a_buffer)
    );
    bigint_multiply(
        to_decimal_multiplier,
        ARRAY_SIZE(to_decimal_multiplier),
        &v,
        1,
        b_buffer,
        ARRAY_SIZE(b_buffer)
    );
    bigint_multiply(
        to_decimal_multiplier,
        ARRAY_SIZE(to_decimal_multiplier),
        &w,
        1,
        c_buffer,
        ARRAY_SIZE(c_buffer)
    );

i32 decimal_exponent = 0;
    if (binary_exponent < 0) {
        decimal_exponent += binary_exponent;
    }

bool all_a_zero = true, all_b_zero = true;
    u32 b_last_digit = 0;

if (config.accept_larger) {
        bigint_subtract_digit(c_buffer, ARRAY_SIZE(c_buffer), 1, c_buffer);
    }

u32 *a = a_buffer, *a_swap = a_swap_buffer;
    u32 *c = c_buffer, *c_swap = c_swap_buffer;

while (true) {
        u32 a_digit = bigint_divide_by_digit(a, ARRAY_SIZE(a_buffer), 10, a_swap);
        u32 c_digit = bigint_divide_by_digit(c, ARRAY_SIZE(c_buffer), 10, c_swap);

if (bigint_compare(a_swap, ARRAY_SIZE(a_buffer), c_swap, ARRAY_SIZE(c_buffer)) >= 0) {
            break;
        }

all_a_zero = all_a_zero && a_digit == 0;

b_last_digit = bigint_divide_by_digit(b_buffer, ARRAY_SIZE(b_buffer), 10, b_buffer);
        all_b_zero = all_b_zero && b_last_digit == 0;

SWAP(u32 *, a, a_swap);
        SWAP(u32 *, c, c_swap);

decimal_exponent += 1;
    }

if (config.accept_smaller && all_a_zero) {
        while (true) {
            u32 a_digit = bigint_divide_by_digit(a, ARRAY_SIZE(a_buffer), 10, a_swap);

if (a_digit != 0) {
                break;
            }

bigint_divide_by_digit(c, ARRAY_SIZE(c_buffer), 10, c);

b_last_digit = bigint_divide_by_digit(b_buffer, ARRAY_SIZE(b_buffer), 10, b_buffer);
            all_b_zero = all_b_zero && b_last_digit == 0;

SWAP(u32 *, a, a_swap);

decimal_exponent += 1;
        }
    }

bool is_tie = b_last_digit == 5 && all_b_zero;
    bool want_round_down = b_last_digit < 5 || is_tie && config.break_tie_down;

u32 b_plus_one[F32_DECIMAL_MANTISSA_BIGINT_SIZE];
    bigint_add_digit(b_buffer, ARRAY_SIZE(b_buffer), 1, b_plus_one);

bool round_down =
        want_round_down &&
            (bigint_compare(a, ARRAY_SIZE(a_buffer), b_buffer, ARRAY_SIZE(b_buffer)) != 0 ||
             all_b_zero) ||
        bigint_compare(b_plus_one, ARRAY_SIZE(b_plus_one), c, ARRAY_SIZE(c_buffer)) > 0;
    if (!round_down) {
        memcpy(b_buffer, b_plus_one, ARRAY_SIZE(b_buffer) * sizeof(u32));
    }

char decimal_mantissa_buffer[F32_DECIMAL_MANTISSA_MAX_DIGITS + 1] = {0};
    StringView decimal_mantissa = bigint_string_print(
        b_buffer,
        ARRAY_SIZE(b_buffer),
        decimal_mantissa_buffer,
        ARRAY_SIZE(decimal_mantissa_buffer) - 1
    );

if (decimal_exponent < -9) {
        sprintf(
            result,
            "%s1.%se%+lld",
            is_negative ? "-" : "",
            *(decimal_mantissa.data + 1) == '\0' ? "0" : decimal_mantissa.data + 1,
            decimal_exponent + decimal_mantissa.size - 1
        );
    } else if (decimal_exponent < 9) {
        if (decimal_exponent >= 0) {
            char *result_iter = result;
            result_iter +=
                sprintf(result_iter, "%s%s", is_negative ? "-" : "", decimal_mantissa.data);

for (isize i = 0; i < decimal_exponent; i += 1) {
                *result_iter = '0';
                result_iter += 1;
            }
        } else {
            char *result_iter = result;
            result_iter += sprintf(result_iter, "%s", is_negative ? "-" : "");

char *decimal_mantissa_iter = decimal_mantissa.data;
            for (isize i = 0; i < decimal_exponent + decimal_mantissa.size; i += 1) {
                *result_iter = *decimal_mantissa_iter;
                result_iter += 1;
                decimal_mantissa_iter += 1;
            }
            if (decimal_exponent + decimal_mantissa.size <= 0) {
                result_iter += sprintf(result_iter, "0.");
                for (isize i = 0; i < -decimal_exponent - decimal_mantissa.size; i += 1) {
                    *result_iter = '0';
                    result_iter += 1;
                }
            } else {
                result_iter += sprintf(result_iter, ".");
            }

if (decimal_mantissa_iter == decimal_mantissa.data + decimal_mantissa.size) {
                result_iter += sprintf(result_iter, "0");
            } else {
                while (decimal_mantissa_iter < decimal_mantissa.data + decimal_mantissa.size) {
                    *result_iter = *decimal_mantissa_iter;
                    result_iter += 1;
                    decimal_mantissa_iter += 1;
                }
            }
        }
    } else {
        sprintf(
            result, "%s%s.0e%+d", is_negative ? "-" : "", decimal_mantissa.data, decimal_exponent
        );
    }
}

i32 main(void) {
    FloatPrintConfig config = {
        .accept_smaller = true,
        .accept_larger = false,
        .break_tie_down = true,
    };

char const *test_cases[] = {
        "1e-46",
        "1e-45",
        "1.23461234e-40",
        "1e-30",
        "1e-20",
        "1e-10",
        "0.0",
        "-0.0",
        "1.23456789e-25",
        "1e-13",
        "0.0001",
        "0.123123123",
        "1.0",
        "12.12",
        "620.0",
        "1357.934142e-5",
        "-1291.0001",
        "-1392.123",
        "1e5",
        "21541251.1",
        "123456789.123456789",
        "2e+10",
        "9e+13",
        "1e+24",
        "1e30",
        "1.1234567e35",
        "3.19999e38",
        "1e40",
    };

char const **test_cases_iter = test_cases;
    while (test_cases_iter < test_cases + ARRAY_SIZE(test_cases)) {
        char result[32] = {0};
        f32_print(strtof(*test_cases_iter, NULL), config, result);

printf(
            "%s\t%20s => \"%s\"\n",
            strtof(result, NULL) == strtof(*test_cases_iter, NULL) ? "OK" : "FAIL",
            *test_cases_iter,
            result
        );

test_cases_iter += 1;
    }

return 0;
}