在高精度计算中，连分式展开相对泰勒展开式究竟有多少优点？

zeroieme 发表于 2016-4-30 13:38:23

对连分式感兴趣1是数值计算原理教科书上的栗子：Ln2的连分式R44比泰勒展开式8阶精确10万倍；2是Mathematica的实现说明也写很多函数是使用连分式计算的；3是“连分式可以减少乘除法计算量。”

然而动手编码却发现没那么简单：
先说减少乘除法计算量，恐怕是针对浮点机器数的，一次除法大概等于乘法。但高精度除法一次大概等于10次乘法。避免每次做除法就使用递推法分别计算分子分母，最后再除法，那么也是每阶要4次乘法。还是不如秦九韶法计算泰勒展开式。

然后就是“同等阶数”下连分式和泰勒展开式的精度问题。原来函数做连分式展开式没有唯一方式。
使用欧拉公式把泰勒展开式转化为连分式？精度同等的，如前所说，没有计算量优势。
原始方式，一次次计算一次泰勒展开？除非找到系数规律，不然对于高精度计算需要上百阶的情形没意义。
还有就是高斯的椭圆函数连分式了，估计Mathematica的“很多函数”就是用这个了。

liangbch 发表于 2016-5-2 16:30:50

高精度除法一次大概等于10次乘法

我不这么认为，我认为一次高精度除法的运行时间大概等于1.5到2次高精度乘法。

zeroieme 发表于 2016-5-2 17:18:09

liangbch 发表于 2016-5-2 16:30
我不这么认为，我认为一次高精度除法的运行时间大概等于1.5到2次高精度乘法。

a/b,前面逐次倍增精度是n*log依位数不同有差别。单算最后倒数 x*b=1,接着 a*x已经是2次。

推一步讲即使除法优化到1.5到2次高精度乘法时间，一般情况下依然是不如泰勒展开。

liangbch 发表于 2016-5-2 19:42:40

使用牛顿迭代法计算b的倒数并精确到p位数字，当p比较小时（在这种情况下，最快的大数乘法依然是硬乘法），计算r=1/b的复杂度约等于0.5次p位乘以p位的大数乘法，故总的复杂度是1.5次。
“前面逐次倍增精度是n*log依位数不同有差别”，
我曾经做过这样的事情，我记得，
若需精确都p位数，则精度逐次倍增时，每次迭代的运算量是上一次迭代的2倍，最后一次运算量是（0.5 * p）^2, 这是一个公比为0.5的等比数列，故总的迭代过程的运算量位 0.25*p*p /(1-0.5)=0.5*p *p.故有这个结论，你可重新检查一下，看看是否有误。

by the way, 我并没有否定你的结论,即，一般情况下，泰勒展开式更快。
不过对于2次方根，其连分数形式一定是一个循环简单连分数，求其连分数形式也不复杂。若求得其“循环节”，其连分数法不失为一种重要的方法。

初等数论是这样说的。
实2次无理数的无限简单连分数表示式一定是一个循环连分数。
对于整系数方程 a(x^2) + bx+c=0,若 d=b^2-4ac>0, 则其根x0一定是 r+s*sqrt(d)的形式，这样形式的无理数称之为实2次无理数。

Pell 方程的解法就是用到这个性质。

liangbch 发表于 2016-5-2 19:53:05

关于整数平方根的计算，可参阅下面两篇博文
http://blog.csdn.net/liangbch/article/details/7263059
http://blog.csdn.net/liangbch/article/details/7336403

liangbch 发表于 2016-5-2 20:15:48

本站帖子 http://bbs.emath.ac.cn/thread-143-1-10.html 讨论了使用牛顿迭代法求高精度平方根的方法，可惜在论坛转移的过程中，代码无法下载。我在我的电脑上搜索了下，找到一份代码，不知道有没有问题，见下。

// sqrt2.cpp : Defines the entry point for the console application.
//
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <iostream.h>

#include "../../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
#include "../../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeInt.h" // 10进制系统
#include "../../../HugeCalc/HugeCalc_API/CppAPI/Include/HugeIntX.h" // 16进制系统

#pragma message( "automatic link to ../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#pragma comment( lib, "../../../HugeCalc/HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )

#define integer CHugeInt

#define TEN9 1000000000
/*
sqrt(a) 算法：

step1: 得到一个 1/sqrt(a) 近似值 x0
step2: x= 3/2*x + a/2 x^3
另一个形式： ( x= x + x ( 1/2- a/2 *x^2)
不断计算step2,x精度逐步提高，每次迭代精度提高一倍。
最后：计算 x* a 得到 sqrt(a)
*/
double int_sqrt_demo(int n)
{
double x0;
double x1;
double c0,c1;

c0=0.5;
c1=double(n)/2;

x0= 1 /sqrt( n );
x0= floor(x0 * 10000.0) /10000; //取前4位有效数字

while (true)
{
x1= x0 + x0 * (c0- c1 * x0 * x0);
if ( (x1- x0)/x0 <1e-15)
{
break;
}
x0=x1;
}
return double(n) * x1;
}

