Skip to content

Commit 3766b2d

Browse files
tgross35tgross35
committed
f128 division works!!
1 parent a17b3a4 commit 3766b2d

File tree

3 files changed

+102
-75
lines changed

3 files changed

+102
-75
lines changed

‎src/float/div.rs‎

Lines changed: 31 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ Relevant notation:
1111
fixed-point number with 1 integral bit, and 31 decimal bits. A `uqN` variable of type `uM`
1212
will have N bits of integer and M-N bits of fraction.
1313
- `hw`: half width, i.e. for `f64` this will be a `u32`.
14+
- `x` is the best estimate of `1/b`
1415
1516
# Method Overview
1617
@@ -110,6 +111,7 @@ const fn get_iterations<F: Float>() -> (usize, usize) {
110111

111112
// If widening multiply will be efficient (uses word-sized integers), there is no reason
112113
// to use half-sized iterations.
114+
// TODO: use half iterations.
113115
if 2 * size_of::<F>() <= size_of::<*const ()>() {
114116
(0, total_iterations)
115117
} else {
@@ -422,6 +424,7 @@ where
422424
// no overflow occurred earlier: ((rep_t)x_UQ0_hw * b_UQ1_hw >> HW) is
423425
// expected to be strictly positive because b_UQ1_hw has its highest bit set
424426
// and x_UQ0_hw should be rather large (it converges to 1/2 < 1/b_hw <= 1).
427+
425428
// let corr_uq1_hw: HalfRep<F> = zero_hw.wrapping_sub(x_uq0_hw.widen_mul(b_uq1_hw).hi());
426429

427430
// Now, we should multiply UQ0.HW and UQ1.(HW-1) numbers, naturally
@@ -436,6 +439,7 @@ where
436439
// The fact corr_UQ1_hw was virtually round up (due to result of
437440
// multiplication being **first** truncated, then negated - to improve
438441
// error estimations) can increase x_UQ0_hw by up to 2*Ulp of x_UQ0_hw.
442+
439443
// x_uq0_hw = (x_uq0_hw.widen_mul(corr_uq1_hw) >> (hw - 1)).cast();
440444

441445
// Now, either no overflow occurred or x_UQ0_hw is 0 or 1 in its half_rep_t
@@ -632,18 +636,35 @@ where
632636
F::from_repr(abs_result | quotient_sign)
633637
}
634638

635-
/// Perform one iteration at any width.
636-
///
637-
/// Given
638-
fn iter_once<I>(x_uq0: I, b_uq1: I) -> I
639-
where
640-
I: Int + HInt,
641-
<I as HInt>::D: ops::Shr<u32, Output = <I as HInt>::D>,
642-
{
643-
let corr_uq1: I = I::ZERO.wrapping_sub(x_uq0.widen_mul(b_uq1).hi());
644-
(x_uq0.widen_mul(corr_uq1) >> (I::BITS - 1)).lo()
639+
mod implementation {
640+
use crate::int::{DInt, HInt, Int};
641+
use core::ops;
642+
643+
/// Perform one iteration at any width to approach `1/b`, given previous guess `x`. It returns
644+
/// the next `x` as a UQ0 number.
645+
///
646+
/// This is the `x_{n+1} = 2*x_n - b*x_n^2` algorithm, implemented as `x_n * (2 - b*x_n)`.
647+
pub fn iter_once<I>(x_uq0: I, b_uq1: I) -> I
648+
where
649+
I: Int + HInt,
650+
<I as HInt>::D: ops::Shr<u32, Output = <I as HInt>::D>,
651+
{
652+
// `corr = 2 - b*x_n`
653+
//
654+
// This looks like `0 - b*x_n`. However, this works - in `UQ1`, `0.0 - x = 2.0 - x`.
655+
let corr_uq1: I = I::ZERO.wrapping_sub(x_uq0.widen_mul(b_uq1).hi());
656+
657+
// `x_n * corr = x_n * (2 - b*x_n)`
658+
(x_uq0.widen_mul(corr_uq1) >> (I::BITS - 1)).lo()
659+
}
645660
}
646661

662+
#[cfg(not(feature = "public-test-deps"))]
663+
use implementation::*;
664+
665+
#[cfg(feature = "public-test-deps")]
666+
pub use implementation::*;
667+
647668
intrinsics! {
648669
#[avr_skip]
649670
#[arm_aeabi_alias = __aeabi_fdiv]

‎testcrate/tests/div_rem.rs‎

Lines changed: 57 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -147,85 +147,77 @@ mod float_div {
147147
f32, __divsf3, Single, all();
148148
f64, __divdf3, Double, all();
149149
}
150-
}
151-
152-
#[cfg(target_arch = "arm")]
153-
mod float_div_arm {
154-
use super::*;
155150

151+
#[cfg(target_arch = "arm")]
156152
float! {
157153
f32, __divsf3vfp, Single, all();
158154
f64, __divdf3vfp, Double, all();
159155
}
160-
}
161-
162-
#[cfg(not(feature = "no-f16-f128"))]
163-
#[cfg(not(all(target_arch = "x86", not(target_feature = "sse"))))]
164-
mod float_div_f128 {
165-
use super::*;
166156

157+
#[cfg(not(feature = "no-f16-f128"))]
167158
#[cfg(not(any(target_arch = "powerpc", target_arch = "powerpc64")))]
168159
float! {
169160
f128, __divtf3, Quad, not(feature = "no-sys-f128");
170161
}
171162

163+
#[cfg(not(feature = "no-f16-f128"))]
172164
#[cfg(any(target_arch = "powerpc", target_arch = "powerpc64"))]
173165
float! {
174166
f128, __divkf3, Quad, not(feature = "no-sys-f128");
175167
}
176168
}
177169

178-
#[test]
179-
fn problem_f128() {
180-
use compiler_builtins::float::div::__divtf3;
181-
182-
let a = f128::from_bits(0x00000000000000000000000000000001);
183-
let b = f128::from_bits(0x0001FFFFFFFFFFFFFFFFFFFFFFFFFFFF);
184-
let res = __divtf3(a, b);
185-
println!(
186-
"{:#036x} / {:#036x} = {:#036x}",
187-
a.to_bits(),
188-
b.to_bits(),
189-
res.to_bits()
190-
);
191-
// got 0x3f8f0000000000000000000000000001
192-
// exp 0x3f8e0000000000000000000000000001
193-
assert_eq!(res.to_bits(), 0x3F8E0000000000000000000000000001);
194-
panic!();
195-
}
196-
197-
#[test]
198-
fn not_problem_f64() {
199-
use compiler_builtins::float::div::__divdf3;
200-
201-
let a = f64::from_bits(0x0000000000000001);
202-
let b = f64::from_bits(0x001FFFFFFFFFFFFF);
203-
let res = __divdf3(a, b);
204-
println!(
205-
"{:#018x} / {:#018x} = {:#018x}",
206-
a.to_bits(),
207-
b.to_bits(),
208-
res.to_bits()
209-
);
210-
// 0x3CA0000000000001
211-
assert_eq!(res.to_bits(), 0x3CA0000000000001);
212-
panic!();
213-
}
214-
215-
#[test]
216-
fn not_problem_f32() {
217-
use compiler_builtins::float::div::__divsf3;
218-
219-
let a = f32::from_bits(0x00000001);
220-
let b = f32::from_bits(0x00FFFFFF);
221-
let res = __divsf3(a, b);
222-
println!(
223-
"{:#010x} / {:#010x} = {:#010x}",
224-
a.to_bits(),
225-
b.to_bits(),
226-
res.to_bits()
227-
);
228-
// 0x33800001
229-
assert_eq!(res.to_bits(), 0x33800001);
230-
panic!();
231-
}
170+
// #[test]
171+
// fn problem_f128() {
172+
// use compiler_builtins::float::div::__divtf3;
173+
174+
// let a = f128::from_bits(0x00000000000000000000000000000001);
175+
// let b = f128::from_bits(0x0001FFFFFFFFFFFFFFFFFFFFFFFFFFFF);
176+
// let res = __divtf3(a, b);
177+
// println!(
178+
// "{:#036x} / {:#036x} = {:#036x}",
179+
// a.to_bits(),
180+
// b.to_bits(),
181+
// res.to_bits()
182+
// );
183+
// // got 0x3f8f0000000000000000000000000001
184+
// // exp 0x3f8e0000000000000000000000000001
185+
// assert_eq!(res.to_bits(), 0x3F8E0000000000000000000000000001);
186+
// panic!();
187+
// }
188+
189+
// #[test]
190+
// fn not_problem_f64() {
191+
// use compiler_builtins::float::div::__divdf3;
192+
193+
// let a = f64::from_bits(0x0000000000000001);
194+
// let b = f64::from_bits(0x001FFFFFFFFFFFFF);
195+
// let res = __divdf3(a, b);
196+
// println!(
197+
// "{:#018x} / {:#018x} = {:#018x}",
198+
// a.to_bits(),
199+
// b.to_bits(),
200+
// res.to_bits()
201+
// );
202+
// // 0x3CA0000000000001
203+
// assert_eq!(res.to_bits(), 0x3CA0000000000001);
204+
// panic!();
205+
// }
206+
207+
// #[test]
208+
// fn not_problem_f32() {
209+
// use compiler_builtins::float::div::__divsf3;
210+
211+
// let a = f32::from_bits(0x00000001);
212+
// let b = f32::from_bits(0x00FFFFFF);
213+
// let res = __divsf3(a, b);
214+
// println!(
215+
// "{:#010x} / {:#010x} = {:#010x}",
216+
// a.to_bits(),
217+
// b.to_bits(),
218+
// res.to_bits()
219+
// );
220+
// // 0x33800001
221+
// assert_eq!(res.to_bits(), 0x33800001);
222+
// panic!();
223+
// }

‎testcrate/tests/div_unit.rs‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
//! Unit tests for division since they can't be in-crate.
2+
3+
// use compiler_builtins::float::div::iter_once;
4+
5+
// #[test]
6+
// fn test_iter_once() {
7+
// dbg!(iter_once(0u32, 0u32));
8+
// dbg!(iter_once(1u32, 0u32));
9+
// dbg!(iter_once(0u32, 1u32));
10+
// dbg!(iter_once(0u32, 1u32));
11+
// dbg!(iter_once(u32::MAX, u16::MAX as u32));
12+
// dbg!(u32::MAX, u16::MAX);
13+
// panic!();
14+
// }

0 commit comments

Comments
 (0)