Форум умных людей

#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;
}