一个计算圆周率任意精度的spigot算法研究 - 第4页 - 算法交流 - 数学研发论坛

liangbch 发表于 2011-4-26 11:44:26

以开平方为例，解释一下我程序的功能

楼主的本质是想找到这样一个表达式，对于给定的数N，求正整数p,q和a. 使得p,q,a 满足 sqrt(N)=(q/p)sqrt(1-1/a), 显然p,q是sqrt(N)的一个有理数逼近。下面的代码中的函数print_pq打印出p小于limit且满足以上公式的表达式#include <stdlib.h>
#include <stdio.h>
#include <math.h>

typedef unsigned long DWORD;
typedef__int64 INT64;

void print_pq(INT64 N,INT64 limit)
{
INT64 low_p,low_q,high_p,high_q,mid_p,mid_q;
double f=sqrt(N);

//初始化low_q,low_q,high_p,high_q
low_q=DWORD(floor(f));
high_q=DWORD(ceil(f));
low_p=high_p=1;

mid_p=low_p+high_p;
mid_q=low_q+high_q;

while ( mid_p<limit )
{
INT64 q=mid_q;
INT64 p=mid_p;

if ( double(q)/double(p)>f)
{
INT64 t1=q*q-N*p*p;
INT64 t2=q*q;
if ( t2 % N ==0)
{
t1/=N;
t2/=N;
}
printf("\$sqrt(%I64d)=(%I64d/%I64d)sqrt(1-%I64d/%I64d)\$\n",N,q,p,t1,t2);
}

if ( double(mid_q)/double(mid_p) > f)
{
high_p=mid_p; high_q=mid_q;
}
else
{
low_p=mid_p; low_q=mid_q;
}
mid_p=low_p+high_p;
mid_q=low_q+high_q;
}
}

int main(int argc, char* argv[])
{
print_pq(2,1000000I64);
print_pq(3,1000000I64);
return 0;
}

plp626 发表于 2011-4-26 11:46:01

31# liangbch

XX？
没看明白。。

plp626 发表于 2011-4-26 11:48:59

$sqrt(n)=(q/p)sqrt(1-1/a)$
liangbch 发表于 2011-4-26 11:44 http://bbs.emath.ac.cn/images/common/back.gif

就是啊，我怎么受启发选了那个级数呢，这个更简单些。。

我再研究下。。

liangbch 发表于 2011-4-26 11:50:03

贴出31楼的运行结果$sqrt(2)=(3/2)sqrt(1-1/9)$
$sqrt(2)=(10/7)sqrt(1-1/50)$
$sqrt(2)=(17/12)sqrt(1-1/289)$
$sqrt(2)=(58/41)sqrt(1-1/1682)$
$sqrt(2)=(99/70)sqrt(1-1/9801)$
$sqrt(2)=(338/239)sqrt(1-1/57122)$
$sqrt(2)=(577/408)sqrt(1-1/332929)$
$sqrt(2)=(1970/1393)sqrt(1-1/1940450)$
$sqrt(2)=(3363/2378)sqrt(1-1/11309769)$
$sqrt(2)=(11482/8119)sqrt(1-1/65918162)$
$sqrt(2)=(19601/13860)sqrt(1-1/384199201)$
$sqrt(2)=(66922/47321)sqrt(1-1/2239277042)$
$sqrt(2)=(114243/80782)sqrt(1-1/13051463049)$
$sqrt(2)=(390050/275807)sqrt(1-1/76069501250)$
$sqrt(2)=(665857/470832)sqrt(1-1/443365544449)$
$sqrt(3)=(7/4)sqrt(1-1/49)$
$sqrt(3)=(26/15)sqrt(1-1/676)$
$sqrt(3)=(97/56)sqrt(1-1/9409)$
$sqrt(3)=(362/209)sqrt(1-1/131044)$
$sqrt(3)=(1351/780)sqrt(1-1/1825201)$
$sqrt(3)=(5042/2911)sqrt(1-1/25421764)$
$sqrt(3)=(18817/10864)sqrt(1-1/354079489)$
$sqrt(3)=(70226/40545)sqrt(1-1/4931691076)$
$sqrt(3)=(262087/151316)sqrt(1-1/68689595569)$
$sqrt(3)=(978122/564719)sqrt(1-1/956722646884)$

plp626 发表于 2011-4-26 11:54:29

31# liangbch

不过，(q/p)sqrt(1-b/a)相比较于(q/p)*1/sqrt(1-b/a)那个级数通项稍复杂点点.

liangbch 发表于 2011-4-26 11:59:42

贴出 (1-x) n次方根的展开式：

plp626 发表于 2011-4-26 12:47:16

这个是1/root{m}{1-x}的展示

1/root{m}{1-x}=1+x/m+(1+m)/{m^2*2!}*x^2+{(1+m)(1+2m)}/{m^3*3!}*x^3+......
=(1+1/mx(1+{1+m}/{2m}x(1+{1+2m}/{3m}x(1+...))))
|x|<1

------------------------------------------------

