一道google面试题--自然数e中出现的连续的第一个10个数字组成的质数

4 篇文章 0 订阅
4 篇文章 0 订阅


博客内容移到 http://www.linuxyu.com/

此CSDN博客将不再更新,欢迎大家访问新的网站~~


Google早几年在美国很多地铁的出站口都有大幅招聘广告,它的第一题如下了:{first 10-digit prime found in consecutive digits e}.com,也就是自然数e中出现的连续的第一个10个数字组成的质数。据说当时只要正确解答了这道题,在浏览器的地址栏中输入这个答案,就可以进入下一轮的测试,整个测试过程如同一个数学迷宫,直到你成为google的一员。

对这题比较感兴趣,就暂时开题,这两天动笔完成这篇博客。

难点分析

初步分析,主要两个难点:

(1)e的小数点后n位有效数字,提供素数判断的数据源;

(2)素数的高效判断,尤其时针对10位数的指数判断。

自然数e计算

计算e可使用泰勒公式推出来的

e^x1 + x + x^2/2! + x^3/3! +……+ x^n/n!  1

当x=1时,

e1 + 1 + 1/2! + 1/3! +……+ 1/n!   2

取n=10,即可算出近似值e≈2.7182818。 当n取较大值时,可获得高精度的e。

版本1

最初步的想法,直接利用泰勒公式计算,printf出来。

//compute e 
//program 1

#include <stdio.h>
#define N 20

int main()
{
    double e = 1.0, t = 1.0;
    int i = 0;

    for(i = 1; i <= N; i++){
        t = t / i;
        e = e + t;
    }

    printf ("e=%.20lf\n",e );
    return 0;
}

版本2

对(2)进行叠乘合并起来,改写成如下(3):

e=2+1/2(1+1/3(1+1/4(1+1/(n-1)(1+1/n))))   3

从最里面的括号往外算,共做n次除法和加法得一段结果,运算效率也是O(N*M),但是由于收敛速度快些。

//compute e 
//program 2

#include <stdio.h>
#define N 20

int main()
{
    double e = 0.0;
    int i = 0;

    for(i = N; i > 0; i--){
        e = (e + 1.0) / i;
    }

    printf ("e=%.20lf\n",e );
    return 0;
}

版本3

版本1和2只能达到小数点后15位精度,假设我们需要的数据源是e达到e后面1000位数,那么就必须要考虑到高精度的问题。自定义数据结构时必须的,但(3)中主要计算加法、减法和乘法不会丢失精度,而除法会丢失精度,因为它始终存在余数,余数是丢失精度的原因,所以要保存精度就是保存余数,在这里的除法将产生两个整数,一个是商,一个是余数,形如:

a / b = (d,r)

如果以10为基,那高精度就是对这样的除法继续做下去,将r*10恢复精度,再除,直到满足精度为止。

//compute e 
//program 3

#include <stdio.h>

//umber of significant digit=(DN-1)*4
#define DN 2504

int euler[DN], data[DN];
int step,n;

void precise_division()
{
    int i=0;
    long y=0,r=0;

    for(i=step;i<DN;i++)
    {
        y = 10000*r + data[i];
        data[i] = y/n;
        euler[i] += data[i];
        r = y%n;
    }
}

void revise()
{
    int c=0,i=0;

    for(i=DN-1;i>=0;i--)
    {
        euler[i] += c;
        if(euler[i]>9999)
        {
            c = euler[i]/10000;
            euler[i] = euler[i]%10000;
        }
        else
            c = 0;
    }
}

void euler_print()
{
    int i=0;

    printf("\n\nE=%d.\n",euler[0]);
    for(i=1;i<DN;i++)
    {
        printf("%.4d ",euler[i]);
        if((i%250)==0)
            printf("\n\tnow it has %d bit\n",i*4);
    }
    printf("\n\neulur number is ok\n");
}

void main()
{
    int i=0;
    step=0;
    n=1;

    for(i=0;i<DN;i++)
    {
        euler[i]=0;
        data[i]=0;
    }
    euler[0]=1;
    data[0]=1;

    while(1)
    {
        i=step;
        while(data[i]==0)
        {
            i++;
            if(i>=DN)
            goto _EXT;
        }

        step=i;
        precise_division();
        revise();
        n++;
    }

_EXT:
    euler_print();
}

版本4

参考文献2,实现本版本,唯一没有想通的部分在于GetLength函数实现的数学原理,这个后续再补充。

//compute e 
//program 4

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

//umber of significant digit
#define N 10000

typedef unsigned int UINT;

