素数個数関数

 投稿者：しばっち  投稿日：2015年10月 1日(木)23時29分30秒
  https://ja.wikipedia.org/wiki/リーマンの素数公式

OPTION ARITHMETIC COMPLEX
PUBLIC NUMERIC P(100)
FOR I=1 TO 100
   READ K
NEXT I
DATA 14.134725142
DATA 21.022039639
DATA 25.010857580
DATA 30.424876126
DATA 32.935061588
DATA 37.586178159
DATA 40.918719012
DATA 43.327073281
DATA 48.005150881
DATA 49.773832478
DATA 52.970321478
DATA 56.446247697
DATA 59.347044003
DATA 60.831778525
DATA 65.112544048
DATA 67.079810529
DATA 69.546401711
DATA 72.067157674
DATA 75.704690699
DATA 77.144840069
DATA 79.337375020
DATA 82.910380854
DATA 84.735492981
DATA 87.425274613
DATA 88.809111208
DATA 92.491899271
DATA 94.651344041
DATA 95.870634228
DATA 98.831194218
DATA 101.317851006
DATA 103.725538040
DATA 105.446623052
DATA 107.168611184
DATA 111.029535543
DATA 111.874659177
DATA 114.320220915
DATA 116.226680321
DATA 118.790782866
DATA 121.370125002
DATA 122.946829294
DATA 124.256818554
DATA 127.516683880
DATA 129.578704200
DATA 131.087688531
DATA 133.497737203
DATA 134.756509753
DATA 138.116042055
DATA 139.736208952
DATA 141.123707404
DATA 143.111845808
DATA 146.000982487
DATA 147.422765343
DATA 150.053520421
DATA 150.925257612
DATA 153.024693811
DATA 156.112909294
DATA 157.597591818
DATA 158.849988171
DATA 161.188964138
DATA 163.030709687
DATA 165.537069188
DATA 167.184439978
DATA 169.094515416
DATA 169.911976479
DATA 173.411536520
DATA 174.754191523
DATA 176.441434298
DATA 178.377407776
DATA 179.916484020
DATA 182.207078484
DATA 184.874467848
DATA 185.598783678
DATA 187.228922584
DATA 189.416158656
DATA 192.026656361
DATA 193.079726604
DATA 195.265396680
DATA 196.876481841
DATA 198.015309676
DATA 201.264751944
DATA 202.493594514
DATA 204.189671803
DATA 205.394697202
DATA 207.906258888
DATA 209.576509717
DATA 211.690862595
DATA 213.347919360
DATA 214.547044783
DATA 216.169538508
DATA 219.067596349
DATA 220.714918839
DATA 221.430705555
DATA 224.007000255
DATA 224.983324670
DATA 227.421444280
DATA 229.337413306
DATA 231.250188700
DATA 231.987235253
DATA 233.693404179
DATA 236.524229666
FOR X=2 TO 100
   PRINT X;":";π(X)
NEXT X
END

EXTERNAL  FUNCTION π(X) !'素数個数関数
OPTION ARITHMETIC COMPLEX
LET S=R(X)-T(X)
!'FOR M=1 TO INT(LOG2(X)) !'第3項
!'   IF U(M)<>0 THEN
!'      LET A=X^(1/M)
!'      RESTORE
!'      LET SS=0
!'      FOR I=1 TO 7 !'ガウス・ラゲール求積
!'         READ TT,W
!'         LET SS=SS+F2(TT)*EXP(TT-A)*W
!'      NEXT I
!'      LET S=S+U(M)/M*(SS-LOG(2))
!'   END IF
!'NEXT M
LET π=S
DATA     .1930436765603624138382479,4.0931895170127390213043288E-01
DATA    1.0266648953391919503451994,4.2183127786171977992928101E-01
DATA    2.5678767449507462069077862,1.4712634865750527839537418E-01
DATA    4.9003530845264845681017144,2.0633514468716939865705615E-02
DATA    8.1821534445628607910818276,1.0740101432807455221319596E-03
DATA   12.7341802917978137580126425,1.5865464348564201268732622E-05
DATA   19.3957278622625403117125821,3.1703154789955805622713222E-08
END FUNCTION

EXTERNAL  FUNCTION F2(T)
OPTION ARITHMETIC COMPLEX
LET F2=1/T/(T*T-1)/LOG(T)
END FUNCTION

EXTERNAL  FUNCTION R(X) !'リーマン関数(主要項)
OPTION ARITHMETIC COMPLEX
FOR M=1 TO INT(LOG2(X))
   IF U(M)<>0 THEN LET S=S+U(M)/M*LI(X^(1/M))
NEXT M
LET R=S
END FUNCTION

EXTERNAL  FUNCTION EI(Z) !'指数積分
OPTION ARITHMETIC COMPLEX
LET EULER=.57721566490153286
LET S=LOG(Z)+EULER
LET A=1
FOR K=1 TO 1000
   LET A=A*Z/K
   LET S=S+A/K
   IF ABS(A)<1E-15 THEN EXIT FOR
NEXT K
LET EI=S
END FUNCTION

EXTERNAL  FUNCTION LI(Z) !'対数積分
OPTION ARITHMETIC COMPLEX
LET LI=EI(LOG(Z))
END FUNCTION

!'EXTERNAL  FUNCTION T(X) !'振動項(うまく収束しない)
!'OPTION ARITHMETIC COMPLEX
!'FOR M=1 TO INT(LOG2(X))
!'   IF U(M)<>0 THEN
!'      FOR K=1 TO 100
!'         LET PP=P(K)/M
!'         LET S=S+U(M)/M*(LI(X^PP)+LI(X^CONJ(PP)))
!'      NEXT K
!'   END IF
!'NEXT M
!'LET T=S
!'END FUNCTION

