Skip to content

Commit bfdb026

Browse files
folkertdevfolkertdev
committed
ppc: update frexp 
1 parent eea6a55 commit bfdb026

2 files changed

Lines changed: 252 additions & 5 deletions

File tree

‎src/ppc.rs‎

Lines changed: 246 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -165,6 +165,30 @@ fn is_power_of_two<F: Float>(x: F) -> bool {
165165
x.is_finite_non_zero() && x.abs().scalbn(-x.ilogb()) == F::from_u128(1).value
166166
}
167167

168+
/// Compute the ULP of the input using a definition from:
169+
/// Jean-Michel Muller. On the definition of ulp(x). [Research Report] RR-5504,
170+
/// LIP RR-2005-09, INRIA, LIP. 2005, pp.16. inria-00070503
171+
fn harrison_ulp<F: Float>(x: F) -> F {
172+
match x.category() {
173+
Category::NaN => return F::qnan(None),
174+
Category::Infinity => return F::INFINITY,
175+
Category::Zero => return F::SMALLEST,
176+
Category::Normal => { /* fall through */ }
177+
}
178+
179+
if x.is_denormal() || x.is_smallest_normalized() {
180+
return F::SMALLEST;
181+
}
182+
183+
// Match LLVM in not considering negative powers of two.
184+
let mut exp = x.ilogb();
185+
if !x.is_negative() && is_power_of_two(x) {
186+
exp -= 1;
187+
}
188+
189+
F::from_u128(1).value.scalbn(exp - (F::PRECISION as i32 - 1))
190+
}
191+
168192
impl<F: FloatConvert<Fallback<F>>> Float for DoubleFloat<F>
169193
where
170194
Self: From<Fallback<F>>,
@@ -565,12 +589,111 @@ where
565589
}
566590

