天天看点

第七章 数学

说明:本文翻译自《TangoRefMan_Sep_1_2008》

由于本人是编程初学者,对很多程序设计概念不是非常熟悉,编程经验不多,再加上英语水平不高,翻译纯属一个D语言爱好者实验之作,很多错误在所难免,还请读者见谅。另外,如果你发现本文有不当和错误之处,还请多提宝贵意见。  

第七章 数学(Doing the Math)

介绍(Introduction)

计算机最初的构想是作为进行数学运算的设备。早期的电脑花费大量时间用于求解方程。尽管现在工程和科学计算只是计算世界的一小部分,然而计算技术的发展留下来一个宝贵的财富:那就是现在几乎所有的电脑硬件具有超高的计算准确性和非常快的计算速度。遗憾的是,大多数编程语言很难让程序员充分利用硬件的这些特性。一个显著的问题是缺乏文档;甚至对许多数学程式员来说,浮点运算仍然是神秘的事物。

在大多数的编程语言中,对数学地支持仅限于单个文件。例如,在C和c++,所有的数学函数都包含在“math.h”中。考虑数学领域的范围,以及数学与计算科学的历史联系,这个有限的存在是几乎是一种侮辱。这种安排也意味著深奥的功能与简单的和广泛应用的功能混杂起来,使得该领域的新人更加感到扑朔迷离。在tango数学库中,基本的“高等学校”数学函数从低级数值分析函数和更少熟悉的特殊函数中被分离出来。

IEEE浮点运算(IEEE floating-point arithmetic)

tango的数学库更多地涉及浮点运算。要理解它,必须对IEEE浮点算法有一个基本的了解。从根本上说,它是一个无穷多的实数在一部分字节上的映射。作为一个32位浮点数,仅有40亿个不同的数字可以被正确表示。正是在这样一个让人沮丧的狭小的表示空间中,IEEE浮点算法做了一个不同寻常的维护幻想的好工作,使数学实数可以正常工作,不过重要的是要明白什么时候打破了幻想。

第一个缺点来自这些可表示数字的分布。IEEE数轴与数学上的数轴迥然不同。

+  +-----------+------------+.. +.. +----------+----------+ +     #

-inf -real.max  -1    -real.min   0   real.min   1   real.max inf  NAN

注意到一半数量IEEE数轴在-1和+ 1之间。当使用80位实数,有超过0.5到1之间一千倍的从0到0.5间的数。(这不是一个纯粹的对数关系,那就是:在1.5和2之间存在与1和1.5之间同样多的可正确表示的数)。这对于精度具有重要意义:有效的精度令人难以相信地接近零。利用这一点在-1和+1间的数看作分开的会让一些例子被覆盖。

同时也要注意这些特殊的数:+ / -无穷大,所谓的“denormals”或"subnormals"在+/-real.min 和 0之间,以减少的精度来表示,最后神秘而令人难以置信的有用的“NaN”(“Not-a-Number”)将在下面讨论。

同时也要注意这些特殊的数:+ / -无穷大,所谓的“denormals”或"subnormals"在+/-real.min 和 0之间,以减少的精度来表示,最后神秘而令人难以置信的有用的“NaN”(“Not-a-Number”)将在下面讨论。

 

 第二个重要概念是需要注意取整误差(round-off error)。由于大多数数学上的数字不能完美表示,只好用最接近的可描述数字来代替。因此,几乎所有浮点运算会遭受一些取整(舍入)误差。不幸的是,取整误差会在漫长的计算中积累。在Tango数学库执行浮点计算的一些函数采用了不同寻常的方式,以最小化取整误差。基于同样的原因,几乎所有函数均只有real的精度。

  

有必要知道哪些运算可以避免取整误差。任何可以存储在int型中的数可以在double中精确表示,任何可以存储在long型中的数可以在real中精确表示。任何与int进行运算时精确的算法采用浮点时也是精确的(如,27.0/3.0 = = 9,精确)。与任意的2的幂进行乘、除也是精确的。