root{m}{1-x}=1-x/m-(m-1)/{m^2*2!}*x^2-{(m-1)(2m-1)}/{m^3*3!}*x^3+......
=(1-x/m(1+{m-1}/{2m}x(1-{2m-1}/{3m}x(1+...(1-... )))))

liangbch 发表于 2011-4-26 12:53:10

OK, 换成 $p/q*1/sqrt(1-1/a)$的形式#include <stdlib.h>
#include <stdio.h>
#include <math.h>

typedef unsigned long DWORD;
typedef__int64 INT64;

void print_pq(INT64 N,INT64 limit)
{
INT64 low_p,low_q,high_p,high_q,mid_p,mid_q;
double f=sqrt(N);

//初始化low_q,low_q,high_p,high_q
low_q=DWORD(floor(f));
high_q=DWORD(ceil(f));
low_p=high_p=1;

mid_p=low_p+high_p;
mid_q=low_q+high_q;

while ( mid_p<limit )
{
INT64 q=mid_q;
INT64 p=mid_p;

if ( double(q)/double(p)<f)
{
INT64 t1=p*p*N-q*q;
INT64 t2=p*p*N;
if ( t1 % N ==0)
{
t1/=N;
t2/=N;
}
printf("sqrt(%I64d)=\$(%I64d/%I64d)1/sqrt(1-%I64d/%I64d)\$\n",N,q,p,t1,t2);
}

if ( double(mid_q)/double(mid_p) > f)
{
high_p=mid_p; high_q=mid_q;
}
else
{
low_p=mid_p; low_q=mid_q;
}
mid_p=low_p+high_p;
mid_q=low_q+high_q;
}
}

int main(int argc, char* argv[])
{
print_pq(2,1000000I64);
print_pq(3,1000000I64);
return 0;
}

liangbch 发表于 2011-4-26 13:05:16

贴出楼上代码的运行结果

sqrt(2)=$(4/3)1/sqrt(1-1/9)$
sqrt(2)=$(7/5)1/sqrt(1-1/50)$
sqrt(2)=$(24/17)1/sqrt(1-1/289)$
sqrt(2)=$(41/29)1/sqrt(1-1/1682)$
sqrt(2)=$(140/99)1/sqrt(1-1/9801)$
sqrt(2)=$(239/169)1/sqrt(1-1/57122)$
sqrt(2)=$(816/577)1/sqrt(1-1/332929)$
sqrt(2)=$(1393/985)1/sqrt(1-1/1940450)$
sqrt(2)=$(4756/3363)1/sqrt(1-1/11309769)$
sqrt(2)=$(8119/5741)1/sqrt(1-1/65918162)$
sqrt(2)=$(27720/19601)1/sqrt(1-1/384199201)$
sqrt(2)=$(47321/33461)1/sqrt(1-1/2239277042)$
sqrt(2)=$(161564/114243)1/sqrt(1-1/13051463049)$
sqrt(2)=$(275807/195025)1/sqrt(1-1/76069501250)$
sqrt(2)=$(941664/665857)1/sqrt(1-1/443365544449)$
sqrt(3)=$(3/2)1/sqrt(1-1/4)$
sqrt(3)=$(5/3)1/sqrt(1-2/27)$
sqrt(3)=$(12/7)1/sqrt(1-1/49)$
sqrt(3)=$(19/11)1/sqrt(1-2/363)$
sqrt(3)=$(45/26)1/sqrt(1-1/676)$
sqrt(3)=$(71/41)1/sqrt(1-2/5043)$
sqrt(3)=$(168/97)1/sqrt(1-1/9409)$
sqrt(3)=$(265/153)1/sqrt(1-2/70227)$
sqrt(3)=$(627/362)1/sqrt(1-1/131044)$
sqrt(3)=$(989/571)1/sqrt(1-2/978123)$
sqrt(3)=$(2340/1351)1/sqrt(1-1/1825201)$
sqrt(3)=$(3691/2131)1/sqrt(1-2/13623483)$
sqrt(3)=$(8733/5042)1/sqrt(1-1/25421764)$
sqrt(3)=$(13775/7953)1/sqrt(1-2/189750627)$
sqrt(3)=$(32592/18817)1/sqrt(1-1/354079489)$
sqrt(3)=$(51409/29681)1/sqrt(1-2/2642885283)$
sqrt(3)=$(121635/70226)1/sqrt(1-1/4931691076)$
sqrt(3)=$(191861/110771)1/sqrt(1-2/36810643323)$
sqrt(3)=$(453948/262087)1/sqrt(1-1/68689595569)$
sqrt(3)=$(716035/413403)1/sqrt(1-2/512706121227)$
sqrt(3)=$(1694157/978122)1/sqrt(1-1/956722646884)$

plp626 发表于 2011-4-26 13:13:33

本帖最后由 plp626 于 2011-4-26 14:07 编辑

我现在发现把b/a强制换成1/a会使得p,q数值过大，p超过(2^32/10000)，做分母的时候，得到的商就不够4位了，这样不好统一输出。。。

当然b的绝对值越小越好，还是假定b的绝对值小于10把，这样不会使得p太大，又适合于32位整数的直接除法

页: 1 2 3 [4] 5 6 7 8

数学研发论坛's Archiver