bcd8.rs

use std::iter;

#[derive(Copy, Clone)]
pub struct Bcd8<const N: usize>(pub(super) [u8; N]);

impl<const N: usize> Bcd8<N> {
    pub const fn zeros() -> Self {
        Self([0; N])
    }

    pub fn digits(&self) -> impl Iterator<Item = u8> {
        self.0.into_iter()
    }
}

impl<const N: usize> std::ops::AddAssign<&Bcd8<N>> for Bcd8<N> {
    fn add_assign(&mut self, rhs: &Bcd8<N>) {
        let mut carry = 0;
        for (a, b) in iter::zip(&mut self.0, &rhs.0).rev() {
            let sum = *a + *b + carry;
            carry = sum / 10;
            *a = sum % 10;
        }
        // overflow is impossible in this problem domain
        debug_assert_eq!(carry, 0);
    }
}

const_friendly_u256.rs

#[derive(Debug, Default, Copy, Clone, PartialEq, Eq)]
#[allow(non_camel_case_types)]
pub struct u256 {
    pub hi: u128,
    pub lo: u128,
}

impl u256 {
    pub const ONE: Self = Self { hi: 0, lo: 1 };

    pub const fn add(self, rhs: Self) -> Self {
        let (lo, carry) = self.lo.overflowing_add(rhs.lo);
        let hi = self.hi + rhs.hi + carry as u128;
        Self { lo, hi }
    }

    pub const fn mul_5(self) -> Self {
        let two = self.add(self);
        let four = two.add(two);
        four.add(self)
    }

    pub const fn div_rem_10(self) -> (Self, u8) {
        let hi = self.hi;
        let lo_a = self.lo >> 64;
        let lo_b = (self.lo << 64) >> 64;

        let (q_hi, r) = dr10(hi, 0);
        let (q_lo_a, r) = dr10(lo_a, r);
        let (q_lo_b, r) = dr10(lo_b, r);

        let q_lo = (q_lo_a << 64) | q_lo_b;
        (Self { hi: q_hi, lo: q_lo }, r as u8)
    }
}

const fn dr10(x: u128, rem: u128) -> (u128, u128) {
    let y = (rem << 64) | x;
    (y / 10, y % 10)
}

main.rs

use rand::Rng;
use std::fmt;

pub mod bcd8;
use bcd8::Bcd8;

pub mod const_friendly_u256;
use const_friendly_u256::u256;

fn main() {
    let mut n = 1 << 63;
    for _ in 0..64 {
        print_n(n);
        n >>= 1;
    }

    print_n(0);
    print_n(u64::MAX);

    print_n(rand::thread_rng().gen());
}

fn print_n(n: u64) {
    println!("{n:064b}");
    println!("{}", Frac64(n));
    println!();
}

#[derive(Copy, Clone, Default, Hash, PartialEq, Eq, PartialOrd, Ord)]
pub struct Frac64(u64);

impl fmt::Display for Frac64 {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        if self.0 == 0 {
            return f.write_str("0.0");
        }

        f.write_str("0.")?;
        let bcd = self.to_bcd();
        let len = 64 - self.0.trailing_zeros(); // it just works
        for digit in bcd.digits().take(len as usize) {
            write!(f, "{digit}")?;
        }

        Ok(())
    }
}

impl Frac64 {
    fn to_bcd(&self) -> Bcd8<64> {
        let mut n = self.0;

        let mut bcd = Bcd8::zeros();

        while n != 0 {
            let i = n.leading_zeros(); // position of the least significant 1
            n ^= 1 << (63 - i);
            bcd += &BCD_TABLE[i as usize];
        }

        bcd
    }
}

// in decimal, the negative powers of two look like the positive powers of five,
// but shifted right according to the magnitude of the exponent
//
// we can take advantage of this to generate the table in a const context

static BCD_TABLE: [Bcd8<64>; 64] = gen_bcd_table();

const fn gen_bcd_table() -> [Bcd8<64>; 64] {
    let mut bufs = [Bcd8::zeros(); 64];
    let mut n = u256::ONE;

    let mut i = 0;
    while i < 64 {
        n = n.mul_5();
        bufs[i] = gen_bcd_row(n, i);
        i += 1;
    }

    bufs
}

const fn gen_bcd_row(mut n: u256, i: usize) -> Bcd8<64> {
    let mut buf = [0; 64];

    let mut j = 0;
    while j <= i {
        let (q, r) = n.div_rem_10();
        buf = rotate_right(buf);
        buf[0] = r;
        n = q;
        j += 1;
    }

    Bcd8(buf)
}

const fn rotate_right(x: [u8; 64]) -> [u8; 64] {
    let mut y = [0; 64];
    y[0] = x[63];
    let mut i = 0;
    while i < 63 {
        y[i + 1] = x[i];
        i += 1;
    }
    y
}