int dig2Unit(int n)
{
return (n+8)/9;
}

void sqrt_proc( double n, integer *x,int lmt_pre[],double limitLen)
{
int i=0;
char *pBuff=NULL;
const char *pTmp;
double esp=1e-15;
UINT32 digLen;//数字个数

integer x0;
integer x1;
integer t1;
integer t2;

integer c0;
integer c1;

{
unsigned int n_01,n_02;
n/=2;
c0= (TEN9/2); // c0=0.5;
n_01= UINT32(n);
n_02= UINT32((n-n_01+esp)*TEN9);

c1 = n_01;
c1 *= TEN9;
c1 += n_02; // c1=double(n)/2;
#ifdef _DEBUG
if ( pBuff!=NULL)
{ delete[] pBuff; }

pBuff=new char;
pTmp= c1.GetStr(FS_NORMAL);
strcpy(pBuff,pTmp);
printf("c1=%s\n",pBuff);

#endif

}

x0= *x;
i=0;
digLen=x0.GetDigits();

while (true)
{

#ifdef _DEBUG
printf("第%d次迭代\n",i+1);
#endif

if (i>0 && digLen -lmt_pre>=9)
{
int shift= (digLen-lmt_pre)/9*9;
x0.DecRShift( shift ); //丢弃多于的数字
digLen=x0.GetDigits();
}

#ifdef _DEBUG
if ( pBuff!=NULL)
{ delete[] pBuff; }

pBuff=new char;
pTmp= x0.GetStr(FS_NORMAL);
strcpy(pBuff,pTmp);
printf("x0=%s\n",pBuff);
#endif

t1.Pow( x0,2); //t1= x0 * x0

#ifdef _DEBUG
if ( pBuff!=NULL)
{ delete[] pBuff; }
pBuff=new char;
pTmp= t1.GetStr(FS_NORMAL);
strcpy(pBuff,pTmp);
printf("t1=%s\n",pBuff);
#endif

// x1= x0 + x0 * (c0- c1 * x0 * x0);
// c0=0.5;
// c1=double(n)/2;

t1 *= c1; // t1= c1* x0 * x0
t1.DecRShift(9);
#ifdef _DEBUG
if ( pBuff!=NULL)
{ delete[] pBuff; }

pBuff=new char;
pTmp= t1.GetStr(FS_NORMAL);
strcpy(pBuff,pTmp);
printf("t1=%s\n",pBuff);
#endif

t2=c0;
t2.DecLShift( dig2Unit(digLen)*2*9-9);//扩展数位

#ifdef _DEBUG
if ( pBuff!=NULL)
{ delete[] pBuff; }

pBuff=new char;
pTmp= t2.GetStr(FS_NORMAL);
strcpy(pBuff,pTmp);
printf("t2=%s\n",pBuff);
#endif

t2-=t1; // t2=(c0- c1* x0 * x0);
t2 *= x0; // t2=x0 * (c0- c1* x0 * x0);
t2.DecRShift(dig2Unit(digLen)*9);

#ifdef _DEBUG
if ( pBuff!=NULL)
{ delete[] pBuff; }
pBuff=new char;
pTmp= t2.GetStr(FS_NORMAL);
strcpy(pBuff,pTmp);
printf("t2=%s\n",pBuff);

#endif
x0.DecLShift( dig2Unit(digLen)*9 );//扩展x0长度

#ifdef _DEBUG
if ( pBuff!=NULL)
{ delete[] pBuff; }
pBuff=new char;
pTmp= x0.GetStr(FS_NORMAL);
strcpy(pBuff,pTmp);
printf("x0=%s\n",pBuff);
#endif
x1= x0 + t2; // x1= x0 + x0 * (c0- c1 * x0 * x0);

#ifdef _DEBUG
if ( pBuff!=NULL)
{ delete[] pBuff; }
pBuff=new char;

pTmp= x1.GetStr(FS_NORMAL);
strcpy(pBuff,pTmp);
printf("x1=%s\n",pBuff);
#endif

x0=x1;
digLen= x0.GetDigits();
if ( digLen >= limitLen)
{
break;
}
i++;
}

*x= x0;
if (pBuff!=NULL)
{
delete[] pBuff;
}
}

