Commit 31efbc49 authored by Andrey Zgarbul's avatar Andrey Zgarbul

softposit-math

parent 18c82dba
......@@ -1076,6 +1076,15 @@ macro_rules! quire_add_sub {
}
}
impl ops::AddAssign<($posit, ($posit, $posit, $posit))> for $quire {
#[inline]
fn add_assign(&mut self, rhs: ($posit, ($posit, $posit, $posit))) {
*self += (rhs.0, (rhs.1).0);
*self += (rhs.0, (rhs.1).1);
*self += (rhs.0, (rhs.1).2);
}
}
impl ops::AddAssign<$posit> for $quire {
#[inline]
fn add_assign(&mut self, rhs: $posit) {
......
......@@ -307,6 +307,22 @@ impl Q16E1 {
pub fn neg(&mut self) {
self.0 = self.0.wrapping_neg();
}
#[inline]
fn into_two_posits(mut self) -> (P16E1, P16E1) {
let p1 = self.to_posit();
self -= p1;
(p1, self.to_posit())
}
#[inline]
fn into_three_posits(mut self) -> (P16E1, P16E1, P16E1) {
let p1 = self.to_posit();
self -= p1;
let p2 = self.to_posit();
self -= p2;
(p1, p2, self.to_posit())
}
}
impl crate::Quire<P16E1> for Q16E1 {
......
use super::P16E1;
use super::{P16E1, Q16E1};
use crate::{MulAddType, WithSign};
const HALF: P16E1 = P16E1::new(0x_3000);
......@@ -35,11 +35,11 @@ impl P16E1 {
}
#[inline]
pub fn floor(self) -> Self {
(self - HALF).round()
floor(self)
}
#[inline]
pub fn ceil(self) -> Self {
(self + HALF).round()
ceil(self)
}
#[inline]
pub fn round(self) -> Self {
......@@ -93,7 +93,7 @@ impl P16E1 {
}
#[inline]
pub fn exp(self) -> Self {
unimplemented!()
exp(self)
}
#[inline]
pub fn exp2(self) -> Self {
......@@ -392,26 +392,25 @@ fn round(p_a: P16E1) -> P16E1 {
let mut mask = 0x2000_u16;
let mut scale = 0_u16;
let mut u_a = p_a.to_bits();
let mut ui_a = u_a; // Copy of the input.
let mut ui_a = p_a.to_bits();
let sign = ui_a > 0x8000;
// sign is True if p_a > NaR.
if sign {
ui_a = ui_a.wrapping_neg() // A is now |A|.
};
if ui_a <= 0x3000 {
let u_a = if ui_a <= 0x3000 {
// 0 <= |p_a| <= 1/2 rounds to zero.
return P16E1::ZERO;
} else if ui_a < 0x4800 {
// 1/2 < x < 3/2 rounds to 1.
u_a = 0x4000;
0x4000
} else if ui_a <= 0x5400 {
// 3/2 <= x <= 5/2 rounds to 2.
u_a = 0x5000;
0x5000
} else if ui_a >= 0x7C00 {
// If |A| is 256 or greater, leave it unchanged.
return P16E1::from_bits(u_a); // This also takes care of the NaR case, 0x8000.
return p_a; // This also takes care of the NaR case, 0x8000.
} else {
// 34% of the cases, we have to decode the posit.
while (mask & ui_a) != 0 {
......@@ -439,8 +438,8 @@ fn round(p_a: P16E1) -> P16E1 {
ui_a += mask << 1;
}
}
u_a = ui_a;
}
ui_a
};
P16E1::from_bits(u_a.with_sign(sign))
}
......@@ -532,6 +531,228 @@ fn sqrt(p_a: P16E1) -> P16E1 {
P16E1::from_bits(ui_z | ((frac_z >> 4) as u16))
}
fn ceil(p_a: P16E1) -> P16E1 {
let mut mask = 0x2000_u16;
let mut scale = 0_u16;
let mut ui_a = p_a.to_bits();
let sign = ui_a > 0x8000;
// sign is True if p_a > NaR.
if sign {
ui_a = ui_a.wrapping_neg() // A is now |A|.
};
let u_a = if ui_a == 0 {
return p_a;
} else if ui_a <= 0x4000 {
// 0 <= |pA| < 1 ceiling to zero.(if not negative and whole number)
if sign && (ui_a != 0x4000) {
0x0
} else {
0x4000
}
} else if ui_a <= 0x5000 {
// 1 <= x < 2 ceiling to 1 (if not negative and whole number)
if sign && (ui_a != 0x5000) {
0x4000
} else {
0x5000
}
} else if ui_a <= 0x5800 {
// 2 <= x < 3 ceiling to 2 (if not negative and whole number)
if sign & (ui_a != 0x5800) {
0x5000
} else {
0x5800
}
} else if ui_a >= 0x7C00 {
// If |A| is 256 or greater, leave it unchanged.
return p_a; // This also takes care of the NaR case, 0x8000.
} else {
// 34% of the cases, we have to decode the posit.
while (mask & ui_a) != 0 {
// Increment scale by 2 for each regime sign bit.
scale += 2; // Regime sign bit is always 1 in this range.
mask >>= 1; // Move the mask right, to the next bit.
}
mask >>= 1; // Skip over termination bit.
if (mask & ui_a) != 0 {
scale += 1; // If exponent is 1, increment the scale.
}
mask >>= scale; // Point to the last bit of the integer part.
mask >>= 1;
let mut tmp = ui_a & mask;
let bit_n_plus_one = tmp; // "True" if nonzero.
ui_a ^= tmp; // Erase the bit, if it was set.
tmp = ui_a & (mask - 1); // tmp has any remaining bits = bitsMore
ui_a ^= tmp; // Erase those bits, if any were set.
if !sign && (bit_n_plus_one | tmp) != 0 {
ui_a += mask << 1;
}
ui_a
};
P16E1::from_bits(u_a.with_sign(sign))
}
fn floor(p_a: P16E1) -> P16E1 {
let mut mask = 0x2000_u16;
let mut scale = 0_u16;
let mut ui_a = p_a.to_bits();
let sign = ui_a > 0x8000;
// sign is True if p_a > NaR.
if sign {
ui_a = ui_a.wrapping_neg() // A is now |A|.
};
let u_a = if ui_a < 0x4000 {
// 0 <= |pA| < 1 floor to zero.(if not negative and whole number)
if sign && (ui_a != 0x0) {
0x4000
} else {
0x0
}
} else if ui_a < 0x5000 {
// 1 <= x < 2 floor to 1 (if not negative and whole number)
if sign && (ui_a != 0x4000) {
0x5000
} else {
0x4000
}
} else if ui_a < 0x5800 {
// 2 <= x < 3 floor to 2 (if not negative and whole number)
if sign & (ui_a != 0x5000) {
0x5800
} else {
0x5000
}
} else if ui_a >= 0x7C00 {
// If |A| is 256 or greater, leave it unchanged.
return p_a; // This also takes care of the NaR case, 0x8000.
} else {
// 34% of the cases, we have to decode the posit.
while (mask & ui_a) != 0 {
// Increment scale by 2 for each regime sign bit.
scale += 2; // Regime sign bit is always 1 in this range.
mask >>= 1; // Move the mask right, to the next bit.
}
mask >>= 1; // Skip over termination bit.
if (mask & ui_a) != 0 {
scale += 1; // If exponent is 1, increment the scale.
}
mask >>= scale; // Point to the last bit of the integer part.
mask >>= 1;
let mut tmp = ui_a & mask;
let bit_n_plus_one = tmp; // "True" if nonzero.
ui_a ^= tmp; // Erase the bit, if it was set.
tmp = ui_a & (mask - 1); // tmp has any remaining bits = bitsMore
ui_a ^= tmp; // Erase those bits, if any were set.
if sign && ((bit_n_plus_one | tmp) != 0) {
ui_a += mask << 1;
}
ui_a
};
P16E1::from_bits(u_a.with_sign(sign))
}
fn exp(p_a: P16E1) -> P16E1 {
// Use names for commonly-used posits.
const LARGE: P16E1 = P16E1::new(0x7FFE);
const SMALL: P16E1 = P16E1::new(0x0002);
// If input is NaR, return NaR.
if p_a.is_nar() {
P16E1::NAR
}
// If input is large enough that result is maxpos, return maxpos.
else if P16E1::new(0x70AE) <= p_a {
P16E1::MAX
}
// If input is negative enough that result is minpos, return minpos.
else if p_a <= P16E1::new(-0x_70ae) {
P16E1::MIN_POSITIVE
}
// This range rounds to maxpos/4, the next-to-largest posit16.
else if P16E1::new(0x706F) < p_a {
LARGE
}
// This range rounds to minpos*4, the next-to-smallest posit16.
else if p_a < P16E1::new(-0x_7067) {
SMALL
} else {
// Scale input by 1/log(2) to trio precision, using the quire.
let mut q = Q16E1::init();
q += (
p_a,
(P16E1::new(0x4715), P16E1::new(0x0087), P16E1::new(0x000A)),
);
let mut p_n = q.to_posit().round();
q -= p_n; // Reduce the argument range.
let (t1, t2, t3) = q.into_three_posits();
// Evaluate 6th-degree polynomial in t using Horner's rule.
let mut q = Q16E1::init();
q += (t1, P16E1::new(0x00D1));
q += (P16E1::new(0x0A20), P16E1::new(0x0F28));
let p1 = q.to_posit(); // c6 * t + c5, solo precision
q.clear();
q += (p1, t1);
q += (P16E1::new(0x0A60), P16E1::new(0x28B8));
let (p1, p2) = q.into_two_posits(); // (...) * t + c4, duo precision
let mut q = Q16E1::init();
q += (t1, (p1, p2));
q += (p1, t2);
q += (P16E1::new(0x0B80), P16E1::new(0x4E50));
let (p1, p2) = q.into_two_posits(); // (...) * t + c3, duo precision
let mut q = Q16E1::init();
q += (t1, (p1, p2));
q += (p1, t2);
q += (P16E1::new(0x11F0), P16E1::new(0x58C1));
let (p1, p2, p3) = q.into_three_posits(); // (...) * t + c2, trio precision
let mut q = Q16E1::init();
q += (t1, (p1, p2, p3));
q += (p1, t2);
q += (P16E1::new(0x2D78), P16E1::new(0x4816));
q += (P16E1::new(0x0180), P16E1::new(0x0180));
let (p1, p2, p3) = q.into_three_posits(); // c1 term, trio precision
let mut q = Q16E1::init();
q += (t1, (p1, p2, p3));
q += (p1, t2);
q += (p1, t3);
q += P16E1::ONE; // (...) * t + c0, (where c0 = 1).
let (p1, p2, p3) = q.into_three_posits(); // polynomial for exp(x), trio precision
// Convert `p_n` to an integer and use to compute `2` to the power `n`.
let mut n = i32::from(p_n);
n = if n < 0 {
((n & 0x1) | 0x0002) << (12 - ((-n) >> 1))
} else {
0x7FFF & (((n & 0x1) | 0x7FFC) << (12 - (n >> 1)))
};
p_n = P16E1::new(n as i16);
// Scale the result by that power.
let mut q = Q16E1::init();
q += (p_n, (p1, p2, p3)); // trio precision; round it and we're done.
q.to_posit()
}
}
#[test]
fn test_mul_add() {
use rand::Rng;
......
......@@ -220,11 +220,11 @@ impl P32E2 {
}
#[inline]
pub fn floor(self) -> Self {
(self - HALF).round()
floor(self)
}
#[inline]
pub fn ceil(self) -> Self {
(self + HALF).round()
ceil(self)
}
#[inline]
pub fn round(self) -> Self {
......@@ -386,7 +386,7 @@ impl P32E2 {
}
#[allow(clippy::cognitive_complexity)]
pub(super) fn mul_add(mut ui_a: u32, mut ui_b: u32, mut ui_c: u32, op: MulAddType) -> P32E2 {
fn mul_add(mut ui_a: u32, mut ui_b: u32, mut ui_c: u32, op: MulAddType) -> P32E2 {
let mut bits_more = false;
//NaR
if (ui_a == 0x8000_0000) || (ui_b == 0x8000_0000) || (ui_c == 0x8000_0000) {
......@@ -573,7 +573,7 @@ pub(super) fn mul_add(mut ui_a: u32, mut ui_b: u32, mut ui_c: u32, op: MulAddTyp
P32E2::from_bits(u_z.with_sign(sign_z))
}
pub(super) fn round(p_a: P32E2) -> P32E2 {
pub fn round(p_a: P32E2) -> P32E2 {
let mut mask = 0x2000_0000_u32;
let mut scale = 0_u32;
......@@ -636,7 +636,7 @@ pub(super) fn round(p_a: P32E2) -> P32E2 {
}
#[inline]
pub(super) fn sqrt(p_a: P32E2) -> P32E2 {
pub fn sqrt(p_a: P32E2) -> P32E2 {
let mut ui_a = p_a.to_bits();
// If NaR or a negative number, return NaR.
......@@ -726,6 +726,148 @@ pub(super) fn sqrt(p_a: P32E2) -> P32E2 {
P32E2::from_bits(ui_z | (exp_z << (27 - shift)) | (frac_z >> (5 + shift)) as u32)
}
pub fn ceil(p_a: P32E2) -> P32E2 {
let mut mask = 0x2000_0000_u32;
let mut scale = 0_u32;
let mut ui_a = p_a.to_bits();
let sign = (ui_a & 0x8000_0000) != 0;
// sign is True if pA > NaR.
if sign {
ui_a = ui_a.wrapping_neg();
} // A is now |A|.
let u_a = if ui_a <= 0x_4000_0000 {
// 0 <= |pA| < 1 floor to zero.(if not negative and whole number)
if sign && (ui_a != 0x0) {
0x0
} else {
0x_4000_0000
}
} else if ui_a <= 0x_4800_0000 {
// 0 <= |pA| < 1 floor to 1.(if not negative and whole number)
if sign && (ui_a != 0x_4800_0000) {
0x_4000_0000
} else {
0x_4800_0000
}
} else if ui_a <= 0x_4C00_0000 {
// 0 <= |pA| < 2 floor to zero.(if not negative and whole number)
if sign && (ui_a != 0x_4C00_0000) {
0x_4800_0000
} else {
0x_4C00_0000
}
} else if ui_a >= 0x7E80_0000 {
// If |A| is 0x7E80_0000 (posit is pure integer value), leave it unchanged.
return p_a; // This also takes care of the NaR case, 0x8000_0000.
} else {
// 34% of the cases, we have to decode the posit.
while (mask & ui_a) != 0 {
scale += 4;
mask >>= 1;
}
mask >>= 1;
//Exponential (2 bits)
if (mask & ui_a) != 0 {
scale += 2;
}
mask >>= 1;
if (mask & ui_a) != 0 {
scale += 1;
}
mask >>= scale;
//the rest of the bits
mask >>= 1;
let mut tmp = ui_a & mask;
let bit_n_plus_one = tmp;
ui_a ^= tmp; // Erase the bit, if it was set.
tmp = ui_a & (mask - 1); // this is actually bits_more
ui_a ^= tmp;
if !sign && (bit_n_plus_one | tmp) != 0 {
ui_a += mask << 1;
}
ui_a
};
P32E2::from_bits(u_a.with_sign(sign))
}
pub fn floor(p_a: P32E2) -> P32E2 {
let mut mask = 0x2000_0000_u32;
let mut scale = 0_u32;
let mut ui_a = p_a.to_bits();
let sign = (ui_a & 0x8000_0000) != 0;
// sign is True if pA > NaR.
if sign {
ui_a = ui_a.wrapping_neg();
} // A is now |A|.
let u_a = if ui_a < 0x_4000_0000 {
// 0 <= |pA| < 1 floor to zero.(if not negative and whole number)
if sign && (ui_a != 0x0) {
0x0
} else {
0x_4000_0000
}
} else if ui_a < 0x_4800_0000 {
// 0 <= |pA| < 1 floor to 1.(if not negative and whole number)
if sign && (ui_a != 0x_4000_0000) {
0x_4800_0000
} else {
0x_4000_0000
}
} else if ui_a <= 0x_4C00_0000 {
// 0 <= |pA| < 2 floor to zero.(if not negative and whole number)
if sign && (ui_a != 0x_4800_0000) {
0x_4C00_0000
} else {
0x_4800_0000
}
} else if ui_a >= 0x7E80_0000 {
// If |A| is 0x7E80_0000 (posit is pure integer value), leave it unchanged.
return p_a; // This also takes care of the NaR case, 0x8000_0000.
} else {
// 34% of the cases, we have to decode the posit.
while (mask & ui_a) != 0 {
scale += 4;
mask >>= 1;
}
mask >>= 1;
//Exponential (2 bits)
if (mask & ui_a) != 0 {
scale += 2;
}
mask >>= 1;
if (mask & ui_a) != 0 {
scale += 1;
}
mask >>= scale;
//the rest of the bits
mask >>= 1;
let mut tmp = ui_a & mask;
let bit_n_plus_one = tmp;
ui_a ^= tmp; // Erase the bit, if it was set.
tmp = ui_a & (mask - 1); // this is actually bits_more
ui_a ^= tmp;
if sign && (bit_n_plus_one | tmp) != 0 {
ui_a += mask << 1;
}
ui_a
};
P32E2::from_bits(u_a.with_sign(sign))
}
#[test]
fn test_mul_add() {
use rand::Rng;
......
......@@ -35,11 +35,11 @@ impl P8E0 {
}
#[inline]
pub fn floor(self) -> Self {
(self - HALF).round()
floor(self)
}
#[inline]
pub fn ceil(self) -> Self {
(self + HALF).round()
ceil(self)
}
#[inline]
pub fn round(self) -> Self {
......@@ -93,7 +93,7 @@ impl P8E0 {
}
#[inline]
pub fn exp(self) -> Self {
unimplemented!()
exp(self)
}
#[inline]
pub fn exp2(self) -> Self {
......@@ -101,7 +101,7 @@ impl P8E0 {
}
#[inline]
pub fn ln(self) -> Self {
unimplemented!()
log(self)
}
#[inline]
pub fn log(self, _base: Self) -> Self {
......@@ -200,26 +200,25 @@ fn round(p_a: P8E0) -> P8E0 {
let mut mask = 0x20_u8;
let mut scale = 0_u8;
let mut u_a = p_a.to_bits();
let mut ui_a = u_a;
let mut ui_a = p_a.to_bits();
let sign = ui_a > 0x80;
// sign is True if p_a > NaR.
if sign {
ui_a = ui_a.wrapping_neg();
}
if ui_a <= 0x20 {
let u_a = if ui_a <= 0x20 {
// 0 <= |p_a| <= 1/2 rounds to zero.
return P8E0::ZERO;
} else if ui_a < 0x50 {
// 1/2 < x < 3/2 rounds to 1.
u_a = 0x40;
0x40
} else if ui_a <= 0x64 {
// 3/2 <= x <= 5/2 rounds to 2.
u_a = 0x60;
0x60
} else if ui_a >= 0x78 {
// If |A| is 8 or greater, leave it unchanged.
return P8E0::from_bits(u_a); // This also takes care of the NaR case, 0x80.
return p_a; // This also takes care of the NaR case, 0x80.
} else {
while (mask & ui_a) != 0 {
scale += 1;
......@@ -239,8 +238,8 @@ fn round(p_a: P8E0) -> P8E0 {
if bit_n_plus_one && (((bit_last as u8) | tmp) != 0) {
ui_a += mask << 1;
}
u_a = ui_a;
}
ui_a
};
P8E0::from_bits(u_a.with_sign(sign))
}
......@@ -414,6 +413,163 @@ fn mul_add(mut ui_a: u8, mut ui_b: u8, mut ui_c: u8, op: MulAddType) -> P8E0 {
P8E0::from_bits(u_z.with_sign(sign_z))
}
fn ceil(p_a: P8E0) -> P8E0 {
let mut mask = 0x20_u8;
let mut scale = 0_u8;
let mut ui_a = p_a.to_bits();
let sign = ui_a > 0x80;
// sign is True if p_a > NaR.
if sign {
ui_a = ui_a.wrapping_neg();
}
let u_a = if ui_a == 0 {
return p_a;
} else if ui_a <= 0x40 {
// 0 <= |pA| < 1 floor to zero.(if not negative and whole number)
if sign && (ui_a != 0x40) {
0x0
} else {
0x40
}
} else if ui_a <= 0x60 {
// 1 <= x < 2 floor to 1 (if not negative and whole number)
if sign && (ui_a != 0x60) {
0x40
} else {
0x60
}
} else if ui_a <= 0x68 {
// 2 <= x < 3 floor to 2 (if not negative and whole number)
if sign && (ui_a != 0x68) {
0x60
} else {