|
投稿者：たろさ
投稿日：2015年10月 2日(金)02時34分40秒
|
|
|
|
! ζ(x) X=500 を求めるプログラム
100 OPTION ARITHMETIC RATIONAL
101 LET t0=TIME
110 LET format$="#." & REPEAT$("#",1140)
120 LET p=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199
130 LET b=-16596380640568557229852123088077134206658664302806671892352650993155331641220960084014956088135770921465025323942809207851857992860213463783252745409096420932509953165466735675485979034817619983727209844291081908145597829674980159889976244240633746601120703300698329029710482600069717866917229113749797632930033559794717838407415772796504419464932337498642714226081743688706971990010734262076881238322867559275748219588404488023034528296023051638858467185173202483888794342720837413737644410765563213220043477396887812891242952336301344808165757942109887803692579439427973561487863524556256869403384306433922049078300720480361757680714198044230522015775475287075315668886299978958150756677417180004362981454396613646612327019784141740499835461/8365830
140 LET q=1
150 FOR n=1 TO 500 ! 階乗数
160 LET k=n*q
170 LET q=k
180 NEXT n
190 LET z=(-1*b*2^499)/k*p^500
200 PRINT USING format$: z
210 PRINT TIME-t0;"秒で計算しました"
220 END
実行結果
1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000305493636349960468205197939321361769978940274057232666389361390928129162652472045770185751013331084857067694824349847835941470546838642267279113733370156013676833730923434530872620333347455442425607860198374600438389341457088354275417213227146995205688275457893692917327448071083856710385116593545333980085810590902040565924433365890972440295472134967177946250193272901910685580918026995490249850797834908996937559304168159737646719914815706639565554808346422994189258230903331312648052123784001572222418899043993491675713595743325898062662150739575848900411992098649088533778468543823487152960410139299824092427165312585873162835834785390156964316151479380817898770583201397753046035925243913393144514022226717256820109064643891961466904626664792337893058039537925892858730037171249346830463620064564420926132913346370610578490178139276012234619998000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
19027/100 秒で計算しました
----------------------------
LET z=(-1*b*2^499)/k*p^500
----------------------------
数式について
P=π
B=ベルヌーイ数 -1と1が出ます。
K=階乗
----------------------------
なお、ζ(x)求められるのは偶数のみで奇数はできせん。
質問です。
REPEAT$("#",1140)
999以上は全て0になるのでしょうか?
上記の数式はWikipedia 「リーマンゼータ関数」ゼータ関数の特殊値を参照しました。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2015年10月22日(木)20時01分35秒
|
|
|
|
> No.3869[元記事へ]
たろささんへのお返事です。
副プログラムの勉強中です。少し改良しました。
! ζ(x)1000桁 X=偶数を求めるプログラム
OPTION ARITHMETIC RATIONAL
LET t0=TIME
LET format$="#." & REPEAT$("#",999)
LET p=PI
FOR K=2 TO 100 STEP 2
LET q=1
FOR n=1 TO k ! 階乗数
LET ka=n*q
LET q=ka
NEXT n
LET b=BERNOULLI(K)
IF b<0 THEN LET L=-1 ELSE LET L=1
LET z=(L*b*2^(k-1))/q*p^k
PRINT "zeta(";k;")";
PRINT USING format$: z
NEXT k
PRINT TIME-t0;"秒で計算しました"
END
EXTERNAL FUNCTION BERNOULLI(K) !'ベルヌーイ定数
OPTION ARITHMETIC RATIONAL
LET C=1
LET D=1
LET N=K/2
IF K=1 THEN
LET BERNOULLI=-1/2
ELSEIF MOD(K,2)=1 THEN
LET BERNOULLI=0
ELSE
FOR M=N TO 1 STEP-1
LET T=T+(-1)^M*D*M^(K-1)
LET C=C*(N+M+1)/(N-M+1)
LET D=D+C
NEXT M
LET BERNOULLI=-T*K/(2^K*(2^K-1))
END IF
END FUNCTION
zeta( 3000 )まで、zeta( 3100 )計算中です。
この辺りで、1000桁 1に収束?
できたら、zata(奇数)1000桁も出して見たいと思っています。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2015年10月25日(日)12時30分34秒
|
|
|
|
> No.3896[元記事へ]
たろささんへのお返事です。
副プログラムの勉強中です。ζ(x)の奇数を求めるプログラム
今のところ、ζ(2n+1) のn=奇数限定です。偶数が上手く出ません。
! Ramanujan ζ(2n+1)
OPTION ARITHMETIC DECIMAL_HIGH
LET t0=TIME
FOR n=1 TO 10 STEP 2
LET f=2^(2*n)*pi^(2*n+1)
LET S=0
FOR k1=0 TO n+1
LET S=S+((-1)^(k1+1))*BERNOULLI(2*k1)/FACT(2*k1)*BERNOULLI(2*n+2-2*k1)/FACT(2*n+2-2*k1)
NEXT k1
LET S1=0
FOR k2=1 TO 330
LET q=2*PI*k2
LET S1=S1+(k2^(-2*n-1))/(EX(q)-1)
NEXT k2
LET zeta=f*S-2*S1
PRINT 2*n+1;zeta
NEXT n
PRINT TIME-t0;"秒"
END
EXTERNAL FUNCTION BERNOULLI(n) !'ベルヌーイ定数
OPTION ARITHMETIC DECIMAL_HIGH
DIM b(0 TO n)
LET b(0)=1
IF n=0 THEN LET BERNOULLI=1
FOR k=2 TO n+1
LET s = 0
FOR i=0 TO k-2
LET s = s + comb(k,i)*b(i)
NEXT i
LET b(k-1) = - s / k
LET BERNOULLI=b(k-1)
NEXT k
END FUNCTION
!1000桁モードで利用する指数関数
EXTERNAL FUNCTION EX(x)
OPTION ARITHMETIC DECIMAL_HIGH
FUNCTION s(y,n)
LET t=y*x/n
IF ABS(t)<=EPS(0) THEN
LET s=y+t
ELSE
LET s=y+s(t,n+1)
END IF
END FUNCTION
LET EX=s(1,1)
END FUNCTION
上限
FACT(452)226まで
FOR n=223 TO 225 STEP 2
447
451
4329.45 秒
1時間12分09.45秒
AMD 2.4GHz
かなり重いです。
数式の展開方法がわからないので、手探り状態です。
出来れば、有理数モードで、と思っていますが、今のところ成功しません。
宜しくお願いします。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2015年11月 5日(木)19時12分52秒
|
|
|
|
> No.3897[元記事へ]
たろささんへのお返事です。
暫定版ですが、出来ました。
#3880
! Cohen 2000 MOD(x,4)=1(-)Ramanujan ζ(2n+1) zeta(x)1000桁精度
OPTION ARITHMETIC DECIMAL_HIGH
PRINT DATE$;"/"; TIME$
LET t0=TIME
LET Na=0
FOR n=5 TO 9 STEP 4
LET Na=Na+2
LET f=((2*PI)^n)/(FACT(n+1)*(n-1))
LET S=0
LET w= INT((n+1)/4)
FOR k1=0 TO w
LET S=S+((-1)^(k1))*(n+1-4*k1)*comb(n+1,2*k1)*BERNOULLI(n+1-2*k1)*BERNOULLI(2*k1)
NEXT k1
LET S1=0
FOR k2=1 TO 365 !数値溢れ370まで
LET q=2*PI*k2
LET S1=S1+(EX(q)*((1+(4*PI*k2)/(n-1))-1))/((k2^n)*(EX(q)-1)^2)
NEXT k2
LET zeta=f*S-2*S1
PRINT n;zeta+Ramanujan(Na)
NEXT n
PRINT TIME-t0;"秒"
PRINT DATE$;"/"; TIME$
END
EXTERNAL FUNCTION Ramanujan(n) !zeta(x) MOD(x,4)=3
OPTION ARITHMETIC DECIMAL_HIGH
FOR n=1 TO n
LET f=2^(2*n)*pi^(2*n+1)
LET S=0
FOR k1=0 TO n+1
LET S=S+((-1)^(k1+1))*BERNOULLI(2*k1)/FACT(2*k1)*BERNOULLI(2*n+2-2*k1)/FACT(2*n+2-2*k1)
NEXT k1
LET S1=0
FOR k2=1 TO 360
LET q=2*PI*k2
LET S1=S1+(k2^(-2*n-1))/(EX(q)-1)
NEXT k2
LET Ramanujan=f*S-2*S1
NEXT n
END FUNCTION
EXTERNAL FUNCTION BERNOULLI(n) !'ベルヌーイ定数
OPTION ARITHMETIC DECIMAL_HIGH
DIM b(0 TO n)
LET b(0)=1
IF n=0 THEN LET BERNOULLI=1
FOR k=2 TO n+1
LET s=0
FOR i=0 TO k-2
LET s=s+comb(k,i)*b(i)
NEXT i
LET b(k-1)=-s/k
LET BERNOULLI=b(k-1)
NEXT k
END FUNCTION
!1000桁モードで利用する指数関数
EXTERNAL FUNCTION EX(x)
OPTION ARITHMETIC DECIMAL_HIGH
FUNCTION s(y,n)
LET t=y*x/n
IF ABS(t)<=EPS(0) THEN
LET s=y+t
ELSE
LET s=y+s(t,n+1)
END IF
END FUNCTION
LET EX=s(1,1)
END FUNCTION
似たような計算を二度繰り返しているので、しばっちさんなら半分以下になるのでは?
宜しくお願いします。
1000桁精度確認画像を添付しました。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2015年11月 5日(木)19時17分42秒
|
|
|
|
> No.3945[元記事へ]
たろささんへのお返事です。
zeta(x)1000桁プロジェクトの副産物です。
OPTION ARITHMETIC DECIMAL_HIGH
LET t0=TIME
FOR n=5 TO 101 STEP 4 !1000桁なので、1つ求めるのに3分前後
PRINT n;Cohen(n)
NEXT n
PRINT TIME-t0;"秒"
END
EXTERNAL FUNCTION Cohen(n)
OPTION ARITHMETIC DECIMAL_HIGH
LET n=n
LET f=((2*PI)^n)/(FACT(n+1)*(n-1))
LET S=0
LET w= INT((n+1)/4)
FOR k1=0 TO w
LET S=S+((-1)^(k1))*(n+1-4*k1)*comb(n+1,2*k1)*BERNOULLI(n+1-2*k1)*BERNOULLI(2*k1)
NEXT k1
LET S1=0
FOR k2=1 TO 365 !数値溢れ370まで
LET q=2*PI*k2
LET S1=S1+(EX(q)*((1+(4*PI*k2)/(n-1))-1))/((k2^n)*(EX(q)-1)^2)
NEXT k2
LET Cohen=f*S-2*S1
END FUNCTION
EXTERNAL FUNCTION BERNOULLI(n) !'ベルヌーイ定数
OPTION ARITHMETIC DECIMAL_HIGH
DIM b(0 TO n)
LET b(0)=1
IF n=0 THEN LET BERNOULLI=1
FOR k=2 TO n+1
LET s=0
FOR i=0 TO k-2
LET s=s+comb(k,i)*b(i)
NEXT i
LET b(k-1)=- s / k
LET BERNOULLI=b(k-1)
NEXT k
END FUNCTION
!1000桁モードで利用する指数関数
EXTERNAL FUNCTION EX(x)
OPTION ARITHMETIC DECIMAL_HIGH
FUNCTION s(y,n)
LET t=y*x/n
IF ABS(t)<=EPS(0) THEN
LET s=y+t
ELSE
LET s=y+s(t,n+1)
END IF
END FUNCTION
LET EX=s(1,1)
END FUNCTION
算出DATA
1.00374187319732128820155269119480001746245242299590766340483002846700467633058071968805856226611181108206128491019325061879605107651261114766526275503361003853072936960008291766340390487115974168534515598148720703074119581432829005076897960080729783313913675902021836616743235177078620582793116121623385243650737831989816409452105664540813213633855367667620147214800141227538339590493683121477293087676562027398949714485347712152335504957823492370068640259966669520659878561730890939369943486268450302072478377125499645534741877271067064330638764179343923989667027533560460548667576913476231269740949181859624014613840494879775926675918637565864354272939964758378784384814853739817998063520249485509081842409908405798966658043605867324226611559761338616788511916966279166969797912586225650116869534812022674061393607939807598695611881201024178316343676062941988614835197888129877691966657602424720720262267094658156963690202633514058638342845140431580867749428636844298782565307497805992730890615134629252844(100)
1.00374187319732128820155269119480001746245242299590766340483002846700467633058071968805856226611181108206128491019325061879605107651261114766526275503361003853072936960008291766340390487115974168534515598148720703074119581432829005076897960080729783313913675902021836616743235177085693653641094489900811945233516538437429962038403426889569939588603957843064739604657825503776859020334288039116530140662302979890324065405785969952934198127439339431268429059901727756268085278992170400641101103299862273911962461857405265315974823371680328100039732348448548943089086839961125441567327555121034642545598469936536206394523268336888456540698826104558575132586428240664077946490845779685218987559700141870359282083557544144299982311690543722569036691797965281872287465008189463863731384876762920380669760295584723060960795665958530142437960984475131578652442900631580728938025896365269527757322803875438206014238908280368546174254091769956430862124958753577254358357695958594014017186970858451472433850067958828382(200)
1.00374187319732128820155269119480001746245242299590766340483002846700467633058071968805856226611181108206128491019325061879605107651261114766526275503361003853072936960008291766340390487115974168534515598148720703074119581432829005076897960080729783313913675902021836616743235177085693653641094489900811945233516538437429962038403426889569939588603957843064739604657825503776859020334288039116530140662302979890324065405785969952934198127439339431268429059901727756268085278992170400641101103299862273911962461857405265315974823371680328100039732348448548943089086839961125441567327555121034642545598469936536206394523268336888456540698826104558575132586428240664077946490845839188975785892257991619023020637406929243458317296814619761392196084064095930968817137316805100037294689503745488706124666899413964895532933209653306488702079477466498993257038263735203242190432270423835610676865358171473549302841722395495015391276330099513338593481922278385751205360336352009745201657984970966051516246476846265266(225)
もしや、と思い。
! 1/tanh(pi)
OPTION ARITHMETIC DECIMAL_HIGH
LET t0=TIME
LET x=PI
LET SINH1=0
LET COSH1=0
FOR n=0 TO 225 !225まで
LET SINH1=SINH1+x^(2*n+1)/FACT(2*n+1)
LET COSH1=COSH1+x^(2*n)/FACT(2*n)
PRINT n;1/(SINH1/COSH1)
NEXT n
PRINT SINH1/COSH1
PRINT 1/(SINH1/COSH1)
PRINT TIME-t0;"秒"
END
収束の仕方が違うようでした。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2015年11月 5日(木)23時56分44秒
|
|
|
|
> No.3946[元記事へ]
たろささんへのお返事です。
zeta(100)までは、速いと思いました。
! ζ(2n) zeta(x)1000桁
OPTION ARITHMETIC DECIMAL_HIGH
LET t0=TIME
FOR n=2 TO 100 STEP 2 !452まで
LET z=2^(n-1)*ABS(BERNOULLI(n))*PI^n/FACT(n)
PRINT n;z
NEXT n
PRINT TIME-t0;"秒"
END
EXTERNAL FUNCTION BERNOULLI(n) !'ベルヌーイ定数
OPTION ARITHMETIC DECIMAL_HIGH
DIM b(0 TO n)
LET b(0)=1
IF n=0 THEN LET BERNOULLI=1
FOR k=2 TO n+1
LET s=0
FOR i=0 TO k-2
LET s=s+comb(k,i)*b(i)
NEXT i
LET b(k-1)=- s / k
LET BERNOULLI=b(k-1)
NEXT k
END FUNCTION
算出DATA
zeta( 452 )1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000085988814176858736022924359538562907258501404151846766798778241932419782413925225831568060979382126233853707992486822907198580350290928486629026407372314803053733508131169217272750860805236981540665396616343561869002941657513484764112563318094346661647313403371674485076997284643252546537393229349439553279900027671289000817570677825614962942645759830524130307477731757211593662413791300985722530966557407104362199638978980114140200660955834824412210737939838949275010203805787030414698753154772220076211141972692525502559001072360551847142097111239758839507089745929847362853730452865224133201726011390589648358909320665895291010852042959240548123664377734867171422477568010977337488166197910909842401283979521056535019491488660676170103552866931249441209527690484943998422899043066150077961072636452156938609302069463626185009750807888631085481518813682245507032
452 1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000085988814176858736022924359538562907258501404151846766798778241932419782413925225831568060979382126233853707992486822907198580350290928486629026407372314803053733508131169217272750860805236981540665396616343561869002941657513484764112563318094346661647313403371674485076997284643252546537393229349439553279900027671289000817570677825614962942645759830524130307477731757211593662413791300985722530966557407104362199638978980114140200660955834824412210737939838949275010203805787030414698753154772220076211141972692525502559001072360551847142097111239758839507089745929847362853730452865224133201726011390589648358909320665895291010852042959240548123664377734867171422477568010977337488166197910909842401283979521056535019491488660676170103552866931249441209527690484943998422899043066150077961072636452156938609302069463626185009750807888631085481518813682245507048
末尾2桁の誤差ありました。どちらが正しいのか?
452で数値溢れです。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2015年11月 6日(金)13時47分14秒
|
|
|
|
> No.3945[元記事へ]
たろささんへのお返事です。
> 暫定版ですが、出来ました。
手直しして、途中からの計算をしました。
! Cohen 2000 MOD(x,4)=1(-)Ramanujan ζ(2n+1) zeta(x)1000桁精度
OPTION ARITHMETIC DECIMAL_HIGH
PRINT DATE$;"/"; TIME$
LET t0=TIME
FOR n=101 TO 105 STEP 4
LET Na=(N-1)/2
LET f=((2*PI)^n)/(FACT(n+1)*(n-1))
LET S=0
LET w= INT((n+1)/4)
FOR k1=0 TO w
LET S=S+((-1)^(k1))*(n+1-4*k1)*comb(n+1,2*k1)*BERNOULLI(n+1-2*k1)*BERNOULLI(2*k1)
NEXT k1
LET S1=0
FOR k2=1 TO 365 !数値溢れ370まで
LET q=2*PI*k2
LET S1=S1+(EX(q)*((1+(4*PI*k2)/(n-1))-1))/((k2^n)*(EX(q)-1)^2)
NEXT k2
LET zeta=f*S-2*S1
PRINT n;zeta+Ramanujan(Na)
NEXT n
PRINT TIME-t0;"秒"
PRINT DATE$;"/"; TIME$
END
EXTERNAL FUNCTION Ramanujan(n) !zeta(x) MOD(x,4)=3
OPTION ARITHMETIC DECIMAL_HIGH
FOR n=1 TO n
LET f=2^(2*n)*pi^(2*n+1)
LET S=0
FOR k1=0 TO n+1
LET S=S+((-1)^(k1+1))*BERNOULLI(2*k1)/FACT(2*k1)*BERNOULLI(2*n+2-2*k1)/FACT(2*n+2-2*k1)
NEXT k1
LET S1=0
FOR k2=1 TO 360
LET q=2*PI*k2
LET S1=S1+(k2^(-2*n-1))/(EX(q)-1)
NEXT k2
LET Ramanujan=f*S-2*S1
NEXT n
END FUNCTION
EXTERNAL FUNCTION BERNOULLI(n) !'ベルヌーイ定数
OPTION ARITHMETIC DECIMAL_HIGH
DIM b(0 TO n)
LET b(0)=1
IF n=0 THEN LET BERNOULLI=1
FOR k=2 TO n+1
LET s=0
FOR i=0 TO k-2
LET s=s+comb(k,i)*b(i)
NEXT i
LET b(k-1)=-s/k
LET BERNOULLI=b(k-1)
NEXT k
END FUNCTION
!1000桁モードで利用する指数関数
EXTERNAL FUNCTION EX(x)
OPTION ARITHMETIC DECIMAL_HIGH
FUNCTION s(y,n)
LET t=y*x/n
IF ABS(t)<=EPS(0) THEN
LET s=y+t
ELSE
LET s=y+s(t,n+1)
END IF
END FUNCTION
LET EX=s(1,1)
END FUNCTION
計算時間
20151106/10:49:06
50 101 1.00000000000000000000000000000039443045261050590335263935513575963608141044204433323680796761625237491928861445447104761428674299745438441574818957622047157087121949588291127548790096524591503013583366154091609592555734040093489167401537877730680292119177244345062720079530048987294893773720825124823557829467375629675004037826824689502078351153082020190034556854466943862945893574932257594460178334935822767226090976835684592392707422365634881380727204483236981418091615024437298781245092204000055752472644607980496146919583719462572020728942861303262753442706897595473814372098686898267645378670075647357998375763074992690066417983655645433158089100683705978144522084297107303957312401644138041378689506763888443128330237105606859250300409176278264467936209253937799101337731282246825242721017474195558896835960631640886967598186269848060428661060515337152476174048005987058003306858744049169321678254362379586923838911912007463446576063920040030151386246756337944355544528686758422059473828545744162746759
52 105 1.00000000000000000000000000000002465190328815661892710139510328781252773254859693427874476585837833141988834731051880378175626243532251405340305439542034477846314545378655567978140112785348324983780886855766308850997331315595325849404274875554200033271876099547506493928559569661347749040920220284617698459756285299907235960323489037631022207883909373233074479485346380846526453164589760838258265649982040133795192689428043960352287449872071163012854161869583666131500588715228284752482681170211212356737159502128684897403886992986286509036909235009575710743607587868873632989523716970514420613758929899183152917522064922535459889570084454194640071087546111610176471433038740037583760393532387745911806779739726065271018823999809581623114128124221327426547497395047222095610853791743910319009075257991685492529861618941593722311219357538689767376851814030062298069692165312911207185810141620298109690603253051928062249868752743843314262982434832062833658466460497038837578288576444294145687146867275876456996
7775.49 秒
20151106/12:58:41
G3258@4.6GHz(2時間09分35.49秒)
このPCの場合は、Ramanujan ζ(2n+1)と個別に、Cohen 2000 MOD(x,4)=1計算した方が速い。
Ramanujan ζ(2n+1)の探究には、役に立つのでは?
参照
双曲線関数にまつわる重要な公式まとめ
OPTION ARITHMETIC COMPLEX
PRINT 1/TANH(PI)
LET t=(EXP(PI)-EXP(-PI))/(EXP(PI)+EXP(-PI))
PRINT 1/t
END
素人思考では、桁合わせが難しい。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2015年11月 6日(金)17時00分1秒
|
|
|
|
> No.3949[元記事へ]
たろささんへのお返事です。
zeta(x)1000桁プロジェクトからの報告です。
既にご存知の方も多いと思います。
!総和式 ゼータ関数 1000桁 精度確認用
OPTION ARITHMETIC RATIONAL
PRINT DATE$;"/"; TIME$
LET format$="#." & REPEAT$("#",999)
LET t0=TIME
LET f=0
FOR s=1000 TO 1100 STEP 1
LET z=ZETA(s)
PRINT s;
PRINT USING format$:z
NEXT s
PRINT TIME-t0;"秒で計算しました"
PRINT DATE$;"/"; TIME$
END
EXTERNAL FUNCTION ZETA(s)
OPTION ARITHMETIC RATIONAL
LET f=0
FOR n=1 TO 10
LET f=f+1/n^s
NEXT n
LET ZETA=f
END FUNCTION
計算時間は、AMD 2.4GHz 191/25 秒
ゼータ関数1000桁の収束
!総和式 ゼータ関数 1000桁 精度確認用
OPTION ARITHMETIC RATIONAL
PRINT DATE$;"/"; TIME$
LET format$="#." & REPEAT$("#",999)
LET t0=TIME
LET f=0
LET s=3318
LET f=0
FOR n=1 TO 5
LET f=f+1/n^s
PRINT USING format$:f
NEXT n
PRINT s;
PRINT USING format$:f
PRINT TIME-t0;"秒で計算しました"
PRINT DATE$;"/"; TIME$
END
計算結果
zeta(3319)から1になりました。
zeta(350)から、総和式の計算速度と比較しています。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2015年11月 8日(日)18時57分10秒
|
|
|
|
> No.3950[元記事へ]
たろささんへのお返事です。
!zeta(-n)=bernoulli(n+1)/(n+1)*(-1)^n
OPTION ARITHMETIC RATIONAL
LET format$="#." & REPEAT$("#",999)
FOR n=0 TO 18
LET ZETA=BERNOULLI(n+1)/(n+1)*(-1)^n
! PRINT n;ROUND(ZETA,999) !1000桁モード
PRINT n;
PRINT USING format$:ZETA
NEXT n
END
EXTERNAL FUNCTION BERNOULLI(n) !'ベルヌーイ定数 ゼロ有り
OPTION ARITHMETIC RATIONAL
DIM b(0 TO n)
LET b(0)=1
IF n=0 THEN LET BERNOULLI=1
FOR k=2 TO n+1
LET s=0
FOR i=0 TO k-2
LET s=s+comb(k,i)*b(i)
NEXT i
LET b(k-1)=- s / k
LET BERNOULLI=b(k-1)
NEXT k
END FUNCTION
1000桁モードにすると誤差が出てます。
高精度計算サイトで50桁の精度確認済です。以外の精度は、わかりません。
1000桁モード用
!zeta(-n)=bernoulli(n+1)/(n+1)*(-1)^n
OPTION ARITHMETIC DECIMAL_HIGH
FOR n=1 TO 100 !最小値1から
LET ZETA=BERNOULLI(n+1)/(n+1)*(-1)^n
PRINT n;ZETA
NEXT n
END
EXTERNAL FUNCTION BERNOULLI(K) !'ベルヌーイ定数
OPTION ARITHMETIC DECIMAL_HIGH
LET C=1
LET D=1
LET N=K/2
IF K=1 THEN
LET BERNOULLI=-1/2
ELSEIF MOD(K,2)=1 THEN
LET BERNOULLI=0
ELSE
FOR M=N TO 1 STEP-1
LET T=T+(-1)^M*D*M^(K-1)
LET C=C*(N+M+1)/(N-M+1)
LET D=D+C
NEXT M
LET BERNOULLI=-T*K/(2^K*(2^K-1))
END IF
END FUNCTION
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：たろさ
投稿日：2015年11月11日(水)00時47分26秒
|
|
|
|
たろささんへのお返事です。
zeta(x)1000桁プロジェクトからの報告です。
既にご存知の方も多いと思います。
zeta(200)まで
総乗=2.294856591673313794183515831344311288713163799441668673275812601908633963673377415179943750789169776454334944288675243342793042798419633146404801344525793271284796225256839225654003934627958980729155934995390141275548002788589007176859578781539342572938835883667622453728503419451999574796782569303658867125848611781958268232787898845074076388755051637266756596760765141784403430197066202675040752290644360507316565126064063431092179824388570169404898711312512187369473331963755179694075163351070114063217051756089090819479358332277505453595993861256395431731457357236915748084669331619317638954916401036782077744414439695124454751474164135077198630609514253282298056699777616960120676770252945053189318783388261904168119241544869222072832636588448489595602568733878779459397450797700894740462069556116532204037831237781647722200884216608486752406491789678587758579198399992074435417660376057312804250425694945138191204552687543664799283291549121666678014588144991108347459115803455784320701350852558
総和=199.999999999999999999999999999999999999999999999999999999999999377698472213885829285593594621987574058543808328370076748177264429074834467343736155613002853049954766406987052052530699466477416942672141087298454572998606328027081453315387356526592868265612471298697007253574231113158123920717933406539845146257229058667849180175269075401329849405512134222752230480868071869034747095044978953975003151190224566538791889836415993430302583144456807832044456465871197249228702470575669556682067140415657905779244091333102034680748258929166707529973064948232758990230726614295844102598989409688437513498293802313965070119460407689940739443188441342829136217851241471406780908062457198439087010034311408649490079755139835199548450511037806632476709522722977641466786259849958138104540946819393655881743052163157343326828173197330149682511837562017283090894077880173675805111183131022676353309074313737377840775361809823507344793108724630141488759810038894766120522057492686896792557623604100748074735874379376
zeta(300)まで
総乗=2.294856591673313794183515831344311288713163799441668673275814030001397012011323157501796803396705941015955181691237315319390932797243110565830702827099507914024450527801450793132935598259843737981895953078506059390785828004842855373554434809833090449162058662083204292673372850877901452196699665738879321749255317003887024502520201446831253054276001628250050303173842547515521791246154856830518308368669172530927305545904074409145174053822784883642775906084242797127615999416688123988591823015642446736156803746066015789626116554756387591572874384703086548196079388279376439075207165933729489128464021841838019551327176009408059504073925569124556783742636557715159756021867310261343882644058740340946445081521492157578096375615674503612809616525079468313071377230215313367519517712608888878631396436335414335158172172690788219956317961925905430615318618613638377775452659931381515345066103772046288641818412200552264798463496562270245827867857999347010816060417035652319070519781553690649101100975757
総和 =299.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999509090653470227344690422804501372435702478448750055039836316797707379471021325023423532610159426365602811084234619402361713037471941250870170806071108662419575027591424840141243522596107804516555388610315116269376309990821412423586357697032007977981824874443433641898471623311005124749627266421309008813383031153787461156051164222916662657809213404186335212722660996446660515979880630977830114965874590456732667170036137742232592133683555064757833964025894304573699433590278195723719694070302217854614603857679073635128245298300290657042214573768811921077597571096946395315271654018257188947533051076349829222715575919956780251045173586291326795538051276617623210028549840587468254830226521579985571362228999743988241867094673920424539800616992013554907457645383304424242504381050703514486857678284464928458885365777539206808702319482214405581735529567500077632629426864668590098266079648831693265713881710923
zeta(400)まで
総乗=2.294856591673313794183515831344311288713163799441668673275814030001397012011323157501796804523272490813842972172103285178872716733698850029879956992175951803354849702580526171406664612088196522080460021713807500658076502971502998123846098374271874441083854204881413666994952980913778697348334663584695676208482379985317759484894728691987194112408151413797255491784323115645746137230512566990318352441162903500137061395109584770807765602501009953086913765172910320827818187460774978612512697346432703217064458982336811950113382825521023319883524511174058041472944490907604654008010118607524939389511101083423493784729204318364359952993375740411003205098158031913677628729971661439146622708771293425292073472244338062608046201482198269544572732854714216916719926809761395556336171872947495895549704446368252461156498656218153718745073707383018333553844406453845974668268199645750272098840282564416854120022278483900576130934343375599621581238321789462794880065498089534550744352674493236257265793092684
総和=399.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999612740808515068172718196936671364815242978080795120913451223705865558358103138653744440701228090124833504690974267373311679642185932606936232198725087917590704101490254872050546367953573243476756008548762011589827953020990140071713614709858030888206146799517992461058182136047569824211836561612221251501556000430861949423598442065625103325915731139217359428617574963440022768536186789789225343000247621196672163401745704597634676683081301230943053360859940617727454289965042375118012589473131286778236329785713087096438468358950368118243961744781247724964375775223825484220136211418843039804787251459280345748536269699626931830345620910646038049590617221520964936691415618865986894762549962894270603873727686073283435517659142750832941114229753979080938457764023024827977232348243264688770792361173933221703589098088276749361075061951003644286267102463888689509011274545296508787
zeta(500)まで
総乗=2.294856591673313794183515831344311288713163799441668673275814030001397012011323157501796804523272490813842972172103285179761421041964123866454436501348506219322701162204307913788296404169597676409351918353183577082712723968818479540668299936497308576996739386896226092306851913661910023030115409143065452753559830897686760857944312378308010923769793746708738584691938250769850935381571561573850639362654960942845930137881376074420143905641812345775856109987119704129127275522313622206103034229295000824502637393537135641136187105721375749028325893149565390922667069120511449608243932675061197200399549574616386440478160999087556626907742500486581897665164079359080031829651768400223488216619065559912141780771091012539866676122000384491989068732205153590201279972341338689414827283271441231784729903181800379294182696643969747145969240281966806752864090771938259661875428623999950586004058790700093587519478719280386285381905820662675166577936148114859086981723441433866538661703981747791441046284327
総和=499.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999694506363650039531794802060678638230021059725942767333610638609071870837347527954229814262737933696158622394953145552724263372055190464667873686477604891653073688633826243926778333450838122636126261447784636405785760920583449475232711932323202030517513084057738040683011645528293349124597188537731080658195318894528639879393918554480688908125177630133418151052985636753678524911639986676859671652184790688582924602168977768023097235929419181491546347903836200394409640270252814725520702898445326252050062064845170225642561784621496618539552808307023966458469089321824121223947551995989550638623731169350453835230113165714670337705946079180782060593132302088963404280257834230053691103948575320171042456633577963887590632261366985458262818299367832446863599708216428500556014747311538712679768321162456121815199614024634953306944738252398901699444867
zeta(600)まで
総乗=2.294856591673313794183515831344311288713163799441668673275814030001397012011323157501796804523272490813842972172103285179761421041964123866454436501349207283407793119246747769333049514157291089213694156204030549107737221017381583609737389583416962074236594260680740123692022740876352061924221918258588989829729945598383539557342895394265881193029867444328902140230642809303882781214774465047538193558068079680385477442426556727265544515254203033312996658260406524824678094275470195611719432943092720444193205971298147144523802427958458694391534341254893773778936352393323779196951396582640742638096083163269342710317945998177448393309240433190835013453027807889180713626393780731753526107459448507586156495688869599242224616696348641475282673796289755703577358692931806446147698520142597383304415728755854604387235534238022705769200935620099633400982537846389217090847328214524488807858129126381134322322116039983484177888594968272855870193639332817520865072195925159757695001271605889054733186538626
総和=599.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999759008013489711588225924996528749106356899504549010202816970002284119614713973896252431840230458492109586536085192921868675641756743967674476654276974757524698953588202249546150870055684971418374392177836853552462809076239753387205398493097588668626448356720023133776037511569227350956286650584099184280186575362228300435675593698995637106254038031396465826683622209059900615116427784283276379617890544933697489671956281881508949475813505921740540178141037534753078077076484376559466607893160286054607299538474224332120756591487134254163586668488844257448538723319873952021324355088395009459356684610755784424406990233228755117243517062861156669421152096893856802781734850625154036679409798182826009063381442278688900791432133067549963176758289742545589630426144728528081610507402986525993135411929761306222885955810445
zeta(700)まで
総乗=2.294856591673313794183515831344311288713163799441668673275814030001397012011323157501796804523272490813842972172103285179761421041964123866454436501349207283407793119246747769333050067199339872795374031927010510608797070381312717040492764750353609688327882844299973351061520894663314644660172723760116545381160804553890999584913408969793951655778882847375165054923885828339244666717197251415245961933741372455389208201749534605830906770462786114202962198897423994610043949731485889353012610952138732951608862299298899283065376655959894912529278664966735884637490358342855881544281410326686011246787755230091901342765676668499946365281107959883862193744946265861799461403336367618031436482951778217287057220732370114639921150568254210829767977493091108844316077494975687292484552887810584907431641040750796979620175677568250869358510681075474897883013155499954399175278186121017298080409165063271373647364049779087243523757341984704071456594222600319854046679448196253154746134797464586057263663945375
総和=699.99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999980989084337048401764849275941648968907351287936264809675582539424341457507225277964427985029473781367973893704967018113252212198400760474248543592111209692917171142905745622645007786535919203609876055447642302009159138973721637811460640469887694353219438396914949910672717455639435561889260045611336835942649276405562177200138618098306967198883776575966891115032929761655760975473420649280252675798964461282817559553763405218508394462926674733394619804893850196902200828289613660448796382496328109870549243176554461371601216232529471336073371544723652181868529442145562325957420476541996606370343793989168168532091531694758969498044288896097627751304851940091405097628943139630976898382341933392494547799443027417055258999895259510140383668644337560938050895785517857891683532258817535092
zeta(800)まで
総乗=2.294856591673313794183515831344311288713163799441668673275814030001397012011323157501796804523272490813842972172103285179761421041964123866454436501349207283407793119246747769333050067199339872795374031927010511045070321610011931838936532973062668842601824348895677245701698531343947365225018889429546000927949415058336221462209025315033642448811440846431928312725821455594172065304077171866852006660176164036160581341595623605467149996952423045553321000222523806973764162452504872438908742392680504920922357778209509476764217384709276965555584498661001089834298803467124875095568392122960152435969852357391037599946536198418338398000180764019204231576044157028304223153621928105037945775320471642371702069115938195323596982286153632382974666746711658313330513015410691607535335241499661768025600837458473393862555564152252388668119116167401144398950027668818067391910055249468615116578870077509826892817868423187110266485988241066521739139479969649595643189084583736928431001550475282566882791898278
総和=799.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999850030318610436904518235556237193464646003830376089170206266555238228914414780969722903187160258516375751035237199413900461635849736852977083761965695190557083973424631012045519742907135061467558345162092532101523346575860045557112857479499240062255800439256955096012619774555912496017969298235291220005630049957261915659845718722041329281744999433069828208579337473025630045797260323661024389819776131149162104014041524244474753639989323702384588910228093783980474484171778313728822880185123057684261653724538067936663617897393967979229974607077280282311980159522039056675915250526003281227587983477016730289770288460805074958966854142879905888870626145281970740953314460855619440256787219314592941074723535341721102075923106950002101249219702630951070140718
zeta(900)まで
総乗=2.294856591673313794183515831344311288713163799441668673275814030001397012011323157501796804523272490813842972172103285179761421041964123866454436501349207283407793119246747769333050067199339872795374031927010511045070321610011931838936532973406827754489809912361266430357810692053814457136526473793329384572458371141029217368988361154920331547169876886057992535733219125036683402416519807986753945834568763376490570969889766016801054213951929538866782538656566669315189774039711425097809978242719950552174551627401350290714936080592407940871290985863516406308162541310582779723908338018253818736748367352409518589882966100990800234158933303609008189547665506551570642704304233171343798755279961917973614051185023591767033846783894345207927548473642421066739999004686525994753593569306238821178904432997582116173618249270285966280865201344951359026072298151819420896679207520476559274678073159506366129801672839832253609401492753106989353348708785390788521357003190629371687145944087117078872443974855
総和=899.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999881694781383322528902724840248197346834781781609522782259683729835579217136708079878174604560934189070008394739144500924872154059681718499589074521955059393023883779986407862394371109688675192533297902304294334784560256806705113004202170228172509488803865271982266603708820713401517915438321843524658406593789420170921022498967946801378755859222636724485297349956005944569821161870075111832909013790307339460548275795752874308205055721298623020940431709875078123329956378139827680726327002398939083806402778469712498381866521613844178883649095033891405547677236342063559096572753474460936499306393301708464061053560923997692311761857731094081920789691793896032414603593890799515335414402831152349806195237996977512243143080505919
zeta(1000)まで
総乗=2.294856591673313794183515831344311288713163799441668673275814030001397012011323157501796804523272490813842972172103285179761421041964123866454436501349207283407793119246747769333050067199339872795374031927010511045070321610011931838936532973406827754489809912361266430358082185564586291869213773116111814384314428154465151563602245330921298004823817633581043446569877118787427222065741493757711392788592692964234419668491537172604036817707964701339953365441587081866824255401603645839591466345733362109320128683767324075237098565117279966184443893216837906350328692282051006130524712407630013856574275302800677040281427120379848143226398041526305819296421892086695519652402309109445460022483589191461306063945446412179376406334836423961115474925425870380879933165507302580783750188735363181555231501184606759476912180171123424220091429858903011491736999729204077588580738417533302576740015806161283933867828766292505659417938616444955083743219122570680808848179190723966682239016272360161862681255581
総和=999.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999906673638149678112100991045527618283038290855362829197537828566020403308902422436554555967290211889764040501006967575737578451247864596760515847918279606924376558933386167484972222297835254611858652078425503202106355736504239930448643478856303903602811798138270099608350776882949990518153631727249453450171310106247975508602109606369950736137520811883145010247961056762466133830425616845900330648793735024943932369490012367613897114196888734207352480031983501509272659129938788056149418719540328070753781423571822439439227681257796179649445158356509625118283763599892930621748994657326721776669516255338337628845212049808084760695428372823165021135274582890815192137916819203955986697097690999697139
精度は1000桁末尾数桁は色々でした。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
戻る
