エラトステネスの篩

 投稿者:しばっち  投稿日:2016年 2月28日(日)17時56分29秒
  エラトステネスの篩をC++で書いてみました
64bit整数(long long)を使用しています
sqrt関数がdouble型のため、isqrt関数を定義しています
また、ライブラリーのatoll関数がunsignedに対応していないので
atoull関数を定義しています。これにより、2^64-1まで対応できます
bool型(1byte)を使用し、余計な処理はしていません)
ASSIGN文でdllファイルのパスを指定してください

OPTION ARITHMETIC DECIMAL_HIGH
LET S=1E+19+1 !'奇数 (1E+19 < 2^64 < 1E+20))
LET E=S+1000  !'奇数
LET S$=STR$(S)
LET E$=STR$(E)
PRINT ERATOS(S$,E$)
END

EXTERNAL FUNCTION ERATOS(S$,E$)
OPTION ARITHMETIC DECIMAL_HIGH
ASSIGN ".\DLL\eratos.dll","eratos",FPU
END FUNCTION

-----------------------------------------------------------------------
                                    eratos.cpp

using namespace std;

unsigned long long isqrt(unsigned long long num)
{
unsigned long long n,n1;
if(num==0) return 0;
n=num/2+1;
n1=(n+(num/n))/2;
while (n1<n) {
    n=n1;
    n1=(n+(num/n))/2;
        }
return n;
}

unsigned long long atoull(char *str)
{
unsigned long long result = 0;
while (*str>='0' && *str<='9'){
    result=(result*10)+(*str++ - '0');
                  }
return result;
}

extern "C" __declspec(dllexport) double eratos(char *s,char *e)
{
unsigned long long n,nn,m,i,j,count=0;
bool *x;
n=atoull(s);
m=atoull(e);
if (n>=m) return -1.0;
x=new bool [m-n+2];
for(i=0;i<=m-n+1;i++) x[i]=true;
for(i=3;i<=isqrt(m);i+=2){
   if (n % i==0)  nn=n;
        else
                  nn=((n/i)+1)*i;
   if (n<=i && x[nn-n]==true) nn+=i;
for(j=nn;j<=m;j+=i) x[j-n]=false;
                   }
for(j=n;j<=m;j++)
    if (j % 2==1 && x[j-n]==true) count++;
delete [] x;
return (double)count;
}
 

Re: エラトステネスの篩

 投稿者:たろさ  投稿日:2016年 3月 7日(月)17時56分20秒
  > No.3985[元記事へ]

しばっちさんへのお返事です。


> ASSIGN文でdllファイルのパスを指定してください

パスが通らないのでルートで動作報告です。

Intel i7-4790K@4.7GHz  Win7 32bt

100億    455052511
200億    882206716
300億   1300005926
400億  1711955433
500億  2119654578
600億  2524038155
700億  2925699539
800億  3325059246
900億  3722428991
1000億 4118054813


100億から1000億 1億刻みで901回  3596.68 秒で計算しました

素数定理

100億と1000億の素数の個数を確認しました。

十進BASICプログラムの計算時間です。

http://blogs.yahoo.co.jp/donald_stinger

 

Re: エラトステネスの篩

 投稿者:たろさ  投稿日:2016年 3月 8日(火)08時56分24秒
  > No.3989[元記事へ]

しばっちさんへのお返事です。

動作報告です。

LET C=4118054813
FOR i=1E+11 TO 1E+12 STEP 1E+8

(1E+12-1E+11)/1E+8=9000回

Intel i7-4790K@4.7GHz Win7 32bt

nまで   素数の個数
1000億     4118054813個
2000億     8007105059個
3000億    11818439135個
4000億   15581005657個
5000億   19308136142個
6000億   23007501786個
7000億   26684074310個
8000億   30341383527個
9000億   33981987586個
1兆    37607912018個

計算開始時刻2016 03 07/17:23:40
計算終了時刻2016 03 08/04:51:15
11時間27分35秒で計算しました
秒換算:41255 秒

LET C=37607912018
(1E+13-1E+12)/1E+8=90000回

計算開始しました。まだ200回程度です。
何時かは終わるでしょう。

予想
4日間と 18時間35分50秒  追記3月8日(1E+13-1E+12)/1E+8=90000回 約30000回を超えた
ところで、「回文数」programをBASIC Acc で爆走中落ちました。また再起動
かなり負荷かかってます。10兆の計算結果は日曜日以降のお楽しみです。

追記
オーバークロックしてますので念のため別のテスト機
計算結果は同様でした。

100000000000 : 4118054813
1000000000000 : 37607912018
37607912018個 確認しました。

Intel G3258@4.6GHz Win7 32bt
計算開始 2016 03 07/22:51:55
計算終了 2016 03 08/11:13:51

経過時間:0日と12時間21分56秒
日換算:0.51523148148148 日
時間換算:12.365555555556 時間
分換算:741.93333333333 分
秒換算:44516 秒

計算コスト
G3258 CPは高い様です。

FOR i=1E+11 TO 1E+12 STEP 1E+9  10億stepも計算中

http://blogs.yahoo.co.jp/donald_stinger

 