int GetLength(int m)
{
    int n=1;
    double x=2.0/2.72;
    while(m)
    {
        while(x<1.0)
        {
            n++;
            x*=(double)(n+1);
        }
        while(x>=1.0 && m)
        {
            x/=10.0;
            --m;
        }
    }
    while(x<1.0)
    {
        n++;
        x*=(double)(n+1);
    }
    return n;
}

void main()
{
    const UINT base=1000000;
    int i,j,k,m;
    UINT *r=NULL;
    UINT *euler=NULL;
    UINT y=0;

    m=GetLength(N);

    r=(UINT *)calloc(m+1,sizeof(UINT));
    euler=(UINT *)calloc(N,sizeof(UINT));

    assert(r!=NULL);
    assert(euler!=NULL);

    for(i=0;i<=m;++i)
    {
        r[i]=1;
        euler[i]=0;
    }

    j=1;
    euler[0] = 2;

    for(k=N;k>0;k-=6)
    {
        y=0;
        for(i=m;i>=2;--i)
        {
            y = y + r[i]*base;
            r[i] = y%i;
            y /= i;
        }

        if(k<6)
        {
            euler[j++] = y%base;
        }
        else
        {
            if(y<base)
                euler[j++] = y;
            else
            {
                if(r)
                    free(r);
                return;
            }
        }
    }
    if(r)
        free(r);

    printf("\n\nE=%d.\n",euler[0]);
    for(i=1;i<j;i++)
    {
        printf("%.6d ",euler[i]);
        if((i%100)==0)
            printf("\n\tnow it has %d bit\n",i*6);
    }
    printf("\n\neulur number is ok\n");

    if(euler)
        free(euler);

    return;
}

素数判断

版本1

素数,即除1以外,只能被1和其本身整除的数。因此根据定义,最初的判断素数的方法为版本1。

//judge prime 
//program 1
bool prime_1(int n)
{
    int i;
    for(i=2;i<n;i++)
    {
        if(n%i==0)
        {
            return false;
        }
    }
    return true;
}

版本2

由于一个数的因式分解可表示为x = a * b, a<=b,因此在判断整数n是不是素数时并不需要整数从2 ~ n,只需要整除2 ~ sqrt(n) 即可。

//judge prime 
//program 2
bool prime_2(int n)
{
    int i,n2;

    n2=sqrt(n);
    for(i=2;i<=n2;i++)
    {
        if(n%i==0)
        {
            return false;
        }
    }
    return true;
}

但是需要注意的是,int类型只能表示2^-31 ~ 2^31 范围,并不能完整表示10位大小的整数n,因此需要使用能够long long int类型。其次,求平方根的函数原型为double sqrt(double),进行sqrt(n)计算时并不确定能保证其精度。最后,对e小数点后,每个10位数字大小的整数n进行判断,需要对每个2 ~ sqrt(n)的整数进行判断,单个判断耗时间太久,所有判断重复运算太多。因此次方法必然不行。

版本3

任何一个合数都可以表现为适当个素数的乘积的形式,所以我们只用素数去除要判断的数即可,比如要判断100以内的素数,只用10以内的2,3,5,7就够了。那么首先就要构建一个素数表,以空间换时间。

可以使用一种被称为“埃拉托色尼筛”的算法。该算法一开始初始化一个2~n的连续整数序列。从2开始,2的倍数肯定不是素数,去除,剩下的数中下一个是3。再把3的倍数去除,再下一个是5(4已经去除了),以此类推,最后剩下的数必然是素数,因为它没有被筛,不是任何一个数的倍数(除了1和自己)。这样只需要很多次加法就可以了。

//create prime less then n 
//program 3

#include <stdio.h>
#include <stdlib.h>
#include <math.h>  //Linux, gcc -lm
#include <stdbool.h>


#define Long64 long long
#define MAX_N 8000000

Long64 prime_table[MAX_N/2]={0};

void create_prime_table(Long64 n)
{
    Long64 i,j,n2;
    bool prime[MAX_N+1]={0};

    n2=sqrt(n);

    i=2;
    for(i=2;i<=n2;i++)
    {
        if(prime[i]==0)
        {
            for(j=i+i;j<=n;j+=i)
            {
                prime[j]=1;
            }
        }
    }

    j=0;
    for(i=2;i<=n;i++)
    {
        if(prime[i]==0)
        {
            prime_table[j]=i;
            j++;
            printf("%lld : %lld \n",j,i);
        }
    }
}

void main()
{
    create_prime_table(MAX_N);
}

