昔ＨＳＰコンテストで円周率計算機なるものを投稿したと言いましたよね。
いまさらですが公開したいと思います。
また素数計算プログラムもおまけでついてます。
またこのプログラムを実行するにはlongint.dllが別途必要です
「多倍長 hsp」でググってダウンロードしてきてください
ではのっけます。
#ifndef __longint__
#define __longint__
#regcmd "_init@4","longint.dll",1
#cmd ExEuclid 0x100
#cmd StoreLongInt 0x101
#cmd GetPowMod 0x300
#cmd GetGCD 0x301
#cmd IsPrime 0x302
#cmd LoadLongInt 0x304
#cmd GetLintSize 0x305
#cmd LongInt 0x800
#endif
screen 0,120,64,0
objsize 120,32
button "円周率を計算する",*enpir
button "素数を検索する",*sosuu
stop
*sosuu
keta=200000000
screen 0,260,80,0
boxf
color 255,255,255
sdim hi,1200
repeat 50,1
hi+=str(cnt)+"\n"
loop
mes "処理優先率（数が大きいほど早い）"
combox de,200,hi
sdim hi,100000
hi=""
exist "素数表.txt"
gt=strsize
if gt<30{
hi="2\n"
wait 5
bsave "素数表.txt",hi,3
f=1
sdim hi,100000
hi=""
gt=3
}else{
sdim r,20
bload "素数表.txt",r,20,gt-20
f=0
repeat 18
if peek(r,cnt)=10:f=cnt
loop
f=int(strmid(r,f+1,18-f-1))
}
wait 10
dim r,1
r=0
onexit*enf
repeat keta
f+=2
if IsPrime(Longint(f))=1:hi+=str(f)+"\n":r+=10
if r>99980{
bsave "素数表.txt",hi,strlen(hi),gt
gt+=strlen(hi)
a+
title ""+a+""
r=0
sdim hi,100000
hi=""
wait 2
}
if cnt\(de*16+1)=0:await 1:color 0,0,0:boxf 0,20,260,80:pos 20,50:color 255,255,255:mes f
loop
onexit 0
wait 60
if hi!"":bsave "素数表.txt",hi,strlen(hi),gt
wait 60
end
*enf
dialog "保存する?",3:if stat=6:if hi!"":bsave "素数表.txt",hi,strlen(hi),gt
end
*enpir
π=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986
buffer 1,400,60
repeat 60
color -cnt*3/2+110,-cnt*3/2+110,-cnt*3/2+110
line 200,cnt,400,cnt
loop
buffer 2,12,12,0:boxf
repeat 6
color cnt*28+50,cnt*28+50,cnt*28+50
circle cnt,cnt,12-cnt,12-cnt,1
loop
screen 0,200,120,0:boxf
title "高速電子式計算機"
font "ＭＳ 明朝",16,16
pos 10,60
objsize 64,32
combox けたけた,200,"5120桁\n10240桁\n20480桁\n40960桁\n81920桁\n163840桁\n327680桁\n655360桁\n1310720桁"
pos 100,60
button "計算実行",*start
虚=0.0
無=0.0
*main
redraw 1:redraw 0
color 0,0,0:boxf
color 160,160,160:gmode 3,12,12,150
if 無=0:無=(rnd(2)*2-1)*(rnd(200)+30):正=36+rnd(30)
if 正=0{
無-=無/abs(無)
負+=無/10
}
if 正!0:正-
repeat 8
虚=π/4.0*cnt+0.01*負
circle sin(虚)*40+92,cos(虚)*40-3,sin(虚)*40+108,cos(虚)*40+13,0
pos sin(虚)*40+94,cos(虚)*40-1:gcopy 2,0,0,12,12
loop
repeat 8
虚=π/4.0*cnt+0.01*負
circle sin(-虚)*40+92,cos(-虚)*40+107,sin(-虚)*40+108,cos(-虚)*40+123,0
pos sin(-虚)*40+94,cos(-虚)*40+109:gcopy 2,0,0,12,12
loop
gsel 1:pos 0,0:gcopy 0,0,60,200,60:gsel 0
color 0,0,0:boxf 0,60,200,120
gmode 7,200,60,156:pos 0,60:gcopy 1,0,0,200,60
color 255,255,255
pos 1,1
mes "求めるπの精度(桁数)を\n 選択してください。"
wait 2
goto*main
*start
gsel 0
redraw 1
clrobj
color 0,0,0:boxf
color 255,255,255:pos 0,0
mes "計算中....."
けたけた+=9
時間=gettime(4)*60000*60+gettime(5)*60000+gettime(6)*1000+gettime(7)
a=Longint(1)
r=longint(10)
dimtype g,longint,けたけた
repeat 10
a=a*10
loop
repeat けたけた
r*=2
g.cnt=a
a*=a
loop
a=g
けた=int(r)
title ""+けた+"桁を計算中"
font "ＭＳ 明朝",12,16
color 155,155,155:line 10,100,190,100
line 10,96,10,104
line 100,96,100,104
line 190,96,190,104
pos 0,108:mes "0%"
pos 86,108:mes "50%"
pos 170,108:mes "100%"
color 255,255,255
chuiyutxdihubyxty=けた/14-1
rtuytrbtyryuvbyts=けたけた
gosub*sq
t=longint("1")
t=y:t=t*426880:t*=a
n=けた/14-1
n=Longint(n)
x=Longint("1")
repeat 3
x*=32
x*=5
x*=23
x*=29
loop
x*=-9
m=Longint("1")
m=2*n-1
y=longint(6*n-1)
y*=m:y*=(6*n-5)
m=n*n*n
c=Longint("545140134")
c*=n:c+=Longint("13591409")
c*=a
p=Longint("545140134"):p*=a
g=longint("0")
;gは多長倍の乱数
;cは０の多長倍（１０桁）
;ｍは短い多長倍（けたけたに比例）０に収束
;xは短い多長倍（２０桁）不変
;yは短い多長倍（１０桁）０に収束
;ｎは整数、＝CNT
;ｐは０の多長倍、不変（１０桁）
repeat int(n)
g+=c
g/=x*m
g*=y
m-=n*3*(n-1)+1
y-=longint(2)*(n*108*(n-2)+113)
n-
c-=p
boxf 10,77,100+90*(cnt+1)/chuiyutxdihubyxty,87
loop
g+=c
x=t/g
sdim k,けた+100
k=str(x)
k="π=3.\n"+strmid(k,1,int(r))
bsave "出力π.txt",k,int(r)+7
時間=gettime(4)*60000*60+gettime(5)*60000+gettime(6)*1000+gettime(7)-時間
color 0,0,0:boxf
color 255,255,255:pos 0,0:sdim estg,2
if 時間\1000<100:estg="0":if 時間\1000<10:estg="00"
mes "πの計算が終わりました。"
mes ""+時間/1000+"."+estg+""+時間\1000+"秒かかりました\n\n「π.txt」に出力しました"
stop
*sq
p=Longint("10005"):p=p*a
x=longint(0)
y=longint(0)
x2=longint("1000249968758")
x=a*longint("10000000000")/x2
y=a*x2/longint("10000000000")
repeat けたけた
a*=g.cnt:p*=g.cnt:x*=g.cnt:y*=g.cnt
boxf 10,77,10+30*(3*cnt+1)/rtuytrbtyryuvbyts,87
x+=(x*(a-x*y/a))/a
boxf 10,77,10+30*(3*cnt+2)/rtuytrbtyryuvbyts,87
y+=(x*(p-y*y/a)/2)/a
boxf 10,77,10+30*(3*cnt+3)/rtuytrbtyryuvbyts,87
loop
return
