オイラーのゼータ関数

 投稿者：たろさ  投稿日：2022年 1月16日(日)21時36分47秒
  ! 6n+1,6n-1 篩　Ver.2
!5億までの素数を配列DIMに登録する。
!オイラーのゼータ関数
OPTION ARITHMETIC NATIVE       !2進モード
LET t0=TIME
LET Pk=26355867
LET k=5E8    !1E7  9999991 素数 (664579th)
LET k1=IP(k/6)
LET k2=IP(SQR(k1))

DIM A5(k1),A7(k1)

MAT A5=ZER
MAT A7=ZER

LET P1=5
LET C1=1

DO
   FOR n=1 TO k1
      LET P6=P1*n+C1
      IF P6 > k1 THEN EXIT FOR
      LET A5(P6)=1
   NEXT n
   FOR n=1 TO k1
      LET P6=P1*n-C1
      IF P6 > k1 THEN EXIT FOR
      LET A7(P6)=1
   NEXT n

   LET P1=P1+6
   LET C1=C1+1
   IF P1 >k1 THEN EXIT DO
LOOP


LET P1=7
LET C1=1
DO
   FOR n=1 TO k1
      LET P6=P1*n+C1
      IF P6 > k1 THEN EXIT FOR
      LET A7(P6)=1
   NEXT n

   FOR n=1 TO k1
      LET P6=P1*n-C1
      IF P6 > k1 THEN EXIT FOR
      LET A5(P6)=1
   NEXT n

   LET P1=P1+6
   LET C1=C1+1
   IF P1 >k1 THEN EXIT DO
LOOP
DIM B(Pk)       !1E8 99999989 素数 (5761455th)
LET B(1)=2      !2E8 199999991 素数 (11078937th)
LET B(2)=3
LET C=3
FOR n=1 TO k1
   LET P5=6*n-1
   LET P7=6*n+1
   IF A5(n)=0 THEN
      LET B(C)=P5
      LET C=C+1
   END IF
   IF A7(n)=0 THEN
      LET B(C)=P7
      LET C=C+1
   END IF
NEXT n

LET Q=1
FOR n=1 TO Pk
   LET P=B(n)^2
   LET Q=Q*(P/(P-1))
NEXT n

PRINT Q
PRINT PI^2/6

LET TM=TIME-t0
PRINT TM;"秒"

END

動作確認　Microsoft　surface　pro6  メモリーDDR3 8GB   Win11 Pro 64-bit
Paract BASIC Ver. 2.1.3.1 (2022.01.09)
64ビット版Lazarus3.2.2

計算結果
1.644934066247
1.644934066848
9.474000000002 秒

5億まで　では?
15億までは確認しました。1.644934066

とある八雲の科学解説 『素数から生まれる円周率』
https://www.youtube.com/watch?v=aq_-cIv92t8
 

戻る