基于泰勒展开式的高精三角函数实现

liangbch

问题，对于任意正数x,求其任意精度的自然对数log(x).

第一步是找到一个分数P/Q, 使得
  1) P和Q的素因子都是小素数，
  2) x 可表示为 P/Q * (1 + d)的形式，并且d足够小。

若Pn表示第n个素数，例，P1=2，P2=3,P3=5
我现在的目标是，
1. P和Q的所有素因子都不大于P48
2. d小于 \( \frac{1}{2^{96}} \)

zeroieme

liangbch 发表于 2016-5-11 10:58
问题，对于任意正数x,求其任意精度的自然对数log(x).

第一步是找到一个分数P/Q, 使得

2进制处理
1、容易求a，$1<=2^(-a)x<2$
2、预计算 (1+2^(-1))^(-1)， (1+2^(-2))^(-1)……， (1+2^(-96))^(-1)的P/Q贴近形式以及所需精度近似值b1,b2,b3……。
3、令$y0=2^(-a)x$，假如y0的$2^(-1)$位为1,y1=b1 y0，否则,y1= y0
……

liangbch

此算法不失为一种好的方法，我之前没想到。但 P/Q= (p1/q1)* (p2/q2) * ... (p96/q96)中的P和Q偏大，不如我的方案好。

落叶

因为上次代码调用的开平方程序是我自已写的，因不知道高效的开平方算法，效率很低，现改调用GMP开平方函数mpz_sqrt（）
这段代码
        For i = 1 To 45                                '针对一万精度大约进行45次分解
                   y = myKpf(y)                         ‘开平方
        Next
        对y进行泰勒级数运算(没能象sinx一样找到简化算法，现只能老老实实硬算)，针对一万精度采用（精度除以13）次泰勒级数运算。
改为
        For i = 1 To 240                              '针对一万精度大约进行240次分解，因为开方函数快了，让它多做点，能者多劳。
                   y = myKpf(y)                         ‘开平方
        Next
        对y进行泰勒级数运算，针对一万精度采用（精度除以70）次泰勒级数运算。
现Ln(9),一万精度运行时间7.5m,疯狂计算器运行36m

不过mpz_sqrt（）只支持正整数，开整数次方，前期传参，和后期指数对位有点麻烦，现把我的方法贴出来，只是简单测试了一下，不一定完全正确，主要是介绍我的思路
前期传参：
       所有的数全部换算成这种形式：2.51E2，可以用类似我前面介绍的高精度头表示
       strlenZ = (jd * 2) * 1.2       ‘jd是你所需的开方精度，strlenZ是你传给mpz_sqrt的参数长度，按说传20位的长度（函数返回10位数）应该获得10位的精度，不过传回的数的后百分之十位数精度不对，所以加乘1.2倍长度
        strlenZ = (strlenZ \ 2) * 2
        If (string1.eE Mod 2) = 0 Then          ’ 这里很重要，例：251.123 ，你需要传奇数长的数（后面补零），内存情况为2511230，例：25.1123你需要参偶数长的数 ，内存情况为251123            
            strlenZ = strlenZ + 1
        End If
        ReDim op1(-1 To strlenZ + 2) As Byte               'op1就是你传的参数，这里定义它的长度

后期指数对位：
       If string1.eE >= 0 Then                                                ‘string1.eE为被开方数用科学计数法整理后的指数
            jg.eE = (string1.eE + 1) \ 2 + ((string1.eE + 1) Mod 2) - 1             ’jg.eE 为最终结果的指数
        Else
            jg.eE = string1.eE \ 2 + string1.eE Mod 2
        End If




补充内容 (2016-5-26 17:41):
现Ln(9),一万精度运行时间7.5m,疯狂计算器运行36m
改为现Ln(9),一万精度运行时间7.5秒,疯狂计算器运行36秒

落叶