了解这些对于进行相等测试时具有重要意义:不要使用==进行精确的相等测试。该测试通常只能得到近似相等,它带有一些舍入误差的修正值。Tango库中包括一个函数feqrel()可用于该目的的,将在下面讨论。

错误处理

当一个提供的参数超出函数可表示的领域时,函数会怎样做呢?例如一个最简单的情况sqrt(-1)。有如下一些可能性,各有优缺点:

使用一个输入的契约。创建调试版本时一个断言被产生,从而帮助人们快速发现程序错误。不幸的是,这种行为在创建发行版本时是未定义的,除非另一个解决方案被很好的使用。另一个问题是,有效的定义域会变得非常复杂,它会不切实际的期望用户代码避免所有的违反契约情况。

抛出一个异常。这会导致代码膨胀,并强迫用户把代码弄乱放入try..catch语句块中。

返回一个NaN。一个Not-a-Number意味着它不可能提供一种明智的答案;例如,当要求CPU计算0.0/0.0,它会创造一个NaN。NaNs被硬件非常好地支持,并且调试时非常有用,不过会被大多数编程语言忽略。

在Tango.math中,所有的函数在表示无效输入时均返回NaN。此外,契约仅用于在明确知道有一个程序的逻辑错误时使用(例如,对于统计函数“statistical functions”来说,所有可能性必须在0 . .1之间)。不过它会尽可能的做得比这更好一些。IEEE NaN可以包含一个有效负载——一个整数存储在NaN内,该整数可以用来解释NaN来自哪里和它为什么会产生。D也许是第一个使用这个性能的流行编程语言,因此必须制定一个ABI。如果可能的话,在可能的地方,tango库NaN包含一个有效负载解释它在哪里被产生。这个有效载荷被I/O函数显示,也可以被库函数取回。注意,在许多处理器上包含NaN的计算通常都比较慢,因此更可取的是要限制他们使用到异常情况下。

math.Math

在math.Math中,最常见和通用的处理浮点数的数学函数都可以找到。几乎所有的函数都会被学过高等数学的人所熟悉。这些函数可分为三类:乘方、三角、舍入。

(常量)Constants

与同等的数学常量有一些微妙的差异。注意,PI不是数学数字pi;它只是在可表示的范围中最接近PI的那个数。所以,例如,sin(pi)= 0,但罪(PI)= 2e-20。这个取整误差的原因不在正弦函数sine上,而在于PI的表示精度。(实际上,情况更复杂的:在x86平台,sin(x)在abs(x) > pi/2的无论何时都会返回不正确的值)。

乘方、根和绝对值(Powers, roots, and absolute value)

abs: 可以是integer, real, complex,和 imaginary 参数

sqrt: real 或 complex,。 real版本决不会返回一个 complex 结果。

cbrt :立方根(cube root)。 仅有real 型。

exp, log, expm1, log1p.

exp,log是指数和自然对数函数.对于大的正数和负数exp()是精确的,不过当任意参数非常接近0时当然会返回1.0。像先前讨论的一样,IEEE数轴表示接近0的数要好于表示接近1的数。expm1(x)等价于exp(x)- 1.0,但它保存精度几乎降为零。同样的,log1p(x)是一个相当于log(x)+ 1当x的值接近零时的值。在适当的时候,用这两个函数代替exp()和log()是一种提高精度的简单的方法。当然,x是一个任意的大负数x时,expm1(x)恒为- 1。

pow -- 乘方, 高速计算整数的乘方