567591
fn frexp_r(self, exp: &mut ExpInt, round: Round) -> Self {
568-
let a = self.0.frexp_r(exp, round);
569-
let mut b = self.1;
570-
if self.category() == Category::Normal {
571-
b = b.scalbn_r(-*exp, round);
592+
// Get the unbiased exponent e of the number, where |self| = m * 2^e for m in [1.0, 2.0).
593+
*exp = self.ilogb();
594+
595+
// For NaNs, quiet any signaling NaN and return the result, as per standard practice.
596+
if *exp == crate::IEK_NAN {
597+
let mut quiet = self;
598+
quiet.0 = quiet.0.add_r(F::ZERO, Round::NearestTiesToEven).value;
599+
return quiet;
572600
}
573-
DoubleFloat(a, b)
601+
602+
// For infinity, return it unchanged. The exponent remains IEK_Inf.
603+
if *exp == crate::IEK_INF {
604+
return self;
605+
}
606+
607+
// For zero, the fraction is zero and the standard requires the exponent be 0.
608+
if *exp == crate::IEK_ZERO {
609+
*exp = 0;
610+
return self;
611+
}
612+
613+
let DoubleFloat(hi, lo) = self;
614+
615+
// frexp requires the fraction's absolute value to be in [0.5, 1.0).
616+
// ilogb provides an exponent for an absolute value in [1.0, 2.0).
617+
// Increment the exponent to ensure the fraction is in the correct range.
618+
*exp += 1;
619+
620+
let signs_disagree = hi.is_negative() != lo.is_negative();
621+
let mut second = lo;
622+
if self.category() == Category::Normal && lo.is_finite_non_zero() {
623+
// The interpretation of Round::TowardZero depends on the sign of the combined
624+
// self rather than the sign of the component.
625+
let lo_rounding_mode = if round == Round::TowardZero {
626+
if self.is_negative() {
627+
Round::TowardPositive
628+
} else {
629+
Round::TowardNegative
630+
}
631+
} else if round == Round::NearestTiesToAway && signs_disagree && *exp > 0 {
632+
// For Round::NearestTiesToAway, we face a similar problem. If signs disagree,
633+
// Lo is a correction *toward* zero relative to Hi. Rounding Lo
634+
// "away from zero" based on its own sign would move the value in the
635+
// wrong direction. As a safe proxy, we use Round::NearestTiesToEven, which is
636+
// direction-agnostic. We only need to bother with this if Lo is scaled
637+
// down.
638+
Round::NearestTiesToEven
639+
} else {
640+
round
641+
};
642+
643+
second = lo.scalbn_r(-*exp, lo_rounding_mode);
644+
645+
// The Round::NearestTiesToEven proxy is correct most of the time, but it
646+
// differs from Round::NearestTiesToAway when the scaled value of Lo is an
647+
// exact midpoint.
648+
// NOTE: This is morally equivalent to roundTiesTowardZero.
649+
if round == Round::NearestTiesToAway && lo_rounding_mode == Round::NearestTiesToEven {
650+
// Re-scale the result back to check if rounding occurred.
651+
let recomposed_lo = second.scalbn_r(*exp, Round::NearestTiesToEven);
652+
if recomposed_lo != lo {
653+
// RoundingError tells us which direction we rounded:
654+
// - RoundingError > 0: we rounded up.
655+
// - RoundingError < 0: we down up.
656+
let rounding_error = (recomposed_lo - lo).value;
657+
// Determine if scalbn(Lo, -Exp) landed exactly on a midpoint.
658+
// We do this by checking if the absolute rounding error is exactly
659+
// half a ULP of the result.
660+
let ulp_of_second = harrison_ulp(second);
661+
let scaled_ulp_of_second = ulp_of_second.scalbn_r(*exp - 1, Round::NearestTiesToEven);
662+
let is_midpoint = rounding_error.abs() == scaled_ulp_of_second;
663+
let rounded_lo_away = second.is_negative() == rounding_error.is_negative();
664+
// The sign of Hi and Lo disagree and we rounded Lo away: we must
665+
// decrease the magnitude of Second to increase the magnitude
666+
// First+Second.
667+
if is_midpoint && rounded_lo_away {
668+
second = if second.is_negative() {
669+
second.next_up().value
670+
} else {
671+
second.next_down().value
672+
};
673+
}
674+
}
675+
}
676+
677+
// Handle a tricky edge case where self is slightly less than a power of two
678+
// (e.g., self = 2^k - epsilon). In this situation:
679+
// 1. Hi is 2^k, and Lo is a small negative value -epsilon.
680+
// 2. ilogb(self) correctly returns k-1.
681+
// 3. Our initial Exp becomes (k-1) + 1 = k.
682+
// 4. Scaling Hi (2^k) by 2^-k would yield a magnitude of 1.0 and
683+
// scaling Lo by 2^-k would yield zero. This would make the result 1.0
684+
// which is an invalid fraction, as the required interval is [0.5, 1.0).
685+
// We detect this specific case by checking if Hi is a power of two and if
686+
// the scaled Lo underflowed to zero. The fix: Increment Exp to k+1. This
687+
// adjusts the scale factor, causing Hi to be scaled to 0.5, which is a
688+
// valid fraction.
689+
if second.is_zero() && signs_disagree && is_power_of_two(hi) {
690+
*exp += 1;
691+
}
692+
}
693+
694+
let first = hi.scalbn_r(-*exp, round);
695+
696+
DoubleFloat(first, second)
574697
}
575698
}
576699