这算是基本的筛法了,但它还有很多优化的空间。首先内存中存的表不只是素数,还有合数的空间,合数比素数多得多。 其次对于合数可能会被筛选多次,比如6,在2的时候筛掉一次。在3的时候又筛掉一次。

版本4

采用线性筛选方式,即对一个合数,只进行一次筛选,极大程度的减少重复计算。 同时,采用堆结构存放素数表。

//create prime less then n
//program 4

#include <stdio.h>
#include <stdlib.h>
#include <math.h>  //Linux, gcc -lm
#include <stdbool.h>


#define Long64 long long
#define MAX_N 1000000000

void create_prime_table2(Long64 n)
{
    Long64 i,j,k;
    bool *prime;
    Long64 *prime_table;

    prime=(bool *)calloc(MAX_N+1,sizeof(bool));
    prime_table=(Long64 *)calloc(MAX_N/100,sizeof(Long64));

    i=2;
    k=0;
    for(i=2;i<=n;i++)
    {
        if(prime[i]==0)
        {
            prime_table[k++]=i;
            //printf("%lld : %lld \n",k,i);
        }
        for(j=0;j<k && prime_table[j]*i <= n;j++)
        {
            prime[prime_table[j]*i]=1;
            if(i%prime_table[j]==0)
                break;
        }
    }
    printf("k = %lld\n",k);

    if(prime_table)
        free(prime_table);
    if(prime)
        free(prime);
}

void main()
{
    create_prime_table2(MAX_N);
}

这样做之后,程序运行速度快多了,但是只能计算到1x10^9....再大时,运行就提示段错误了。还要再在其他地方想办法。。。 暂时先记录到这里。

自然数e

上上一篇博客中,虽然实现了计算n为自然数e的程序,但是GetLength 函数实现的数学原理一直没有想清楚。其实是这样的。 计算自然数e时,我们使用的时泰勒公式

e^x = 1 + x + x^2/2! +... + x^n/n!...

当考虑到泰勒公式的余项时,余项可以表示为

R(x) = (e^esp * x^(n+1))/(n+1)!

其中esp在0~x之间。再取x=1,得到

R(1) = e^esp/(n+1)! <= e/(n+1)! < 10^(-M)

10^(-M)就是我们要达到的误差限,也就是小数点后的后M位。 因此,若要计算到10000位的精度,只需计算(n+1)!/e > 10^10000即可。 程序中GetLength的实现原理就是这样。

素数判断

上一篇博客中,只计算到1x10^9内的素数,并没有达到要求,主要时卡在了数据结构上,再大的话在我2G内存32位系统上就已经提示段错误了。现在回过头了想想,其实还可以在算法上再改进。 我们知道a*b = c,一个M位的数乘上另外一个M位的数,得到结果位数范围在M+1~2M的范围内。我们只要求出5位数以内的素数,就可以用来判断10位数的整数是否为素数了。 那么素数判断的程序就可以修改成这样,也就是版本2和版本4的结合。

版本5

#include <stdio.h>
#include <stdlib.h>
#include <math.h>  //Linux, gcc -lm
#include <stdbool.h>

#define Long64 long long
#define MAX_N 100000

//create prime less then n 
//program 5

Long64 create_prime_table3(bool *prime, Long64 *prime_table, Long64 n)
{
    Long64 i,j,k;

    i=2;
    k=0;
    for(i=2;i<=n;i++)
    {
        if(prime[i]==0)
        {
            prime_table[k++]=i;
            //printf("%lld : %lld \n",k,i);
        }
        for(j=0;j<k && prime_table[j]*i <= n;j++)
        {
            prime[prime_table[j]*i]=1;
            if(i%prime_table[j]==0)
                break;
        }
    }

    return k;
}

bool prime3(Long64 prime_table[], Long64 length, Long64 n)
{
    Long64 i;

    for(i=0;i<length;i++)
    {
        if(n%(prime_table[i])==0)
        {
            return false;
        }
    }
    return true;
}

void main()
{   
    bool *prime;
    Long64 *prime_table,i=0,j=0,length=0;

    prime=(bool *)calloc(MAX_N+1,sizeof(bool));
    prime_table=(Long64 *)calloc(MAX_N,sizeof(Long64));

    length = create_prime_table3(prime,prime_table,MAX_N);

    for(i=1000000000; i<10000000000; i++)
        if(prime3(prime_table,length,i)!=false)
            printf("%lld : %lld\n",j++,i);
}

收尾:判断e中首个连续10位素数

