программа на С.
число вычисляемых знаков пропорционально N:(2.4*N).
#include <stdio.h>
#include <stdlib.h>
#define N 500
#define u16 unsigned short
#define u32 unsigned
typedef struct{u16 q[N];short nz;}u16N;
typedef struct{u16N u0,u1;}u16N2;
u16N new16N(){u16N a;for(int i=0;i<N;i++)a.q[i]=0;a.nz=0;return a;}
u16N toU16N(u32 n){u16N a=new16N();a.q[0]=n;a.q[1]=(n>>16);a.nz=a.q[1]?2:(a.q[0]?1:0);return a;}
u16N sum(u16N a,u16N b){
int m=a.nz>=b.nz?a.nz:b.nz;
u16N s=new16N();
u32 c=0;
for(int i=0;i<m;i++){u32 t=c+a.q[i]+b.q[i];s.q[i]=t;c=t>>16;}
if(c&&m<N)s.q[m]=1;
s.nz=nZ(s.q);
return s;
}
u16N inv(u16N a){
u16N b;
for(int i=0;i<N;i++)
b.q[i]=~a.q[i];
b.nz=nZ(b.q);
return b;
}
u16N subt(u16N a, u16N b){return sum(a,sum(inv(b),toU16N(1)));}
int nZ(u16 q[N]){int cnt=N;for(int i=N-1;i>=0;i--)if(q[i])break;else cnt--;return cnt;}
u16N shL0(u16N a,int n){
u16N b=new16N();
u16 c=0;
for(int i=0;i<a.nz;i++){u32 r=c|((u32)a.q[i]<<n);b.q[i]=r;c=r>>16;}
if(a.nz<N)b.q[a.nz]=c;
b.nz=nZ(b.q);
return b;
}
u16N shL(u16N a,int n){
u16N b=shL0(a,n&0xf);
int n1=n>>4;
if(n1)for(int i=a.nz;i>=0;i--){int j=i+n1;if(j<N)b.q[j]=b.q[i];if(i<N)b.q[i]=0;}
b.nz=nZ(b.q);
return b;
}
u16N shR(u16N a){
u16N b=new16N();
for(int i=a.nz-1;i>=0;i--){u32 u=(i==a.nz-1?0:a.q[i+1]);u32 w=(u<<16)|a.q[i];b.q[i]=(w>>1);}
b.nz=nZ(b.q);
return b;
}
int max1(u16N a){
int k=-1;
if(a.nz){for(int i=15;i>=0;i--)if((a.q[a.nz-1]>>i)&1){k=i;break;}k+=(a.nz-1)<<4;}
return k;
}
int cmp(u16N a, u16N b){
int s=a.nz-b.nz;
if(s)return s;
for(int i=a.nz-1;i>=0;i--)if(s=a.q[i]-b.q[i])return s;
return 0;
}
void show10(u16N a, u16N b){
u16 m[N];
for(int i=N-1;i>=0;i--){
u16 r=0;
for(int j=0;j<16;j++){
b=shR(b);
r<<=1;
if(cmp(a,b)>=0){a=subt(a,b);r|=1;}
}
m[i]=r;
}
int n=(int)(N*0.6);
u16 d[n];
int cnt=0;
int k=0;
while(k<N&&cnt<n){
u32 c=0;
for(int i=k;i<N;i++){
u32 t=(u32)10000*m[i]+c;
m[i]=t;
c=t>>16;
}
d[++cnt-1]=c;
for(int i=k;i<N;i++)
if(!m[i]){k=i;break;}
}
for(int i=0;i<cnt;i++){printf(" %04d",d[i]);if((i+1)%15==0)printf("\n");}
printf("\n");
}
u16N fn(int n,u16N*p){
if(n>1){
p[n]=shL(p[n-1],1);
if(!(n&1))p[n]=subt(p[n],p[n>>1]);
}
return p[n];
}
int main()
{
int size=N<<3;
u16N* a=(u16N*)malloc(size*sizeof(u16N));
a[0]=toU16N(0);a[1]=toU16N(1);
u16N numin=toU16N(0);
int i=0;
while(i<size-1&&(N<<4)-max1(numin)>4)
numin=sum(shL0(numin,2),fn(++i,a));
u16N denom=shL(toU16N(1),i+i-1);
show10(numin,denom);
return 0;
}