EXTERNAL  FUNCTION T(X) !'振動項
OPTION ARITHMETIC COMPLEX
FOR M=1 TO INT(LOG2(X))
   FOR K=1 TO 100
      IF U(M)<>0 THEN LET S=S+U(M)/M*2*SQR(X)/ABS(P(K))/LOG(X)*COS(IM(P(K))*LOG(X)-ANGLE(RE(P(K)),IM(P(K))))
   NEXT K
NEXT M
LET T=S
END FUNCTION

EXTERNAL  FUNCTION U(N) !'メビウス関数
OPTION ARITHMETIC COMPLEX
FOR K=1 TO N
   IF GCD(K,N)=1 THEN
      LET S=S+COS(2*PI*K/N)
   END IF
NEXT K
LET U=INT(S+.1)
END FUNCTION

EXTERNAL FUNCTION GCD(M,N)
OPTION ARITHMETIC COMPLEX
DO WHILE N<>0
   LET TT=MOD(M,N)
   LET M=N
   LET N=TT
LOOP
LET GCD=M
END FUNCTION

----------------------------------------------------------------------------------

OPTION ARITHMETIC RATIONAL
FOR J=2 TO 100
   PRINT J;":";π(J)
NEXT J
END

EXTERNAL  FUNCTION π(M) !'素数個数関数
OPTION ARITHMETIC RATIONAL
LET S=0
FOR J=2 TO M
   LET S=S+F(J)
NEXT J
LET π=S
END FUNCTION

EXTERNAL  FUNCTION F(J)
OPTION ARITHMETIC RATIONAL
LET F=INT(((FACT(J-1)+1)/J)-INT(FACT(J-1)/J))
END FUNCTION

----------------------------------------------------------------------------------

OPTION ARITHMETIC NATIVE
FOR I=2 TO 100
   PRINT I;":";PRIMECOUNT(I)
NEXT I
END

EXTERNAL  FUNCTION PRIMECOUNT(N) !'素数個数関数
OPTION ARITHMETIC NATIVE
FOR I=2 TO N
   LET FL=0
   FOR J=2 TO INT(SQR(I))
      IF MOD(I,J)=0 THEN
         LET FL=1
         EXIT FOR
      END IF
   NEXT J
   IF FL=0 THEN LET C=C+1
NEXT I
LET PRIMECOUNT=C
END FUNCTION
 		 

Re: 素数個数関数

 投稿者：たろさ  投稿日：2015年10月 2日(金)01時44分34秒
  > No.3866[元記事へ]

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

今、私が一番見たいプログラムです。感激しました。書店で売られている本よりも

興奮しました。

私の個人的予想では、プログラムで書きました。

FOR n=1 TO 2
   PRINT n;
   PRINT n-1
NEXT n

FOR n=3 TO 48
   PRINT n;
   LET  p=MOD(n,2)/2+MOD(n,3)/3-MOD(n,2*3)/(2*3)
   PRINT p
NEXT n

FOR n=49 TO 120
   PRINT n;
   LET  p=MOD(n,2)/2+MOD(n,3)/3+MOD(n,5)/5+MOD(n,7)/7-MOD(n,6)/6-MOD(n,10)/10-MOD(n,14)/14-MOD(n,15)/15-MOD(n,21)/21-MOD(n,35)/35+MOD(n,30)/30+MOD(n,42)/42+MOD(n,70)/70+MOD(n,105)/105-MOD(n,210)/210
   PRINT p
NEXT n

FOR n=121 TO 168
   PRINT n;
   LET  p=MOD(n,2)/2+MOD(n,3)/3+MOD(n,5)/5+MOD(n,7)/7+MOD(n,11)/11-MOD(n,6)/6-MOD(n,10)/10-MOD(n,14)/14-MOD(n,22)/22-MOD(n,15)/15-MOD(n,21)/21-MOD(n,33)/33-MOD(n,35)/35-MOD(n,55)/55-MOD(n,77)/77+MOD(n,30)/30+MOD(n,42)/42+MOD(n,66)/66+MOD(n,70)/70+MOD(n,105)/105+MOD(n,110)/110+MOD(n,165)/165+MOD(n,385)/385+MOD(n,154)/154+MOD(n,231)/231-MOD(n,210)/210-MOD(n,330)/330-MOD(n,462)/462-MOD(n,770)/770-MOD(n,1155)/1155+MOD(n,2310)/2310
   PRINT p
NEXT n


END



数式を解読していた時、数式の部分を抜き取ると

実行結果
---------------------
18  0
19  .666666666666667
20  .333333333333333
21  0
22 -.333333333333333
23  .333333333333333
24  0
25  .666666666666667
26  .333333333333333
27  0
--------------------

一次篩6n±1の探究時見つけた性質

(6n-1)/3    小数点以下.66666666

(6n+1)/3    小数点以下.333333333

私の場合は、判別法に使用しています。


この度々登場するゼロと、非自明なゼロ点には、何らかの関係があるのでは?

そのように思っていました。

ちなみにこのゼロは　K=4 (2,3,5,7)　以降は出ません。105  0 　確認

素数の性質探究をしていると、ある一定領域にだけある性質が見られます。

今回enlong chiouさんの数式もどこまでなのか?

検証中です。K=25 でも10000　　せめて10万までは、と思ってていますが

どうでしょうか?

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

  

Re: 素数個数関数

 投稿者：しばっち  投稿日：2015年10月 2日(金)20時41分31秒
  > No.3868[元記事へ]

たろささんへのお返事です。

