@@ -32,31 +32,87 @@ fn power_of_ten(e: i16) -> Fp {
3232 Fp { f : sig, e : exp }
3333}
3434
35+ // Most architectures floating point operations with explicit bit size, therefore the precision of
36+ // the computation is determined on a per-operation basis.
3537#[ cfg( any( not( target_arch="x86" ) , target_feature="sse2" ) ) ]
3638mod fpu_precision {
3739 pub fn set_precision < T > ( ) { }
3840}
3941
42+ // On x86, the x87 FPU is used for float operations if the SSE[2] extensions are not available.
43+ // The x87 FPU operates with 80 bits of precision by default, which means that operations will
44+ // round to 80 bits causing double rounding to happen when values are eventually represented as
45+ // 32/64 bit float values. To overcome this, the FPU control word can be set so that the
46+ // computations are performed in the desired precision.
4047#[ cfg( all( target_arch="x86" , not( target_feature="sse2" ) ) ) ]
4148mod fpu_precision {
4249 use mem:: size_of;
4350 use ops:: Drop ;
4451
52+ /// A structure used to preserve the original value of the FPU control word, so that it can be
53+ /// restored when the structure is dropped.
54+ ///
55+ /// The x87 FPU is a 16-bits register whose fields are as follows:
56+ ///
57+ /// 1111 11
58+ /// 5432 10 98 76 5 4 3 2 1 0
59+ /// +----+--+--+--+-+-+-+-+-+-+
60+ /// | |RC|PC| |P|U|O|Z|D|I|
61+ /// | | | | |M|M|M|M|M|M|
62+ /// +----+--+--+--+-+-+-+-+-+-+
63+ /// The fields are:
64+ /// - Invalid operation Mask
65+ /// - Denormal operand Mask
66+ /// - Zero divide Mask
67+ /// - Overflow Mask
68+ /// - Underflow Mask
69+ /// - Precision Mask
70+ /// - Precision Control
71+ /// - Rounding Control
72+ ///
73+ /// The fields with no name are unused (on FPUs more modern than 287).
74+ ///
75+ /// The 6 LSBs (bits 0-5) are the exception mask bits; each blocks a specific type of floating
76+ /// point exceptions from being raised.
77+ ///
78+ /// The Precision Control field determines the precision of the operations performed by the
79+ /// FPU. It can set to:
80+ /// - 0b00, single precision i.e. 32-bits
81+ /// - 0b10, double precision i.e. 64-bits
82+ /// - 0b11, double extended precision i.e. 80-bits (default state)
83+ /// The 0b01 value is reserved and should not be used.
84+ ///
85+ /// The Rounding Control field determines how values which cannot be represented exactly are
86+ /// rounded. It can be set to:
87+ /// - 0b00, round to nearest even (default state)
88+ /// - 0b01, round down (toward -inf)
89+ /// - 0b10, round up (toward +inf)
90+ /// - 0b11, round toward 0 (truncate)
4591 pub struct FPUControlWord ( u16 ) ;
4692
4793 fn set_cw ( cw : u16 ) {
4894 unsafe { asm ! ( "fldcw $0" :: "m" ( cw) ) :: "volatile" }
4995 }
5096
97+ /// Set the precision field of the FPU to `T` and return a `FPUControlWord`
5198 pub fn set_precision < T > ( ) -> FPUControlWord {
5299 let cw = 0u16 ;
100+
101+ // Compute the value for the Precision Control field that is appropriate for `T`.
53102 let cw_precision = match size_of :: < T > ( ) {
54103 4 => 0x0000 , // 32 bits
55104 8 => 0x0200 , // 64 bits
56105 _ => 0x0300 , // default, 80 bits
57106 } ;
107+
108+ // Get the original value of the control word to restore it later, when the
109+ // `FPUControlWord` structure is dropped
58110 unsafe { asm ! ( "fnstcw $0" : "=*m" ( & cw) ) :: : "volatile" }
111+
112+ // Set the control word to the desired precision. This is achieved by masking away the old
113+ // precision (bits 8 and 9, 0x300) and replacing it with the precision flag computed above.
59114 set_cw ( ( cw & 0xFCFF ) | cw_precision) ;
115+
60116 FPUControlWord ( cw)
61117 }
62118
@@ -71,10 +127,6 @@ mod fpu_precision {
71127///
72128/// This is extracted into a separate function so that it can be attempted before constructing
73129/// a bignum.
74- ///
75- /// The fast path crucially depends on arithmetic being correctly rounded, so on x86
76- /// without SSE or SSE2 it requires the precision of the x87 FPU stack to be changed
77- /// so that it directly rounds to 64/32 bit.
78130pub fn fast_path < T : RawFloat > ( integral : & [ u8 ] , fractional : & [ u8 ] , e : i64 ) -> Option < T > {
79131 let num_digits = integral. len ( ) + fractional. len ( ) ;
80132 // log_10(f64::max_sig) ~ 15.95. We compare the exact value to max_sig near the end,
@@ -91,11 +143,16 @@ pub fn fast_path<T: RawFloat>(integral: &[u8], fractional: &[u8], e: i64) -> Opt
91143 return None ;
92144 }
93145
146+ // The fast path crucially depends on arithmetic being rounded to the correct number of bits
147+ // without any intermediate rounding. On x86 (without SSE or SSE2) this requires the precision
148+ // of the x87 FPU stack to be changed so that it directly rounds to 64/32 bit.
149+ // The `set_precision` function takes care of setting the precision on architectures which
150+ // require setting it by changing the global state (like the control word of the x87 FPU).
94151 let _cw = fpu_precision:: set_precision :: < T > ( ) ;
95152
96153 // The case e < 0 cannot be folded into the other branch. Negative powers result in
97154 // a repeating fractional part in binary, which are rounded, which causes real
98- // (and occasioally quite significant!) errors in the final result.
155+ // (and occasionally quite significant!) errors in the final result.
99156 if e >= 0 {
100157 Some ( T :: from_int ( f) * T :: short_fast_pow10 ( e as usize ) )
101158 } else {
0 commit comments