n阶六边形果园是在边长为n的正六边形中的三角形网格。如下所示是5阶六边形果园:
1

站在中间向外看,图中用绿色标出的点都会被更近的点挡住。可以看出,对于5阶六边形果园,从中间向外看时有30个点会被挡住。

对于n阶六边形果园,记从中间向外看时会被挡住的点的数目为H(n)。

已知H(5) = 30,H(10) = 138,H(1 000) = 1177848。

求H(100 000 000).

分析
下面列出1-10阶六边形挡住点的个数。

阶数 点数 增量
1: 0 0   
2: 6 6   
3: 12 6   
4: 24 12  
5: 30 6  
6: 54 24   
7: 60 6   
8: 84 24   
9: 102 18   
10: 138 36

可以看出当阶数为质数时,增量为6,即只有六边形对角线被挡住。k*n阶上的点会被n阶挡住。 其他情况见下图:
2

在oeis.org上查找的结果:http://oeis.org/A216453

测试2-12的程序,用到了容斥原理

#include <cstdio>
const int NN=12;
int a[NN+1][10]=
{ //prime factor count, each prime factor ...
    {0,0}, //0
    {0,0}, //1
    {0,0}, //2
    {0,0}, //3
    {1,2}, //4
    {0,0}, //5
    {2,2,3}, //6
    {0,0}, //7
    {1,2}, //8
    {1,3}, //9
    {2,2,5},  //10
    {0,0},  //11
    {2,2,3}  //12
};
int f(int n)
{
    int t=0;
    for(int i=1; i<=a[n][0]; i++)
        t+=n/a[n][i]-1;
    for(int j=1; j<a[n][0]; j++) //exclude prime1*prime2 's times
        for(int k=j+1; k<=a[n][0]; k++)
            t-=n/(a[n][j]*a[n][k])-1;
    return t;
}
int main()
{
    int sum=0;
    int x=0;
    for(int i=2; i<=NN; i++)
    {
        x=f(i);
        sum+=x+1;
        printf("%d: %d %d\n",i,x,sum*6);
    }
}

输出

2: 0 6
3: 0 12
4: 1 24
5: 0 30
6: 3 54
7: 0 60
8: 3 84
9: 2 102
10: 5 138
11: 0 144
12: 7 192

最后的程序,需要64位系统,当NN取1e7时,可以在win32上编译。中间犯了一个错误。把p=1; for(int lv=1; lv<=path[cnt][0]; lv++)写成了for(int lv=1, p=1; lv<=path[cnt][0]; lv++),结果是又声明了一个局部变量p,另外,在老版本g++中,不支持int t[n+1]= {0};,要改写成int t[MM+1]= {0};

#include <cstdio>
typedef long long ll;
const int NN=1e8;
const int MM=10; //100000000 <2*3*5*7*11*13*17*19*23 ,MM = max prime factor count
static int path[MM+1][1<<MM][MM+1]= {0}; //set pattern
static int a[NN+1][MM+1]= {0}; //use path
/*
int a[NN+1][10]=
{ //prime factor count, each prime factor ...
    {0,0}, //0
    {0,0}, //1
    {0,0}, //2
    {0,0}, //3
    {1,2}, //4
    {0,0}, //5
    {2,2,3}, //6
    {0,0}, //7
    {1,2}, //8
    {1,3}, //9
    {2,2,5},  //10
    {0,0},  //11
    {2,2,3}  //12
};
*/

void my_subset(const int n,int path[][MM+1])//int *S, const int n)
{
    int t[n+1]= {0};

    for(int i=1; i<=n; i++) t[i]=i; //init t
    int cnt=1;

    for(int lv=1,begin=0,end=0; lv<=n; lv++,begin=end+1,end=cnt-1) //level 1- n
    {
        for(int idx=begin; idx<=end; idx++) //prior level
            for(int i=1; i<=n; i++) //t[i]
                if(t[i]>path[idx][lv-1])//last elem
                {
                    path[cnt][0]=lv;
                    for(int j=1; j<lv; j++)
                        path[cnt][j]=path[idx][j];
                    path[cnt][lv]=t[i];
                    cnt++;
                }
    }
}

static int p[NN+1]= {0};

void fact()
{
    for(int i=2; i<=NN; i++)
        if(p[i]==0)
            for(int j=i*2; j<=NN; j+=i)
            {
                a[j][0]++; //count ,primer's count ==0
                a[j][a[j][0]]=i; //prime
                if(p[j]==0)
                    p[j]=1;
            }

}
int f(int n,int path[][MM+1])
{
    int t=0,p=1,sign=1;
    for(int cnt=1; path[cnt][0]>0; cnt++)
    {
        p=1;
        for(int lv=1; lv<=path[cnt][0]; lv++)
        {
            p*=a[n][path[cnt][lv]];
        }
        t+=(n/p-1)*sign;

        if(path[cnt+1][0]!=path[cnt][0])
            sign=-sign;
    }
    return t;
}

int main()
{
    ll sum=0;
    int x=0;
    for(int i=1; i<=MM; i++)
        my_subset(i,path[i]); //calc for exclude

    fact();

    for(int i=2; i<=NN; i++)
    {
        if(p[i]>0)
            x=f(i,path[a[i][0]]);
        else
            x=0;
        sum+=x+1;
    }
    printf("%lld \n",sum*6);
}

评论

推荐 0
改用short类型存储质因数,因为大于6万的因数只可能是最后一个,所以用short足够
再配合分步计算,在win32上30秒出结果

推荐 0
保持NN不变。
把static int a[NN+1][MM+1]= {0}; //use path
改为int *a;
然后, a=(int*)malloc((NN+1)*(MM+1)*sizeof(int));
memset(a,0,(NN+1)*(MM+1)*sizeof(int));
用a[j*(MM+1)+0]代替a[j][0]
可以编译,但执行过程出错

我要评论

需要登录后才能发言
登录未成功,请修改提交。

× 4644
× 2