> 私の個人的予想では、プログラムで書きました。
>
> FOR n=1 TO 2
>    PRINT n;
>    PRINT n-1
> NEXT n
>
> FOR n=3 TO 48
>    PRINT n;
>    LET  p=MOD(n,2)/2+MOD(n,3)/3-MOD(n,2*3)/(2*3)
>    PRINT p
> NEXT n
>
> FOR n=49 TO 120
>    PRINT n;
>    LET  p=MOD(n,2)/2+MOD(n,3)/3+MOD(n,5)/5+MOD(n,7)/7-MOD(n,6)/6-MOD(n,10)/10-MOD(n,14)/14-MOD(n,15)/15-MOD(n,21)/21-MOD(n,35)/35+MOD(n,30)/30+MOD(n,42)/42+MOD(n,70)/70+MOD(n,105)/105-MOD(n,210)/210
>    PRINT p
> NEXT n
>
> FOR n=121 TO 168
>    PRINT n;
>    LET  p=MOD(n,2)/2+MOD(n,3)/3+MOD(n,5)/5+MOD(n,7)/7+MOD(n,11)/11-MOD(n,6)/6-MOD(n,10)/10-MOD(n,14)/14-MOD(n,22)/22-MOD(n,15)/15-MOD(n,21)/21-MOD(n,33)/33-MOD(n,35)/35-MOD(n,55)/55-MOD(n,77)/77+MOD(n,30)/30+MOD(n,42)/42+MOD(n,66)/66+MOD(n,70)/70+MOD(n,105)/105+MOD(n,110)/110+MOD(n,165)/165+MOD(n,385)/385+MOD(n,154)/154+MOD(n,231)/231-MOD(n,210)/210-MOD(n,330)/330-MOD(n,462)/462-MOD(n,770)/770-MOD(n,1155)/1155+MOD(n,2310)/2310
>    PRINT p
> NEXT n
>
> END

山中氏のプログラムを少しいじってみました。


OPTION ARITHMETIC RATIONAL
DIM P(25)
MAT READ P
DATA 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97
PRINT "SELECT CASE N"
FOR K=1 TO 10
   IF K=1 THEN
      PRINT "CASE 2 TO 8"
   ELSE
      PRINT "CASE";P(K)^2;"TO";P(K+1)^2-1
   END IF
   LET X=1
   LET Y=1
   FOR I=1 TO K
      LET X=X*(P(I)-1)
      LET Y=Y*P(I)
   NEXT I
   PRINT "P=N*";STR$(X/Y);
   LET L=0
   FOR J=1 TO 2^K-1
      LET T=J
      LET S=1
      LET B=-1
      FOR I=1 TO K
         IF MOD(T,2)=1 THEN
            LET S=S*P(I)
            LET B=-B
         END IF
         LET T=INT(T/2)
      NEXT I
      IF B<0 THEN PRINT "-"; ELSE PRINT "+";
      LET S$="MOD(N,"&STR$(S)&")/"&STR$(S)
      PRINT S$;
      LET L=L+LEN(S$)
      IF L>1000 THEN
         LET L=0
         PRINT
         PRINT "P=P";
      END IF
   NEXT J
   IF K>1 THEN PRINT "+";STR$(K-1) ELSE PRINT
NEXT K
PRINT "CASE ELSE"
PRINT "P=0"
PRINT "END SELECT"
END
 		 

Re: 素数個数関数

 投稿者：たろさ  投稿日：2015年10月 2日(金)21時15分29秒
  > No.3871[元記事へ]

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

ありがとうございます。最初は電卓で計算してました。

今年の4月から10進BASIC を始めました。

まだ、副プログラムの勉強中です。年齢的には3回目の二十歳代に爆走中なので

中々進みません。


ゼータ関数のζ(x)奇数を求める為

OPTION ARITHMETIC COMPLEX

LET P1=3.141592653589793
LET  E=2.718281828459045
FOR X= 3 TO 12
   FOR N= 1 TO 500
      LET P=(2*P1*n)
      LET Z=(1/(n^x*(E^P-1)))+R
      LET R=Z
   NEXT N
NEXT X
PRINT Z
END



エラー
EXTYP 1002
数値演算の桁あふれ

どのようにしたら良いかわかりません。


カシオ高精度計算サイト

zeta(3)-(7/180)*pi^3=-0.0037427455187320547576740910980487115437749296661

zeta(7)-(19/56700)*pi^7=-0.0037419276931926807928289645835153181149727186343

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

  

Re: 素数個数関数

 投稿者：しばっち  投稿日：2015年10月 2日(金)21時53分11秒
  > No.3872[元記事へ]

たろささんへのお返事です。

> カシオ高精度計算サイト
>
> zeta(3)-(7/180)*pi^3=-0.0037427455187320547576740910980487115437749296661
>
> zeta(7)-(19/56700)*pi^7=-0.0037419276931926807928289645835153181149727186343

OPTION ARITHMETIC NATIVE
LET S=0
FOR N=1 TO 100
   LET S=S+1/N^3/(EXP(2*PI*N)-1)
NEXT N
PRINT 7/180*PI^3-2*S

LET S=0
FOR N=1 TO 100
   LET S=S+1/N^7/(EXP(2*PI*N)-1)
NEXT N
PRINT 19/56700*PI^7-2*S
END

とりあえず下記のプログラムでゼータ関数計算できます

OPTION ARITHMETIC COMPLEX
LET X=COMPLEX(.5,14)
DO
   LET X=X-ZETA(X)/DF(X)
LOOP UNTIL ABS(ZETA(X))<1E-12
PRINT X
END

EXTERNAL  FUNCTION DF(X)
OPTION ARITHMETIC COMPLEX
LET H=1/128
LET DF=(ZETA(X+H)-ZETA(X))/H
END FUNCTION

EXTERNAL  FUNCTION ZETA(S)
OPTION ARITHMETIC COMPLEX
FOR M=1 TO 300
   LET SS=0
   FOR J=1 TO M
      LET SS=SS+(-1)^(J-1)*COMB(M-1,J-1)*J^(-S)
   NEXT J
   LET SUM=SUM+SS*2^(-M)
   IF ABS(SUM-S0)<1E-12 THEN
      LET ZETA=SUM/(1-2^(1-S))
      EXIT FUNCTION
   END IF
   LET S0=SUM
