数学研发论坛

 找回密码
 欢迎注册
楼主: 无心人

[讨论] 10秒内计算出1亿素数表, 不输出

[复制链接]
发表于 2008-7-5 10:19:20 | 显示全部楼层
无心人#52的的程序,我改写了一下,基本上是一样的。gcc 编译1.6双核CPU 运行时间是>1000秒。

#include <stdio.h>
#include <math.h>
#include <time.h>

#define C17       510510
#define SIZE      255255

char BaseCache[SIZE],WorkCache[SIZE];
int PRIMES [43000];
int LOWSQRT, HIGHSQRT; //起始位置和结束位置的平方根


void initCache(void)
{
        int p,i;
        for ( i = 0; i < SIZE ; i++ )
                BaseCache[i] = 1;
        for ( p = 3; p <= 17; p+=2 )
        {
                if ( p == 9 || p == 15 )
                        continue;
                for ( i = p >> 1; i < SIZE; i+=p )
                        BaseCache[i] = 0;
        }
}

void firstSieve(void)
{
    int p,i,j,delta = 36,cigema = 180;
        for ( i = 0; i < SIZE; i++ )
                WorkCache[i] = BaseCache[i];
    LOWSQRT = (int)sqrt(C17-1);
        for ( p = 19; p <= LOWSQRT; p+=2 )
    {
                if (WorkCache[p >> 1])
                        for ( j = cigema; j < SIZE; j+=p )
                                WorkCache[j] = 0;
                delta += 4;
                cigema += delta;
        }
    WorkCache[1] = 1;
    WorkCache[2] = 1;
        WorkCache[3] = 1;
        WorkCache[5] = 1;
        WorkCache[6] = 1;
        WorkCache[8] = 1;
        j = 1; PRIMES[j] = 2;
        for (i = 1; i < SIZE; i++)
                if (WorkCache[i])
                PRIMES[++j] = (i << 1) + 1;
         PRIMES[0] = j;
}

void __fastcall nextSieve(int beg, int end)
{
        int p,i,j,k;
        HIGHSQRT = (int)sqrt(end);
        for ( i = 0; i < SIZE; i++ )
                WorkCache[i] = BaseCache[i];
        i = 8; p = 19;
        while ( p <= LOWSQRT )
        {
                j = beg % p;
                if (j != 0)
                {
                        j = beg + p - j;
                        if (j % 2 == 0) j+=p;
                } else j = beg;
                j = (j - beg + 1) >> 1;
                for ( k = j; k < SIZE; k+=p )
                        WorkCache[k] = 0;
                p = PRIMES[++i];
        }
        while ( p <= HIGHSQRT )
        {
                j = (p * p - beg + 1) >> 1;
                for ( k = j; k < SIZE; k+=p )
                        WorkCache[k] = 0;
            p = PRIMES[++i];
        }
        LOWSQRT = HIGHSQRT;
}

int __fastcall PrimeCount(int len)
{
        int c = 0,i;
        for ( i = 0; i < len; i++ )
                if (WorkCache[i]) c++;
    return c;
}

int main()
{
        int c,i,j,k;
        clock_t t1 = clock(),t2;

        initCache();
        firstSieve();
        c = PRIMES[0];
        int beg = C17 + 1,end = (C17 << 1) - 1;
        while ( end < 100000000 )
        {
                nextSieve(beg,end);
                c += PrimeCount(SIZE);
            beg += C17; end += C17;
        }
    end = 99999999;
        nextSieve(beg,end);
        c += PrimeCount((end - beg)>>1);
        t2 = clock();
        printf("--> time : %d\n",t2-t1);
        printf("Prime Count: %d\n",c);
        return 0;
}

/********************************************************************************************************
--> time : 1023
Prime Count: 5761455
********************************************************************************************************/
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-5 10:21:16 | 显示全部楼层
//为什么您的程序只有500多毫秒,而我的却要1000多毫秒?

E:\worker>gcc -o shai shai.c -O3

E:\worker>shai
--> time : 390
Prime Count: 5761455

不好意思,刚才忘了打-O3了哈哈!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-7-5 11:03:46 | 显示全部楼层
呵呵

我的估计输在内存访问
用位数组,内存访问快
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-5 11:35:25 | 显示全部楼层
同样是普通筛法,我的程序太慢了,没法和P1比阿。liangbch能发下你的p1源吗可以么?

62#的附件包含了多个算法的源代码,你可以下载下来看
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-7-5 20:38:40 | 显示全部楼层
liang什么时候回来阿
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-7 10:13:50 | 显示全部楼层
按照这边给我定的日程表,八月底才能回来,时间太长了。

[ 本帖最后由 liangbch 于 2008-7-7 12:39 编辑 ]
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-7-7 10:41:31 | 显示全部楼层


少了一份交流
呵呵
希望早点回来
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-7 11:27:53 | 显示全部楼层
原帖由 liangbch 于 2008-7-7 10:13 发表
按照这边给我定的日程表,八月底能回来,时间太长了。


应为:可能回不来
还是:才能回得来

确实:时间太长了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-7 14:37:57 | 显示全部楼层
看看附件这个ecprime的威力吧
100亿的筛法只要4s(AMD 3600+ 单线程)
还能计算素数间距, K生素数等(这个不如我的程序)
有时间我将它改成多进程的,
四核机器百亿只要1s
以下是作者的help

23.10.2003, e0325944@student.tuwien.ac.at

        help file, ecprime 1.4 by Kim Walisch

#--------------------------------------------------
version 1.4, 2003

ecprime [nr] [-n[nr]] [-s[s;u]] [-r[nr]] [-g[nr]] [-t[k]] [-p[k]] [-c] [-d]

#--------------------------------------------------

Usage:

# Count Functions

ecprime 1000
counts primes (0 - 1000)

ecprime 1000 -n100
counts primes (100 - 1000)

ecprime 1000 -n100 -t2
counts 2-tuplets (100 - 1000)

i.e -t3,-t4,-t5,-t6,-t7

# Print Functions

ecprime 1000 -p1
prints primes (0 - 1000) to primes.txt

ecprime 1000 -p2
prints 2-tuplets (0 - 1000) to 2-tuplets.txt

i.e -p3,-p4,-p5,-p6,-p7

# Settings

ecprime 1000 -s3
uses a sieve of 2^3KB, -s(0 - 24)

max speed for sieve = L1-cache

ecprime 1000 -s7;4
uses a sieve of 2^7 KB, divided into subsieves of 2^4 KB, fast for numbers > 10^12

ecprime 1000 -r17
uses a repetion field of 17#/30 bytes, -r(7 - 23)

see limits file for details

# Detailed Time Output, -d

time1: initialization repetition field
time2: initialization small primes

time3: sieving time (std)

# Prime Gaps

ecprime 1000 -g
prints first occurrence gaps to gaps.txt

ecprime 1000 -g100
prints first occurrence gaps (>100) to gaps.txt

ecprime 1000 -g100 -c
prints all prime gaps (>100) to gaps.txt

ecprime.zip

172.28 KB, 下载次数: 30, 下载积分: 金币 -1 枚, 经验 1 点, 下载 1 次

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2009-1-1 22:38:35 | 显示全部楼层

收藏

一下理解不了这么多,
以后慢慢来看,
以前自己写的筛法,用的是24的模的,
其他部分就跟死算一样,
慢的要死,
筛法确实是很讲究的东东,
期待新颖的素数判断法出现。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2020-12-5 19:05 , Processed in 0.065514 second(s), 17 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表