POJ 3090 Visible Lattice Points 法雷级数

来源:互联网 发布:mac版微信开发工具 编辑:程序博客网 时间:2024/06/06 21:02

题意:从原点看第一象限里的所有点,能直接看到的点的数目是多少。(不包含原点)

 

法雷级数定义  R.亨斯贝尔格著李忠翻译的《数学中的智巧》一书,介绍了法雷级数。这里每一行从0/1开始,以1/1结尾,其它数自左至右将所有的真分数按增加顺序排列;第n行是由所有分母小于或等于n的真分数组成,我们称为n阶法雷级数。如下表:

  F1: 0/1 1/1

  F2: 0/1 1/2 1/1

  F3: 0/1 1/3 1/2 2/3 1/1

  F4: 0/1 1/4 1/3 1/2 2/3 3/4 1/1

  F5: 0/1 1/5 1/4 1/3 2/5 1/2 3/5 2/3 3/4 4/5 1/1

  F6:0/1 1/6 1/5 1/4 1/3 2/5 1/2 3/5 2/3 3/4 4/5 5/6 1/1

  …… ………………………………

  这里我们想问的是第n行Fn的真分数的个数有多少个呢?

  我们设Fn的个数为ψ(n), ψ(n)比 ψ(n-1)增加的个数是分母是n,分子比n小且与n互质的数的个数,这正是欧拉函数φ(n)。即

  ψ(n)=ψ(n-1)+ φ(n)

  ψ(1)=1+φ(1)

  ψ(2)=ψ(1)+φ(2)

  ψ(3)=ψ(2)+φ(3)

  ………………

  ψ(n)= ψ(n-1)+ φ(n)

  所以 ψ(n)=1+φ(1)+φ(2)+φ(3)+……+φ(n)很容易证明,当n≥3时,欧拉函数φ(n)是个偶数。由此我们得到除ψ(1)=2是偶数外,法雷级数其它各级的个数都是奇数,并且许多是素数。ψ(1)=2,ψ(2)=3,ψ(3)=5,ψ(4)=7,ψ(5)=11,ψ(6)=13,ψ(7)=19,ψ(8)=23,ψ(9)=29,……。

性质

  法雷级数Fn具有很多美妙的性质,下面是一些常见的性质:

  1.如果a/b,c/d是相邻的两项,则abs(a*d-b*c)=1。

  2.如果a/b,c/d,e/f是相邻的三项,则 (a+e)/(b+f)=c/d,特别的,如果c/d是新添加的,即c/d不属于F(n-1),则c=a+e;d=b+f。

  性质2对于这个问题至关重要,它的证明可以参见哈代(Hardy)写的数论导引第三章

  关于Farey级数的介绍。根据这条性质可以知道,丛F(n−1)到F(n)的构造过程中,F(n)的新项的分母一定是其相领两项的分母和。另一方面,如果F(n−1)中的相邻两项 a/b,c/d, b+d=n,则(a+c)/n一定会被添加到F(n)中。

 

#include<cstdio>#include<cstring>using namespace std;#define MAX 1009int p[MAX], a[MAX], pn;int eul[MAX];void prime(){    pn = 0;    memset(a,0,sizeof(a));    int i, j;    for ( i = 2; i < MAX; i++ )    {        if ( a[i] == 0 ) p[pn++] = i;        for ( j = 0; j < pn && i * p[j] < MAX && (p[j] <= a[i] || a[i] == 0); j++ )            a[i*p[j]] = p[j];    }}void Euler_Farey (){    for ( int i = 2; i < MAX; i++ )    {        if ( a[i] == 0 )            eul[i] = i - 1;        else        {            int k = i / a[i];            if ( k % a[i] == 0 ) eul[i] = eul[k] * a[i];            else eul[i] = eul[k] * ( a[i] - 1 );        }    }    eul[1] = 1;    for ( int i = 2; i < MAX; i++ )        eul[i] += eul[i-1];}int main(){    int n, t;    prime();    Euler_Farey();    scanf("%d",&t);    for ( int i = 1; i <= t; i++ )    {        scanf("%d",&n);        printf("%d %d %d\n",i,n,eul[n]*2+1);    }    return 0;}