NEXT M
STOP
END FUNCTION

 		 

Re: 素数個数関数

 投稿者：たろさ  投稿日：2015年10月 3日(土)17時41分24秒
  しばっちさんへのお返事です。

ありがとうございます。自分の不勉強を棚に上げて、ろくに計算もしないで聞いてしまいました。

> たろささんへのお返事です。
>
> > カシオ高精度計算サイト
> >
> > zeta(3)-(7/180)*pi^3=-0.0037427455187320547576740910980487115437749296661
> >
> > zeta(7)-(19/56700)*pi^7=-0.0037419276931926807928289645835153181149727186343

(1/(n^3*e^(2*pi*n)-1))
(1/(n^7*e^(2*pi*n)-1))
カシオ高精度計算サイトで計算すると

計算回数　約300回　1000桁近くに
n^3---300買い目　8.7664534439000242640980007802080033699960889978E-827
n^7---300買い目 1.08227820295062027951827170126024732962914678985E-836

これをファイルにして読み込んで足しました。

!ゼータ関数　無限級数
OPTION ARITHMETIC RATIONAL
LET t0=TIME
LET A1=300
DIM A(A1)  !(1/(n^3*e^(2*pi*n)-1))

ASK DIRECTORY d$
LET f_name3$="ZATA_300_7"  ! ファイル名
LET f3$=d$&"\"&f_name3$&".txt"
!PRINT "ファイル保存場所[3] = ";f3$
OPEN #3: NAME f3$ ,ACCESS INPUT
FOR i2=1 TO A1
   INPUT #3: A(i2)
NEXT i2
CLOSE #3

FOR i2=1 TO A1
   LET Z=Z+A(i2)
NEXT i2
PRINT z

PRINT TIME-t0;"秒で計算しました"
END

------------------------------------------
有理数で出力すと

n^3

LET A=935686379017925191250796553673194052656087227132291731219603613975701972323372795070695707341051694333041461416617798292371870729138153420944855330166120553231062822792396969033035360742846471276962070515879208592530422237728427074787476956806821438649529241985758293460776861047916507563818694858789693813033935181462985075817986091401949987781387743696302049044433168708633425757646360787940379289393388500551855962450087280840910692903919579357382237509178105019962393425014909447675389857239761822966852476224426083709633194368749215953739803792878253703381223654493097348243698510213924057817098878283764085832228580612853735483995801712437547770824140608563623408511336820377837558294827506581439784956022278968540875346486789721414494706684449204326539775370091752190464631981938297646713105054959842789605364543275838448245221799947599825254618548523356213744489
LET B=500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
LET C=A/B
PRINT C*2
END
---------------------------------------
n^7

LET A=3741927693004140332863323157154235280467492756639334868567736653531814410941068770273590457531817355482005567002520408237720864547747721915916632587337776871933554687440080221071837310612327048769869521368660532830662354618883340779056093909806287487307684737302990106378585033070689464667036609902257636263307834325932679257576365681308595684841223259003559811275552937939579468842247542992272639058147977458713961480879276330002345365041104429191683348220167165522183005997757054728185792689086254331004637345347886072827617812926614305566197373185614245504655532125170388501140465714366951122966623111461256342374237003511786934661579718335126518422632496452302601740010396932535438346433415789501967033598982822413606269794652725843767089198596350512947505152492439682962193409574018904035316878264253458871562563344052079407237129934768119930808096245605085658529547003491997
LET B=2000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
LET C=A/B
PRINT C*2
END

--------------------------------------
この他、最初は20個のDATA読み込みで試しましたところ、n^3の場合は20個で20桁近くになり収束は速いと思いました。

この有理数を1000桁モードで小数にすると

              300 　 - .003742745516071700765003186214692776210624348908529166924878414455902807889293491180282782829364206777332165845666471193169487482916552613683779421320664482212924251291169587876132141442971385885107848282063516834370121688950913708299149907827227285754598116967943033173843107444191666030255274779435158775252135740725851940303271944365607799951125550974785208196177732674834533703030585443151761517157573554002207423849800349123363642771615678317429528950036712420079849573700059637790701559428959047291867409904897704334838532777474996863814959215171513014813524894617972389392974794040855696231268395513135056343328914322451414941935983206849750191083296562434254493634045347281511350233179310026325759139824089115874163501385947158885657978826737796817306159101480367008761858527927753190586852420219839371158421458173103353792980887199790399301018474194093424854977956
              100    - .0037427455160717007650031862146927762106243489085291669248784144559028078892934911802827828293642067773321658456664711931694874829165526136837794213206644822129242512911695878761321414429713858851078482820635168343701216889509137082991499078272272857545981169679430331738431074441868261957336010182715750390027869414412865137
zeta(3)-(7/180)*pi^3=-0.0037427455187320547576740910980487115437749296661

----------------------------------------------
zeta(7)-(19/56700)*pi^7=-0.0037419276931926807928289645835153181149727186343
              300       - .003741927693004140332863323157154235280467492756639334868567736653531814410941068770273590457531817355482005567002520408237720864547747721915916632587337776871933554687440080221071837310612327048769869521368660532830662354618883340779056093909806287487307684737302990106378585033070689464667036609902257636263307834325932679257576365681308595684841223259003559811275552937939579468842247542992272639058147977458713961480879276330002345365041104429191683348220167165522183005997757054728185792689086254331004637345347886072827617812926614305566197373185614245504655532125170388501140465714366951122966623111461256342374237003511786934661579718335126518422632496452302601740010396932535438346433415789501967033598982822413606269794652725843767089198596350512947505152492439682962193409574018904035316878264253458871562563344052079407237129934768119930808096245605085658529547003491997


他サイトの計算結果を持ち込み大変失礼だとお詫びします。

素人の興味本位の探究なので、お許しください。

自分の計算の仕方が間違っているのか不安でもあるし、

少々困惑しています。

別件で恐縮ですが


