|
投稿者：しばっち
投稿日：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
|
|
|
|
投稿者：たろさ
投稿日：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
|
|
|
|
投稿者：しばっち
投稿日：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
|
|
|
|
投稿者：たろさ
投稿日：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
|
|
|
|
投稿者：しばっち
投稿日：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
|
|
|
|
投稿者：たろさ
投稿日：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
|
|
|
|
投稿者：たろさ
投稿日：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
|
|
|
|
投稿者：しばっち
投稿日：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
|
|
|
|
投稿者：たろさ
投稿日：2015年10月 5日(月)06時11分4秒
|
|
|
|
投稿者：たろさ
投稿日：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
|
|
|
|
投稿者：たろさ
投稿日：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
|
|
|
|
投稿者：しばっち
投稿日：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
|
|
|
戻る
