@@ -32,31 +32,87 @@ fn power_of_ten(e: i16) -> Fp {
32
32
Fp { f : sig, e : exp }
33
33
}
34
34
35
+ // Most architectures floating point operations with explicit bit size, therefore the precision of
36
+ // the computation is determined on a per-operation basis.
35
37
#[ cfg( any( not( target_arch="x86" ) , target_feature="sse2" ) ) ]
36
38
mod fpu_precision {
37
39
pub fn set_precision < T > ( ) { }
38
40
}
39
41
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.
40
47
#[ cfg( all( target_arch="x86" , not( target_feature="sse2" ) ) ) ]
41
48
mod fpu_precision {
42
49
use mem:: size_of;
43
50
use ops:: Drop ;
44
51
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)
45
91
pub struct FPUControlWord ( u16 ) ;
46
92
47
93
fn set_cw ( cw : u16 ) {
48
94
unsafe { asm ! ( "fldcw $0" :: "m" ( cw) ) :: "volatile" }
49
95
}
50
96
97
+ /// Set the precision field of the FPU to `T` and return a `FPUControlWord`
51
98
pub fn set_precision < T > ( ) -> FPUControlWord {
52
99
let cw = 0u16 ;
100
+
101
+ // Compute the value for the Precision Control field that is appropriate for `T`.
53
102
let cw_precision = match size_of :: < T > ( ) {
54
103
4 => 0x0000 , // 32 bits
55
104
8 => 0x0200 , // 64 bits
56
105
_ => 0x0300 , // default, 80 bits
57
106
} ;
107
+
108
+ // Get the original value of the control word to restore it later, when the
109
+ // `FPUControlWord` structure is dropped
58
110
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.
59
114
set_cw ( ( cw & 0xFCFF ) | cw_precision) ;
115
+
60
116
FPUControlWord ( cw)
61
117
}
62
118
@@ -71,10 +127,6 @@ mod fpu_precision {
71
127
///
72
128
/// This is extracted into a separate function so that it can be attempted before constructing
73
129
/// 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.
78
130
pub fn fast_path < T : RawFloat > ( integral : & [ u8 ] , fractional : & [ u8 ] , e : i64 ) -> Option < T > {
79
131
let num_digits = integral. len ( ) + fractional. len ( ) ;
80
132
// 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
91
143
return None ;
92
144
}
93
145
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).
94
151
let _cw = fpu_precision:: set_precision :: < T > ( ) ;
95
152
96
153
// The case e < 0 cannot be folded into the other branch. Negative powers result in
97
154
// 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.
99
156
if e >= 0 {
100
157
Some ( T :: from_int ( f) * T :: short_fast_pow10 ( e as usize ) )
101
158
} else {
0 commit comments