Java™ 語言規范第 5 版向 java.lang.Math 和 java.lang.StrictMath 添加了 10 種新方法,Java 6 又添加了 10 種。這個共兩部分的系列文章的 第 1 部分 介紹了很有意 義的新的數學方法。它提供了在還未出現計算機的時代中數學家比較熟悉的函數。在第 2 部 分中,我主要關注這樣一些函數,它們的目的是操作浮點數,而不是抽象實數。
就像 我在 第 1 部分中 提到的一樣,實數(比如 e 或 0.2)和它的 計算機表示(比如 Java double)之間的區別是非常重要的。最理想的數字應該是無限精確 的,然而 Java 表示的位數是固定的(float 為 32 位,double 為 64 位)。float 的最大 值約為 3.4*1038。這個值還不足以表示某些東西,比如宇宙中的電子數目。
double 的最大值為 1.8*10308,幾乎能夠表示任何物理量。不過涉及到 抽象數學量的計算時,可能超出這些值的范圍。例如,光是 171! (171 * 170 * 169 * 168 * ... * 1) 就超出了 double 最大值。float 只能表示 35! 以內的數字。非常小的數(值 接近於 0 的數字)也會帶來麻煩,同時涉及到非常大的數和非常小的數的計算是非常危險的 。
為了處理這個問題,浮點數學 IEEE 754 標准(參見 參考資料)添加了特殊值 Inf 和 NaN,它們分別表示無窮大(Infinity)和非數字(Not a Number)。IEEE 754 還定 義了正 0 和負 0(在一般的數學中,0 是不分正負的,但在計算機數學中,它們可以是正的 ,也可以是負的)。這些值給傳統的原理帶來了混亂。例如,當使用 NaN 時,排中律就不成 立了。x == y 或 x != y 都有可能是不正確的。當 x 或 y 為 NaN 時,這兩個式子都不成 立。
除了數字大小問題外,精度是一個更加實際的問題。看看這個常見的循環,將 1.0 相加 10 次之後等到的結果不是 10,而是 9.99999999999998:
for (double x = 0.0; x <= 10.0; x += 0.1) {
System.err.println(x);
}
對於簡單的應用程序,您通常讓 java.text.DecimalFormat 將最終的輸出格式化為與其 值最接近的整數,這樣就可以了。不過,在科學和工程應用方面(您不能確定計算的結果是 否為整數),則需要加倍小心。如果需要在特別大的數字之間執行減法以得到較小的數字, 則需要萬分 小心。如果以特別小的數字作為除數,也需要加以注意。這些操作能夠將很小的 錯誤變成大錯誤,並給現實應用帶來巨大的影響。由有限精度浮點數字引起的很小的捨入錯 誤就會嚴重歪曲數學精度計算。
浮點數和雙精度數字的二進制表示
由 Java 實現的 IEEE 754 浮點數有 32 位。第一位是符號位,0 表示正,1 表示負。接 下來的 8 位表示指數,其值的范圍是 -125 到 +127。最後的 23 位表示尾數(有時稱為有 效數字),其值的范圍是 0 到 33,554,431。綜合起來,浮點數是這樣表示的: sign * mantissa * 2exponent 。
敏銳的讀者可能已經注意到這些數字有些不對勁。首先,表示指數的 8 位應該是從 -128 到 127,就像帶符號的字節一樣。但是這些指數的偏差是 126,即用不帶符號的值(0 到 255)減去 126 獲得真正的指數(現在是從 -126 到 128)。但是 128 和 -126 是特殊值。 當指數都是 1 位(128)時,則表明這個數字是 Inf、-Inf 或 NaN。要確定具體情況,必須 查看它的尾數。當指數都是 0 位(-126)時,則表明這個數字是不正常的(稍後將詳細介紹 ),但是指數仍然是 -125。
尾數一般是一個 23 位的不帶符號的整數 — 它非常簡單。23 位可以容納 0 到 224-1,即 16,777,215。等一下,我剛才是不是說尾數的范圍是從 0 到 33,554,431?即 225-1。多出的一位是從哪裡來的?
因此,可以通過指數表示第 1 位是什麼。如果指數都是 0 位,則第 1 位為 0。否則第 1 位為 1。因為我們通常知道第 1 位是什麼,所以沒有必要包含在數字中。您 “免費” 得 到一個額外的位。是不是有些離奇?
尾數的第 1 位為 1 的浮點數是正常的。即尾數的值通常在 1 到 2 之間。尾數的第 1 位為 0 的浮點數是不正常的,盡管指數通常為 -125,但它通常能夠表示更小的數字。
雙精度數是以類似的方式編碼的,但是它使用 52 位的尾數和 11 位的指數來獲得更高的 精度。雙精度數的指數的偏差是 1023。
尾數和指數
在 Java 6 中添加的兩個 getExponent() 方法在表示浮點數或雙精度數時返回無偏差 指 數。對於浮點數,這個數字的范圍是 -125 到 +127,對於雙精度數,這個數字的范圍是 - 1022 到 +1023(Inf 和 NaN 為 +128/+1024)。例如,清單 1 根據更常見的以 2 為底數的 對數比較了 getExponent() 方法的結果:
清單 1. Math.log(x)/Math.log(2) 和 Math.getExponent()
public class ExponentTest {
public static void main(String[] args) {
System.out.println("x\tlg(x)\tMath.getExponent(x)");
for (int i = -255; i < 256; i++) {
double x = Math.pow(2, i);
System.out.println(
x + "\t" +
lg(x) + "\t" +
Math.getExponent(x));
}
}
public static double lg(double x) {
return Math.log(x)/Math.log(2);
}
}
對於使用捨入的一些值,Math.getExponent() 比一般的計算要准確一些:
x lg(x) Math.getExponent(x)
...
2.68435456E8 28.0 28
5.36870912E8 29.000000000000004 29
1.073741824E9 30.0 30
2.147483648E9 31.000000000000004 31
4.294967296E9 32.0 32
如果要執行大量此類計算,Math.getExponent() 會更快。不過需要注意,它僅適用於計 算 2 的冪次方。例如,當改為 3 的冪次方時,結果如下:
x lg(x) Math.getExponent(x)
...
1.0 0.0 0
3.0 1.584962500721156 1
9.0 3.1699250014423126 3
27.0 4.754887502163469 4
81.0 6.339850002884625 6
getExponent() 不處理尾數,尾數由 Math.log() 處理。通過一些步驟,就可以找到尾數 、取尾數的對數並將該值添加到指數,但這有些費勁。如果想要快速估計數量級(而不是精 確值),Math.getExponent() 是非常有用的。
與 Math.log() 不同,Math.getExponent() 從不返回 NaN 或 Inf。如果參數為 NaN 或 Inf,則對應的浮點數和雙精度數的結果分別是 128 和 1024。如果參數為 0,則對應的浮點 數和雙精度數的結果分別是 -127 和 -1023。如果參數為負數,則數字的指數與該數字的絕 對值的指數相同。例如,-8 的指數為 3,這與 8 的指數相同。
沒有對應的 getMantissa() 方法,但是使用簡單的數學知識就能構造一個:
public static double getMantissa(double x) {
int exponent = Math.getExponent(x);
return x / Math.pow(2, exponent);
}
盡管算法不是很明顯,但還是可以通過位屏蔽來查找尾數。要提取位,僅需計算 Double.doubleToLongBits(x) & 0x000FFFFFFFFFFFFFL。不過,隨後則需要考慮正常數 字中多出的 1 位,然後再轉換回范圍在 1 到 2 之間的浮點數。
最小的精度單位
實數是非常密集的。任意兩個不同的實數中間都可以出現其他實數。但浮點數則不是這樣 。對於浮點數和雙精度數,也存在下一個浮點數;連續的浮點數和雙精度數之間存在最小的 有限距離。nextUp() 方法返回比第一個參數大的最近浮點數。例如,清單 2 打印出所有在 1.0 和 2.0 之間的浮點數:
清單 2. 計算浮點數數量
public class FloatCounter {
public static void main(String[] args) {
float x = 1.0F;
int numFloats = 0;
while (x <= 2.0) {
numFloats++;
System.out.println(x);
x = Math.nextUp(x);
}
System.out.println(numFloats);
}
}
結果是 1.0 和 2.0 之間包含 8,388,609 個浮點數;雖然很多,但還不至於是無窮多的 實數。相鄰數字的距離為 0.0000001。這個距離稱為 ULP,它是最小精度單位(unit of least precision) 或最後位置單位(unit in the last place)的縮略。
如果需要向後查找小於指定數字的最近浮點數,則可以改用 nextAfter() 方法。第二個 參數指定是否查找在第一個參數之上或之下的最近數字:
public static double nextAfter(float start, float direction)
public static double nextAfter(double start, double direction)
如果 direction 大於 start,則 nextAfter() 返回在 start 之上的下一個數字。如果 direction 小於 start,則 nextAfter() 返回在 start 之下的下一個數字。如果 direction 等於 start,則 nextAfter() 返回 start 本身。
這些方法在某些建模或圖形工具中是非常有用的。從數字上來說,您可能需要在 a 和 b 之間的 10,000 個位置上提取樣例值,但如果您具備的精度僅能識別 a 和 b 之間的 1,000 個獨立的點,那麼有十分之九的工作是重復的。您可以只做十分之一的工作,但又獲得相同 的結果。
當然,如果一定需要額外的精度,則可以選擇具有高精度的數據類型,比如 double 或 BigDecimal。例如,我曾經在 Mandelbrot 集合管理器看見過這種情況。在其中可以放大曲 線圖,讓其落在最近的兩個雙精度數之間。Mandelbrot 集合在各個級別上都是非常細微和復 雜的,但是 float 或 double 可以在失去區分相鄰點的能力之前達到這個細微的級別。
Math.ulp() 返回一個數字和距其最近的數字之間的距離。清單 3 列出了 2 的各種冪次 方的 ULP:
清單 3. 浮點數 2 的冪次方的 ULP
public class UlpPrinter {
public static void main(String[] args) {
for (float x = 1.0f; x <= Float.MAX_VALUE; x *= 2.0f) {
System.out.println(Math.getExponent(x) + "\t" + x + "\t" + Math.ulp(x));
}
}
}
下面給出了一些輸出:
0 1.0 1.1920929E-7
1 2.0 2.3841858E-7
2 4.0 4.7683716E-7
3 8.0 9.536743E-7
4 16.0 1.9073486E-6
...
20 1048576.0 0.125
21 2097152.0 0.25
22 4194304.0 0.5
23 8388608.0 1.0
24 1.6777216E7 2.0
25 3.3554432E7 4.0
...
125 4.2535296E37 5.0706024E30
126 8.507059E37 1.0141205E31
127 1.7014118E38 2.028241E31
可以看到,對於比較小的 2 的冪次方,浮點數是非常精確的。但是在許多應用程序中, 在數值約為 220 時,這一精度將出現問題。在接近浮點數的最大極限時,相鄰 的值將被 千的七乘方(sextillions)隔開(事實上可能更大一點,但我找不到詞匯來表達 )。
如清單 3 所示,ULP 的大小並不是固定的。隨著數字變大,它們之間的浮點數就會越來 越少。例如,10,000 和 10,001 之間只有 1,025 個浮點數;它們的距離是 0.001。在 1,000,000 和 1,000,001 之間僅有 17 個浮點數,它們的距離是 0.05。精度與數量級成反 比關系。對於浮點數 10,000,000,ULP 的精確度變為 1.0,超過這個數之後,將有多個整數 值映射到同一個浮點數。對於雙精度數,只有達到 4.5E15 時才會出現這種情況,但這也是 個問題。
Math.ulp() 為測試提供一個實用的用途。很明顯,我們一般不會比較兩個浮點數是否完 全相等。相反,我們檢查它們是否在一定的容錯范圍內相等。例如,在 JUnit 中,像以下這 樣比較預期的實際浮點值:
assertEquals(expectedValue, actualValue, 0.02);
這表明實際值與預期值的偏差在 0.02 之內。但是,0.02 是合理的容錯范圍嗎?如果預 期值是 10.5 或 -107.82,則 0.02 是完全可以接受的。但當預期值為幾十億時,0.02 則與 0 沒有什麼區別。通常,就 ULP 進行測試時考慮的是相對錯誤。一般選擇的容錯范圍在 1 至 10 ULP 之間,具體情況取決於計算所需的精度。例如,下面指定實際結果必須在真實值 的 5 個 ULP 之內:
assertEquals(expectedValue, actualValue, 5*Math.ulp (expectedValue));
根據期望值不同,這個值可以是萬億分之一,也可以是數百萬。
浮點數的有限精度會導致一個難以預料的結果:超過某個點時,x+1 == x 便是真的。例 如,下面這個簡單的循環實際上是無限的:
for (float x = 16777213f; x < 16777218f; x += 1.0f) {
System.out.println(x);
}
實際上,這個循環將在一個固定的點上停下來,准確的數字是 16,777,216。這個數字等 於 224,在這個點上,ULP 比增量大。
scalb
Math.scalb(x, y) 用 2y 乘以 x,scalb 是 “scale binary(二進法)” 的縮寫。
public static double scalb(float f, int scaleFactor)
public static double scalb(double d, int scaleFactor)
例如,Math.scalb(3, 4) 返回 3 * 24,即 3*16,結果是 48.0。也可以使 用 Math.scalb() 來實現 getMantissa():
public static double getMantissa(double x) {
int exponent = Math.getExponent(x);
return x / Math.scalb(1.0, exponent);
}
Math.scalb() 和 x*Math.pow(2, scaleFactor) 的區別是什麼?實際上,最終的結果是 一樣的。任何輸入返回的值都是完全一樣的。不過,性能方面則存在差別。Math.pow() 的性 能是非常糟糕的。它必須能夠真正處理一些非常少見的情況,比如對 3.14 采用冪 -0.078。 對於小的整數冪,比如 2 和 3(或以 2 為基數,這比較特殊),通常會選擇完全錯誤的算 法。
我擔心這會對總體性能產生影響。一些編譯器和 VM 的智能程度比較高。一些優化器會將 x*Math.pow(2, y) 識別為特殊情況並將其轉換為 Math.scalb(x, y) 或類似的東西。因此性 能上的影響體現不出來。不過,我敢保證有些 VM 是沒有這麼智能的。例如,使用 Apple 的 Java 6 VM 進行測試時,Math.scalb() 幾乎總是比 x*Math.pow(2, y) 快兩個數量級。當然 ,這通常不會造成影響。但是在特殊情況下,比如執行數百萬次求冪運算時,則需要考慮能 否轉換它們以使用 Math.scalb()。
Copysign
Math.copySign() 方法將第一個參數的標記設置為第二個參數的標記。最簡單的實現如清 單 4 所示:
清單 4. 可能實現的 copysign 算法
public static double copySign(double magnitude, double sign) {
if (magnitude == 0.0) return 0.0;
else if (sign < 0) {
if (magnitude < 0) return magnitude;
else return -magnitude;
}
else if (sign > 0) {
if (magnitude < 0) return -magnitude;
else return magnitude;
}
return magnitude;
}
不過,真正的實現如清單 5 所示:
清單 5. 來自 sun.misc.FpUtils 的真正算法
public static double rawCopySign(double magnitude, double sign) {
return Double.longBitsToDouble((Double.doubleToRawLongBits(sign) &
(DoubleConsts.SIGN_BIT_MASK)) |
(Double.doubleToRawLongBits(magnitude) &
(DoubleConsts.EXP_BIT_MASK |
DoubleConsts.SIGNIF_BIT_MASK)));
}
仔細觀察這些位就會看到,NaN 標記被視為正的。嚴格來說,Math.copySign() 並不對此 提供保證,而是由 StrictMath.copySign() 負責,但在現實中,它們都調用相同的位處理代 碼。
清單 5 可能會比清單 4 快一些,但它的主要目的是正確處理負 0。Math.copySign(10, -0.0) 返回 -10,而 Math.copySign(10, 0.0) 返回 10.0。清單 4 中最簡單形式的算法在 兩種情況中都返回 10.0。當執行敏感的操作時,比如用極小的負雙精度數除以極大的正雙精 度數,就可能出現負 0。例如,-1.0E-147/2.1E189 返回負 0,而 1.0E-147/2.1E189 返回 正 0。不過,使用 == 比較這兩個值時,它們是相等的。因此,如果要區分它們,必須使用 Math.copySign(10, -0.0) 或 Math.signum()(調用 Math.copySign(10, -0.0))來執行比 較。
對數和指數
指數函數是一個很好的例子,它表明處理有限精度浮點數(而不是無限精度實數)時是需 要非常小心的。在很多等式中都出現 ex(Math.exp())。例如,它可用於定義 cosh 函數,這已經在 第 1 部分中 討論:
cosh(x) = (ex + e-x)/2
不過,對於負值的 x,一般是 -4 以下的數字,用於計算 Math.exp() 的算法表現很差, 並且容易出現捨入錯誤。使用另一個算法計算 ex - 1 會更加精確,然後在最終 結果上加 1。Math.expm1() 能夠實現這個不同的算法(m1 表示 “減 1”)。例如,清單 6 給出的 cosh 函數根據 x 的大小在兩個算法之間進行切換:
清單 6. cosh 函數
public static double cosh(double x) {
if (x < 0) x = -x;
double term1 = Math.exp(x);
double term2 = Math.expm1(-x) + 1;
return (term1 + term2)/2;
}
這個例子有些呆板,因為在 Math.exp() 與 Math.expm1() + 1 之間的差別很明顯的情況 下,常常使用 ex,而不是 e-x。不過,Math.expm1() 在帶有多種 利率的金融計算中是非常實用的,比如短期國庫券的日利率。
Math.log1p() 與 Math.expm1() 剛好相反,就像 Math.log() 與 Math.exp() 的關系一 樣。它計算 1 的對數和參數(1p 表示 “加 1”)。在值接近 1 的數字中使用這個函數。 例如,應該使用它計算 Math.log1p(0.0002),而不是 Math.log(1.0002)。
現在舉一個例子,假設您需要知道在日利率為 0.03 的情況下,需要多少天投資才能使 $1,000 增長到 $1,100。清單 7 完成了這個計算任務:
清單 7. 計算從當前投資額增長到未來特定值所需的時間
public static double calculateNumberOfPeriods(
double presentValue, double futureValue, double rate) {
return (Math.log(futureValue) - Math.log(presentValue))/Math.log1p (rate);
}
在這個例子中,1p 的含義是很容易理解的,因為在計算類似數據的一般公式中通常出現 1+r。換句話說,盡管投資方很希望獲得初始投資成本的 (1+r)n ,貸方通常將 利率作為附加的百分比(+r 部分)。實際上,以 3% 的利率貸款的投資者如果僅能取回投資 成本的 3% 的話,那就非常糟糕了。
結束語
浮點數並不是實數。它們的數量是有限的。它們能夠表示最大和最小的值。更值得注意的 是,它們的精度雖然很高,但范圍很窄,並且容易出現捨入錯誤。相反,浮點數和雙精度數 處理整數時獲得的精度遠比整型數和長型數差。您必須仔細考慮這些限制,尤其是在科研和 工程應用方面,以生產出健壯、可靠的代碼。對於財務應用程序(尤其是需要精確到最後一 位的會計應用程序),處理浮點數和雙精度數時也需要格外小心。
java.lang.Math 和 java.lang.StrictMath 類經過了精心設計,可以解決這些問題。適 當地使用這些類及其包含的方法能夠改善程序。本文特別展示了良好的浮點算法有多麼巧妙 !最好使用專家提供的算法,而不是自己獨創算法。如果適合使用 java.lang.Math 和 java.lang.StrictMath 中提供的方法,最好繼續使用。它們通常是最佳的選擇。