diff --git a/vlib/math/bits/bits.v b/vlib/math/bits/bits.v index 184e7a35bc..5b37855750 100644 --- a/vlib/math/bits/bits.v +++ b/vlib/math/bits/bits.v @@ -96,6 +96,7 @@ pub fn trailing_zeros_64(x u64) int { } // --- OnesCount --- + // ones_count_8 returns the number of one bits ("population count") in x. pub fn ones_count_8(x byte) int { return int(pop_8_tab[x]) @@ -142,6 +143,7 @@ pub fn ones_count_64(x u64) int { } // --- RotateLeft --- + // rotate_left_8 returns the value of x rotated left by (k mod 8) bits. // To rotate x right by k bits, call rotate_left_8(x, -k). // @@ -187,6 +189,7 @@ pub fn rotate_left_64(x u64, k int) u64 { } // --- Reverse --- + // reverse_8 returns the value of x with its bits in reversed order. [inline] pub fn reverse_8(x byte) byte { @@ -218,6 +221,7 @@ pub fn reverse_64(x u64) u64 { } // --- ReverseBytes --- + // reverse_bytes_16 returns the value of x with its bytes in reversed order. // // This function's execution time does not depend on the inputs. @@ -246,6 +250,7 @@ pub fn reverse_bytes_64(x u64) u64 { } // --- Len --- + // len_8 returns the minimum number of bits required to represent x; the result is 0 for x == 0. pub fn len_8(x byte) int { return int(len_8_tab[x]) @@ -296,3 +301,195 @@ pub fn len_64(x u64) int { return n + int(len_8_tab[y]) } +// --- Add with carry --- + +// Add returns the sum with carry of x, y and carry: sum = x + y + carry. +// The carry input must be 0 or 1; otherwise the behavior is undefined. +// The carryOut output is guaranteed to be 0 or 1. +// + +// add_32 returns the sum with carry of x, y and carry: sum = x + y + carry. +// The carry input must be 0 or 1; otherwise the behavior is undefined. +// The carryOut output is guaranteed to be 0 or 1. +// +// This function's execution time does not depend on the inputs. +fn add_32(x u32, y u32, carry u32) (u32, u32) { + sum64 := u64(x) + u64(y) + u64(carry) + sum := u32(sum64) + carry_out := u32(sum64>>32) + return sum, carry_out +} + +// add_64 returns the sum with carry of x, y and carry: sum = x + y + carry. +// The carry input must be 0 or 1; otherwise the behavior is undefined. +// The carryOut output is guaranteed to be 0 or 1. +// +// This function's execution time does not depend on the inputs. +fn add_64(x u64, y u64, carry u64) (u64, u64) { + sum := x + y + carry + // The sum will overflow if both top bits are set (x & y) or if one of them + // is (x | y), and a carry from the lower place happened. If such a carry + // happens, the top bit will be 1 + 0 + 1 = 0 (&^ sum). + carry_out := ((x & y) | ((x | y) & ~sum ))>>63 + return sum, carry_out +} + +// --- Subtract with borrow --- + +// Sub returns the difference of x, y and borrow: diff = x - y - borrow. +// The borrow input must be 0 or 1; otherwise the behavior is undefined. +// The borrowOut output is guaranteed to be 0 or 1. +// + +// sub_32 returns the difference of x, y and borrow, diff = x - y - borrow. +// The borrow input must be 0 or 1; otherwise the behavior is undefined. +// The borrowOut output is guaranteed to be 0 or 1. +// +// This function's execution time does not depend on the inputs. +fn sub_32(x u32, y u32, borrow u32) (u32, u32) { + diff := x - y - borrow + // The difference will underflow if the top bit of x is not set and the top + // bit of y is set (^x & y) or if they are the same (^(x ^ y)) and a borrow + // from the lower place happens. If that borrow happens, the result will be + // 1 - 1 - 1 = 0 - 0 - 1 = 1 (& diff). + borrow_out := ((~x & y) | (~(x ^ y) & diff))>>31 + return diff, borrow_out +} + +// sub_64 returns the difference of x, y and borrow: diff = x - y - borrow. +// The borrow input must be 0 or 1; otherwise the behavior is undefined. +// The borrowOut output is guaranteed to be 0 or 1. +// +// This function's execution time does not depend on the inputs. +fn sub_64(x u64, y u64, borrow u64) (u64, u64) { + diff := x - y - borrow + // See Sub32 for the bit logic. + borrow_out := ((~x & y) | (~(x ^ y) & diff))>>63 + return diff, borrow_out +} + +// --- Full-width multiply --- + +const ( + two32 = u64(0x1_0000_0000) + mask32 = two32 - 1 + overflow_error = "Overflow Error" + divide_error = "Divide Error" +) + +// mul_32 returns the 64-bit product of x and y: (hi, lo) = x * y +// with the product bits' upper half returned in hi and the lower +// half returned in lo. +// +// This function's execution time does not depend on the inputs. +fn mul_32(x u32, y u32) (u32, u32) { + tmp := u64(x) * u64(y) + hi := u32(tmp>>32) + lo := u32(tmp) + return hi, lo +} + +// mul_64 returns the 128-bit product of x and y: (hi, lo) = x * y +// with the product bits' upper half returned in hi and the lower +// half returned in lo. +// +// This function's execution time does not depend on the inputs. +fn mul_64(x u64, y u64) (u64, u64) { + x0 := x & mask32 + x1 := x>>32 + y0 := y & mask32 + y1 := y>>32 + w0 := x0 * y0 + t := x1*y0 + (w0>>32) + mut w1 := t & mask32 + w2 := t>>32 + w1 += x0 * y1 + hi := x1*y1 + w2 + (w1>>32) + lo := x * y + return hi, lo +} + +// --- Full-width divide --- + +// div_32 returns the quotient and remainder of (hi, lo) divided by y: +// quo = (hi, lo)/y, rem = (hi, lo)%y with the dividend bits' upper +// half in parameter hi and the lower half in parameter lo. +// div_32 panics for y == 0 (division by zero) or y <= hi (quotient overflow). +fn div_32(hi u32, lo u32, y u32) (u32, u32) { + if y != 0 && y <= hi { + panic(overflow_error) + } + z := (u64(hi)<<32) | u64(lo) + quo := u32(z/u64(y)) + rem := u32(z%u64(y)) + return quo, rem +} + +// div_64 returns the quotient and remainder of (hi, lo) divided by y: +// quo = (hi, lo)/y, rem = (hi, lo)%y with the dividend bits' upper +// half in parameter hi and the lower half in parameter lo. +// div_64 panics for y == 0 (division by zero) or y <= hi (quotient overflow). +fn div_64(hi u64, lo u64, y1 u64) (u64, u64) { + mut y := y1 + if y == 0 { + panic(overflow_error) + } + if y <= hi { + panic(overflow_error) + } + + s := u32(leading_zeros_64(y)) + y <<= s + + yn1 := y>>32 + yn0 := y & mask32 + un32 := (hi<>(64-s)) + un10 := lo<>32 + un0 := un10 & mask32 + mut q1 := un32 / yn1 + mut rhat := un32 - q1*yn1 + + for q1 >= two32 || q1*yn0 > two32*rhat+un1 { + q1-- + rhat += yn1 + if rhat >= two32 { + break + } + } + + un21 := un32*two32 + un1 - q1*y + mut q0 := un21 / yn1 + rhat = un21 - q0*yn1 + + for q0 >= two32 || q0*yn0 > two32*rhat+un0 { + q0-- + rhat += yn1 + if rhat >= two32 { + break + } + } + + return q1*two32 + q0, (un21*two32 + un0 - q0*y)>>s +} + +// rem_32 returns the remainder of (hi, lo) divided by y. Rem32 panics +// for y == 0 (division by zero) but, unlike Div32, it doesn't panic +// on a quotient overflow. +fn rem_32(hi u32, lo u32, y u32) u32 { + return u32((u64(hi)<<32 | u64(lo)) % u64(y)) +} + +// rem_64 returns the remainder of (hi, lo) divided by y. Rem64 panics +// for y == 0 (division by zero) but, unlike div_64, it doesn't panic +// on a quotient overflow. +fn rem_64(hi, lo, y u64) u64 { + // We scale down hi so that hi < y, then use div_64 to compute the + // rem with the guarantee that it won't panic on quotient overflow. + // Given that + // hi ≡ hi%y (mod y) + // we have + // hi<<64 + lo ≡ (hi%y)<<64 + lo (mod y) + _, rem := div_64(hi%y, lo, y) + return rem +} diff --git a/vlib/math/bits/bits_test.v b/vlib/math/bits/bits_test.v index d808920e94..79eb059d63 100644 --- a/vlib/math/bits/bits_test.v +++ b/vlib/math/bits/bits_test.v @@ -1,8 +1,12 @@ +// +// test suite for bits and bits math functions +// module bits fn test_bits(){ mut i := 0 mut i1:= u64(0) + // // --- LeadingZeros --- // @@ -10,26 +14,29 @@ fn test_bits(){ // 8 bit i = 1 for x in 0..8 { - //C.printf("x:%02x lz: %d cmp: %d\n",i<> 1 } - //C.printf("x:%02x lz: %llu cmp: %d\n",byte(i),reverse_8(byte(i)), rv) + //C.printf("x:%02x lz: %llu cmp: %d\n", byte(i), reverse_8(byte(i)), rv) assert reverse_8(byte(i)) == rv i = (i << 1) + 1 } + // 16 bit i = 0 for x in 0..17 { @@ -103,10 +114,11 @@ fn test_bits(){ bc++ n = n >> 1 } - //C.printf("x:%04x lz: %llu cmp: %d\n",u16(i),reverse_16(u16(i)), rv) + //C.printf("x:%04x lz: %llu cmp: %d\n", u16(i), reverse_16(u16(i)), rv) assert reverse_16(u16(i)) == rv i = (i << 1) + 1 } + // 32 bit i = 0 for x in 0..33 { @@ -118,10 +130,11 @@ fn test_bits(){ bc++ n = n >> 1 } - //C.printf("x:%08x lz: %llu cmp: %d\n",u32(i),reverse_32(u32(i)), rv) + //C.printf("x:%08x lz: %llu cmp: %d\n", u32(i), reverse_32(u32(i)), rv) assert reverse_32(u32(i)) == rv i = (i << 1) + 1 } + // 64 bit i1 = 0 for x in 0..64 { @@ -133,8 +146,144 @@ fn test_bits(){ bc++ n = n >> 1 } - //C.printf("x:%016llx lz: %016llx cmp: %016llx\n",u64(i1),reverse_64(u64(i1)), rv) - assert reverse_64(u64(i1)) == rv + //C.printf("x:%016llx lz: %016llx cmp: %016llx\n", u64(i1), reverse_64(u64(i1)), rv) + assert reverse_64(i1) == rv i1 = (i1 << 1) + 1 } -} \ No newline at end of file + + // + // --- add --- + // + + // 32 bit + i = 1 + for x in 0..32 { + v := u32(i) << x + sum,carry := add_32(v, v, u32(0)) + //C.printf("x:%08x [%llu,%llu] %llu\n", u32(i) << x, sum, carry, u64(v) + u64(v)) + assert ((u64(carry) << 32) | u64(sum)) == u64(v) + u64(v) + } + mut sum_32t, mut carry_32t := add_32(0x8000_0000, 0x8000_0000, u32(0)) + assert sum_32t == u32(0) + assert carry_32t == u32(1) + + sum_32t, carry_32t = add_32(0xFFFF_FFFF, 0xFFFF_FFFF, u32(1)) + assert sum_32t == 0xFFFF_FFFF + assert carry_32t == u32(1) + + // 64 bit + i = 1 + for x in 0..63 { + v := u64(i) << x + sum,carry := add_64(v, v, u64(0)) + //C.printf("x:%16x [%llu,%llu] %llu\n", u64(i) << x, sum, carry, u64(v >> 32) + u64(v >> 32)) + assert ((carry << 32) | sum) == v + v + } + mut sum_64t, mut carry_64t := add_64(0x8000_0000_0000_0000, 0x8000_0000_0000_0000, u64(0)) + assert sum_64t == u64(0) + assert carry_64t == u64(1) + + sum_64t, carry_64t = add_64(0xFFFF_FFFF_FFFF_FFFF, 0xFFFF_FFFF_FFFF_FFFF, u64(1)) + assert sum_64t == 0xFFFF_FFFF_FFFF_FFFF + assert carry_64t == u64(1) + + // + // --- sub --- + // + + // 32 bit + i = 1 + for x in 1..32 { + v0 := u32(i) << x + v1 := v0 >> 1 + mut diff, mut borrow_out := sub_32(v0, v1, u32(0)) + //C.printf("x:%08x [%llu,%llu] %08x\n", u32(i) << x, diff, borrow_out, v0 - v1) + assert diff == v1 + + diff, borrow_out = sub_32(v0, v1, u32(1)) + //C.printf("x:%08x [%llu,%llu] %08x\n", u32(i) << x, diff, borrow_out, v0 - v1) + assert diff == (v1 - 1) + assert borrow_out == u32(0) + + diff, borrow_out = sub_32(v1, v0, u32(1)) + //C.printf("x:%08x [%llu,%llu] %08x\n", u32(i) << x, diff, borrow_out, v1 - v0) + assert borrow_out == u32(1) + } + + // 64 bit + i = 1 + for x in 1..64 { + v0 := u64(i) << x + v1 := v0 >> 1 + mut diff, mut borrow_out := sub_64(v0, v1, u64(0)) + //C.printf("x:%08x [%llu,%llu] %08x\n", u64(i) << x, diff, borrow_out, v0 - v1) + assert diff == v1 + + diff, borrow_out = sub_64(v0, v1, u64(1)) + //C.printf("x:%08x [%llu,%llu] %08x\n", u64(i) << x, diff, borrow_out, v0 - v1) + assert diff == (v1 - 1) + assert borrow_out == u64(0) + + diff, borrow_out = sub_64(v1, v0, u64(1)) + //C.printf("x:%08x [%llu,%llu] %08x\n",u64(i) << x, diff, borrow_out, v1 - v0) + assert borrow_out == u64(1) + } + + // + // --- mul --- + // + + // 32 bit + i = 1 + for x in 0..32 { + v0 := u32(i) << x + v1 := v0 - 1 + hi, lo := mul_32(v0, v1) + //C.printf("x:%08x [%llu,%llu] %llu\n", v0, hi, lo, u64(v0 * v1)) + assert (u64(hi) << 32) | (u64(lo)) == u64(v0 * v1) + } + + // 64 bit + i = 1 + for x in 0..64 { + v0 := u64(i) << x + v1 := v0 - 1 + hi, lo := mul_64(v0, v1) + //C.printf("v0: %llu v1: %llu [%llu,%llu] tt: %llu\n", v0, v1, hi, lo, (v0 >> 32) * (v1 >> 32)) + assert (hi & 0xFFFF_FFFF_0000_0000) == (((v0 >> 32)*(v1 >> 32)) & 0xFFFF_FFFF_0000_0000) + assert (lo & 0x0000_0000_FFFF_FFFF) == (((v0 & 0x0000_0000_FFFF_FFFF) * (v1 & 0x0000_0000_FFFF_FFFF)) & 0x0000_0000_FFFF_FFFF) + } + + // + // --- div --- + // + + // 32 bit + i = 1 + for x in 0..31 { + hi := u32(i) << x + lo := hi - 1 + y := u32(3) << x + quo, rem := div_32(hi, lo, y) + //C.printf("[%08x_%08x] %08x (%08x,%08x)\n", hi, lo, y, quo, rem) + tst := ((u64(hi) << 32) | u64(lo)) + assert quo == (tst / u64(y)) + assert rem == (tst % u64(y)) + assert rem == rem_32(hi, lo, y) + } + + // 64 bit + i = 1 + for x in 0..62 { + hi := u64(i) << x + lo := u64(2) //hi - 1 + y := 0x4000_0000_0000_0000 + quo, rem := div_64(hi, lo, y) + //C.printf("[%016llx_%016llx] %016llx (%016llx,%016llx)\n", hi, lo, y, quo, rem) + assert quo == u64(2)<<(x+1) + _, rem1 := div_64(hi%y, lo, y) + assert rem == rem1 + assert rem == rem_64(hi, lo, y) + } + +}