素数個数関数のプログラムで1000万刻みの最後の11を算出しています。

毎度説明の仕方が悪いので画像を添付させて頂きます。

現在5000万のラスト11個を計算中です。

どこまで、計算できるかわかりませんが、興味があります。

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

  

Re: 素数個数関数

 投稿者：たろさ  投稿日：2015年10月 4日(日)01時53分50秒
  > No.3874[元記事へ]

たろささんへのお返事です。

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


ありがとうございます。


こんな計算もしてみました。


高精度計算サイト[実数 DEG]

(1/tanh(pi*n))/n^7　　1000個の総和(50桁)
------------------------------------------------------------------------
1.0120912050751155074664592644340120758815100958282601246804669728404888


-------------------------------------------------------------------------------------
1.0120912050751155074664592644340120758815100958282601246804669728404888/pi^7*56700
=18.999999999999999996880540277244697295616850505925



(1/tanh(pi*n))/n^7　　2000個の総和(50桁)
--------------------------------------------------------------------------
1.012091205075115507635321564404603510287059971500990852791159096798337985


1.012091205075115507635321564404603510287059971500990852791159096798337985/pi^7*56700
=19.00000000000000000005059420455169038793700513113


(1/tanh(pi*n))/n^7　　3000個の総和(50桁)
--------------------------------------------------------------------------
1.01209120507511550763769343194927773249232817807833932916910770970633823112

1.01209120507511550763769343194927773249232817807833932916910770970633823112/pi^7*56700
=19.0000000000000000000951213016382277088487439393


(1/tanh(pi*n))/n^7　　4000個の総和(50桁)
--------------------------------------------------------------------------
1.012091205075115507637881167504425825408204474419642945099075257796859948129

1.012091205075115507637881167504425825408204474419642945099075257796859948129/pi^7*56700
=19.000000000000000000098645663405158616205502952777


(1/tanh(pi*n))/n^7　　5000個の総和(50桁)
--------------------------------------------------------------------------
1.012091205075115507637911166831755326765437554701091235620833422619592271006

1.012091205075115507637911166831755326765437554701091235620833422619592271006/pi^7*56700
=19.00000000000000000009920884112708953883115534931


(1/tanh(pi*n))/n^7　　6000個の総和(50桁)
--------------------------------------------------------------------------
1.0120912050751155076379182566406059766316666495830741765231805968595393737569

1.0120912050751155076379182566406059766316666495830741765231805968595393737569/pi^7*56700
=19.00000000000000000009934193819135298118050986422


(1/tanh(pi*n))/n^7　　7000個の総和(50桁)
--------------------------------------------------------------------------
1.0120912050751155076379204110636549249928637609364754523867038275935567615542

1.0120912050751155076379204110636549249928637609364754523867038275935567615542/pi^7*56700
=19.00000000000000000009938238320038453463162448429


(1/tanh(pi*n))/n^7　　8000個の総和(50桁)
--------------------------------------------------------------------------
1.01209120507511550763792119155542153225318634345472572283306849628438330899987

1.01209120507511550763792119155542153225318634345472572283306849628438330899987/pi^7*56700
=19.000000000000000000099397035381424322780077309755


(1/tanh(pi*n))/n^7　　9000個の総和(50桁)
--------------------------------------------------------------------------
1.01209120507511550763792151359170216143722870003985016630132115095218768678188

1.01209120507511550763792151359170216143722870003985016630132115095218768678188/pi^7*56700
=19.000000000000000000099403080972277471611835907716


(1/tanh(pi*n))/n^7　　10000個の総和(50桁)
--------------------------------------------------------------------------
1.0120912050751155076379216604832428266938275622085900642365819094261613209145

1.0120912050751155076379216604832428266938275622085900642365819094261613209145/pi^7*56700
=19.00000000000000000009940583856888402718121212196


ラストDATA
10000 1E-28
-------------------------------------------------------------------------------------
1.0120912050751155076379533102494968961760606292754/pi^7*56700
=19.000000000000000000099999999999999999999999999998

1.0120912050751155076379533102494968961760606292755/pi^7*56700
=19.0000000000000000001

1.01209120507511550763795331024949689617606062927545/pi^7*56700
=19.000000000000000000099999999999999999999999999999

1.01209120507511550763795331024949689617606062927546/pi^7*56700
=19.0000000000000000001

1.012091205075115507637953310249496896176060629275455/pi^7*56700
=19.000000000000000000099999999999999999999999999999

1.012091205075115507637953310249496896176060629275456/pi^7*56700
=19.0000000000000000001

1.0120912050751155076379533102494968961760606292754556/pi^7*56700
=19.000000000000000000099999999999999999999999999999

1.0120912050751155076379533102494968961760606292754557/pi^7*56700
=19.0000000000000000001

1.01209120507511550763795331024949689617606062927545567/pi^7*56700
=19.000000000000000000099999999999999999999999999999

1.01209120507511550763795331024949689617606062927545568/pi^7*56700
=19.0000000000000000001

1.01209120507511550763795331024949689617606062927545567/pi^7*56700
=19.000000000000000000099999999999999999999999999999

1.01209120507511550763795331024949689617606062927545568/pi^7*56700
=19.0000000000000000001

1.012091205075115507637953310249496896176060629275455679/pi^7*56700
=19.000000000000000000099999999999999999999999999999

どこかで、DATA飛ばしてしまったのだろうか?

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

  

Re: 素数個数関数

 投稿者：しばっち  投稿日：2015年10月 5日(月)00時15分29秒
  > No.3880[元記事へ]

たろささんへのお返事です。

下記のようにするとムダな計算をしないで済みます

