尤拉計劃433題:歐幾里德演算法的步數

記E(x₀, y₀)為運用歐幾里德演算法求x₀和y₀的最大公約數的步數。更加正式的說：
x₁ = y₀, y₁ = x₀ mod y₀
x_n = y_n-1, y_n = x_n-1 mod y_n-1
E(x₀, y₀)是使得y_n = 0的最小的n。

已知E(1,1) = 1，E(10,6) = 3以及E(6,10) = 4。

對於1 ≤ x,y ≤ N，定義S(N)為所有E(x,y)的和。
已知S(1) = 1，S(10) = 221以及S(100) = 39826。

求S(5·10⁶)。

求gcd步數方法1，迴圈，方法2，遞迴，方法3，遞迴+快取，方法4，迴圈+快取

#include <cstdio>

#include <cmath>
#define NN 100

static char gs[NN+1][NN+1]= {0};
int a,b;

int srv=0,s=0;


int gcdstep4sieve(int n)
{
    for(a=1; a<=NN; a++)
    {
        gs[a][1]=1;

        for(b=a; b<=NN; b+=a)
        {
            gs[b][a]=1;
        }
    }
    for(a=1; a<=NN; a++)
        for(b=1; b<=a-1; b++)
            for(int n=b; n<=NN; n+=a)
                gs[n][a]=gs[a][b]+1;

    for(a=1; a<=NN; a++)
        for(b=1; b<=NN; b++)
        {
            //if(gs[a][b]==0)printf("(%d,%d)",a,b);
            s+=gs[a][b];
        }
    printf("s=%d\n",s);
}

int gcdstep3(int a,int b)
{
    int r1=a%b;
    if(r1!=0)
    {
        int rv=gs[b][r1];
        if (rv>0)
        {
            srv+=rv;
            s++;
            gs[a][b]=1+rv;
            return 1+rv;
        }
        else
            return 1+gcdstep3(b,r1);

    }
    gs[a][b]=1;
    return 1;
}

int gcdstep2(int a,int b)
{
    int r1=a%b;
    if(r1!=0)
    {
        return 1+gcdstep2(b,r1);
    }
    return 1;
}

int gcdstep(int a,int b)
{
    int a1,b1,r1;
    int c=1;
    a1=a;
    b1=b;
    r1=a1%b1;
    while(r1!=0)
    {
        a1=b1;
        b1=r1;
        r1=a1%b1;
        c++;
    }
    return c;
}
int main()
{

    gcdstep4sieve(NN);
    return 1;
    int cnt=0,e=0;
    for(a=1; a<=NN; a++)
        for(b=1; b<=NN; b++)
        {
            e=gcdstep3(a,b);
            cnt+=e;
            //printf("E(%d,%d)=%d\n",a,b,e);
        }

    printf("S(%d)=%d,s=%d,srv=%d",NN,cnt,s,srv);
    return 0;

}

--結果
S(10000)=782828208
1 11.151
2 11.273
3 2.802
4 2.213

逐個計算太慢,考慮:
1. 所有E(n,1)=1
2. 所有E(n*m,m)=1
3. 所有E(n,m)=1+E(m,(n%m)),即E(r+m*r,m)=1+E(m,r)
如下圖所示，如果能計算出每個返回值各有多少個出現，就能快速合計。
顯然返回值1的總數是N/1+N/2...N/N