有了以上问题解决后,这道题终于可以解决了。不仅可以找到第一个10位长度的素数,还可以找到前10个!前10个结果如下,格式为“第几个:在e中的起始位置:要找的素数”

1:99:7427466391 
2:123:7413596629 
3:149:6059563073 
4:171:3490763233 
5:182:2988075319 
6:201:1573834187 
7:214:7021540891 
8:218:5408914993 
9:254:6480016847 
10:295:9920695517

完整实现代码如下

#include <stdio.h>
#include <stdlib.h>
#include <math.h>  //Linux, gcc -lm
#include <stdbool.h>
#include <assert.h>

#define MAX_E_BITE_N 1000
#define PRIME_N 1000000

#define Long64 long long

typedef unsigned int UINT;


int get_length(int nbite)
{
    int m = 1;
    double x = 2.0/2.72;
    while(nbite)
    {
        while(x<1.0)
        {
            m++;
            x *= (double)(m+1);
        }
        while(x >= 1.0 && nbite)
        {
            x /= 10.0;
            --nbite;
        }
    }
    while(x < 1.0)
    {
        m++;
        x *= (double)(m+1);
    }
    return m;
}


int create_e(Long64 *euler, const UINT base, int m)
{
    UINT *r = NULL;
    UINT y = 0;
    int i = 0,j = 0,k = 0;

    r=(UINT *)calloc(m+1,sizeof(UINT));

    assert(r != NULL);
    assert(euler != NULL);

    for(i = 0;i <= m; ++i)
    {
        r[i] = 1;
        euler[i] = 0;
    }

    j = 1;
    euler[0] = 2;

    for(k = MAX_E_BITE_N; k > 0; k -= 5)
    {
        y = 0;
        for(i = m; i >= 2; --i)
        {
            y = y + r[i]*base;
            r[i] = y%i;
            y /= i;
        }

        if(k < 5)
        {
            euler[j++] = y%base;
        }
        else
        {
            if(y < base)
                euler[j++] = y;
            else
            {
                if(r)
                    free(r);
                return 0;
            }
        }
    }
    if(r)
        free(r);

    printf("\n\nE=%lld.\n",euler[0]);
    for(i = 1; i < j; i++)
    {
        printf("%.5lld ",euler[i]);
    }
    printf("\n\neulur number is ok j=%d\n",j);

    return j;
}

Long64 create_prime_table(bool *prime, Long64 *prime_table, Long64 n)
{
    Long64 i,j,k;

    i = 2;
    k = 0;
    for(i = 2;i <= n;i++)
    {
        if(prime[i] == 0)
        {
            prime_table[k++] = i;
//            printf("%lld : %lld \n",k,i);
        }
        for(j = 0; j < k && (prime_table[j]*i) <= n; j++)
        {
            prime[prime_table[j] * i] = 1;
            if(i%prime_table[j] == 0)
                break;
        }
    }

    return k;
}

bool is_prime(Long64 prime_table[], Long64 length, Long64 n)
{
    Long64 i;

    for(i = 0; i < length; i++)
    {
        if(n%(prime_table[i]) == 0)
        {
            return false;
        }
    }
    return true;
}


void main()
{
    const UINT base = 100000;
    int m = 0,elength = 0,j = 0,find_number = 10;

    Long64 *euler = NULL;
    bool *prime = NULL;
    Long64 *prime_table = NULL;

    Long64 i = 0,length = 0,data = 0;

    Long64 ten[10] = {1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000};

    m = get_length(MAX_E_BITE_N);

    euler = (Long64 *)calloc(MAX_E_BITE_N,sizeof(Long64));

    elength = create_e(euler,base,m);


    prime = (bool *)calloc(PRIME_N+1,sizeof(bool));
    prime_table = (Long64 *)calloc(PRIME_N,sizeof(Long64));

    length = create_prime_table(prime,prime_table,PRIME_N);

    if(prime)
        free(prime);

    m = 0;
    for(i = 1; (i+2) < elength; i++)
    {
        for(j = 0; j < 5; j++)
        {
            data = euler[i]%ten[5-j]*ten[5+j] + euler[i+1] * ten[j] + euler[i+2]/ten[5-j];

//          printf("data : %lld \n",data);

            if(is_prime(prime_table,length,data) == true)
            {
                printf("find it! %d:%lld:%lld\n",++m,5*(i-1)+j+1,data);

                if(m == find_number)
                    break;    
            }               
        }

        if(m == find_number)
            break;    
    }


    if(euler)
        free(euler);
    if(prime_table)
        free(prime_table);

    return;
}

评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值