精确几何计算的非代数方法与超几何函数的实现

53- × ÷《理论计算机科学电子札记》66卷第1期（2002年）网址：http://www.elsevier.nl/locate/entcs/volume66.html12页超几何函数精确几何计算ZilinDua，1，MariaEleftheriou b，Jos'eE. Moreirab，CheeYapa，1a 计算机科学系，柯朗研究所，纽约大学，251默瑟街，纽约，NY 10012bIBM T.J.沃森研究中心，邮政编码地址：Box 218，Yorktown Heights，NY10598摘要计算几何中的大多数问题都是代数问题。精确几何计算（EGC）是解决此类问题的非鲁棒性的一般方法。现在有一些通用的库为通用程序员支持EGC（例如，核心库，LEDA Real）。许多应用程序也需要非代数函数。在本文中，我们描述了如何在其他EGC功能的背景下提供非代数函数我们实现了一个多精度超几何级数包，它可以用来评估常见的初等数学函数的任意精度。使用CoreLibrary可以相对容易地实现这一点，它支持有保证的精度水平。我们解决了几个问题的效率在这样一个超几何包：自动错误分析，参数减少，超几何参数的预处理，和预先计算的常数。本文报道了一些初步的实验结果1介绍数值计算有两种基本模式：一种是假定固定精度，另一种是要求任意高精度。数值计算中的第二个二分法是基于将计算分类为代数或超越代数的计算根据定义，代数计算只能使用代数函数作为原语（例如，+，代数函数，如sin或log。代数问题是可以使用代数计算解决的问题[22]。在计算几何[7，1]中出现的大多数问题以及几何建模的很大一部分都是代数问题。我们可以把非代数计算的算法问题1由NSF/ITR资助#CCR-0082056支持。本文在西班牙马拉加举行的2002年共同国家评估会议上发表（2002年6月12日至13日）。2002年由ElsevierScienceB出版。 诉 操作访问根据C CB Y-NC-N D许可证进行。DU，ELEFTHERIOU，MOREIRA和YAP54使用可计算分析方法[19]。相反，我们试图generalize更强大的技术可用于数值代数计算非代数设置。本文讨论的数值计算是非代数的，这需要任意高的精度。 更确切地说，我们希望以保证精度的方式计算超几何函数。代数设置的保证精度首先在[20]中提出，并体现在核心库精度API [10]中。问题是，当我们将非代数函数纳入我们的系统时，是否可以做出类似的保证是一个主要的开放问题。保证精度的概念源于几何算法[18，21]中普遍存在的数值非鲁棒性问题的一般解决方案。 这个众所周知的问题困扰着许多科学和工程计算。原则上，可以使用我们称之为精确几何计算（EGC）的方法[22]来消除大类问题中的此类问题基本上，EGC相当于保证准确性的计算。重要的是要认识到，保证的准确性不需要像点过滤器[2，6，12]、构造性根边界[11，5，4，12]和低次谓词[13]这样的技术都是由EGC计算模式驱动的技术目前有两种通用的数值库支持EGC模式：LEDA Real[3]和CoreLibrary[10]。许多数值计算需要代数和非代数计算。目前，核心库和LEDA Real不支持非代数函数. 本文将为填补这一空白迈出第一步 本文描述了一个超几何软件包的设计与实现. 最常见的非代数函数，如exp（x），log（x），erf（x）和三角函数是超几何的。事实上，标准库（如math. h）中的所有非代数数学函数都是超几何的。Jeandel [9]描述了一个最近的尝试，以提供hypergeomet，Ric函数在通用多精度数字包中（在本例中为gmp）。虽然多精度算法是EGC的先决条件，但它仍然缺乏关键的EGC精确比较能力。EGC是否可能用于非代数计算？核心库的设计旨在使任何程序员都能轻松地构建健壮的程序。 为此，我们在[20]中定义了一个自然而简单的数值精度API，具有四个精度级别：• 一级：机器精度（即，IEEE 754标准）• 等级II：任意精度（例如，计算到1000位）• III级：保证准确性（例如，保证100位）• IV级：混合精度（即，3个先前级别的组合目标是允许单个程序在这些级别中的任何一个中运行，只需调用库即可。第二层没有保证：计算到1000位不会DU，ELEFTHERIOU，MOREIRA和YAP55∈∈- ×∈k=0保证1000比特的任何部分由于误差传播而正确第二级和第三级之间存在一个可能不明显的根本差距：第三级不仅仅是简单地以越来越高的精度迭代第二级计算。虽然我们知道如何为所有代数计算提供Level III能力 [21]，但在非代数情况下Level III是否可能是一个悬而未决的问题我们称之为EGC的基本问题设f是一组实数或复数函数或常数。在实际中，它的运算符中包含+，nZE上的常数的集合记为E（E）。表达式e E（n）的值是一个实数（或复数），但它也可以是未定义的。常零问题，记作CZP（），是决定给定e E（），e的值是否被定义并等于0。当矩阵包含至少一个非代数函数（如sinx或logx）时，矩阵的常零问题是不可判定的，并且与不可判定的问题密切相关[14]。然而，这正是我们所需要的保证准确性。因此EGC的基本问题实际上是一个问题族，每个问题都有CZP（n）本文的贡献。• 我们描述了一个多精度超几何库的设计和实现。 这很容易实现，具有核心库的III级精度。• 本文介绍了超几何函数有效求值的超几何参数处理问题• 我们介绍了超几何函数的自动误差分析技术。• 我们解决了减少论点的问题• 该库的一个自然应用是计算任意精度的数学这既是超几何包的应用，也是对包的应用我们还定义了存储和访问这些常量的格式。2超几何级数一个超几何级数其中t= 1和2的比值连续项是求和指数k的有理函数（一）tk+1tkP（k）=x，Q（k）其中P（k）和Q（k）是k中的monic多项式，x是称为自变量的常数。通过分解多项式P（k）和Q（k），我们可以写为：（二）tk+1 =P（k）x =（a1+ k）（a2+ k）···（ap+ k） x。Q（k）（b1+k）（b2+k）·· ·（bq+k）（k+1）DU，ELEFTHERIOU，MOREIRA和YAP56≤||12Q K12Q K升 阶 乘 或 Pochhammer 符 号 （ a ） k 由 （ a ） k=a （ a+ 1 ）（a+2）···（a+k−1）给出，其中k≥1且（a）0= 1。使用这个符号，项tk的一般表达式变为（a1）k（a2）k···（ap）kxk（三）tk=（b）第（1）款（b）第（1）款.···（b） k！因此，超几何级数完全由序列a =（a1，a2， . ） 定 义 。 ，ap）andb=（b1，b2， . ，bq）。Bi's可能不会被零或零整数分割（否则，我们将有一个由0分割的整数）。另一方面，一方面，把k≥0tk变成一个有限级数，因此是x中的多项式。注意传统的K！在（3）中的tk的分母中，隐式Lower参数r0fb0=1。给出了与这些参数对应的超几何级数∞（a1）k（a2）k·· ·（ap）kxk（四）pFq（a1，a2，. 、ap;b1、b1、b2、. ，bq;x）=（b）第（1）款（b）第（1）款、···（b） k！上述级数收敛于任何复杂的x时，p q，并为x1时，p=q+1。相应的复变函数是超几何函数。 在本文中，我们只对x以及s参数sai，bj是实数的情况感兴趣。基本功能。 最常见的基本函数[16]是超几何函数例如，exp（x）的通常级数是∞xkk=0k！初始项t0= 1。两个连续任期之间的比例是tk+1/tk=x/（k+1）。从等式（2），我们得出p=q=0，因此exp（x）=0。表1列出了一些初等函数的超几何级数表示 对于每个函数，我们列出（i）该函数的常用幂级数表示，（ii）该幂级数的两个召唤项之间的比率k + 1/tk，以及（iii）相应的超几何级数。在某些情况下，幂级数的第一项t0不是1，必须被分解。表1某些初等函数的超几何级数表示。初等函数幂级数比率tk+1/tk超几何级数exp（x）∞xkk=0k！1 X（k+1）0F0（; ;x）.Σ电流变率2xπk =0（2 k +1）k！X-2（k+1）x32（k+2）（k+1）102x1F1（1;3;−x2）π2 2正弦（x）∞x2kxk=0（−1）k（2k+1）！121（−4x）（k+3）（k+1）223−xx·0F1（; 2;四、cos（x）∞x2kk=0（−1）k（2k）！121-4x（k+1）（k+1）220F1（;1;−x）2 4反正弦（x）x ∞（2k）！x2 kk=022k（2k +1）（k！）22（k+1）2x32（k+2）（k+1）x·2F1（1，1;3;x2）22 2KKk=0KKDU，ELEFTHERIOU，MOREIRA和YAP57arctan（x）x∞（−1）kx2kk=02k+1-2（k+1）X2（k+3）2x·2F1（1，1;3;−x2）2 2log（1+x）∞（−x）kxk=0k+1（k+1）（k+1）（−x）（k+2）（k+1）x·2F1（1， 1; 2;−x）表1中log（1+x）的标准级数收敛性较差通过从log（1 −x）中减去log（1+x）的标准级数，我们DU，ELEFTHERIOU，MOREIRA和YAP58∞√∞i=0时||≤|≤|| ≤|≤ ≥|≤ ≥||获得日志1−x1+x=− 2x2x41++3 5Σ+···=−2x·2F1（1，1; 3; x2）。2 2改变变量，我们得到.- 是的.Σ2Σlogy =− 21−y1+y2F1131， ;;2 21−y，1+y0 a且n≥2p（a≠）2。p（b−a）p(ii) 当p ≤ q时，n ≥ max {2，2a}.DU，ELEFTHERIOU，MOREIRA和YAP59|| ≤n| | ≥Q我(B) 几何案例当序列不交替时，我们使用以下观察：引理3.2存在单调递减函数g（n）使得f（n）≤g（n）和limn→∞g（n）= limn→∞f（n）.证据观察到对于任何a≤0，b≤ 0，如果a≤b，则（n+a）/（n+b）≤1，如果a > b，n > m，则（n+a）/（n+b）（m+a）/（m+b） bi − 1。 ，r，并且对于ri=r+1， . . . ，p.Thenpf（n）=1 （n+ai）=n+ai1Qj=0（n+bi）i=1 n+bi−1j=p n+birn+aq≤n+b1n+B =g（n）.i=1i−1j=pi所以g（n）是单调递减的。注意limn→∞g（n）= limn→∞f（n）（这个极限是0或1，取决于p q+ 1与否）。Q.E.D.引理3.3对于引理3.2中给出的g（n），如果xg（n）<1，则Rtn。1−xg（n）证据注意，对于i≥0，i−1i−1tn+i=tnxij=0f（n+j）≤tnxij=0g（n + j）≤ t n x i g（n）i.总结，Rn= Σi≥0n+i≤nΣi≥0x i g（n）i= t n/（1-xg（n））。Q.E.D.对于g（n）的一个更简单的选择是g（n）=np−q−1。实际上，在选择尽可能小的g（n）和计算g（n）的计算量之间存在权衡。实施. 上述界限直接应用于表1中的基本函数。下面是一个动态计算任何超几何级数的算法：(1) 首先，我们检查p> q +1或（p = q +1和x1）。如果是这样，我们报告一个(2) 否则，我们检查是否x0。如果是这样，我们可以将情形（A）应用于交错级数。(3) 否则，我们计算g（n）并将情形（B）应用于几何情形。一般来说，为了确定n，使得|R n| ≤ε/2，我们需要上界t n的值。一个简单的方法是累积Sn中的项，同时检查下一项tn。从表面上看，计算tn需要迭代，因为tn中的因子数随着n而增长。但在许多情况下，tn中的因子数与n无关，因为pDU，ELEFTHERIOU，MOREIRA和YAP60≤∈| − |≤≤取消（例如，以log（1+x）为单位。）在这种情况下，我们可以通过直接求值在常数时间内确定n这产生了更快的实现。这种情况在参数处理中被利用（见下文）。4变元缩减超几何函数的有效求值中的一个问题是著名的变元缩减问题每一个超几何级数在一个有界的范围内是一般有效的，问题是把一般的参数降到这个范围内。即使当参数在有效范围内时，仍然可以应用参数减少以实现更快的收敛。正如在[16，p.145-147]中所指出的每当我们执行参数归约时，都会在修改后的参数中引入一个错误。我们需要约束这个错误的后果例如，三角函数的自变量约简使用了它们具有2π周期的事实。通过利用其他性质，参数可以减少到大小为π/2的范围。如果r是对应于x的原始自变量的约化自变量，我们有、、、π2r= x − 2πx。但我们只能计算出r的近似值r。使用一个非常精确的π，我们可以将rr=y，并且是期望的误差界εJ。标准函数的εJ的选择由以下引理导出。引理4.1对于ε>0，我们有如下的界：|≤ε| ≤ ε|≤ε| ≤ ε|≤4ε，0≤x≤π/4，επ/12| ≤ 4 ε,0 ≤ x ≤ π/4,ε <π/12|≤2ε，π/4≤x≤π/2，ε≤π/4| ≤ 2 ε,π/4 ≤ x ≤ π/2,ε ≤ π/4|≤2ε，|X|0的情况。|<0. 5，ε≤ 1/ 4|≤2ε，|X|0的情况。|<0. 5，ε≤ 1/ 4|≤ε，|X| 1|<1|≤ε/x，x>0| ≤ ε/x,x > 0|≤2εexp（x），ε≤log（2）|≤ 2 ε exp(x),ε ≤ log(2)证明使用泰勒展开式的余数形式，f（x+h）=f（x）+hfJ（x + θ），0θ h。自然对数函数。如果x > 2，我们设x = 2kr，其中kN和一<对二。则log（x）=k log（2）+ log（r）。 下面是近似这个表达式的步骤：1. 首先计算k =[log 2xk]。DU，ELEFTHERIOU，MOREIRA和YAP61˜˜≤[−|≤|≤K2. 计算log（2）作为log（2）的近似值，绝对误差≤ ε/（2k）。3. Computerérésuc hthat |r−r|≤ε/4，其中r=x2−k。4. 计算log（r）近似log（r）的绝对误差≤ε/4。5. 返回z=klog（2）+log（r）.使用引理4.1（除其他外），我们可以证明这个过程的正确性。也就是说，zlog（x）ε。此外，每个步骤都可以使用Core Library轻松计算。第二步需要对常数log（2）进行近似，这是我们预先计算的（见第6节）。指数函数。设k = x/log（2），r = xk log（2）。则exp（x）= 2kexp（r）。注意1r 2。1. 首先，我们计算k（需要对log 2进行适当的近似2. 计算ε2−k−2e−2。r=x−klog2的近似值，绝对误差˜3. 计算e ex p（r）作为 ex p（r）对绝对误差ε2的近似值-k-14. 返回z˜=2ex p（r）。三角函数 计算arcsin（x）当0。5 1、我们使用. Σπ1arctan（x）= 2 − arctan x。同样，我们需要计算1/x到ε/2的绝对误差界sin， cos， tan， cot的情况5超几何参数预处理在我们的实现中，预处理是通过在第一次调用超几何函数之前为它们引入“注册设施”来实现的预处理包括计算截断误差的边界函数（例如，（3）函数g（n）但是我们现在讨论预处理的另一个为了得到这些级数的标准形式，有时会人为地引入额外的超几何参数。例如，x·2F1（1，1; 2;-x）（= log（1 + x））中的上参数等于取消了fb0=1的隐式下参数r。 这导致了一个系数k！/k！在第k项中，虽然在数学上无害，但这会对性能产生重大影响核心图书馆评价机制。log（1 + x）的例子还说明了另一种可能的改进：1的上部参数具有DU，ELEFTHERIOU，MOREIRA和YAP62R-≥比tk+1/tk中的因子1/（k+1）的2个更小的参数。同样，重要的是不要将该因子评估为（1）k/（2）k。更一般地，每当上参数和下参数相差一个整数时，就会发生抵消，并且可以通过认识到这一点来提高效率。我们概述了处理超几何参数的一般算法。 设 1， 2， ， . . . ，ap和b0 ，b1， . ，bq是fpFq的参数s。 我们没有将b0=1添加到较低参数的标准列表中。(1) 我们先对a设a1≤a2≤···≤ap，b0≤b1≤·· ·≤bq是排序后的结果。(2) 通过一个类似合并的算法，我们从两个列表中删除共同的术语请注意，我们仍然保留单独的列表。(3) Nextweformthemaximmnumberof（ai，bj）ofanupper peranda lowerparameter其中a i-bi是整数。让我们称两个实数x，y等价，如果x-y∈Z。设（A i，B i）（i= 1，.，r）是这样的等价对的集合;这些被称为ab-对，因为A i是上参数，B i是下参数。从原始参数列表中删除其对应值A i、B i。很容易看出，ab对的最大数目r是唯一的。然而，这些对的集合不是唯一的。 为了确保最大效率，我们需要将Ai与Bi进行匹配，以使和最小化i=1|A i− B i|.这个(4) 我们计算连续项tn如下：令sn是从剩余的上下参数列表计算的项，如下所示：sn=sn−1×fn其中fn=（a1+n）（a2+n）·· ·（ap+n）/（b0+n）·· ·（bq+n）。然后我们初始化tn到sn。然后对于每个对（A，B），其中B−A=k≥ 1，我们更新tn=tnA（A+1）···（A+k−1）（A+n）···（A+k+n−1）如果A B1，则有一个类似的因子. 有一种特殊类型的对可以进一步利用：当A，B是一半的倍数时（这也可以推广如果A=α/ 2且B−A=k≥ 1，则（A，B）对tn贡献以下因子：α（α+2）···（α+2（k−1））（α+2n）···（α+2（k+n−1））在核心库中，这个公式将再次导致更小深度的表达式和更有效的评估。下表显示了当我们利用参数减少时的加速（在log（1+x）的标准系列∗DU，ELEFTHERIOU，MOREIRA和YAP63{−|−|}√位数100200300400500无预处理（秒）1.5.016.998.8910.46参数预处理（秒）0.220.520.881.351.83加速：4.59.67.96.85.7最小匹配问题首先，我们需要解决以下问题2。 给定两个实数排序列表，（A1 α（j）。 很容易看出，如果我们“不交叉”这样一对，我们得到其和S不大于原解的解。因此，我们只考虑非交叉对。考虑包括输入列表（A1，.，A1）和（B1，...，B j）。设S（i，j）为最小值对于子问题P（i，j）。当i=j时，解是唯一的。否则，对于i<，j，S（i，j）= min S（i，j 1），S（i1，j1）+A iBJ. 这样，我们就可以在O（mn）的时间内解决这个问题。6数学常数：求值、文件夹和访问超几何求值算法需要任意精确的常数。当对三角函数进行参数缩减时，我们需要π。为了减少exp（x）和log（1+x）的参数，我们需要log 2。对于误差函数erf（x），我们需要1/π。我们可以计算这些常数在dummy上，但性能提高了这些常数的预计算，存储在文件中，并根据需要访问它们下表比较了将π计算到一定数量的位所需的时间（使用Machin公式）与从文本文件中读取相同值所比特10010003000500070009000 1000020000上图（秒）0.040.502.495.8810.6617.1921.38107.61预先计算（秒）0.010.010.010.010.020.020.030.11加速比450249588533859713978我们现在描述计算、存储和读取文件格式的常量的设施。一个基本的决定是使用文本文件而不是二进制文件，因为前者是人类可读的。该格式支持整数、浮点数和有理数表示。接下来，数字的基数2这似乎是一个已知的问题，但我们没有找到参考。DU，ELEFTHERIOU，MOREIRA和YAP64可以是二进制、十六进制或十进制。二进制/十六进制的优点是转换成CoreLibrary的内部格式需要线性时间。正式规范随Core Library版本分发1.3或更高3.7结论本文描述了一种将非代数函数纳入支持保证精度的框架中的方法。这样一个包的实现自然是使用核心库的基本功能来实现的。有几个问题的效率和自动化，我们已经解决：自动误差界计算，超几何参数处理，参数减少，最终常数预先计算，存储和检索。全文将描述实现细节，如注册超几何函数的设施今后工作的主题包括：系 列 的 加 速 超 几 何 包 与 核 心 库 一 起 分 发 ， 在 我 们 的 网 站www.example.com上http://cs.nyu.edu/exact/。引用[1] J. - D. Boissonnat和M. Yvinec 几何学剑桥大学出版社，1997年。 把她和布鲁尼曼的电话换了。[2] H. 布伦尼曼角 Burnikel，andS. 皮昂在计算几何中，输入向量算术产生了高效的动态滤波器。ACM Symp.计算几何，14：165 - 174，1998。[3] C.伯尼克尔河Fleischer，K. Mehlhorn和S.席拉精确的几何计算变得简单。第15届ACM Symp.对比几何图形第341-450页[4] C.伯尼克尔河Fleischer，K. Mehlhorn和S.席拉含根式算术表达式的一个强而易计算的分离界。2000年，《自然》，27：87-99。[5] C. Burnikel，S. Funke，K. Mehlhorn和S.席拉A separation bound for real algebraicexpressions，2001.预印本（提交给期刊）。[6] C. Burnikel，S.Funke和M.西尔使用级联的精确几何计算集团有限公司之 J. Computational Geometry and Applications，11（3）：245 -266，2001. 特刊。[7] M. de Berg，M. van Kreveld，M. Overmars和O.施瓦茨科普夫计算几何：算法与应用。Springer-Verlag，柏林，1997年。[8] A. Fabri，G.-吉泽曼湖Kettner，S. Schirra和S. 舍恩赫尔 CGAL内核：几何计算的基础。In M. C. Lin和D. Manocha，编辑，应用计算几何：走向几何工程，第191-202页斯普林格。计算机科学讲义第1148号;第一届ACM应用计算几何研讨会（WACG）论文集，1996年联邦计算研究会议，美国费城。3在progs/fileIO目录下查找。DU，ELEFTHERIOU，MOREIRA和YAP65[9] E. 让德尔你的评价很快，你的作品也很好。RapportTechnique242，INRIA，2000年。 17页[10] V. Karamcheti，C.李岛Pechtchanski和C.是啊一个Core库，用于强大的数值和几何库。 第15届ACM Symp. 计算几何学，第351-359页，1999年。[11] C. Li和C.是啊代数式的一个新的构造性根界。 在Proc.12th ACM-SIAM Symposiumon Discrete Algorithms，第496-505页中。ACM和SIAM，1月。2001年[12] C. Li和C.是啊精确几何计算的新进展In S.巴苏和L. Gonzalez-Vega，编辑，Proc. DIMACS研讨会关于数学和计算机科学中实代数几何的数学和定量方面，3月12日-2001年16DIMACS系列图书美国数学社会，2001年。已提交。 论文下载http://cs.nyu.edu/exact/doc/。[13] G. Liotta，F.和R.塔马西亚鲁棒的邻近查询：度驱动算法设计的例证。SIAMJ.Computing，2001年。出现。[14] Y. V. Matiyasevich.希尔伯特麻省理工学院出版社，马萨诸塞州剑桥，1994年。[15] K. 梅尔霍尔和S。不是她。 LEDA：一个用于组合和几何计算的平台. CACM，38：96[16] J. - M. 监听器手机监听器基本函数：一个出租和实现。 波士顿，1997年。[17] M.奥弗马斯 设计了计算几何算法库CGAL。In M. C. Lin和D. Manocha，编辑，应用计算几何：走向几何工程，第53-58页，柏林，1996年。斯普林格。计算机科学讲义 No.1148;Proc.1stACMWorkshoponAppliedComputationalGeometry（WACG）。[18] S.席拉几何计算的精确性和鲁棒性。 In M. 作者：J. Roos和P. Widmayer，编辑，地理信息系统的数学基础，1340卷讲义Comp. Science，第255-287页。Springer，1997年。第9[19] K.威劳奇可计算分析。施普林格，柏林，2000年。[20] C. 是啊一个新的数字核心，用于强大的数字和几何库。 第三届CGC几何计算研讨会，1998 年。邀请谈话。 布朗大学，1998 年10 月11 日至12 日。摘要见http://www.cs.brown.edu/cgc/cgc98/home.html。[21] C. K. 是啊强大的几何计算。 在j. E. Goodman和J. O'Rourke，编辑，离散和计算几何手册，第35章，第653-668页。CRCPress LLC，Boca Raton，FL，1997.[22] C. K.是啊精确的几何计算。计算几何：理论与应用，7：3-23，1997。邀请谈话，继续。1993年8月5日至9日在滑铁卢举行的第五届加拿大比较几何会议