使用電腦去進行浮點數運算的時候，在某些條件下會出現顯而易見的誤差，而這樣的誤差，就會導致與浮點數有關的計算或是邏輯判斷出現問題。
在這篇文章開始前，您可以先看看以下的程式：
RustJavaTypeScript
fn main() {
println!("0.1 + 1.1 = {}", 0.1 + 1.1);
println!("0.1 * 1.1 = {}", 0.1 * 1.1);
}
public static void main(final String[] args) {
System.out.println("0.1 + 1.1 = " + (0.1 + 1.1));
System.out.println("0.1 * 1.1 = " + (0.1 * 1.1));
}
console.log(`0.1 + 1.1 = ${0.1 + 1.1}`);
console.log(`0.1 * 1.1 = ${0.1 * 1.1}`);
以上程式的輸出結果如下：
0.1 + 1.1 = 1.2000000000000002
如何？很奇怪吧！
嘗試寫個程式判斷兩個浮點數是否能整除
當我們要判斷兩個數是否能整除的時候，只要將這兩個數相除，再看看除出來的商是否為整數，如果是整數，表示被除數是除數的整數倍。讓我們直接根據這個作法來寫個程式判斷兩個浮點數是否能整除，程式碼如下：
RustJavaTypeScript
pub fn is_divisible(dividend: f64, divisor: f64) -> bool {
let quotient = dividend / divisor;
quotient == quotient.round()
}
public static boolean isDivisible(final double dividend, final double divisor) {
final double quotient = dividend / divisor;
return quotient == Math.round(quotient);
}
function isDivisible(dividend: number, divisor: number) {
const quotient = dividend / divisor;
return quotient == Math.round(quotient);
}
以上的is_divisible或isDivisible函數，當傳入的引數是0.1和0.000001時，就會回傳false(不能整除)。但0.1很顯然是可以被0.000001整除的。問題出現的原因是在於電腦計算0.1/0.000001之後得到了100000.00000000001這樣的非整數結果。
這裡要特別注意的是，在判斷一個浮點數是否為整數的時候，不要使用ceil或是floor函數來消除小數，因為有誤差的值可能會比原值大或是小。好比1/0.00001會算出99999.99999999999，如果我們是用floor直接消除小數，那麼99999.99999999999和99999將會差距巨大，使我們無法正確判斷1/0.00001的結果是否為整數。
判斷兩個浮點數是否相等
在大部分的情況下，要寫程式去判斷兩個浮點數型別的值(#{{a}}#、#{{b}}#)是否相等的時候，要先將兩數相減所得出的差(difference)取絕對值，再判斷差的絕對值是否小於某個誤差值(epsilon)，而不是直接去進行相等的邏輯運算。即：
至於這個#{{\epsilon}}#的值該怎麼定，就看我們要算得多精準。首先要先決定好是要用單倍精準浮點數，還是雙倍精準浮點數。
浮點數是以二進制的科學記號(scientific notation，或稱科學計數法)來表示數值，即#{{\pm a \times 2^n}}#。浮點數的第一個位元用來表示正負號(sign)；中間段的幾個位元用來表示科學記號的指數#{{n}}#，尾段的幾個位元用來表示科學記號的實數#{{a}}#的小數部份(尾數)。
單倍精準浮點數 有32個位元，尾段的23個位元用來表示實數#{{a}}#的尾數，共可表示#{{2^{23} = 8388608}}#個尾數，所以32位元的浮點數的精確度為#{{\lfloor log_{10}{2^{23}} \rfloor = 6}}#個位數 。它可以保證在不計算開頭那些連續的0的位數下(像是0.00106只能算後面的3個位數)，至少能精確地從前面開始表示6個位數。
您可能會有看過32位元的浮點數的有效位數有7個或是6至7個的說法，這跟指數#{{n}}#是有關係的。
當指數#{{n}}#是0時，32位元的浮點數可表示的數值範圍為[1.0, 2.0)(先不管正負號，不重要)。32位元的浮點數可以將這個範圍的數值平均切成8_388_608段，此時的精準度就是7個位數，因為0.000000、0.000001、……、0.999999也不過999_999個數字，離「32位元的浮點數可以表示出8_388_608個尾數」還差得很遠呢！把精確的小數點後6個位數再加上小數點前的1個位數(也就是1)，不就是7位了嗎？
當指數#{{n}}#是23時，浮點數可能的數值範圍為[223 , 224 )，也就是[8_388_608, 16_777_216)，32位元的浮點數可以將這個範圍的數值平均切成8_388_608段，若以精準7位數來看的話，[9_000_000, 10_000_000)共有1_000_000個數字要表示，且因為這個範圍的每個32位元的浮點數相差1，剛好可以對應這1_000_000個數字，所以此時是有7位數的精確度的。
當指數#{{n}}#是33時，浮點數可能的數值範圍為[233 , 234 )，也就是[8_589_934_592, 17_179_869_184)，32位元的浮點數可以將這個範圍的數值平均切成8_388_608段，若以精準7位數來看的話，[9_000_000_000, 10_000_000_000)共有1_000_000個數字要表示(9_000_000_000, 9_000_001_000, 9_000_002_000, ...)，但這個範圍的每個32位元的浮點數相差1024，只能夠表示其中的976_562個數字，所以精準度就不足7位數了！
至於小於1的數字，也可以找到有效位數非7位數的例子，像是0.0009962818，32位元的浮點數就只能夠使用最接近的0.0009962817來表示它。
雙倍精準浮點數 有64個位元，尾段的52個位元用來表示實數#{{a}}#的尾數，所以64位元的浮點數的精確度為#{{\lfloor log_{10}{2^{52}} \rfloor = 15}}#個位數 。
與32位元的浮點數同理，也有64位元的浮點數的有效位數有15個或是15至16個的說法，它們表示的數值範圍會影響到它們的有效位數。
判斷兩個浮點數是否相等的時候，若要設定誤差能容忍到小數點後第n位，可以將#{{\epsilon}}#設定為10-n (1e-n)。但這樣的設定方式會有個問題是我們的浮點數有效位可能不一定都會靠近小數點，好比我們設定了n = 6，但使用的浮點數的絕對值都小於10-6 ，此時會將所有的浮點數都視為相等。
還有一種判斷兩個浮點數是否相等方式，要利用大於1的最小浮點數與1的差，即「machine epsilon」，程式語言通常會以常數的方式提供這個值。不過也要把判斷兩個浮點數是否相等的方式改為以下的樣子以符合machine epsilon的定義，不然像是0.1234567和0.1234568這兩個32位元的浮點數，就會被判定為是一樣的浮點數。
#{{{
為了避免#{{b = 0}}#導致除以0的情況發生，我們或許可以對不等式做一些改變，如下：
#{{{
以上的不等式看起來可行，但實際使用會發現當我們等比例增大或減少#{{a}}#、#{{b}}#時，#{{\mid {a - b} \mid}}#的誤差值並不會等比例做變化，所以要把#{{\mid b \mid \times \epsilon}}#再乘上某值來擴大容許的誤差範圍才行，筆者覺得變成#{{8}}#倍還不錯。
此外，浮點數也有可能會是無限大的，因此在程式實作時，也要納入考量。
程式實作如下：
RustJavaTypeScript
pub fn is_equal(a: f64, b: f64) -> bool {
if a == b {
true
} else {
let a_abs = a.abs();
let b_abs = b.abs();
let d_abs = (a - b).abs();
if a_abs >= b_abs {
d_abs < b_abs * f64::EPSILON * 8.0
} else {
d_abs < a_abs * f64::EPSILON * 8.0
}
}
}
a == b可以判斷無限大是否相等。
public static boolean isEqual(final double a, final double b) {
if (a == b) {
return true;
} else {
final double aAbs = Math.abs(a);
final double bAbs = Math.abs(b);
final double dAbs = Math.abs(a - b);
if (aAbs >= bAbs) {
return dAbs < bAbs * Math.ulp(1.0) * 8.0;
} else {
return dAbs < aAbs * Math.ulp(1.0) * 8.0;
}
}
}
a == b可以判斷無限大是否相等。Java並沒有提供machine epsilon常數，所以用Math.ulp(1.0)函數來算出。
function isEqual(a: number, b: number) {
if (a === b) {
return true;
} else {
const aAbs = Math.abs(a);
const bAbs = Math.abs(b);
const dAbs = Math.abs(a - b);
if (aAbs >= bAbs) {
return dAbs < bAbs * Number.EPSILON * 8.0;
} else {
return dAbs < aAbs * Number.EPSILON * 8.0;
}
}
}
a === b可以判斷無限大是否相等。
以上程式還可以再化簡如下：
RustJavaTypeScript
pub fn is_equal(a: f64, b: f64) -> bool {
if a == b {
true
} else {
let a_abs = a.abs();
let b_abs = b.abs();
let d_abs = (a - b).abs();
d_abs < a_abs.min(b_abs) * f64::EPSILON * 8.0
}
}
public static boolean isEqual(final double a, final double b) {
if (a == b) {
return true;
} else {
final double aAbs = Math.abs(a);
final double bAbs = Math.abs(b);
final double dAbs = Math.abs(a - b);
return dAbs < Math.min(aAbs, bAbs) * Math.ulp(1.0) * 8.0;
}
}
function isEqual(a: number, b: number) {
if (a === b) {
return true;
} else {
const aAbs = Math.abs(a);
const bAbs = Math.abs(b);
const dAbs = Math.abs(a - b);
return dAbs < Math.min(aAbs, bAbs) * Number.EPSILON * 8.0;
}
}
判斷兩個浮點數是否能整除
稍微修改以上的is_divisible或isDivisible函數，結合我們後來想出的判斷浮點數是否相等的作法，再利用大於0的最小浮點數來判斷被除數和除數是否(近似)等於0，可以實作出如下的程式：
RustJavaTypeScript
pub fn is_divisible(dividend: f64, divisor: f64) -> bool {
if divisor.abs() < f64::MIN_POSITIVE {
dividend.abs() < f64::MIN_POSITIVE
} else {
let quotient = dividend / divisor;
is_equal(quotient, quotient.round())
}
}
Rust的f64::MIN_POSITIVE或f32::MIN_POSITIVE是大於0的最小正規數(normal number)，0到最小正規數之間還有許多次正規數(subnormal number)。這樣做可以避免在某些機器上面，0到最小正規數之間的次正規數被當作0。如果要使用大於0的最小次正規數，可以直接使用浮點數定數4.9406564584124654e-324(64位元)、1.401298e-45(32位元)。
public static boolean isDivisible(final double dividend, final double divisor) {
final double DBL_MIN_POSITIVE = Math.ulp(0.0);
if (Math.abs(divisor) < DBL_MIN_POSITIVE) {
return Math.abs(dividend) < DBL_MIN_POSITIVE;
} else {
final double quotient = dividend / divisor;
return isEqual(quotient, Math.round(quotient));
}
}
Java的Math.ulp函數會找出次正規數(subnormal number)間的差。
function isDivisible(dividend: number, divisor: number) {
const DBL_MIN_POSITIVE = Number.MIN_VALUE;
if (Math.abs(divisor) < DBL_MIN_POSITIVE) {
return Math.abs(dividend) < DBL_MIN_POSITIVE;
} else {
const quotient = dividend / divisor;
return isEqual(quotient, Math.round(quotient));
}
}
如果用相除後餘數(remainder)是否為0來判斷兩個浮點數能不能整除呢？大部份程式語言不是都會有提供餘數運算子%，如果沒有也會有相關的mod函數吧？
您可以先看看以下的程式：
RustJavaTypeScript
fn main() {
println!("1.0 % 0.5 = {}", 1.0 % 0.5);
println!("1.0 % 0.2 = {}", 1.0 % 0.2);
println!("1.0 % 0.1 = {}", 1.0 % 0.1);
}
public static void main(final String[] args) {
System.out.println("1.0 % 0.5 = " + (1.0 % 0.5));
System.out.println("1.0 % 0.2 = " + (1.0 % 0.2));
System.out.println("1.0 % 0.1 = " + (1.0 % 0.1));
}
console.log(`1.0 % 0.5 = ${1.0 % 0.5}`);
console.log(`1.0 % 0.2 = ${1.0 % 0.2}`);
console.log(`1.0 % 0.1 = ${1.0 % 0.1}`);
以上程式的輸出結果如下：
1.0 % 0.5 = 0
如何？很奇怪吧！這已經不是餘數是不是近似於0的問題了，0.19999999999999996和0明顯差很多。所以我們不能 用這樣的方式來找浮點數相除的餘數。
計算兩個浮點數相除後的餘數
注意，這邊我們定義用商取最小整數後乘上除數，再用被除數去減掉所得到的乘積即為餘數，所以餘數正負號會跟除數一樣。(Python的%運算子所得到的餘數也是如此，但其它程式語言就不一定了。)
要計算#{{a}}#、#{{b}}#兩個數相除(#{{a \div b}}#)後的餘數，可以使用以下公式：
#{{{
然而我們知道電腦在運算#{{a \over b}}#的時候在有可能會因為出現誤差而導致整數部份有所變化，所以我們不能直接使用floor函數來計算餘數。
我們可以先將#{{a \over b}}#的商代入ceil函數，再判斷計算結果與商是否(近似)相等。如果是的話就使用ceil函數算出來的結果作為#{{\lfloor {a \over b} \rfloor}}#的結果；如果不是的話才使用floor函數來計算。
程式實作如下：
RustJavaTypeScript
pub fn calculate_remainder(dividend: f64, divisor: f64) -> f64 {
let quotient = dividend / divisor;
let quotient_ceil = quotient.ceil();
let floor = if is_equal(quotient, quotient_ceil) {
quotient_ceil
} else {
quotient.floor()
};
dividend - (floor * divisor)
}
public static double calculateRemainder(final double dividend, final double divisor) {
final double quotient = dividend / divisor;
final double quotientCeil = Math.ceil(quotient);
final double floor;
if (isEqual(quotient, quotientCeil)) {
floor = quotientCeil;
} else {
floor = Math.floor(quotient);
}
return dividend - (floor * divisor);
}
function calculateRemainder(dividend: number, divisor: number) {
const quotient = dividend / divisor;
const quotientCeil = Math.ceil(quotient);
let floor;
if (isEqual(quotient, quotientCeil)) {
floor = quotientCeil;
} else {
floor = Math.floor(quotient);
}
return dividend - floor * divisor;
}
