karatsuba multiplication: doesn't work yet

This commit is contained in:
EvilMuffinHa 2022-08-06 00:08:25 -04:00
parent 3510274fce
commit 4d6017b974
2 changed files with 303 additions and 56 deletions

View File

@ -1,6 +1,6 @@
use crate::SmallUint; use crate::SmallUint;
use crate::smallint::SmallUintType; use crate::smallint::SmallUintType;
use core::ops::Add; use core::ops::{Add, Sub, Mul};
use core::mem::ManuallyDrop; use core::mem::ManuallyDrop;
macro_rules! basic_op { macro_rules! basic_op {
@ -83,12 +83,13 @@ fn add_two_slices(slice1: &[u32], slice2: &[u32]) -> Vec<u32> {
if carry { if carry {
res.push(1); res.push(1);
} }
while res[res.len() - 1] == 0 { while res.len() != 1 && res[res.len() - 1] == 0 {
res.pop(); res.pop();
} }
res res
} }
fn add(a: &SmallUint, b: &SmallUint) -> SmallUint { fn add(a: &SmallUint, b: &SmallUint) -> SmallUint {
match (&a.0, &b.0) { match (&a.0, &b.0) {
@ -153,6 +154,103 @@ fn add(a: &SmallUint, b: &SmallUint) -> SmallUint {
basic_op!(Add, SmallUint, add); basic_op!(Add, SmallUint, add);
fn sub_two_slices(slice1: &[u32], slice2: &[u32]) -> Vec<u32> {
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 // Taken from https://github.com/rust-lang/rust/issues/85532#issuecomment-916309635. Credit to
// AaronKutch. // AaronKutch.
const fn carrying_mul_u128(lhs: u128, rhs: u128, carry: u128) -> (u128, u128) { 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) (sum0, sum1)
} }
// fn mul_two_slices(slice1: &[u32], slice2: &[u32]) -> Vec<u32> { fn mul_two_slices(slice1: &[u32], slice2: &[u32]) -> Vec<u32> {
//
// let l1 = slice1.len(); // https://en.wikipedia.org/wiki/Karatsuba_algorithm
// let l2 = slice2.len();
// let l1 = slice1.len();
// let l2 = slice2.len();
// }
//
// if l1 == 0 || l2 == 0 {
// fn mul(a: &SmallUint, b: &SmallUint) -> SmallUint { return vec![];
// match (&a.0, &b.0) { } else if l1 == 1 {
// (&SmallUintType::Inline(i), &SmallUintType::Inline(j)) => { let mut overflow = 0;
// match carrying_mul_u128(i, j, 0) { let mut res: Vec<u32> = Vec::with_capacity(l2 + 1);
// (t, 0) => SmallUint(SmallUintType::Inline(t)),
// (t, o) => { #[allow(clippy::needless_range_loop)]
// let mut res = Vec::with_capacity(8); for i in 0..l2 {
// let mut r = (slice2[i] as u64) * (slice1[0] as u64);
// let mut v = t; r += overflow as u64;
// #[allow(clippy::needless_range_loop)] let m = r as u32;
// for _ in 0..4 { overflow = (r >> 32) as u32;
// res.push(v as u32); res.push(m);
//
// v >>= 32; }
// }
// if overflow != 0 {
// let mut v = o; res.push(overflow);
// for _ in 4..8 { }
// res.push(v as u32);
// return res;
// v >>= 32; } else if l2 == 1 {
// }
// let mut overflow = 0;
// while res[res.len() - 1] == 0 { let mut res: Vec<u32> = Vec::with_capacity(l2 + 1);
// res.pop();
// } #[allow(clippy::needless_range_loop)]
// for i in 0..l1 {
// let size = res.len(); let mut r = (slice1[i] as u64) * (slice2[0] as u64);
// r += overflow as u64;
// let mut slice = ManuallyDrop::new(res.into_boxed_slice()); let m = r as u32;
// overflow = (r >> 32) as u32;
// SmallUint(SmallUintType::Heap((slice.as_mut_ptr(), size))) res.push(m);
//
// } }
// }
// }, if overflow != 0 {
// (&SmallUintType::Heap((r, s)), &SmallUintType::Inline(i)) | (&SmallUintType::Inline(i), &SmallUintType::Heap((r, s))) => { 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);

View File

@ -45,12 +45,54 @@ fn test_op_add_u_u() {
let q = i + k; let q = i + k;
assert_eq!(BigUint::from(&q), BigUint::new(vec![5, 4, 9, 3, 1, 81]) + u128::MAX); 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 k = SmallUint::from(&BigUint::new(vec![5, 4, 9, 3, 1, 81]));
let q = i + k; 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] #[test]
#[cfg(feature = "num-bigint")] #[cfg(feature = "num-bigint")]
fn test_bigint() { fn test_bigint() {