const enum3 = @import("enum3.zig").enum3; const enum3_data = @import("enum3.zig").enum3_data; const lookup_table = @import("lookup.zig").lookup_table; const HP = @import("lookup.zig").HP; const math = @import("../../math/index.zig"); const mem = @import("../../mem.zig"); const assert = @import("../../debug.zig").assert; pub const FloatDecimal = struct { digits: []u8, exp: i32, }; const u128 = @IntType(false, 128); /// Corrected Errol3 double to ASCII conversion. pub fn errol3(value: f64, buffer: []u8) -> FloatDecimal { const bits = @bitCast(u64, value); const i = tableLowerBound(bits); if (i < enum3.len and enum3[i] == bits) { const data = enum3_data[i]; const digits = buffer[0..data.str.len]; mem.copy(u8, digits, data.str); return FloatDecimal { .digits = digits, .exp = data.exp, }; } return errol3u(value, buffer); } /// Uncorrected Errol3 double to ASCII conversion. fn errol3u(val: f64, buffer: []u8) -> FloatDecimal { // check if in integer or fixed range if (val >= 9.007199254740992e15 and val < 3.40282366920938e+38) { return errolInt(val, buffer); } else if (val >= 16.0 and val < 9.007199254740992e15) { return errolFixed(val, buffer); } // normalize the midpoint var e: i32 = undefined; _ = math.frexp(val, &e); var exp = i16(math.floor(307 + f64(e) * 0.30103)); if (exp < 20) { exp = 20; } else if (usize(exp) >= lookup_table.len) { exp = i16(lookup_table.len - 1); } var mid = lookup_table[usize(exp)]; mid = hpProd(mid, val); const lten = lookup_table[usize(exp)].val; exp -= 307; var ten: f64 = 1.0; while (mid.val > 10.0 or (mid.val == 10.0 and mid.off >= 0.0)) { exp += 1; hpDiv10(&mid); ten /= 10.0; } while (mid.val < 1.0 or (mid.val == 1.0 and mid.off < 0.0)) { exp -= 1; hpMul10(&mid); ten *= 10.0; } // compute boundaries var high = HP { .val = mid.val, .off = mid.off + (fpnext(val) - val) * lten * ten / 2.0, }; var low = HP { .val = mid.val, .off = mid.off + (fpprev(val) - val) * lten * ten / 2.0, }; hpNormalize(&high); hpNormalize(&low); // normalized boundaries while (high.val > 10.0 or (high.val == 10.0 and high.off >= 0.0)) { exp += 1; hpDiv10(&high); hpDiv10(&low); } while (high.val < 1.0 or (high.val == 1.0 and high.off < 0.0)) { exp -= 1; hpMul10(&high); hpMul10(&low); } // digit generation var buf_index: usize = 0; while (true) { var hdig = u8(math.floor(high.val)); if ((high.val == f64(hdig)) and (high.off < 0)) hdig -= 1; var ldig = u8(math.floor(low.val)); if ((low.val == f64(ldig)) and (low.off < 0)) ldig -= 1; if (ldig != hdig) break; buffer[buf_index] = hdig + '0'; buf_index += 1; high.val -= f64(hdig); low.val -= f64(ldig); hpMul10(&high); hpMul10(&low); } const tmp = (high.val + low.val) / 2.0; var mdig = u8(math.floor(tmp + 0.5)); if ((f64(mdig) - tmp) == 0.5 and (mdig & 0x1) != 0) mdig -= 1; buffer[buf_index] = mdig + '0'; buf_index += 1; return FloatDecimal { .digits = buffer[0..buf_index], .exp = exp, }; } fn tableLowerBound(k: u64) -> usize { var i = enum3.len; var j: usize = 0; while (j < enum3.len) { if (enum3[j] < k) { j = 2 * k + 2; } else { i = j; j = 2 * j + 1; } } return i; } /// Compute the product of an HP number and a double. /// @in: The HP number. /// @val: The double. /// &returns: The HP number. fn hpProd(in: &const HP, val: f64) -> HP { var hi: f64 = undefined; var lo: f64 = undefined; split(in.val, &hi, &lo); var hi2: f64 = undefined; var lo2: f64 = undefined; split(val, &hi2, &lo2); const p = in.val * val; const e = ((hi * hi2 - p) + lo * hi2 + hi * lo2) + lo * lo2; return HP { .val = p, .off = in.off * val + e, }; } /// Split a double into two halves. /// @val: The double. /// @hi: The high bits. /// @lo: The low bits. fn split(val: f64, hi: &f64, lo: &f64) { *hi = gethi(val); *lo = val - *hi; } fn gethi(in: f64) { const bits = @bitCast(u64, in); const new_bits = bits & 0xFFFFFFFFF8000000; return @bitCast(f64, new_bits); } /// Normalize the number by factoring in the error. /// @hp: The float pair. fn hpNormalize(hp: &HP) { const val = hp.val; hp.val += hp.off; hp.off += val - hp.val; } /// Divide the high-precision number by ten. /// @hp: The high-precision number fn hpDiv10(hp: &HP) { const val = hp.val; hp.val /= 10.0; hp.off /= 10.0; val -= hp.val * 8.0; val -= hp.val * 2.0; hp.off += val / 10.0; hpNormalize(hp); } /// Multiply the high-precision number by ten. /// @hp: The high-precision number fn hpMul10(hp: &HP) { const val = hp.val; hp.val *= 10.0; hp.off *= 10.0; var off = hp.val; off -= val * 8.0; off -= val * 2.0; hp.off -= off; hpNormalize(hp); } /// Integer conversion algorithm, guaranteed correct, optimal, and best. /// @val: The val. /// @buf: The output buffer. /// &return: The exponent. fn errolInt(val: f64, buffer: []u8) -> FloatDecimal { const pow19 = 1e19; assert((val >= 9.007199254740992e15) and val < (3.40282366920938e38)); var mid = u128(val); var low: u128 = mid - fpeint((fpnext(val) - val) / 2.0); var high: u128 = mid + fpeint((val - fpprev(val)) / 2.0); if (@bitCast(u64, val) & 0x1 != 0) { high -= 1; } else { low -= 1; } const l64 = u64(low % pow19); const lf = u64((low / pow19) % pow19); const h64 = u64(high % pow19); const hf = u64((high / pow19) % pow19); if (lf != hf) { l64 = lf; h64 = hf; mid = mid / (pow19 / 10); } var mi: i32 = mismatch10(l64, h64); var x: u64 = 1; { var i = i32(lf == hf); while (i < mi) : (i += 1) { x *= 10; } } const m64: u64 = mid / x; if (lf != hf) mi += 19; var buf_index = u64toa(m64, buffer) - 1; if (mi != 0) { buffer[buf_index - 1] += (buffer[buf_index] >= '5'); } else { buf_index += 1; } return FloatDecimal { .digits = buffer[0..buf_index], .exp = buf_index + mi, }; } /// Fixed point conversion algorithm, guaranteed correct, optimal, and best. /// @val: The val. /// @buf: The output buffer. /// &return: The exponent. fn errolFixed(val: f64, buffer: []u8) -> FloatDecimal { assert((val >= 16.0) and (val < 9.007199254740992e15)); const u = u64(val); const n = f64(u); var mid = val - n; var lo = ((fpprev(val) - n) + mid) / 2.0; var hi = ((fpnext(val) - n) + mid) / 2.0; var buf_index = u64toa(u, buffer); var exp: i32 = buf_index; var j: i32 = exp; buffer[j] = 0; if (mid != 0.0) { while (mid != 0.0) { lo *= 10.0; var ldig = i32(lo); lo -= ldig; mid *= 10.0; var mdig = i32(mid); mid -= mdig; hi *= 10.0; var hdig = i32(hi); hi -= hdig; buffer[j] = mdig + '0'; j += 1; if(hdig != ldig or j > 50) break; } if (mid > 0.5) { buffer[j-1] += 1; } else if ((mid == 0.5) and (buffer[j-1] & 0x1)) { buffer[j-1] += 1; } } else { while (buffer[j-1] == '0') { buffer[j-1] = 0; j -= 1; } } buffer[j] = 0; return FloatDecimal { .digits = buffer[0..j], .exp = exp, }; } fn fpnext(val: f64) -> f64 { return @bitCast(f64, @bitCast(u64, val) + 1); } fn fpprev(val: f64) -> f64 { return @bitCast(f64, @bitCast(u64, val) - 1); } pub const c_digits_lut = []u8 { '0', '0', '0', '1', '0', '2', '0', '3', '0', '4', '0', '5', '0', '6', '0', '7', '0', '8', '0', '9', '1', '0', '1', '1', '1', '2', '1', '3', '1', '4', '1', '5', '1', '6', '1', '7', '1', '8', '1', '9', '2', '0', '2', '1', '2', '2', '2', '3', '2', '4', '2', '5', '2', '6', '2', '7', '2', '8', '2', '9', '3', '0', '3', '1', '3', '2', '3', '3', '3', '4', '3', '5', '3', '6', '3', '7', '3', '8', '3', '9', '4', '0', '4', '1', '4', '2', '4', '3', '4', '4', '4', '5', '4', '6', '4', '7', '4', '8', '4', '9', '5', '0', '5', '1', '5', '2', '5', '3', '5', '4', '5', '5', '5', '6', '5', '7', '5', '8', '5', '9', '6', '0', '6', '1', '6', '2', '6', '3', '6', '4', '6', '5', '6', '6', '6', '7', '6', '8', '6', '9', '7', '0', '7', '1', '7', '2', '7', '3', '7', '4', '7', '5', '7', '6', '7', '7', '7', '8', '7', '9', '8', '0', '8', '1', '8', '2', '8', '3', '8', '4', '8', '5', '8', '6', '8', '7', '8', '8', '8', '9', '9', '0', '9', '1', '9', '2', '9', '3', '9', '4', '9', '5', '9', '6', '9', '7', '9', '8', '9', '9', }; fn u64toa(value: u64, buffer: []u8) -> usize { const kTen8: u64 = 100000000; const kTen9: u64 = kTen8 * 10; const kTen10: u64 = kTen8 * 100; const kTen11: u64 = kTen8 * 1000; const kTen12: u64 = kTen8 * 10000; const kTen13: u64 = kTen8 * 100000; const kTen14: u64 = kTen8 * 1000000; const kTen15: u64 = kTen8 * 10000000; const kTen16: u64 = kTen8 * kTen8; var buf_index: usize = 0; if (value < kTen8) { const v = u32(value); if (v < 10000) { const d1: u32 = (v / 100) << 1; const d2: u32 = (v % 100) << 1; if (v >= 1000) { buffer[buf_index] = c_digits_lut[d1]; buf_index += 1; } if (v >= 100) { buffer[buf_index] = c_digits_lut[d1 + 1]; buf_index += 1; } if (v >= 10) { buffer[buf_index] = c_digits_lut[d2]; buf_index += 1; } buffer[buf_index] = c_digits_lut[d2 + 1]; buf_index += 1; } else { // value = bbbbcccc const b: u32 = v / 10000; const c: u32 = v % 10000; const d1: u32 = (b / 100) << 1; const d2: u32 = (b % 100) << 1; const d3: u32 = (c / 100) << 1; const d4: u32 = (c % 100) << 1; if (value >= 10000000) { buffer[buf_index] = c_digits_lut[d1]; buf_index += 1; } if (value >= 1000000) { buffer[buf_index] = c_digits_lut[d1 + 1]; buf_index += 1; } if (value >= 100000) { buffer[buf_index] = c_digits_lut[d2]; buf_index += 1; } buffer[buf_index] = c_digits_lut[d2 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d3]; buf_index += 1; buffer[buf_index] = c_digits_lut[d3 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d4]; buf_index += 1; buffer[buf_index] = c_digits_lut[d4 + 1]; buf_index += 1; } } else if (value < kTen16) { const v0: u32 = (uint32_t)(value / kTen8); const v1: u32 = (uint32_t)(value % kTen8); const b0: u32 = v0 / 10000; const c0: u32 = v0 % 10000; const d1: u32 = (b0 / 100) << 1; const d2: u32 = (b0 % 100) << 1; const d3: u32 = (c0 / 100) << 1; const d4: u32 = (c0 % 100) << 1; const b1: u32 = v1 / 10000; const c1: u32 = v1 % 10000; const d5: u32 = (b1 / 100) << 1; const d6: u32 = (b1 % 100) << 1; const d7: u32 = (c1 / 100) << 1; const d8: u32 = (c1 % 100) << 1; if (value >= kTen15) { buffer[buf_index] = c_digits_lut[d1]; buf_index += 1; } if (value >= kTen14) { buffer[buf_index] = c_digits_lut[d1 + 1]; buf_index += 1; } if (value >= kTen13) { buffer[buf_index] = c_digits_lut[d2]; buf_index += 1; } if (value >= kTen12) { buffer[buf_index] = c_digits_lut[d2 + 1]; buf_index += 1; } if (value >= kTen11) { buffer[buf_index] = c_digits_lut[d3]; buf_index += 1; } if (value >= kTen10) { buffer[buf_index] = c_digits_lut[d3 + 1]; buf_index += 1; } if (value >= kTen9) { buffer[buf_index] = c_digits_lut[d4]; buf_index += 1; } if (value >= kTen8) { buffer[buf_index] = c_digits_lut[d4 + 1]; buf_index += 1; } buffer[buf_index] = c_digits_lut[d5]; buf_index += 1; buffer[buf_index] = c_digits_lut[d5 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d6]; buf_index += 1; buffer[buf_index] = c_digits_lut[d6 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d7]; buf_index += 1; buffer[buf_index] = c_digits_lut[d7 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d8]; buf_index += 1; buffer[buf_index] = c_digits_lut[d8]; buf_index += 1; buffer[buf_index] = c_digits_lut[d8]; buf_index += 1; buffer[buf_index] = c_digits_lut[d8 + 1]; buf_index += 1; } else { const a: u32 = (uint32_t)(value / kTen16); // 1 to 1844 value %= kTen16; if (a < 10) { buffer[buf_index] = (char)('0' + (char)(a)); buf_index += 1; } else if (a < 100) { const i: u32 = a << 1; buffer[buf_index] = c_digits_lut[i]; buf_index += 1; buffer[buf_index] = c_digits_lut[i + 1]; buf_index += 1; } else if (a < 1000) { buffer[buf_index] = (char)('0' + (char)(a / 100)); buf_index += 1; const i: u32 = (a % 100) << 1; buffer[buf_index] = c_digits_lut[i]; buf_index += 1; buffer[buf_index] = c_digits_lut[i + 1]; buf_index += 1; } else { const i: u32 = (a / 100) << 1; const j: u32 = (a % 100) << 1; buffer[buf_index] = c_digits_lut[i]; buf_index += 1; buffer[buf_index] = c_digits_lut[i + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[j]; buf_index += 1; buffer[buf_index] = c_digits_lut[j + 1]; buf_index += 1; } const v0: u32 = (uint32_t)(value / kTen8); const v1: u32 = (uint32_t)(value % kTen8); const b0: u32 = v0 / 10000; const c0: u32 = v0 % 10000; const d1: u32 = (b0 / 100) << 1; const d2: u32 = (b0 % 100) << 1; const d3: u32 = (c0 / 100) << 1; const d4: u32 = (c0 % 100) << 1; const b1: u32 = v1 / 10000; const c1: u32 = v1 % 10000; const d5: u32 = (b1 / 100) << 1; const d6: u32 = (b1 % 100) << 1; const d7: u32 = (c1 / 100) << 1; const d8: u32 = (c1 % 100) << 1; buffer[buf_index] = c_digits_lut[d1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d1 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d2]; buf_index += 1; buffer[buf_index] = c_digits_lut[d2 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d3]; buf_index += 1; buffer[buf_index] = c_digits_lut[d3 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d4]; buf_index += 1; buffer[buf_index] = c_digits_lut[d4 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d5]; buf_index += 1; buffer[buf_index] = c_digits_lut[d5 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d6]; buf_index += 1; buffer[buf_index] = c_digits_lut[d6 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d7]; buf_index += 1; buffer[buf_index] = c_digits_lut[d7 + 1]; buf_index += 1; buffer[buf_index] = c_digits_lut[d8]; buf_index += 1; buffer[buf_index] = c_digits_lut[d8 + 1]; buf_index += 1; } return buf_index; } fn fpeint(from: f64) -> u128 { const bits = @bitCast(u64, from); assert((bits & ((1 << 52) - 1)) == 0); return 1 << ((bits >> 52) - 1023); } /// Given two different integers with the same length in terms of the number /// of decimal digits, index the digits from the right-most position starting /// from zero, find the first index where the digits in the two integers /// divergent starting from the highest index. /// @a: Integer a. /// @b: Integer b. /// &returns: An index within [0, 19). fn mismatch10(a: u64, b: u64) -> i32 { const pow10 = 10000000000; const af = a / pow10; const bf = b / pow10; var i: i32 = 0; var a_copy = a; var b_copy = b; if (af != bf) { i = 10; a_copy = af; b_copy = bf; } while (true) : (i += 1) { a_copy /= 10; b_copy /= 10; if (a_copy == b_copy) return i; } }