OPTION ARITHMETIC NATIVE
FOR N=2 TO 10000000
   LET FL=0
   FOR J=2 TO INT(SQR(N))
      IF MOD(N,J)=0 THEN
         LET FL=1
         EXIT FOR
      END IF
   NEXT J
   IF FL=0 THEN LET COUNT=COUNT+1
   IF MOD(N,100000)=0 THEN PRINT N;COUNT
   !' IF N>=9999990 AND N<=10000000 THEN PRINT N;COUNT
NEXT N
END
 		 

Re: 素数個数関数

 投稿者：たろさ  投稿日：2015年10月 5日(月)06時11分4秒
  > No.3882[元記事へ]

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

お世話になります。

素数個数関数Ver.1
FOR I=100000036 TO 100000036
3273.54秒(54分33.54秒)

素数個数関数Ver.2
FOR N=2 TO 100000036
IF MOD(N,10000000)=0
1659.37秒(27分39.37秒)


IF N>=9999990 AND N<=10000000 THEN PRINT N;COUNT
も試しました。計算時間は、同様でした。

今朝は、10億で計算中です。

1億　　5761455個　約28分で表示されました。
------------------------------------------------

何もお礼は出来ませんが、先週買った。インテルのCPUが1個約8千円で

4.3GHz で、自作PC組み立ててから、ずっと素数個数関数Ver.1です。

画像のハンドルネーム間違えてました。目が老眼なのでお許しください。

------------------------------------------------------------------
追記

３億まで算出しました。まだ計算中


99999989(5761455th prime)
199999991(11078937th prime)
299999977(16252325th prime)

素数個数関数Ver.2  算出時間　約2時間(インテルG3258 4.3GHz)


自作プログラムとの比較

自然数 n億からn億までの素数リスト作成所要時間

8億から9億まで　2時間17分48.61秒(インテルi3-4170 3.7GHz)
7億から8億まで　2時間11分23.74秒(以下同じ)
6億から7億まで　2時間04分51.08秒
5億から6億まで　1時間56分54.41秒
4億から5億まで　1時間48分24.34秒
3億から4億まで　1時間27分19.75秒
2億から3億まで　1時間27分04.62秒
1億から2億まで　1時間12分04.92秒
2から1億まで　  0時間49分08.87秒


億単位の素数の個数を数える場合は、実用性があると思います。

------------------------------------------------------------------
追記

10億　45006.83秒(12時間30分06.83秒)

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

  

Re: 素数個数関数

 投稿者：たろさ  投稿日：2015年11月22日(日)02時50分55秒
  > No.3871[元記事へ]

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

K=13 の数式から

OPTION ARITHMETIC RATIONAL     !有理数モード
FOR N=1681 TO 1848 !K=13  19110297600/131710070791
   LET P=MOD(N,2)/2+MOD(N,3)/3-MOD(N,6)/6+MOD(N,5)/5-MOD(N,10)/10-MOD(N,15)/15+MOD(N,30)/30+MOD(N,7)/7-MOD(N,14)/14-MOD(N,21)/21+MOD(N,42)/42-MOD(N,35)/35+MOD(N,70)/70+MOD(N,105)/105-MOD(N,210)/210+MOD(N,11)/11-MOD(N,22)/22-MOD(N,33)/33+MOD(N,66)/66-MOD(N,55)/55+MOD(N,110)/110+MOD(N,165)/165-MOD(N,330)/330-MOD(N,77)/77+MOD(N,154)/154+MOD(N,231)/231-MOD(N,462)/462+MOD(N,385)/385-MOD(N,770)/770-MOD(N,1155)/1155+MOD(N,2310)/2310+MOD(N,13)/13-MOD(N,26)/26-MOD(N,39)/39+MOD(N,78)/78-MOD(N,65)/65+MOD(N,130)/130+MOD(N,195)/195-MOD(N,390)/390-MOD(N,91)/91+MOD(N,182)/182+MOD(N,273)/273-MOD(N,546)/546+MOD(N,455)/455-MOD(N,910)/910-MOD(N,1365)/1365+MOD(N,2730)/2730-MOD(N,143)/143+MOD(N,286)/286+MOD(N,429)/429-MOD(N,858)/858+MOD(N,715)/715-MOD(N,1430)/1430-MOD(N,2145)/2145+MOD(N,4290)/4290+MOD(N,1001)/1001-MOD(N,2002)/2002-MOD(N,3003)/3003+MOD(N,6006)/6006-MOD(N,5005)/5005+MOD(N,10010)/10010+MOD(N,15015)/15015-MOD(N,30030)/30030+MOD(N,17)/17-MOD(N,34)/34-MOD(N,51)/51+MOD(N,102)/102-MOD(N,85)/85+MOD(N,170)/170+MOD(N,255)/255-MOD(N,510)/510-MOD(N,119)/119+MOD(N,238)/238+MOD(N,357)/357

長いので、途中省略

   LET P=P+MOD(N,3951302123730)/3951302123730+MOD(N,921970495537)/921970495537-MOD(N,1843940991074)/1843940991074-MOD(N,2765911486611)/2765911486611+MOD(N,5531822973222)/5531822973222-MOD(N,4609852477685)/4609852477685+MOD(N,9219704955370)/9219704955370+MOD(N,13829557433055)/13829557433055-MOD(N,27659114866110)/27659114866110+MOD(N,1448810778701)/1448810778701-MOD(N,2897621557402)/2897621557402-MOD(N,4346432336103)/4346432336103+MOD(N,8692864672206)/8692864672206-MOD(N,7244053893505)/7244053893505+MOD(N,14488107787010)/14488107787010+MOD(N,21732161680515)/21732161680515-MOD(N,43464323361030)/43464323361030-MOD(N,10141675450907)/10141675450907+MOD(N,20283350901814)/20283350901814+MOD(N,30425026352721)/30425026352721-MOD(N,60850052705442)/60850052705442+MOD(N,50708377254535)/50708377254535-MOD(N,101416754509070)/101416754509070-MOD(N,152125131763605)/152125131763605+MOD(N,304250263527210)/304250263527210
   LET p1=P*131710070791
   PRINT STR$(p1)!&",";  !168個
