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 * 2 exponent
。
敏锐的读者可能已经注意到这些数字有些不对劲。首先,表示指数的 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()
对于使用舍入的一些值,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 之间的浮点数:
结果是 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:
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)); |
根据期望值不同,这个值可以是万亿分之一,也可以是数百万。
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()
。
Math.copySign()
方法将第一个参数的标记设置为第二个参数的标记。最简单的实现如清单 4 所示:
不过,真正的实现如清单 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)
)来执行比较。
指数函数是一个很好的例子,它表明处理有限精度浮点数(而不是无限精度实数)时是需要非常小心的。在很多等式中都出现 e x(Math.exp()
)。例如,它可用于定义 cosh 函数,这已经在 第 1 部分中 讨论:
cosh(x) = (e x + e -x)/2
不过,对于负值的 x,一般是 -4 以下的数字,用于计算 Math.exp()
的算法表现很差,并且容易出现舍入错误。使用另一个算法计算 e x - 1 会更加精确,然后在最终结果上加 1。Math.expm1()
能够实现这个不同的算法(m1
表示 “减 1”)。例如,清单 6 给出的 cosh 函数根据 x
的大小在两个算法之间进行切换:
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
之间的差别很明显的情况下,常常使用 e x,而不是 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 完成了这个计算任务:
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
中提供的方法,最好继续使用。它们通常是最佳的选择。