From 4d6017b97493ebd3ab28b4143184d899a1c84d47 Mon Sep 17 00:00:00 2001 From: EvilMuffinHa Date: Sat, 6 Aug 2022 00:08:25 -0400 Subject: [PATCH] karatsuba multiplication: doesn't work yet --- src/ops.rs | 313 ++++++++++++++++++++++++++++++++++++++++++--------- src/tests.rs | 46 +++++++- 2 files changed, 303 insertions(+), 56 deletions(-) diff --git a/src/ops.rs b/src/ops.rs index cb50c2e..d197bc5 100644 --- a/src/ops.rs +++ b/src/ops.rs @@ -1,6 +1,6 @@ use crate::SmallUint; use crate::smallint::SmallUintType; -use core::ops::Add; +use core::ops::{Add, Sub, Mul}; use core::mem::ManuallyDrop; macro_rules! basic_op { @@ -83,12 +83,13 @@ fn add_two_slices(slice1: &[u32], slice2: &[u32]) -> Vec { if carry { res.push(1); } - while res[res.len() - 1] == 0 { + while res.len() != 1 && res[res.len() - 1] == 0 { res.pop(); } res } + fn add(a: &SmallUint, b: &SmallUint) -> SmallUint { match (&a.0, &b.0) { @@ -153,6 +154,103 @@ fn add(a: &SmallUint, b: &SmallUint) -> SmallUint { basic_op!(Add, SmallUint, add); +fn sub_two_slices(slice1: &[u32], slice2: &[u32]) -> Vec { + let b = slice1.len(); + let s = slice2.len(); + + if b < s { + panic!("First number is smaller than second."); + } + + let mut res = Vec::with_capacity(std::cmp::max(s, b)); + let mut borrow = false; + + for i in 0..b { + let mut value1 = slice1[i]; + + let value2 = if i < s { + slice2[i] + } else { + 0 + }; + + if borrow { + let (temp, b) = value1.overflowing_sub(1); + value1 = temp; + borrow = b; + } + + if value2 > value1 { + borrow = true; + } + + let val = value1.wrapping_sub(value2); + res.push(val); + + } + + if borrow { + panic!("First number is smaller than second. "); + } + + res + +} + +fn sub(a: &SmallUint, b: &SmallUint) -> SmallUint { + + match (&a.0, &b.0) { + (&SmallUintType::Inline(i), &SmallUintType::Inline(j)) => { + if let (t, false) = i.overflowing_sub(j) { + SmallUint(SmallUintType::Inline(t)) + } else { + panic!("First number is smaller than second. "); + } + }, + (&SmallUintType::Heap((r, s)), &SmallUintType::Inline(i)) => { + let slice1 = unsafe { core::slice::from_raw_parts(r, s) }; + + let mut res = [0, 0, 0, 0]; + + let mut v = i; + #[allow(clippy::needless_range_loop)] + for r in 0..4 { + res[r] = v as u32; + + v >>= 32; + } + + let result = sub_two_slices(slice1, &res[..]); + let size = result.len(); + + let mut slice = ManuallyDrop::new(result.into_boxed_slice()); + + SmallUint(SmallUintType::Heap((slice.as_mut_ptr(), size))) + + + }, + (&SmallUintType::Inline(_), &SmallUintType::Heap((_, _))) => { + panic!("First number is smaller than second. "); + + }, + (&SmallUintType::Heap((r, s)), &SmallUintType::Heap((i, j))) => { + + let slice1 = unsafe { core::slice::from_raw_parts(r, s) }; + let slice2 = unsafe { core::slice::from_raw_parts(i, j) }; + + let res = sub_two_slices(slice1, slice2); + let size = res.len(); + + let mut slice = ManuallyDrop::new(res.into_boxed_slice()); + + SmallUint(SmallUintType::Heap((slice.as_mut_ptr(), size))) + + } + } +} + +basic_op!(Sub, SmallUint, sub); + // Taken from https://github.com/rust-lang/rust/issues/85532#issuecomment-916309635. Credit to // AaronKutch. const fn carrying_mul_u128(lhs: u128, rhs: u128, carry: u128) -> (u128, u128) { @@ -189,55 +287,162 @@ const fn carrying_mul_u128(lhs: u128, rhs: u128, carry: u128) -> (u128, u128) { (sum0, sum1) } -// fn mul_two_slices(slice1: &[u32], slice2: &[u32]) -> Vec { -// -// let l1 = slice1.len(); -// let l2 = slice2.len(); -// -// -// } -// -// -// fn mul(a: &SmallUint, b: &SmallUint) -> SmallUint { -// match (&a.0, &b.0) { -// (&SmallUintType::Inline(i), &SmallUintType::Inline(j)) => { -// match carrying_mul_u128(i, j, 0) { -// (t, 0) => SmallUint(SmallUintType::Inline(t)), -// (t, o) => { -// let mut res = Vec::with_capacity(8); -// -// let mut v = t; -// #[allow(clippy::needless_range_loop)] -// for _ in 0..4 { -// res.push(v as u32); -// -// v >>= 32; -// } -// -// let mut v = o; -// for _ in 4..8 { -// res.push(v as u32); -// -// v >>= 32; -// } -// -// while res[res.len() - 1] == 0 { -// res.pop(); -// } -// -// let size = res.len(); -// -// let mut slice = ManuallyDrop::new(res.into_boxed_slice()); -// -// SmallUint(SmallUintType::Heap((slice.as_mut_ptr(), size))) -// -// } -// } -// }, -// (&SmallUintType::Heap((r, s)), &SmallUintType::Inline(i)) | (&SmallUintType::Inline(i), &SmallUintType::Heap((r, s))) => { -// -// } -// } -// } -// -// +fn mul_two_slices(slice1: &[u32], slice2: &[u32]) -> Vec { + + // https://en.wikipedia.org/wiki/Karatsuba_algorithm + + let l1 = slice1.len(); + let l2 = slice2.len(); + + + if l1 == 0 || l2 == 0 { + return vec![]; + } else if l1 == 1 { + let mut overflow = 0; + let mut res: Vec = Vec::with_capacity(l2 + 1); + + #[allow(clippy::needless_range_loop)] + for i in 0..l2 { + let mut r = (slice2[i] as u64) * (slice1[0] as u64); + r += overflow as u64; + let m = r as u32; + overflow = (r >> 32) as u32; + res.push(m); + + } + + if overflow != 0 { + res.push(overflow); + } + + return res; + } else if l2 == 1 { + + let mut overflow = 0; + let mut res: Vec = Vec::with_capacity(l2 + 1); + + #[allow(clippy::needless_range_loop)] + for i in 0..l1 { + let mut r = (slice1[i] as u64) * (slice2[0] as u64); + r += overflow as u64; + let m = r as u32; + overflow = (r >> 32) as u32; + res.push(m); + + } + + if overflow != 0 { + res.push(overflow); + } + + return res; + } + + let m = std::cmp::min(l1, l2); + let m2 = (m as u32) / 2; + + let (low1, high1) = slice1.split_at(m2 as usize); + let (low2, high2) = slice2.split_at(m2 as usize); + + println!("z0: {:?}, {:?}", low1, low2); + println!("z1: {:?}, {:?}", &add_two_slices(low1, high1), &add_two_slices(low2, high2)); + println!("z2: {:?}, {:?}", high1, high2); + + let z0 = mul_two_slices(low1, low2); + let z1 = mul_two_slices(&add_two_slices(low1, high1), &add_two_slices(low2, high2)); + let z2 = mul_two_slices(high1, high2); + + let mut op0 = z0.clone(); + + for _ in 0..(m2 * 2) { + op0.insert(0, 0); + } + + let mut op1 = sub_two_slices(&sub_two_slices(&z1, &z2), &z0); + + for _ in 0..m2 { + op1.insert(0, 0); + } + + add_two_slices(&add_two_slices(&op0, &op1), &z0) + +} + + +fn mul(a: &SmallUint, b: &SmallUint) -> SmallUint { + match (&a.0, &b.0) { + (&SmallUintType::Inline(i), &SmallUintType::Inline(j)) => { + match carrying_mul_u128(i, j, 0) { + (t, 0) => SmallUint(SmallUintType::Inline(t)), + (t, o) => { + let mut res = Vec::with_capacity(8); + + let mut v = t; + #[allow(clippy::needless_range_loop)] + for _ in 0..4 { + res.push(v as u32); + + v >>= 32; + } + + let mut v = o; + for _ in 4..8 { + res.push(v as u32); + + v >>= 32; + } + + while res[res.len() - 1] == 0 { + res.pop(); + } + + let size = res.len(); + + let mut slice = ManuallyDrop::new(res.into_boxed_slice()); + + SmallUint(SmallUintType::Heap((slice.as_mut_ptr(), size))) + + } + } + }, + + (&SmallUintType::Heap((r, s)), &SmallUintType::Inline(i)) | (&SmallUintType::Inline(i), &SmallUintType::Heap((r, s))) => { + let slice1 = unsafe { core::slice::from_raw_parts(r, s) }; + + let mut res = [0, 0, 0, 0]; + + let mut v = i; + #[allow(clippy::needless_range_loop)] + for r in 0..4 { + res[r] = v as u32; + + v >>= 32; + } + + let result = mul_two_slices(slice1, &res[..]); + let size = result.len(); + + let mut slice = ManuallyDrop::new(result.into_boxed_slice()); + + SmallUint(SmallUintType::Heap((slice.as_mut_ptr(), size))) + + + }, + (&SmallUintType::Heap((r, s)), &SmallUintType::Heap((i, j))) => { + + let slice1 = unsafe { core::slice::from_raw_parts(r, s) }; + let slice2 = unsafe { core::slice::from_raw_parts(i, j) }; + + let res = mul_two_slices(slice1, slice2); + let size = res.len(); + + let mut slice = ManuallyDrop::new(res.into_boxed_slice()); + + SmallUint(SmallUintType::Heap((slice.as_mut_ptr(), size))) + + } + } +} + +basic_op!(Mul, SmallUint, mul); + diff --git a/src/tests.rs b/src/tests.rs index fc95bc9..7e4818b 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -45,12 +45,54 @@ fn test_op_add_u_u() { let q = i + k; assert_eq!(BigUint::from(&q), BigUint::new(vec![5, 4, 9, 3, 1, 81]) + u128::MAX); - let i = SmallUint::from(&BigUint::new(vec![3, 9, 8])); + let i = SmallUint::from(&BigUint::new(vec![3, 9, 8, 3, 1])); let k = SmallUint::from(&BigUint::new(vec![5, 4, 9, 3, 1, 81])); let q = i + k; - assert_eq!(BigUint::from(&q), BigUint::new(vec![3, 9, 8]) + BigUint::new(vec![5, 4, 9, 3, 1, 81])); + assert_eq!(BigUint::from(&q), BigUint::new(vec![3, 9, 8, 3, 1]) + BigUint::new(vec![5, 4, 9, 3, 1, 81])); } +#[test] +#[cfg(feature = "num-bigint")] +fn test_op_mul_u_u() { + let i = SmallUint::from(u128::MAX); + let k = SmallUint::from(u128::MAX); + let q = i * k; + assert_eq!(BigUint::from(&q), BigUint::from(u128::MAX) * u128::MAX); + + + let i = SmallUint::from(u32::MAX); + let k = SmallUint::from(&BigUint::new(vec![5, 4, 9, 3, 1, 81])); + let q = i * k; + assert_eq!(BigUint::from(&q), BigUint::new(vec![5, 4, 9, 3, 1, 81]) * u32::MAX); + + let i = SmallUint::from(&BigUint::new(vec![3, 9, 8, 3, 1])); + let k = SmallUint::from(&BigUint::new(vec![5, 4, 9, 3, 1, 81])); + let q = i * k; + assert_eq!(BigUint::from(&q), BigUint::new(vec![3, 9, 8, 3, 1]) * BigUint::new(vec![5, 4, 9, 3, 1, 81])); + +} + +#[test] +#[cfg(feature = "num-bigint")] +fn test_op_sub_u_u() { + let i = SmallUint::from(u128::MAX); + let k = SmallUint::from(u128::MAX); + let q = i - k; + assert_eq!(BigUint::from(&q), BigUint::from(u128::MAX) - u128::MAX); + + let k = SmallUint::from(&BigUint::new(vec![5, 4, 9, 3, 1, 81])); + let i = SmallUint::from(u128::MAX); + let q = k - i; + assert_eq!(BigUint::from(&q), BigUint::new(vec![5, 4, 9, 3, 1, 81]) - u128::MAX); + + let k = SmallUint::from(&BigUint::new(vec![5, 4, 9, 3, 1, 81])); + let i = SmallUint::from(&BigUint::new(vec![3, 9, 8, 3, 1])); + let q = k - i; + assert_eq!(BigUint::from(&q), BigUint::new(vec![5, 4, 9, 3, 1, 81]) - BigUint::new(vec![3, 9, 8, 3, 1])); + +} + + #[test] #[cfg(feature = "num-bigint")] fn test_bigint() {