hypot -- 计算直角三角形斜边的长度,应用毕达哥拉斯规则(Pythagoras' rule),但用了一个避免溢出的方法。

poly -- 计算矩阵多项式(polynomial )。

conj -- 共轭复数的计算 。

max, min, maxNum, minNum

最大值和最小值是出乎意料的复杂。如果一个参数是一个NaN将会发生什么情况呢?根据754R标准,如果一个参数是NaN而另一个是正常数,该数应当被返回(对于检测一个函数的范围来说,该行为是有用的)。因此,这些函数被命名为minNum 和maxNum,以提醒不论何时应尽可能地返回一个数。遗憾的是,它的速度要比朴素的执行要慢一点。max and min 是 maximum和 minimum 函数工作在任意类型 时的版本;若任何参数是一个NaN,则返回结果是未定义的。

三角函数

对于实数(real)、复数( complex)和虚数(imaginary )参数,全部用弧度(而不是度)。

cos, sin, tan

acos, asin, atan

cosh, sinh, tanh

acosh, asinh, atanh

即使对于那些不熟悉的复数的人,也值得知道下面的联立的(同时的)快速计算正弦和余弦的方法。

creal z = expi(real x);

real c = z.re;    // = cos(x)

real s = z.im;    // = sin(x)

在大多数处理器(包括x86),这是一种固有的操作,速度是分别计算正弦和余弦的两倍。

注意就像 exp(x), cos(x) ,当x很小时通常返回1。cos(x) - 1形式的表达式将会用 cos(x) - 1 = -2 * sin(x/2)2 形的结果重写。

sinPi, cosPi, atanPi

要精确地利用度代替弧度进行三角函数运算,可以使用sinPi(x / 180)。自从360可以被浮点数正确地表示,但pi不能,这些函数仍然对大的转数保持准确。(例如,sinPi(1e30)正确地返回0,但sin(1e30 * PI)的返回值是无意义的。

函数sec = 1/cos, cosec = 1/sin、cot = 1/tan,它们曾经在三角学中是重要的,但在对数被介绍之后就变得不那么重要了(http://mathworld.wolfram.com/Secant.html)。他们并不包括在tango.math中,因为它们用手工计算是微不足道的,还因为sec()这个名字和单词seconds(秒)非常相似,会造成不必要的混乱。

atan2(y, x)

这是用来转换直角坐标(x,y)成极坐标(angle, radius(角,半径)),这对于复数运算是基本操作。(它也是适合用在零的符号特别重要的场合上(本句不够准确))。

舍入到整数(Rounding to integers)

将一个实数转换到最接近的整数是非常普遍的需要。有一个小小的困难:当存在一个难以取舍的平局情况怎么办呢?3.5应当舍入到3还是4 ?为了兼容C语言,cast(int)通常向上舍入。不幸的是,在许多处理器上这个操作出乎意料的慢,因为它需要CPU的的舍入的方式被改变,随后又重新恢复。(它是如此缓慢以致于Pentium 4专门包括一个特殊的优化器以缓解速度的惩罚,SSE4指令集进行进一步的优化)。然而在大多数应用程序中,采用何种舍入方式不是很重要的。函数rndint和rndlong是cast(int)和 cast(long)的高速替代者。它们的决定权(tie-breaking behaviour)取决于舍入模式,稍后将进行讨论。

还提供了一些另外有用的从实数到整数转换的模式

floor

ceil

round

trunc

值得注意的是这些函数返回实数(real),而不是整数( integers)。要储存它们为 integral类型,仍然还要使用rndint或rndlong。

math.IEEE

math.IEEE是低级函数的家,它与计算机硬件的密切关系要超过与数学的关系。要让这些函数变得有用,需要对IEEE算法有一定的了解,不过传统上它们通常被给予非常晦涩的名字,在数学库之外很少被使用。人们希望他们会更广泛用于D,因为他们既能大大提高计算精度和健壮性,他还可让一些有趣的技巧成为可能。

检查近似的相等(Checking for approximate equality)

操作符==在大多数相等测试中通常只进行精确地相等比较,因为哪怕是一个微小的取整误差就会导致测试失败。相反的,一些足够接近(close enough)形式的测试必须被使用。不幸的是,很容易就能让写出的这个测试不正确。例如,很明显abs(x - y)<=0.000001是不正确的,比如,它意味着1e-10近似于7.8e-300,但这并不意味着1.000000000000000000000001e60近似等于1e60。取而代之的是,我们需要一个在IEEE数轴附近的测试而不是“在数学实数数轴附近”的测试。

分析师通常所说的对待取整错误的态度—"units in the last place" (ulp)(把误差留在最后的地方),但这仅仅是误差很小的时候才有用。一个更直观的措施来统计舍入错误的位数(或做另一个选择,保留准确的位数)。虽然这听起来有点复杂,它可以充分地利用二进制表示浮点数的优势来有效地执行。在编译时指定允许的最大舍入位数,可以使测试速度非常快——只是稍稍慢于正常的==比较。为确定清爽的舍入误差是多大,使用tango.math.IEEE中的准确函数feqrel() 。

feqrel ( D特有)

Returns the number of significand bits which are equal in both x and y. Formally, real.mant_dig - ceil(log2) of the number of times nextup() would have to be called to move from the smaller number to the larger one.

IEEE 异常状态标志(The IEEE Exception Status Flags)

  

所有IEEE-compiliant处理器包含特殊状态位,用来指示程序可能想了解的当“异常(weird)”的事情发生时状态。例如,ieeeFlags.divideByZero(IEEE标志.被零除)告诉我们一个数被零除时一个无穷大发生了。在计算的末尾仅检查一次,就可以避免相比成千上万次的失败比较。

  

这是一个可以探测到的异常状态的列表:

无效(invalid):设置为当任意NaN形成时。这个状态可能发生于∞-∞,∞×0,∞×0,0÷0,∞÷∞,∞%∞,x%0(x是任意数),以上这些情况几乎总是预示着一种编程错误。

被零除(divisionByZero): 设置为被零除正负无穷大发生时,这通常预示着一个编程错误,但不是经常发生; 某些类型的计算当被零除情况发生时返回正确的结果。 (如, 1/(1+ 1/x) == 0当 x == 0时)。

溢出(overflow):设置为把两个非常大的数字相加或相乘,它们的结果大于real.max时。

下溢(underflow):当把两个非常微小的数相减或相除时,它们的结果接近零时就产生下溢。下溢几乎从不产生问题。

不精确(inexact):几乎所有浮点运算操作都设置这个标记!对精确性问题,已经进行了讨论。

IEEE取整模式(The IEEE Rounding Modes)

取整被函数setIeeeRounding和getIeeeRounding控制 。可以设置四种取整模式。默认模式舍入到最接近的数,在统计学上是最准确的,但缺少直观。在平局的情况下,舍入到偶数。

取整模式             rndint(4.5)             rndint(5.5)               rndint(-4.5)   

舍入到最接近的数  4                              6                          -4                 平局情况下舍入到偶数

向下舍入               4                             5                          -5   

向上舍入                5                             6                         -4   

向0方向舍入           4                             5                         -4   

通常有两种类型的原因要改变取整模式。一是简单的检查数值稳定性:当取整模式被改变时如果计算结果产生难以控制的不同结果,它(译者按:它指改变取整模式)可以作为一个清除符号可以免受取整误差带来的麻烦。第二个原因是为了实现区间四则运算。但最明显的用途是暂时控制rndint的意义,这样它才能准确复制在一个内循环中的cast(int)的工作情况。

操作指数和有效数字的函数(Functions for manipulating exponent and Significant)

由于IEEE采用了基于2的指数,用2乘或除一个数是一种精确的操作(完全不会产生任何舍入错误),当然,这个结果也并不会产生上溢或下溢。利用这一事实,计算的精度有时可以得到提高,或者延展其使用范围。函数ldexp()和frexp()在计算结束时,允许从一个数或它们的恢复中去除所有的因素。

frexp -- 获取指数和有效位数。

ldexp -- N * 2y

ilogb, logb

scalbn -- x * 2n

nextUp, nextDown, nextafter

nextUp(x)和nextDown(x)返回 IEEE数轴上相邻的点。他们已经加入了754R IEEE标准草案。nextafter()的使用是沮丧的,许多IEEE标准委员会成员现在认为这是个错误。

ieeeMean ( D特有)

我们在前面看到从IEEE数轴上存在的0和1间的数是1和2间存在的数的几百或几千倍。试图用二分法(binary chop)获得对数性能的分而克服(divide-and-conquer)算法是有问题的。考虑查找f(x) = 0中的x的值的问题,使用平分法把区间一分为二。如众所周知的0 < x < 2,一个朴素的二等分会把该区间分成[0,1]和[1,2]。糟糕的是这是错误的,因为在区间[0,1]上包含了原区间上超过99%的可表示数字! ieeeMean(x,y)给出了表示范围的中点,在所有的情况下提供良好的性能。注意:这个函数是没有被IEEE标准描述。

signbit

除了要区别正负零外,与x > 0相同。

copysign

可以提高精度和速度的辅助函数(Helper functions that can improve accuracy or speed)

fabs --内在的。被 math.Math.abs使用。

  fma --代替“融合乘加”("fused multiply add")。(x * y)+ z。有些cpu(但不是x87)可以在一条指令中执行,这增加了准确性,因为它仅包含一个舍入(取整)误差。

  expi—sin + cos.。在x86机器上这是一条指令。被exp(ireal)使用。

处理特殊实数类型的函数(Functions for dealing with special types of reals)

isNormal

isSubnormal

isInfinity

isNaN

  这些函数存在的主要为了兼容c。注意isNaN(x)等价于x!<>=0,除了NaNs信号产生时不抛出异常外,isInfinity(x)与abs(x)= = x.infinity一样。

  isIdentical(D特有,判断同一性)

  除了正负零被看作不同外,相同作为相等。如果他们有相同的有效载荷,NaNs被视为相同。

NaN有效载荷(NaN Payloads)

根据IEEE 754标准,有效载荷(payload)可以存储在NaN的尾部。该有效载荷会包含 NaN怎样和为什么产生的信息。从历史上看,几乎没有任何编程语言使用过这个潜在的强大功能。在tango中,该有效载荷由一个正整数组成。有一个困难是每当转换到一个更小类型的情况发生时,某些位会从有效载荷中丢失。IEEE标准没有明确说明此时发生了什么事。

NaN — 带“有效载荷 ”创造一个NaN,该“有效载荷”是一个ulong型数字。

ulong getNaNPayload(real x)——返回整数载荷。注:如果缩小转换已经发生,高序位可能已经发生了变化。决不要存储一个指针作为一个整数载荷在一个NaN中。垃圾收集器将无法找到它!

 math.BigInt

 

对于可靠的应用程序来说,内建的整数类型精度不足,如某些应用高精度的数学计算和密码程序等。tango使用BigInt类型提供任意精度的算法,BigInt能够执行数千位精度的整数计算。一般来说,BigInt可以作为int和long型的插入代替者。要做到这一点,BigInt使用数值语义理论为指导,在D1.0 中用写时复制(Copy-on-Write)实现。不幸的是,这个结果会产生很多临时堆对象(在D2.0中不会发生这个问题)。

BigInt旨在提供优良的性能,用小一些的数字来缓和数位(胜任几千个十进制数位)。高度优化汇编程序是用于最广泛应用的处理器。对极其巨大的数据,或较为不通用的处理器,可以考虑用GMP(http://www.gmplib.org/)。

math.Bracket(支架、括号)

几个重要的数值算法都包含在Tango中。数学。Math.Bracket包含查找根和求单一参数实数函数的极值。他们被设计用于较好地解决高等非线性函数。

FindRoot

经典的一元根查找问题就是:给定a,b和一个实数函数f(real),计算f(x)==0(a<=x<=b)中的x。如果您可以提供一个范围[a..b],这样f(a) 和 f(b)有相反的符号,则要求的根可以在有效的范围之内找到。(“有效”是指以尽可能少的调用f())。

Math.Bracket.findRoot使用一个原本基于TOMS478的算法(Alefeld et. al.)。Math.Bracket.findRoot使用一个原本基于TOMS478的算法(Alefeld最初et.。),尽可能用三次插值方法,在必要的时候退化到较慢收敛方法。显著的速度和稳健性增强了。TOMS478与几乎所有其他公众可用算法遭受极其差的糟糕性能(“行业标准”算法由于P.R.布伦特可能要求为80位实约 40000次调用f());该实现避免了这样的问题,使它比平均速度快1.5倍,在比最糟的情况快200倍。典型的10 - 15次调用f()需求达到全部80位精度。

例如,计算x5 = 0.2中的x,已知x在0到5之间,使用:

real root = findRoot( (real x) { return pow(x, 5) - 0.2; }, 0.0, 5.0);

如果该算法的支架(bracket)没有包括零会工作得更快。通过再次观察IEEE数轴,你可以看到,如果你可以指定一个范围1e-5 ..5.0,你将给该算法比起你给一个包括0的范围少做很多的工作。

注意不是所有的查找根据的问题都能够用这种方法解决;要提供 f(x)上相反符号的两点可能会变得复杂(有时f(x) 仅接触x轴)。

FindMinimum

 

在一个维度的最小值(Minimisation)普遍使用布伦特(Brent)的算法,该算法结合了抛物插值法和平分法。

  在一个维度的最小值(Minimisation)普遍使用布伦特(Brent)的算法,该算法结合了抛物插值法和平分法。如 J. Crenshaw (http://www.embedded.com/story/OEG20010628S0046)说的,布伦特的算法并不是最佳的,将来tango的发行版中可能会有一个结构来取代它。这个函数仅仅找出机器精度的一半的最小值,因为取整误差会让更高的精度毫无意义(虽然不是经常,例如abs(x)有一个非常精确的确定最小值)。

  

如果你需要找到一个函数f(x)的最大值,只要找到-f(x)的最小值。

注意,提供两个点a和b给minimum做为托架,该算法也需要一个第三点c在[a..b]之间并满足f(c)<= f(a) 和 f(c)<=f(b)。依赖于你的函数,查找三个这样的点会极其复杂。

math.Probability(可能性、概率)

统计分布函数被广泛地应用于数理统计理论,经济学,生物学,以及整个实验科学与工程学。Tango提供高精度最重要的概率分布工具。我们之前已注意到,介于0到0.5之间存在着0.5到1之间的数千倍的可表示的数。当表示极端的可能性时这可能会是个问题,因为通常存在同样多的几乎必然发生的事和几乎一定不发生的事。要保留精度,函数被提供使用补助概率q,定义为q = 1 - p;对于接近1的概率使用蒙特卡罗版本(Compl)的分布函数。(这完全类似exp() 和expm1()之间的关系)。

在Tango中包含下面的累积分布函数。在每种情况下,补充(complement,后缀Compl),相反的(inverse,逆,后缀Inv)和补充和逆(complement inverse,后缀ComplInv)被提供。

MATHEMATICAL  TANGO NAME

 normalDistribution

Student's t  studentsTDistribution

f  fDistribution

chi2 chiSqrDistribution

 poissonDistribution

 binomialDistribution

 negativeBinomialDistribution

 betaDistribution

 gammaDistribution

  逆函数是用来统计数据构造置信区间(confidence intervals,也可以是可靠区间)。

  

  几乎所有的这些函数的实现都是通过伽玛(gamma)和β( beta )特殊函数,在下面将进行讨论

  数学的特殊函数(Mathematical Special Functions)

  数学术语的特殊函数(Special Functions)涵盖大范围的超越一般常识的函数。这是罕见的同时用在一个单独的应用程序中查找倍数特殊函数(find multiple special functions)。

math.ErrorFunction(错误函数)

错误函数和正态分布函数的从根本上来说是相同的函数,从比例因子中分离出来。

erf, erfc

错误函数和互补错误函数。它们出现在最常见正态分布的集成中,也被称为“钟形曲线”(bell curve)或“高斯分布”(Gaussian distribution,尽管高斯并没有发明它!)。错误函数和互补错误函数。它们出现在最常见正态分布的集成中,也被称为“钟形曲线”(bell curve)或“高斯分布”(Gaussian distribution,尽管高斯并没有发明它!)。erf()接近0时有相对高的精度,而erfc对大正数有相对高的精度。(这类似于前面提到的expm1() 和 exp()的关系)。利用下面的结果,它给任意x尽可能获得相对高的精度:

erf(x) = 1 - erfc(x)

erf(x) = -erf(-x)

它们与“Q 函数”密切相关,广泛应用于通信工程学和正态分布集成。

// 累计分布(Cumulative distribution)

real normalDistributionFunction(real x) {

       return - 0.5 * erfc(-x * SQRT1_2);

}

real normalDistributionFunctionCompl(real x) {

       return -normalDistributionFunction(-x);

}

alias normalDistributionFunctionCompl Qfunction;

相反的累计正态分布函数,同时提供的还有被熟知的分位点函数(quantile function)。例如:如果考试分数依比例达60的平均分和12.5的标准方差,有学生分数会达到85或更多吗?

import tango.io.Stdout;

import tango.math.ErrorFunction;

void main()

{

    // calculate how many standard deviations

    // we are from the mean.

    real x = (85-60)/12.5;

    Stdout.println(1-normalDistribution(x));

    // prints 0.0227501319481792, which tells us that 2.3% of students

    // will recieve a high distinction.

    // Now for the inverse calculation.

    // What mark is required to be in the top 1% of students?

    Stdout.println(60 + 12.5*normalDistributionInv(0.99));

    // prints: 89.0793484255105

math.GammaFunction

gamma, logGamma, sgnGamma

伽玛函数是阶乘函数包括分数和负数的一种归纳。在数理统计中这是极为重要的。 loggamma(x)与log(abs(gamma(x)))完全相等,但它有效地保留了更大的x的值。在C语言中,由于一些历史原因把它们叫做tgamma和lgamma。sgnGamma(x)给出符号(x>0时gamma 通常是正的,但是符号有规律地变换为负x)。

beta, betaIncomplete, betaIncompleteInv

beta函数与伽玛函数有密切关系。正如伽玛函数gamma(x)类似x的阶乘,beta(a,b)函数类似于二项式系数(binomial coefficient (a choose b)),但它对real型和整型(integers)数同样有效。

不完全的beta函数betaIncomplete(a、b、x),是beta函数的一个部分,是最重要、广泛应用于数理统计概率的计算。绝大多数基本的数理统计测试(如,学生的t,菲舍尔F测试( Fischer F tests))是通过不完全beta函数和它的逆运算(its inverse)来实现的。

gammaIncomplete, gammaIncompleteInv

不完全伽玛函数(incomplete gamma function)在一些统计分布函数中用作低能级函数(low-level function ),其中最引人注目的是Chi-squared累积分布函数。

math.Bessel(贝塞耳函数)

贝塞尔函数是用来解决贝塞尔微分方程的:

x2 d2y/dx2 + x dy/dx + (x2 - a2)y = 0

这类方程最常见于圆柱或球形空间中的电磁波研究(studies of waves)(下一个管道,比如);因此,他们在电磁波的传播学的计算中是很重要的,例如。有几种类型的贝塞尔函数对应圆柱和球面坐标,不管x = 0时的值是finite ("first kind") 还是 infinite ("second kind")。传统的数学中的这些函数的名字相当晦涩(J0 J1,Jn,Y0,Y1,Yn、)。在Tango中,他们被给予了轻松的更具描述性的名字[参考:WalterBrown]。

Value of a  First kind  Second kind

zero  cylBessel_J0  cylBessel_Y0

one  cylBessel_J1  cylBessel_Y1

integer n cylBessel_Jn  cylBessel_Yn

注:其它的贝塞尔( bessel)函数的变化不在目前的Tango实现中,但在将来可能添加到Tango的这个模块里。被修改的贝塞尔函数被重命名为cylBessel_I;第三种被修改的Bessel函数将被命名为cylBessel_K,球形的贝塞尔函数会被命名为sphBessel_j 和sphBessel_y

math.Elliptic(椭圆积分)

最初的椭圆积分出现在计算一个椭圆弧形长时,但通常适用于整个一类积分方程。任何此类方程可以转换成这三种椭圆积分的结合。Tango的命名约定是基于数学应用软件的。不完整的椭圆积分依赖振幅、φ和系数m 。完整的椭圆积分是φ= PI / 2时的特殊情形。

 Incomplete (不完整) Complete (完整)

First (第一) ellipticF  ellipticKComplete

Second (第二) ellipticE  ellipticEComplete

Third(第三)  ellipticPi*  ellipticPiComplete*

当前的实现受限于较低的精度。

在回顾

  对于准确性和快速代码的一些一般性的建议:

  尽可能地利用零附近的相对高的精度。

  不要过分担心除零错误。通常情况下,计算将依然是正确的。异常状态可以用来检测反常行为。

  正确使用NaN会非常有利于调试。

  你想要进行精确相等的测试使用feq忙天文超过= =,除非你想测试确切的平等。

  使用rndint和rndlong比用cast(int) 和 cast(long)好,除非平局决胜行为是很重要的。

  

附录

 

 附录1:Tango NaN结构

  “有效载荷”可以被储存在NaN的尾部。一个位(最有意义的)必须区分NaN是安静的还是在发信号。这就留下float的22位有效载荷,double的51位;80位real的62位;128位quad的111位。不幸的是,每当转换到一个较小的类型的情况发生,位就会从有效载荷中丢失。在在Tango的NaN中, float的22位尾数是有效载荷整数的低22位。22..51位储存在double的尾部,(80位real)位52 ..62储存在real的尾部。

附录2:数值算法目前还没有包含在Tango中

值得注意的是,经典文本“数值分析方法库”("Numerical Recipes",http://www.nrbook.com/a/),虽然是对数值算法进行广泛的概述的优秀资源,但它某些部分包含一些错误的陈述,以及相关的源代码有严重的著作权限制。

快速傅里叶变换利用FFTW库可以有效执行:http://www.fftw.org

 脚注:可移植性和80位实数的使用

为了可移植性,80位实数有时被认为应该牺牲。如果代码是一个目标,轻便,认识到这一点很重要,尽管所有系统都能够存储数据和双精度的浮动,并非所有的都能执行计算那些精度。特别的,x87浮点单元,现在几乎所有的英特尔与AMD的机器的真实硬件支持80-bit浮点计算,在计算精度较低时包括一些仿真。奔腾III是英特尔的第一个原生支持双精度浮点计算的处理器(需要使用SSE2指令集)。Tango库尽可能使用用80位real型数,必要时优雅地退化到使用64位实数。

参考和延伸阅读(References and Further Reading)

"What Every Computer Scientist Should Know About Floating-Point Arithmetic",

《每一个计算机科学家都应该知道什么是浮点算法》

http://docs.sun.com/source/806-3568/ncg_goldberg.html

"A Proposal to Add Mathematical Special Functions to the C++ Standard Library", Walter E. Brown, 《建议增加数学特殊函数到c++标准库》,沃尔特-布朗,

http://std.dkuug.dk/JTC1/SC22/WG21/docs/papers/2003/n1542.pdf. (Note that this proposal has been accepted into the C++0x TR1 standard,注意这个提案已经加入了c++ 0x TR1 标准).

"A Proposal to Add Mathematical Functions for Statistics to the C++ Standard Library", Paul A. Bristow,(《建议增加数理统计函数到c++标准库》, Paul A. Bristow)

http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1069.pdf

"An Interview with the Old Man of Floating-Point: Reminiscences elicited from William Kahan by Charles Severance",

http://www.cs.berkeley.edu/~wkahan/ieee754status/754story.html

继续阅读