//高精度整数平方根
void int_sqrt(int num, int dig,char *fileName)
{
double fNum;
double fLen;
int e; //100^e 表示需要把被开方数缩小 100^e,最终结果需要扩大10^e倍
int i,k;
int tmp;

int lmt_pre;//第i次迭代最低要求精度,单位10进制位

UINT32 i0;
char *pBuff=NULL;
const char *pTmp;

HugeCalc::EnableTimer( TRUE );
HugeCalc::ResetTimer();

integer x;

tmp=num;
e=0;
while (tmp>=100)
{
tmp /= 100;
e++;
}

fNum= (double)(num) / pow(100.00,e);

k=0;
fLen=dig+18;
while (fLen >=9)
{
fLen /= 2.0;
k++;
}
for (i=0;i<=k;i++)
{
lmt_pre=(int)(ceil(fLen));
fLen *=2;
}

i0= (UINT32)floor(1.0/sqrt(fNum) * TEN9+0.5);

if (i0>=TEN9)
i0=TEN9-1;

x=i0;
sqrt_proc(fNum, &x,lmt_pre,dig+18);
x *= integer(num);

if ( pBuff!=NULL)
{ delete[] pBuff; }

pBuff=new char;

pTmp= x.GetStr(FS_NORMAL);
strcpy(pBuff+1,pTmp);

int intLen= (int)ceil(log10(sqrt(num))); //计算整数部分的数字个数

memcpy(pBuff,pBuff+1,intLen);
pBuff='.'; //插入小数点

pBuff=0; //去掉多余的数字

// 记录计算耗时
UINT32 u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
double ftime= u32Timer_Calc;
printf("计算用时 %.6f 秒\n", ftime/1000000.0);

HugeCalc::EnableTimer( TRUE );
HugeCalc::ResetTimer();

FILE *fp=fopen(fileName,"wb");
fwrite(pBuff,dig+1,1,fp);
fclose(fp);

u32Timer_Calc = HugeCalc::GetTimer( TIMER_UNIT_us );
ftime= u32Timer_Calc;
printf("输出到文件用时 %.6f 秒\n", ftime/1000000.0);

delete[] pBuff;
pBuff=NULL;
}

void test_demo()
{
double f1=int_sqrt_demo(2);
double f2=int_sqrt_demo(3);
double f3=int_sqrt_demo(10);
double f4=int_sqrt_demo(1000);
}

void nums_sqrt()
{
int n;
int dig;
char fileName;

printf("calc sqrt(n)\n");

while (true)
{
printf("n=?(0 exit)");
scanf("%d",&n);

if (n==0 )
break;
printf("how many digital is needed?");
scanf("%d",&dig);
if (dig<=0)
break;
sprintf(fileName,"sqrt_%d_%08d.txt",n,dig);
int_sqrt(n,dig,fileName);
}

}

int main(int argc, char* argv[])
{
//test_demo();
nums_sqrt();
return 0;
}

zeroieme 发表于 2016-5-2 23:42:16

liangbch 发表于 2016-5-2 19:42
使用牛顿迭代法计算b的倒数并精确到p位数字，当p比较小时（在这种情况下，最快的大数乘法依然是硬乘法）， ...

我不是想辩论什么，只是想对前段时间对连分式的学习作个小结，虽然结论并不那么美好。

对于2次方根或者若干次方根的连分式有这几个维基百科/(上面再推荐)的地址
https://en.wikipedia.org/wiki/Solving_quadratic_equations_with_continued_fractions
https://en.wikipedia.org/wiki/Generalized_continued_fraction
http://myreckonings.com/Dead_Reckoning/Online/Materials/General%20Method%20for%20Extracting%20Roots.pdf

如你下面说的，对于2次方根我也倾向用牛顿法。
我对高精度计算分了以下层次。是以我觉得实现的难度，而不是按数学上算术、初等函数等区分。请给些意见。

1 加减法，基本上就是延长数位和进位了，线性复杂度。

2 乘法，直乘或Toom-Took，FFT，数论等优化，要依据不同情况采用不同方法。 /*都试写过代码，而不同方法的选择条件反而是我最困惑。*/

3 能够表示成加减乘法方程的函数，除法开方和解多项式方程。 /*看了你对除法的分析，似乎是理论上达到精度就直接算下一步了。而我是代入方程确认达到精度，这里就多做了一次计算。*/

4 就是更复杂的其他函数。/*现在看来依然是使用泰勒展开合适*/

5 个别采用特殊方法的函数：log（/ *算术几何平均，但我对初期开方优化不好*/） ,高斯连分式还有挖掘的地方。

落叶发表于 2016-5-3 00:12:51

我的除法现在也是直接估出迭代次数，可以减少一次迭代。

页: [1]

数学研发论坛's Archiver

在高精度计算中，连分式展开相对泰勒展开式究竟有多少优点？