@@ -788,4 +911,122 @@ mod tests {
788911
assert!(quiet_nan.round_to_integral(Round::NearestTiesToAway).value.bitwise_eq(quiet_nan));
789912
assert!(quiet_nan.round_to_integral(Round::NearestTiesToEven).value.bitwise_eq(quiet_nan));
790913
}
914+
915+
#[test]
916+
fn ppc_double_double_frexp() {
917+
let double_from_f64 = |f: f64| ieee::Double::from_bits(f.to_bits().into());
918+
let dd = |hi, lo| DoubleFloat(double_from_f64(hi), double_from_f64(lo));
919+
920+
let frexp_test_cases = [
921+
// Input: +infinity
922+
TestCase::new(dd(f64::INFINITY, 0.0)).with_all((dd(f64::INFINITY, 0.0), ExpInt::MAX)),
923+
// Input: -infinity
924+
TestCase::new(dd(f64::NEG_INFINITY, 0.0)).with_all((dd(f64::NEG_INFINITY, 0.0), ExpInt::MAX)),
925+
// Input: 2^-1074
926+
TestCase::new(dd(f64::from_bits(1), 0.0)).with_all((dd(0.5, 0.0), -1073)),
927+
// Input: (2^1, -2^-1073 + -2^-1074)
928+
TestCase::new(dd(2.0, -3.0 * f64::from_bits(1)))
929+
.with_all((dd(1.0, -2.0 * f64::from_bits(1)), 1))
930+
.with(Round::NearestTiesToAway, (dd(1.0, -f64::from_bits(1)), 1))
931+
.with(Round::TowardPositive, (dd(1.0, -f64::from_bits(1)), 1)),
932+
// Input: (2^1, -2^-1073)
933+
TestCase::new(dd(2.0, -2.0 * f64::from_bits(1))).with_all((dd(1.0, -f64::from_bits(1)), 1)),
934+
// Input: (2^1, -2^-1074)
935+
TestCase::new(dd(2.0, -f64::from_bits(1)))
936+
.with_all((dd(0.5, -0.0), 2))
937+
.with(Round::TowardNegative, (dd(1.0, -f64::from_bits(1)), 1))
938+
.with(Round::TowardZero, (dd(1.0, -f64::from_bits(1)), 1)),
939+
// Input: (2^2, -2^-1072 + -2^-1073 + -2^-1074)
940+
TestCase::new(dd(4.0, -7.0 * f64::from_bits(1)))
941+
.with_all((dd(1.0, -2.0 * f64::from_bits(1)), 2))
942+
.with(Round::TowardPositive, (dd(1.0, -f64::from_bits(1)), 2)),
943+
// Input: (2^2, -2^-1072 + -2^-1073)
944+
TestCase::new(dd(4.0, -6.0 * f64::from_bits(1)))
945+
.with_all((dd(1.0, -2.0 * f64::from_bits(1)), 2))
946+
.with(Round::NearestTiesToAway, (dd(1.0, -f64::from_bits(1)), 2))
947+
.with(Round::TowardPositive, (dd(1.0, -f64::from_bits(1)), 2)),
948+
// Input: (2^2, -2^-1072 + -2^-1074)
949+
TestCase::new(dd(4.0, -5.0 * f64::from_bits(1)))
950+
.with_all((dd(1.0, -f64::from_bits(1)), 2))
951+
.with(Round::TowardNegative, (dd(1.0, -2.0 * f64::from_bits(1)), 2))
952+
.with(Round::TowardZero, (dd(1.0, -2.0 * f64::from_bits(1)), 2)),
953+
// Input: (2^2, -2^-1072)
954+
TestCase::new(dd(4.0, -4.0 * f64::from_bits(1))).with_all((dd(1.0, -f64::from_bits(1)), 2)),
955+
// Input: (2^2, -2^-1073 + -2^-1074)
956+
TestCase::new(dd(4.0, -3.0 * f64::from_bits(1)))
957+
.with_all((dd(1.0, -f64::from_bits(1)), 2))
958+
.with(Round::TowardPositive, (dd(0.5, -0.0), 3)),
959+
// Input: (2^2, -2^-1073)
960+
TestCase::new(dd(4.0, -2.0 * f64::from_bits(1)))
961+
.with_all((dd(0.5, -0.0), 3))
962+
.with(Round::TowardNegative, (dd(1.0, -f64::from_bits(1)), 2))
963+
.with(Round::TowardZero, (dd(1.0, -f64::from_bits(1)), 2)),
964+
// Input: (2^2, -2^-1074)
965+
TestCase::new(dd(4.0, -f64::from_bits(1)))
966+
.with_all((dd(0.5, -0.0), 3))
967+
.with(Round::TowardNegative, (dd(1.0, -f64::from_bits(1)), 2))
968+
.with(Round::TowardZero, (dd(1.0, -f64::from_bits(1)), 2)),
969+
// Input: 3+3*2^-53 canonicalized to (3+2^-51, -2^-53)
970+
// Output: 0.75+0.75*2^-53 canonicalized to (.75+2^-53, -2^-55)
971+
TestCase::new(dd(f64::from_bits(0x4008_0000_0000_0001), -2.0f64.powi(-53)))
972+
.with_all((dd(f64::from_bits(0x3fe8_0000_0000_0001), -2.0f64.powi(-55)), 2)),
973+
// Input: (2^1021+2^969, 2^968-2^915)
974+
TestCase::new(dd(f64::from_bits(0x7fc0_0000_0000_0001), f64::from_bits(0x7c6f_ffff_ffff_ffff)))
975+
.with_all((dd(f64::from_bits(0x3fe0_0000_0000_0001), f64::from_bits(0x3c8f_ffff_ffff_ffff)), 1022)),
976+
// Input: (2^1023, -2^-1)
977+
TestCase::new(dd(2.0f64.powi(1023), -0.5)).with_all((dd(1.0, -f64::from_bits(1 << 50)), 1023)),
978+
// Input: (2^1023, -2^-51)
979+
TestCase::new(dd(2.0f64.powi(1023), -2.0f64.powi(-51))).with_all((dd(1.0, -f64::from_bits(1)), 1023)),
980+
// Input: (2^1023, -2^-52)
981+
TestCase::new(dd(2.0f64.powi(1023), -2.0f64.powi(-52)))
982+
.with_all((dd(0.5, -0.0), 1024))
983+
.with(Round::TowardNegative, (dd(1.0, -f64::from_bits(1)), 1023))
984+
.with(Round::TowardZero, (dd(1.0, -f64::from_bits(1)), 1023)),
985+
// Input: (2^1023, 2^-1074)
986+
TestCase::new(dd(2.0f64.powi(1023), f64::from_bits(1)))
987+
.with_all((dd(0.5, 0.0), 1024))
988+
.with(Round::TowardPositive, (dd(0.5, f64::from_bits(1)), 1024)),
989+
// Input: (2^1024-2^971, 2^970-2^918)
990+
TestCase::new(DoubleDouble::largest())
991+
.with_all((dd(f64::from_bits(0x3fef_ffff_ffff_ffff), f64::from_bits(0x3c8f_ffff_ffff_fffe)), 1024)),
992+
];
993+
994+
let negate = |test_case: TestCase<DoubleDouble, (DoubleDouble, ExpInt)>| TestCase {
995+
input: -test_case.input,
996+
nearest_ties_to_even: test_case.nearest_ties_to_even.map(|(v, e)| (-v, e)),
997+
toward_positive: test_case.toward_negative.map(|(v, e)| (-v, e)),
998+
toward_negative: test_case.toward_positive.map(|(v, e)| (-v, e)),
999+
toward_zero: test_case.toward_zero.map(|(v, e)| (-v, e)),
1000+
nearest_ties_to_away: test_case.nearest_ties_to_away.map(|(v, e)| (-v, e)),
1001+
};
1002+
1003+
let mut actual_exp = 0;
1004+
1005+
for case in frexp_test_cases.iter().flat_map(|v| [*v, negate(*v)]) {
1006+
if let Some((expected, expected_exp)) = case.nearest_ties_to_even {
1007+
assert_eq!(case.input.frexp_r(&mut actual_exp, Round::NearestTiesToEven), expected);
1008+
assert_eq!(expected_exp, actual_exp);
1009+
}
1010+
1011+
if let Some((expected, expected_exp)) = case.nearest_ties_to_away {
1012+
assert_eq!(case.input.frexp_r(&mut actual_exp, Round::NearestTiesToAway), expected);
1013+
assert_eq!(expected_exp, actual_exp);
1014+
}
1015+
1016+
if let Some((expected, expected_exp)) = case.toward_positive {
1017+
assert_eq!(case.input.frexp_r(&mut actual_exp, Round::TowardPositive), expected);
1018+
assert_eq!(expected_exp, actual_exp);
1019+
}
1020+
1021+
if let Some((expected, expected_exp)) = case.toward_negative {
1022+
assert_eq!(case.input.frexp_r(&mut actual_exp, Round::TowardNegative), expected);
1023+
assert_eq!(expected_exp, actual_exp);
1024+
}
1025+
1026+
if let Some((expected, expected_exp)) = case.toward_zero {
1027+
assert_eq!(case.input.frexp_r(&mut actual_exp, Round::TowardZero), expected);
1028+
assert_eq!(expected_exp, actual_exp);
1029+
}
1030+
}
1031+
}
7911032
}

‎tests/ppc.rs‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -467,4 +467,10 @@ fn ppc_double_double_frexp() {
467467
let result = DoubleDouble::from_bits(input).frexp(&mut exp);
468468
assert_eq!(2, exp);
469469
assert_eq!(0x3c98000000000000_3fe8000000000000, result.to_bits());
470+
471+
// frexp quiets NaN.
472+
let snan = DoubleDouble::snan(None);
473+
let mut exp = 0;
474+
let result = snan.frexp(&mut exp);
475+
assert!(!result.is_signaling());
470476
}

0 commit comments

Comments
 (0)