NEXT N

END

----------------------------------------
DATA を見ると

934817502941-915707205341=19110297600


19110297600-(-112599773191)=131710070791

最終的に

19110297600
±131710070791

三種類の差で、DATAが出来てます。

驚きました。

山中和義様　しばっち様　お二人に感謝します。


追記
K=12の場合も
　477757440
±3212440751
三種類の差で、DATAが出来てました。

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

  

Re: 素数個数関数

 投稿者：たろさ  投稿日：2015年11月25日(水)01時39分39秒
  > No.3960[元記事へ]

たろささんへのお返事です。

与えられた数より小さい素数の個数を求めるプログラム

暫定版です。

K=4

!OPTION ARITHMETIC RATIONAL     !有理数モード
LET N=49  !K=4  8/35
LET P=MOD(N,2)/2+MOD(N,3)/3-MOD(N,6)/6+MOD(N,5)/5-MOD(N,10)/10-MOD(N,15)/15+MOD(N,30)/30+MOD(N,7)/7-MOD(N,14)/14-MOD(N,21)/21+MOD(N,42)/42-MOD(N,35)/35+MOD(N,70)/70+MOD(N,105)/105-MOD(N,210)/210
LET p1=P*35

LET R=72 !DATAの個数(120-49+1)
DIM A(R)
LET A(1)=P1 !p1=28
LET S=P1
!PRINT 49;S
LET R=2
FOR N=49+1 TO 120 STEP 1
   IF ISPRIME(N)=1 THEN
      LET S=S-(-27 ) !35-8
      !  PRINT n;S
      LET A(R)=S
      LET R=R+1
   ELSE
      LET S=S+(-8 ) !等差減少
      !  PRINT n;S
      LET A(R)=S
      LET R=R+1
   END IF
NEXT N

!FOR R=1 TO 72
!   PRINT A(R)   !DATA プリント
!NEXT R

LET R=1
FOR N=49 TO 120                 !K=4  8/35
   LET pn=(n*8+A(R))/35
   PRINT n;pn+3
   LET R=R+1
NEXT n

END

EXTERNAL  FUNCTION ISPRIME(N)
!  OPTION ARITHMETIC RATIONAL
IF N = 2 THEN
   LET ISPRIME=1
   EXIT FUNCTION
END IF
IF N = 1 OR MOD(N , 2) = 0 THEN
   LET ISPRIME=0
   EXIT FUNCTION
END IF
LET D = (N - 1) / 2
LET S = 0
DO WHILE MOD(D , 2) = 0
   LET D = INT(D / 2)
   LET S=S+1
LOOP
FOR I=1 TO 10
   LET ISP=0
   READ A    !' n < 341550071728321 なら a = 2, 3, 5, 7, 11, 13, 17
   DATA 2,3,5,7,11,13,17,23,29,31
   LET ISP = 0
   LET R = POWMOD(A, D, N)
   IF R = 1 OR R = N - 1 THEN
      LET ISP = 1
   END IF
   LET R = POWMOD(R, 2, N)
   FOR J = 0 TO S-1
      IF R = N - 1 THEN
         LET ISP = 1
      END IF
      LET R = POWMOD(R, 2, N)
   NEXT J
   IF ISP=0 THEN
      LET ISPRIME=0
      EXIT FUNCTION
   END IF
NEXT I
LET ISPRIME=1
END FUNCTION

EXTERNAL  FUNCTION POWMOD(B,P,M)
!  OPTION ARITHMETIC RATIONAL
LET RESULT = 1
DO WHILE P > 0
   IF MOD(P , 2)= 1 THEN
      LET RESULT = MOD(RESULT * B , M)
   END IF
   LET B = MOD(B * B, M)
   LET P = INT(P / 2)
LOOP
LET POWMOD=RESULT
END FUNCTION



----------------------------------------------
K=14  から　K=16 まで

!OPTION ARITHMETIC RATIONAL     !有理数モード
LET R=360 !DATAの個数(2208-1849+1)
DIM A(R)
LET S=45086430862710
LET A(1)=S
LET R=2
FOR N=1849+1 TO 2208 STEP 1
   IF ISPRIME(N)=1 THEN
      LET S=S-(-4860900544813 ) !5663533044013-802632499200
      !  PRINT n;S
      LET A(R)=S
      LET R=R+1
   ELSE
      LET S=S+(-802632499200 )
      !  PRINT n;S
      LET A(R)=S
      LET R=R+1
   END IF
NEXT N

LET R=1
FOR N=1849 TO 2208                 !K=14  802632499200/5663533044013
   LET pn=(n*802632499200+A(R))/5663533044013
   PRINT n;pn+13
   LET R=R+1
NEXT n

LET R=600 !DATAの個数(2808-2209+1)
DIM B(R)
LET S=99561214908855
LET B(1)=S
LET R=2
FOR N=2209+1 TO 2808 STEP 1 !CASE 2209 TO 2808
   IF ISPRIME(N)=1 THEN
      LET S=S-(-9968041656757 )!11573306655157-1605264998400
      !  PRINT n;S
      LET B(R)=S
      LET R=R+1
   ELSE
      LET S=S+(-1605264998400 )
      !  PRINT n;S
      LET B(R)=S
      LET R=R+1
   END IF
NEXT N

LET R=1
FOR N=2209 TO 2808                 !K=15  1605264998400/11573306655157
   LET pn=(n*1605264998400+B(R))/11573306655157
   PRINT n;pn+14
   LET R=R+1
NEXT n

LET R=672 !DATAの個数(3480-2809+1)
DIM C(R)
LET S=553533983592098
LET C(1)=S
LET R=2
FOR N=2809+1 TO 3480 STEP 1
   IF ISPRIME(N)=1 THEN
      LET S=S-(-40762420985117) !47183480978717-6421059993600
      !  PRINT n;S
      LET C(R)=S
      LET R=R+1
   ELSE
      LET S=S+(-6421059993600)
      !  PRINT n;S
      LET C(R)=S
      LET R=R+1
   END IF