Re: エラトステネスの篩

 投稿者:たろさ  投稿日:2016年 3月 9日(水)11時41分43秒
  > No.3985[元記事へ]

しばっちさんへのお返事です。

1億STEPと10億STEPの比較DATA

100000000000=1E+11から1000000000000=1E+12までの素数の個数
100000000000 : 4118054813個
1000000000000 : 37607912018個
        37607912018個 確認しました。

【テスト機】              【範囲】         【計算時間】
Intel i7-4790K@4.7GHz Win7 32bt  FOR i=1E+11 TO 1E+12 STEP 1E+8(11時間27分35秒)
Intel G3258@4.6GHz Win7 32bt      FOR i=1E+11 TO 1E+12 STEP 1E+8(12時間21分56秒)
Intel G3258@4.6GHz Win7 32bt      FOR i=1E+11 TO 1E+12 STEP 1E+9(13時間53分44秒)


     nまで  :  π(x)
100000000000 : 4118054813個 【STEP 1E+9】(13時間53分44秒)
100000000000 : 4118054813  【STEP 1E+8】(11時間27分35秒)

200000000000 : 8007105059個
200000000000 : 8007105059

300000000000 : 11818439135個
300000000000 : 11818439135

400000000000 : 15581005657個
400000000000 : 15581005657

500000000000 : 19308136142個
500000000000 : 19308136142

600000000000 : 23007501786個
600000000000 : 23007501786

700000000000 : 26684074310個
700000000000 : 26684074310

800000000000 : 30341383527個
800000000000 : 30341383527

900000000000 : 33981987586個
900000000000 : 33981987586

1000000000000 : 37607912018個
1000000000000 : 37607912018

http://blogs.yahoo.co.jp/donald_stinger

 

Re: エラトステネスの篩

 投稿者:しばっち  投稿日:2016年 7月17日(日)16時12分34秒
  > No.3985[元記事へ]

旧版 eratos.dllをOPEN MP(マルチスレッド化)にとりあえず成功しました。
#3986

VC++2015(x86)にて、ビルド(コンパイル)しています。
こちらの環境では、およそ1.4倍高速化しました。(8スレッドで実行)
但し、OPEN MP上の制約により、unsignedが使えず、2^63-1までになります。

実行にはOPEN MPランタイム(vcomp140.dll)が必要です。(basic.exeと同じフォルダに入れてください)
なお、vcomp140.dllはライセンス上再配布可能らしいので、ビルドしたdllファイルと同梱し
アップローダーにupしておきました。(eratos.zip) 下記 URLの有効期限は現時点より1ヶ月間となります。


ダウンロードパス:shibacchi


OPTION ARITHMETIC NATIVE
LET S=1E+13+1
LET E=S+1E+8
LET S$=STR$(S)
LET E$=STR$(E)
LET T=TIME
PRINT "シングルスレッド"
PRINT ERATOS(S$,E$);"個"
LET T1=TIME-T
PRINT T1;"秒"
LET T=TIME
PRINT "マルチスレッド"
PRINT ERATOS2(S$,E$,0);"個"
LET T2=TIME-T
PRINT T2;"秒"
PRINT T1/T2;"倍"
END

EXTERNAL FUNCTION ERATOS(S$,E$)
OPTION ARITHMETIC NATIVE
ASSIGN "eratos.dll","eratos",FPU
END FUNCTION

EXTERNAL FUNCTION ERATOS2(S$,E$,THREADS) !'THREADS=実行スレッド数(0の時 自動設定)
OPTION ARITHMETIC NATIVE
ASSIGN "eratos_parallel.dll","eratos",FPU
END FUNCTION

-------------------------------------------------------------------------------------------
                                          eratos_parallel.cpp

#include <cstdlib>
#include <omp.h>
using namespace std;

long long isqrt(long long num)
{
long long n,n1;
if(num==0) return 0;
n=num/2+1;
n1=(n+(num/n))/2;
while (n1<n) {
    n=n1;
    n1=(n+(num/n))/2;
        }
return n;
}

extern "C" __declspec(dllexport) double eratos(char *s,char *e,int threads)
{
long long n,nn,m,i,j,count=0;
bool *x;
n=atoll(s);
m=atoll(e);
if (n>=m || n<1) return -1.0;
if (n==1 || n==2) count++;
if (n==1) n+=2;
if (n % 2==0) n++;
x=new bool [m-n+2];
if (x==NULL) return -999.0;
if (threads>0) omp_set_num_threads(threads);
#pragma omp parallel for
for(i=0;i<=m-n+1;i++) x[i]=true;
#pragma omp parallel for private(j,i,nn) schedule(dynamic)
for(i=3;i<=isqrt(m);i+=2){
   if (n % i==0)  nn=n;
        else
                  nn=((n/i)+1)*i;
   if (n<=i && x[nn-n]==true) nn+=i;
for(j=nn;j<=m;j+=i) x[j-n]=false;
                              }
#pragma omp parallel for
for(j=n;j<=m;j+=2)
    if (x[j-n]==true) {
#pragma omp atomic
         count++;
                 }
delete [] x;
return (double)count;
}
 

戻る