针对泰勒展开式：ln(1+x)=x-x^2/2+x^3/3-……+(-1)^(k-1)*(x^k)/k(|x|
转换成：ln(1+x)=x(1-x(1/2-x(1/3-x(1/4-x(1/5..........)))))
这样转换并没有提高效率，但可以预算1/2 1/3 1/4 1/5.....针对一万精度，因为经过了开方优化，这里只需要10000\70=133次泰勒级数运算，所以可以预制133个万位长数的表格供调用，
粗略估计可以获得一秒的速度提升。

落叶

还可转换成：ln(1+x)=(x(N!-x(N!/2-x(N!/3-x(N!/4-x(N!/5..........)))))   )/N!        N是你需要运算的泰勒级数
这个转换不仅不能提高效率，在我的程序万位运算中，速度还有所下降。
但它还是有用的：
现在计算十万精度的ln(9),1-100000的倒数绝大多数都是无理数，我的意思就是除不断，为了获得十万精度，这些倒数你需要计算或预制十万位，100000\70=1428，这里需要1428次泰勒级数运算，
1428！的位长是3886，1428！除以（1-100000）中的任何数都能除断（整除），都只需要3886位除法运算（当然这里需要除法子程序有判断除法除尽的功能），这样，1428个十万位除法就变成了（1428个3886位除法，再加上一个十万除以十万的除法），  为什么在我程序是没有效果是因为我的精度要求过低造成了泰勒级数运算过少。说明一下：1\99999取十万位，虽然1只有一位，但还是要经过十万位运算，1需要扩展到十万位才能进行除法运算。

落叶

今天专门针对ln(1+x)=x(1-x(1/2-x(1/3-x(1/4-x(1/5..........))))) 进行了测试
ln(9)，精度十万，泰勒级数1428，用时1322秒
针对ln(1+x)=(x(N!-x(N!/2-x(N!/3-x(N!/4-x(N!/5..........)))))   )/N!    进行了测试
ln(9)，精度十万，泰勒级数1428，用时28秒，提速47倍
精度十万，1/99999,用时820ms,
精度十万，1428!/99999,用时10-15ms,除法运算的提速和总程序的提速大致吻合。

针对我经常讲的精度控制，可能很多人不理解，如果要进行算式运算或综合运算，没有精度控制，整个程序就和没有监督的权力一样，是不健全的....,错误多多，
我早期的程序就没有精度控制，结果在实现三角函数和开N次方函数时，一直不成功。
N!就是1-N的公倍数，所以能整除，1-N的最少公倍数是多少？我却不知道？
这个公式的提速原理是精度大于N！的位长，大多少公式开始提速？要通过测试确定，当然是越大越好。
第二个提速原理就是尽量把除不尽的除法换算成能整除的除法。

补充内容 (2016-5-27 08:10):
大多少，公式开始提速？

zeroieme

落叶 发表于 2016-5-26 10:24
今天专门针对ln(1+x)=x(1-x(1/2-x(1/3-x(1/4-x(1/5..........))))) 进行了测试
ln(9)，精度十万，泰勒级数 ...

所有小于N的素数pi的乘幂  pi^[Log(N)/Log(pi)]，的乘积，方括号[x]表示小于等于 x 的最大整数

zyqjw

BTCAL2.8专门针对ln函数进行了优化，对于任意数特别是极大数如Ln(1E999999999999999999)，万位精度都已控制在4秒内。

落叶

目前我的程序也解决了在万进制或千进制下的小数对位问题，ln(9)现万位需要时间3.4秒，比你的2.7版首次运算快，因为我预算ln(10)，你的再次运算因已存了ln(2)，只需要2秒，还是比我的快，目前我所能优化的程序运行效率都只有你的一半左右，如三角函数现1.5秒左右，你的0.7左右，大的算法应该和你没有大的区别，目前分析主要还是乘法比你的稍慢，另外vb的效率也还是低了点，如果我还是也vb为主，效率上很难有大的改进了。