NEXT N

LET R=1
FOR N=2809 TO 3480              !K=16  6421059993600/47183480978717
   LET pn=(n*6421059993600+C(R))/47183480978717
   PRINT n;pn+15
   LET R=R+1
NEXT n

END

EXTERNAL  FUNCTION ISPRIME(N)
!  OPTION ARITHMETIC RATIONAL
IF N = 2 THEN
   LET ISPRIME=1
   EXIT FUNCTION
END IF
IF N = 1 OR MOD(N , 2) = 0 THEN
   LET ISPRIME=0
   EXIT FUNCTION
END IF
LET D = (N - 1) / 2
LET S = 0
DO WHILE MOD(D , 2) = 0
   LET D = INT(D / 2)
   LET S=S+1
LOOP
FOR I=1 TO 10
   LET ISP=0
   READ A    !' n < 341550071728321 なら a = 2, 3, 5, 7, 11, 13, 17
   DATA 2,3,5,7,11,13,17,23,29,31
   LET ISP = 0
   LET R = POWMOD(A, D, N)
   IF R = 1 OR R = N - 1 THEN
      LET ISP = 1
   END IF
   LET R = POWMOD(R, 2, N)
   FOR J = 0 TO S-1
      IF R = N - 1 THEN
         LET ISP = 1
      END IF
      LET R = POWMOD(R, 2, N)
   NEXT J
   IF ISP=0 THEN
      LET ISPRIME=0
      EXIT FUNCTION
   END IF
NEXT I
LET ISPRIME=1
END FUNCTION

EXTERNAL  FUNCTION POWMOD(B,P,M)
!  OPTION ARITHMETIC RATIONAL
LET RESULT = 1
DO WHILE P > 0
   IF MOD(P , 2)= 1 THEN
      LET RESULT = MOD(RESULT * B , M)
   END IF
   LET B = MOD(B * B, M)
   LET P = INT(P / 2)
LOOP
LET POWMOD=RESULT
END FUNCTION

-----------------------------------------------

K=19あたりから数式の出力も増大してマス。

簡単にDATAの初期値を求める数式があれば、K=30もできそうです。

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

  

Re: 素数個数関数

 投稿者：しばっち  投稿日：2015年12月19日(土)18時12分4秒
  > No.3882[元記事へ]

素数個数関数をC++で書いてみました
64bit浮動小数(double)ではなく、高速化のため64bit整数(long long)を使用しました
PRIMECOUNT関数の第1引数に求めたい範囲の下限を第2引数に上限を指定します

更に速くするため、ＯＰＥＮ　ＭＰを使って並列化(マルチスレッド化)してみました
第3引数には実行スレッド数です。0(自動)又は、スレッド数を指定して下さい
ASSIGN文でdllファイルのパスを指定してください

LET S$=PACKDBL$(2)
LET E$=PACKDBL$(10000000)
PRINT PRIMECOUNT(S$,E$,0)
END

EXTERNAL FUNCTION PRIMECOUNT(S$,E$,THREADS)
ASSIGN "primecount.dll","primecount",FPU
END FUNCTION

-------------------------------------------------------------
                        primecount.cpp

#include <cmath>
#include <omp.h>

extern "C" __declspec(dllexport) double primecount(double *s,double *e,int threads)
{
long long i,j,count=0;
int flag;
if (threads>0) omp_set_num_threads(threads);
#pragma omp parallel for private(j,flag)
for(i=(long long)*s;i<=(long long)*e;i++){
    flag=0;
        for(j=2;j<=floor(sqrt((double)i));j++){
            if (i % j==0) {
                flag=1;
                break;
                        }
                                   }
        if (i>=2 && flag==0) {
            #pragma omp atomic
            count++;
                      }
                                                    }
return (double)count;
}

ネット上のサンプルを参考に書いてみました
コンパイルにはOPEN MP対応のC/C++コンパイラーを使用してください
64bit整数型を使用し、C++で記述しています
C++にしたのは、私の開発環境では64bit整数がCで機能しなかったからです
unsignedとしていないのは、単にコンパイルエラー(open mp上の制約 ?)が出たからです

1000万以下の素数個数を計測してみました(2～10000000まで 664579個)
十進BASIC(2進モード)  約220秒
BASIC ACC             約110秒

primecount.dll        約28秒(THREADS=0を指定、※2スレッド)
となりました。

こちらはデュアルコアCPU(ハイパースレッディング未対応?)なので、上位CPU(クアッドコアCPU等)なら
更に速くなると思われます(220/28=約7.8倍)

また、実行環境や開発環境、オプションの有無などにより結果は変わると思われます

※1～10000000とするより
1～1000000 , 1000001～2000000 , 2000001～3000000 ... 90000001～10000000
と分割した方が若干早く終わるようです

OPTION ARITHMETIC NATIVE
LET T=TIME
LET ST=1000000
FOR I=1 TO ST*10 STEP ST
   LET L=PRIMECOUNT(I,I+ST-1)
   LET S=S+L
   PRINT I+ST-1;":";L;S
NEXT I
PRINT TIME-T

LET T=TIME
PRINT ST*10;":";PRIMECOUNT(1,ST*10)
PRINT TIME-T
END

EXTERNAL  FUNCTION PRIMECOUNT(S,E)
OPTION ARITHMETIC NATIVE
FOR I=S TO E
   LET FL=0
   FOR J=2 TO INT(SQR(I))
      IF MOD(I,J)=0 THEN
         LET FL=1
         EXIT FOR
      END IF
   NEXT J
   IF I>=2 AND FL=0 THEN LET C=C+1
NEXT I
LET PRIMECOUNT=C
END FUNCTION
 		 

戻る