仮眠プログラマーのつぶやき : １億桁×１億桁

自作ゲームやツール、ソースなどを公開しております。
①ポンコツ自動車シュライシュラー DOWNLOAD	②流体力学ソース付き汚いほうDOWNLOAD 綺麗なほうDOWNLOAD	③ミニスタヲズ DOWNLOAD
④地下鉄でGO DOWNLOAD	⑤ババドン DOWNLOAD	⑥圧縮拳（ツール） DOWNLOAD
⑦複写拳 DOWNLOAD	⑧布シミュレーション DOWNLOAD	⑨minecraft巨大電卓地形データ DOWNLOAD
⑩フリュードランダーデジゲー博頒布α版 DOWNLOAD	⑪パズドラルート解析GPGPU版 DOWNLOAD	⑫ゲーム「流体de月面着陸」 DOWNLOAD

１億桁×１億桁

2013年08月12日02:04

カテゴリ: GPGPU; １億桁×１億桁

円周率プロジェクト途中経過

なんかFFTよりも多倍長計算に向いている高速剰余変換（FMT）というものがあるらしい

で早速やってみた。

回転子であるωを4とし、すべての計算を257の剰余下で行う。
変換する全要素は8つであり、4^8を257で割った余りは1かつ、4^(8/2)を257で割った余りは-1(=256)であることがミソ

以下のソースでは4桁の10進数の数同士を乗算し結果を表示している

FFTを用いた計算よりもFMTが優れている点として、全て整数値同士の計算で完結するということがあげられる。
これは非常に重要だ。
FFTでは計算したい桁数を上げれば上げるほど、浮動小数点の誤差が蓄積してしまう。
doubleの精度を持ってしてでも1億桁×1億桁をやろうとすれば、double１要素につき4桁が限界だ。
これはdoubleつまり8byteに、10進数4桁≒1.66byteしか格納できなく、記憶するメモリ容量的にも無駄が多い。当然桁数あたりの必要な演算量も多くなってくる。
FMTでは4byte整数に約1.5～2byteは敷き詰められるしdoubleと違って演算スピードも早い。

誤差の話に戻るが、結果的にFFTを用いた多倍長の計算では1兆桁付近で精度の限界が生じてしまう。一方FMTは全て整数で計算するので理論上無限大まで誤差の蓄積なく計算することが可能だ。

というわけでFMT試作verのソースをUPする。
~~正直醜いスパゲッティなのであまり聞かないで・・・・~~
2014/2/3　YSRさんのご好意によりきれいなソースとなりました！ありがとうございます

#const beki 3

#const LIST_SIZE 1 << beki

#const LIST_SIZE2 LIST_SIZE >> 1

#const LIST_SIZE3 LIST_SIZE << 1

#const LIST_SIZE8 LIST_SIZE << 2

#const n LIST_SIZE

#const p 257

#const pn p/n

#const omega 4

#module

#deffunc fft array A, int inv

if(inv != 0){

dim c, LIST_SIZE@

repeat LIST_SIZE@ ;ビット逆順

ic = cnt / 65536

ic = (ic & 0x00005555) << 1 | (ic & 0x0000AAAA) >> 1

ic = (ic & 0x00003333) << 2 | (ic & 0x0000CCCC) >> 2

ic = (ic & 0x00000F0F) << 4 | (ic & 0x0000F0F0) >> 4

ic = (ic & 0x000000FF) << 8 | (ic & 0x0000FF00) >> 8

iic = cnt \ 65536

iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1

iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2

iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4

iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8

ic = (iic * 32768 + ic / 2) >> (31 - beki@)

C(cnt) = A(ic)

loop

memcpy a,c,LIST_SIZE8@

}

repeat beki@ ;fft

bekii = 1 << cnt

repeat LIST_SIZE2@

t2 = cnt \ bekii

t0 = (cnt / bekii) * bekii * 2 + t2

t1 = t0 + bekii

w = jouyo(omega@, t2*LIST_SIZE2@ / bekii, p@)

r4 = A(t1) * w \ p@

A(t1) = (A(t0) + p@ - r4) \ p@

A(t0) += r4

A(t0)\ = p@

loop

return

#deffunc ufft array A, int inv

repeat beki@ ;fft

bekii = LIST_SIZE2@ >> cnt

repeat LIST_SIZE2@

t2 = cnt \ bekii

t0 = (cnt / bekii) * bekii * 2 + t2

t1 = t0 + bekii

w = jouyo(omega@, t2 * LIST_SIZE2@ / bekii, p@)

r4 = A(t1)

A(t1) = (A(t0) - r4 + p@) * w \ p@

A(t0) += r4

A(t0) \= p@

loop

if(inv != 0){

dim c, LIST_SIZE@

repeat LIST_SIZE@ ;ビット逆順

ic =cnt / 65536

ic = (ic & 0x00005555) << 1 | (ic & 0x0000AAAA) >> 1

ic = (ic & 0x00003333) << 2 | (ic & 0x0000CCCC) >> 2

ic = (ic & 0x00000F0F) << 4 | (ic & 0x0000F0F0) >> 4

ic = (ic & 0x000000FF) << 8 | (ic & 0x0000FF00) >> 8

iic = cnt \ 65536

iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1

iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2

iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4

iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8

ic=(iic * 32768 + ic / 2) >> (31 - beki@)

C(cnt) = A(ic)

loop

memcpy a, c, LIST_SIZE8@

}

return

#deffunc ifft array A,int inv

if(inv != 0){

dim c, LIST_SIZE@

repeat LIST_SIZE@ ;ビット逆順

ic =cnt / 65536

ic = (ic & 0x00005555) << 1 | (ic & 0x0000AAAA) >> 1

ic = (ic & 0x00003333) << 2 | (ic & 0x0000CCCC) >> 2

ic = (ic & 0x00000F0F) << 4 | (ic & 0x0000F0F0) >> 4

ic = (ic & 0x000000FF) << 8 | (ic & 0x0000FF00) >> 8

iic = cnt \ 65536

iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1

iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2

iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4

iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8

ic=(iic * 32768 + ic / 2) >> (31 - beki@)

C(cnt) = A(ic)

loop

memcpy a, c, LIST_SIZE8@

}

repeat beki@ ;fft

bekii = 1 << cnt

repeat LIST_SIZE2@

t2 = cnt \ bekii

t0 = (cnt / bekii) * bekii * 2 + t2

t1 = t0 + bekii

w = jouyo(omega@, LIST_SIZE@ - t2 * LIST_SIZE2@ / bekii, p@)

r4 = A(t1) * w \ p@

A(t1) = (A(t0) + p@ - r4) \ p@

A(t0) += r4

A(t0) \= p@

loop

return

#deffunc iufft array A,int inv

repeat beki@;fft

bekii=LIST_SIZE2@>>cnt

repeat LIST_SIZE2@

t2 = cnt \ bekii

t0 = (cnt / bekii) * bekii * 2 + t2

t1 = t0 + bekii

w = jouyo(omega@, LIST_SIZE@ - t2 * LIST_SIZE2@ / bekii, p@)

r4 = A(t1)

A(t1) =(A(t0) -r4 + p@) * w \ p@

A(t0) += r4

A(t0) \= p@

loop

if(inv != 0){

dim c,LIST_SIZE@

repeat LIST_SIZE@ ;ビット逆順

ic =cnt / 65536

ic = (ic & 0x00005555) << 1 | (ic & 0x0000AAAA) >> 1

ic = (ic & 0x00003333) << 2 | (ic & 0x0000CCCC) >> 2

ic = (ic & 0x00000F0F) << 4 | (ic & 0x0000F0F0) >> 4

ic = (ic & 0x000000FF) << 8 | (ic & 0x0000FF00) >> 8

iic = cnt \ 65536

iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1

iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2

iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4

iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8

ic=(iic * 32768 + ic / 2) >> (31 - beki@)

C(cnt) = A(ic)

loop

memcpy a, c, LIST_SIZE8@

}

return

#deffunc narasi array A,int inv

repeat LIST_SIZE@

arhgea = a.cnt \ n@

a.cnt /= n@

if(arhgea != 0) :a.cnt += (n@ - arhgea) * pn@ + 1

loop

return

#defcfunc jouyo int a, int inv, int cc

ninv=inv

m = 1.0 * cc

x = (1.0 * a) \ m

jouyoo = 1.0

repeat

/* n&1 はnが奇数(bitの1桁目が1)の時1、偶数(bitの1桁目が0)の時0を返す */

if(int(ninv) \ 2 == 1) :jouyoo = (jouyoo * x) \ m

ninv = ninv / 2

x = (x * x) \ m

if(ninv == 0) :break

loop

if(jouyoo < 0) :jouyoo += m

return int(jouyoo)

#global

dim aa, n

aa = 6,5,9,3,0,0,0,0

dim bb, n

bb = 7,8,6,5,0,0,0,0

dim cccc, n

dim dddd, n

memcpy cccc, aa, n*4

memcpy dddd, bb, n*4

;aa.1 *= -2

;aa.2 *= 4

ufft bb

ufft aa

repeat 8

bb.cnt *= aa.cnt

bb.cnt \= p

loop

ifft bb

narasi bb

repeat 7

bb.(cnt+1) += bb.cnt / 10

bb.cnt \= 10

loop

pos 290, 0

repeat 4

mes cccc.cnt

pos ginfo_cx - 12, ginfo_cy - 18

loop

mes "*"

pos ginfo_cx-12, ginfo_cy-18

repeat 4

mes dddd.cnt

pos ginfo_cx - 12,ginfo_cy - 18

loop

mes "="

pos ginfo_cx - 12, ginfo_cy - 18

repeat 8

mes bb.cnt

pos ginfo_cx - 12, ginfo_cy - 18

loop

pos 0, 40

mes jouyo(400, 1024 * 64, 164 * 1024 * 128 + 1)

stop

ここで、FFTに時間間引きや周波数間引きの理論を取り入れることで、変換→逆変換の際2回あるビット逆順の計算を省くことが可能となる。
これも将来的にOpenCLへ落としこむ時、ランダムアクセスを減らしてパフォーマンス低下を防ぐ大きな意味を持つ。

さらにRADEON HD 7900シリーズ・・いわゆるGCNと言われるタイプのGPUコアは整数演算が非常に早い！
これでさらに面白くなってきそうだ！

もともとdouble精度の計算ができないGPUでも、mad24の機能で高速に整数演算ができる物が多いため、そういった点でもFFTよりFMTが有利だ！

というわけで次なるステップは、HSPCLへの移植だな

2019/10追記
FMTをPythonとCUDAで実装しました。ガウスルジャンドルで3億桁まで求めた記事をQiitaに公開しました。(この量だとlivedoorブログだと文字制限で圧倒的にオーバーする・・・)
https://qiita.com/Red_Black_GPGPU/items/e933e0d846b874b86c32
ソースはこちら
https://github.com/toropippi/FMT_QiitaSample

2012年12月20日03:44

カテゴリ: HSP; １億桁×１億桁

HSPでFFT(高速フーリエ変換)モジュール作った、あと多倍長乗算も

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

#const beki 3;2^bekiこの配列が作られます、例beki=23の時　約400万配列×400万配列が800万配列に格納されます

#const kisuu 10000;1～32768まで、浮動小数配列ひとつに格納できる上限数＝基数、何進数かという話、2^beki*log10底のkisuu が乗算後の値の桁数に値する

#module

#const beki beki@

#deffunc fft初期設定

n=1

repeat beki;桁数計算

n*=2

loop

ddim qcosq,n/2+1;sin cos予め計算格納

ddim qsinq,n/2+1

n2=n/2

repeat n/2+1

itjrad=3.1415926535897932384626*cnt/n2

qcosq.cnt=cos(itjrad)

qsinq.cnt=sin(itjrad)

loop

qcosq.n2=-1.0

qsinq.n2=0.0

itjn=1;バタフライインデックス作成

dim inind1,itjn

dim inind2,itjn*2

repeat beki

repeat itjn:inind2.cnt*=2:loop

dim inind1,itjn

memcpy inind1,inind2,4*itjn,0,0

dim inind2,itjn*2

memcpy inind2,inind1,4*itjn,0,0

repeat itjn:inind1.cnt+:loop

memcpy inind2,inind1,4*itjn,4*itjn,0

itjn*=2

loop

dim inind1,1;2の方に結果が

itjn=0

ddim chasr,n;キャッシュ

ddim chasi,n;キャッシュ

return

#deffunc fft array reitjo,array imitjo

q0q=1

q1q=n/2

repeat beki-1

qct0q=0

qct1q=q0q

qcn2q=0

repeat n/2

dup reitjoqccnt1,reitjo(inind2.qct0q)

dup imitjoqccnt1,imitjo(inind2.qct0q)

dup reitjoqccnt2,reitjo(inind2.qct1q)

dup imitjoqccnt2,imitjo(inind2.qct1q)

ita=reitjoqccnt2*qcosq.qcn2q-imitjoqccnt2*qsinq.qcn2q

itb=imitjoqccnt2*qcosq.qcn2q+reitjoqccnt2*qsinq.qcn2q

reitjoqccnt2=reitjoqccnt1-ita

imitjoqccnt2=imitjoqccnt1-itb

reitjoqccnt1+=ita

imitjoqccnt1+=itb

qct0q+

qct1q+

qcn2q+=q1q

if qct0q\q0q=0:qct0q+=q0q:qct1q+=q0q:qcn2q=0

loop

q0q*=2

q1q/=2

loop

qct0q=0

qct1q=q0q

qcn2q=0

repeat n/2

dup reitjoqccnt1,reitjo(inind2.cnt)

dup imitjoqccnt1,imitjo(inind2.cnt)

dup reitjoqccnt2,reitjo(inind2.qct1q)

dup imitjoqccnt2,imitjo(inind2.qct1q)

ita=reitjoqccnt2*qcosq.cnt-imitjoqccnt2*qsinq.cnt

itb=imitjoqccnt2*qcosq.cnt+reitjoqccnt2*qsinq.cnt

chasr.cnt=reitjoqccnt1+ita

chasi.cnt=imitjoqccnt1+itb

chasr.qct1q=reitjoqccnt1-ita

chasi.qct1q=imitjoqccnt1-itb

qct1q+

loop

memcpy reitjo,chasr,8*n,0,0

memcpy imitjo,chasi,8*n,0,0

return

#deffunc ifft array reitjo,array imitjo

q0q=1

q1q=n/2

repeat beki-1

qct0q=0

qct1q=q0q

qcn2q=n2

repeat n/2

dup reitjoqccnt1,reitjo(inind2.qct0q)

dup imitjoqccnt1,imitjo(inind2.qct0q)

dup reitjoqccnt2,reitjo(inind2.qct1q)

dup imitjoqccnt2,imitjo(inind2.qct1q)

ita=reitjoqccnt2*qcosq.qcn2q-imitjoqccnt2*qsinq.qcn2q

itb=imitjoqccnt2*qcosq.qcn2q+reitjoqccnt2*qsinq.qcn2q

reitjoqccnt2=reitjoqccnt1+ita

imitjoqccnt2=imitjoqccnt1+itb

reitjoqccnt1-=ita

imitjoqccnt1-=itb

qct0q+

qct1q+

qcn2q-=q1q

if qct0q\q0q=0:qct0q+=q0q:qct1q+=q0q:qcn2q=n2

loop

q0q*=2

q1q/=2

loop

qct0q=0

qct1q=q0q

qcn2q=n2

repeat n/2

dup reitjoqccnt1,reitjo(inind2.cnt)

dup imitjoqccnt1,imitjo(inind2.cnt)

dup reitjoqccnt2,reitjo(inind2.qct1q)

dup imitjoqccnt2,imitjo(inind2.qct1q)

ita=reitjoqccnt2*qcosq.qcn2q-imitjoqccnt2*qsinq.qcn2q

itb=imitjoqccnt2*qcosq.qcn2q+reitjoqccnt2*qsinq.qcn2q

chasr.cnt=(reitjoqccnt1-ita)/n

chasi.cnt=(imitjoqccnt1-itb)/n

chasr.qct1q=(reitjoqccnt1+ita)/n

chasi.qct1q=(imitjoqccnt1+itb)/n

qct1q+

qcn2q-

loop

memcpy reitjo,chasr,8*n,0,0

memcpy imitjo,chasi,8*n,0,0

return

#global

screen 0,480,700

fft初期設定

randomize

n=1

repeat beki

n*=2

loop

ddim rex,n

ddim rey,n

ddim reitjo,n;実数代入先

ddim imitjo,n;虚数代入先

ddim chasr,n;キャッシュ

ddim chasi,n;キャッシュ

; gosub*indexsakusei

repeat n/2;「A＊B＝ｃ」のAとBにランダムな値を格納

rex.cnt=1.0*rnd(kisuu);2.7*sin(3.1415926535*490*cnt/n)+2.0*sin(3.1415926535*190.0*cnt/n)

loop

repeat n/2

rey.cnt=1.0*rnd(kisuu);1.0*sin(3.1415926535*50.0*cnt/n)+20.0*sin(3.1415926535*90.0*cnt/n)

loop

; aatim=gettime(4)*3600000+gettime(5)*60000+gettime(6)*1000+gettime(7)

; mes aatim

memcpy reitjo,rex,8*n,0,0

memcpy imitjo,rey,8*n,0,0

fft reitjo,imitjo

memcpy chasr,reitjo,8*n,0,0

memcpy chasi,imitjo,8*n,0,0

ReTm=0.0

imTm=0.0

ReTN=0.0

imTN=0.0

ReTm=chasr

ImTm=chasi

ReTN=chasr

ImTN=chasi

reitjo=(Retm*ImTm+ReTN*ImTN)/2.0

imitjo=( ((ReTN+ImTN)*(ReTN-ImTN)) + ( (ImTm+ReTm)*(ImTm-ReTm) ) )/4.0

itjn=n-1

repeat n-1,1

ReTm=chasr.cnt

ImTm=chasi.cnt

ReTN=chasr.itjn

ImTN=chasi.itjn

reitjo.cnt=0.5*(Retm*ImTm+ReTN*ImTN)

imitjo.cnt=0.25*( (ReTN+ImTN)*(ReTN-ImTN) + (ImTm+ReTm)*(ImTm-ReTm) )

itjn-

loop

ddim chasr,1;キャッシュ

ddim chasi,1;キャッシュ

ifft reitjo,imitjo

; mes gettime(4)*3600000+gettime(5)*60000+gettime(6)*1000+gettime(7)

; mes gettime(4)*3600000+gettime(5)*60000+gettime(6)*1000+gettime(7)-aatim

repeat n

if absf(reitjo.cnt-int(reitjo.cnt))<0.7:if absf(reitjo.cnt-int(reitjo.cnt))>0.3:mes reitjo.cnt

; if absf(imitjo.cnt-int(imitjo.cnt))<0.7:if absf(imitjo.cnt-int(imitjo.cnt))>0.3:mes imitjo.cnt

loop

;: bsave "reitjo",reitjo

; bsave "imitjo",imitjo

; stop

dim kekka,n;計算結果入るint型変数

gosub*繰り上げ計算

goto*hyouji2

*繰り上げ計算

w=n-1

w1=n-2

repeat n-1

kekka.w+=int(reitjo.w+0.5)

if kekka.w>=kisuu:kekka.w1+=kekka.w/kisuu:kekka.w\kisuu

w1-

loop

kekka.w+=int(reitjo.w+0.5)

return

*hyouji;波形表示

color 0,0,0

repeat n

a=reitjo.cnt

if cnt=0:pos 0,int(240.0-a)

line cnt*640/n,int(240.0-a)

loop

color 255,0,0

repeat n

a=imitjo.cnt

if cnt=0:pos 0,int(240.0-a)

line cnt*640/n,int(240.0-a)

loop

color 0,255,0

repeat n

a=(imitjo.cnt*imitjo.cnt+reitjo.cnt*reitjo.cnt)

if cnt=0:pos 0,int(240.0-a)

line cnt*640/n,int(240.0-a)

loop

stop

*hyouji2

font "MS ゴシック",10,1

pos 0,0

mes "掛け算結果配列表示"

repeat limit(n,1,300)

if strlen(str(kekka.cnt))=1:mes "000"+str(kekka.cnt)

if strlen(str(kekka.cnt))=2:mes "00"+str(kekka.cnt)

if strlen(str(kekka.cnt))=3:mes "0"+str(kekka.cnt)

if strlen(str(kekka.cnt))>=4:mes str(kekka.cnt)

loop

pos 150,0

mes "掛ける数A配列表示"

repeat limit(n,1,300)

if strlen(str(int(rex.cnt)))=1:mes "000"+str(int(rex.cnt))

if strlen(str(int(rex.cnt)))=2:mes "00"+str(int(rex.cnt))

if strlen(str(int(rex.cnt)))=3:mes "0"+str(int(rex.cnt))

if strlen(str(int(rex.cnt)))>=4:mes str(int(rex.cnt))

loop

pos 300,0

mes "掛ける数B配列表示"

repeat limit(n,1,300)

if strlen(str(int(rey.cnt)))=1:mes "000"+str(int(rey.cnt))

if strlen(str(int(rey.cnt)))=2:mes "00"+str(int(rey.cnt))

if strlen(str(int(rey.cnt)))=3:mes "0"+str(int(rey.cnt))

if strlen(str(int(rey.cnt)))>=4:mes str(int(rey.cnt))

loop

stop

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

これはhsp 3.31で動きます。他のverでの動作は保証できません

ソースの補足として、
命令として「fft」「ifft」「fft初期設定」の３つを登録した。それぞれ離散フーリエ変換、逆フーリエ、バタフライ演算のインデックス作成、となっている。
「fft」「ifft」ともに、とる引数は実数配列（double）,虚数配列（double）で、変換結果がそれに返される。
配列の要素数は、2^beki個であることが条件、それ以上でも以下でも、１でもずれたらうまく計算されない。

実行すると、
実行画面の左には、A*B＝CのCの配列が表示されている
実行画面の真ん中には、A*B＝CのAの配列が表示されている
実行画面の右には、A*B＝CのBの配列が表示されている
bekiの値を増やすほど、多くの桁数の乗算が行える
仮にbekiを23にすれば、約400万配列×400万配列が800万配列に格納される
bekiが23でkisuuが10000の場合、約1600万桁×1600万桁=3200万桁の乗算となる

fftモジュールはもちろん乗算以外の用途にも使える。

2009年07月08日00:00

カテゴリ: １億桁×１億桁

１億桁×１億桁の開発状況報告③

１億桁×１億桁の開発は・・・ごめんなさい。なんかもう行き詰ってます（泣

だいたい１億桁×１億桁計算する前に１億桁を入力するだけでもすごい時間かかってしまうじゃないですか！！

だれもやりませんよそんなの・・・てなんで自分で自分を否定するか・・・

もう決めた。

これはあきらめよう！

うん。決断は早いほどいい。

せっかくだから作りかけの１００万桁×１００万桁のプログラムを応募しちゃえ！

いやーこれでも画期的ですよ！

本来なら２分３分かかってしまう計算がたった３秒足らずで！

と、いうわけなので、ＨＳＰコンテスト２００９応募作品１

その名も「１００万桁×１００万桁を計算するソフト」

うむ、我ながら良いネーミングセンスだ

でもこれじゃ誰もぎゃふんと言わせられないので、なにか別のプログラムを作らないと。

そうだな、特にショートプログラムで作りたいものだ。

そして去年のレイトレみたいに技術的にすごいのがいいなぁ

となるとやはり視覚的にうったえるもので、なにか派手なものは・・・

やっぱアルゴリズムにかかってくるに違いない。

例えば、ちっこい文字がたくさん飛んでて、どっかに吸い寄せられて形を作るものだったり、

幻想的な模様を再現したり、それこそレイトレみたいに光源の計算だったり、

全部アルゴリズムしだいでインパクトさが決まるような

よし、なにか次回までに考えてくるとするか。

今日の名言

「手で割ればすむものをわざわざ機械を使うなんてねぇ。ああいうものを買う人の気持ちが分かりませんよ。どうせ買うのは卵なんか割ったことの無い関白亭主ですよぉ」　ｂｙ　ノリスケ

2009年06月01日00:00

カテゴリ: １億桁×１億桁

１億桁×１億桁の開発状況報告②

１億桁×１億桁の開発状況報告２

前はkaratsuba法の名前が出てくるところまで話が進んだんですよね？

前も言いましたがlognint.dllとは、ものすごいながい桁数の数値も計算できるようにするためのプラグインです。

これを使います。

ＨＳＰユーザーでこれを使ってる人はなかなか見ないもんで、その技術解説をしているホームページなども見ないので、このブログがＨＳＰユーザーとして一番最初に詳しくlognintについて技術解説することになりそうです！！

なんかうれしい！！

もともと私は四則演算とか（自称）とても得意なんで、一般の人には負けないつもりです！

さてkaratsuba法のプログラムについてですが
ちょっとプログラム的な話ではなく、ほんと数学的な話になりそうですが、まぁ仕方ないですが、許してください

karatsuba法は一言で言うと、とても長い桁数を演算するときに使え、計算手順を効率よくすることで計算時間を短くするという方法の一つです。

実はもっと効率のいい高速フーリエ変換という方法もあるのですが、自分の技術不足でkaratsubaが限界です・・・・

さて実際の手順ですが、

まず四則演算の計算で一番時間がかかるのが掛け算、割り算。

桁数に比例して時間が伸びます。

足し算引き算は、掛け算と比べほとんどゴミのようなものです。

なので、効率化するべきは掛け算です。（割り算は置いておく）

ここでちょっと分かりやすく、１０桁×１０桁の掛け算を考えます。

本当は千桁くらいから使わないと意味がないのですが

例えば１３８４５０１２３９×５６７８２２１０９８の計算

このとき、普通

a=１３８４５０１２３９

b=５６７８２２１０９８

c=a*b

とこうやりますね？

ここで、桁数を最大５桁に制限し変数を分割することを考えます。

a1=１３８４５

a2=１２３９

b1=５６７８２

b2=２１０９８

とし

c=a1*b1*10^10+(a1*b2+a2*b1)*10^5+a2*b2

とこうするとします。

ここで*10^10と*10^5は、変数上でも後ろに0がつくだけなので計算時間がまったくかからないものとします

一つ一つの掛け算が５桁×５桁になったので、１０桁にときと比べると４分の１になったにですが、これでも、掛け算が合計４つになり、結局計算時間的にはなにも変わりません。

ここでやっとkaratsuba法の登場

括弧でくくりまくって

c=(a1*b1)*10^10+(a2*b2)-((a1-a2)*(b1-b2)+(a1*b1)+(a2*b2))*10^5

とこうします。

一見掛け算が増えたように見えますが、a1*b1とa2*b2のところが２箇所ずつあります

だからそこの計算はしなくていいのです。

簡潔に書くと

c1=a1*b1

c2=a2*b2

c=c1*10^10-((a1-a2)*(b1-b2)+c1+c2)*10^5+c2

とこうなります！

よく見てみると*10^10とか省くと、掛け算が３回に減っています！

これが計算の効率化です。

さっきまでずっと、*10^10をの計算時間を省くのを気にしている人がいるかもしれません

が、実際の計算で、計算結果を一つの変数にまとめなければいけないというルールはありません。

１億桁の計算結果を出力するプログラムでは、１万桁ずつ１万個の変数に区切って出力する方法をとっても

なんら問題はおきないのです（当然繰り上がり処理は別途必要になります）

だから、「この変数は10^10が本来かかっている数」として自分の中で考えておけば、わざわざコンピューターに

*10^10の計算をやらせる必要はなくなるのです。

このようにして掛け算が効率化されました。

次のハードルは、karatsuba法の帰納的活用法です。

４回から３回に減った掛け算の中でさらにまたkaratsuba法を適応して計算時間を４分の３に減らします

2009年05月14日00:00

カテゴリ: １億桁×１億桁

１億桁×１億桁の開発状況報告①

１億桁×１億桁の開発状況報告

いまさらですが、１億桁×１億桁の計算をしなければならない人ってこの世にいるのでしょうか・・・

あまり気にしすぎると体に悪いので、早速解説していきたいと思います

まず１億桁×１億桁の計算に絶対必要なのは多倍長の変数

普通の変数は９桁あたりで限界振り切ってしまいますが多倍長の変数

は何桁までいっても限界を振り切りません

そのためにlongint をプラグインとして使います

まず試しに１００万桁×１００万桁の掛け算を行います。

#include "longint.hsp"
a=LongInt(1)
repeat 100000;百万桁の数を生成
a*=LongInt("10000000000")
loop
b=LongInt(a);百万桁の数を複製
c=longint(0)
wait 2
mes "掛け算開始"
wait 1
c=a*b
mes "掛け算終了"

・・・・・・・・・・・・コピペして実行してみれば分かりますが

まず百万桁の変数を作るのに１分

百万桁×百万桁に１分１０秒

長い・・・・長すぎる・・・

こんなんで１億桁×１億桁の計算をやろうとすると、計算時間が尋常じゃないほど長くなります。

一応どのくらいかかるのかと言うと、桁数１００倍より１００^２=１００００だから

計算時間が１万倍になるので、掛け算部分だけでも７０００００秒！

日にちにして８．１日かかる！！

これに１億桁の数を生成＆複製の計算時間を含めるとさらに長く・・・

これではデバッグだけで１年が過ぎてしまう・・・

結論はつまり、「普通の掛け算」ではダメだということです！！！

何か、計算効率のいい掛け算方法（アルゴリズム）を考えなければいけません！

ちなみにわれわれがよく手動で計算する掛け算の手順

一桁一桁ずつかけて足して・・・という方法はじつはとてつもなく非効率な掛け算方法なのです！

ただ人間にとっては手順が単純明快なので、一番楽な方法として認識されてますが

コンピューターにやらせるのには向いてません

ここで、karatsuba法というけ掛け算アルゴリズムを採用します。

これは多倍長の桁を２分割して掛け算４回分の仕事を掛け算３回分と足し算４回に式変形する方法です。

足し算は掛け算と比べて計算時間が非常に短いのでないものと考えると、１回の分割で計算時間が３／４に減る

ことになります。

すると８回分割で計算時間がなんと１０分の１に半減します！

それでもまだ１つの掛け算に４０万桁×４０万桁の計算が必要なのでさらに５回分割して

１万２千桁×１万２千桁の掛け算になるまで分割すれば、計算時間が（３／４）^13≒0.024

だいたい４０分の１くらいに時間が抑えられることになります！

するとなんとたった５時間で、１億桁×１億桁を求められる計算になります！

おおー早い！！

これで、だいたいメドがついてきましたね！

あとはプログラムを打ち込むだけです。

karatsuba法をどうプログラムで再現するか、ここが一番力量が試されるところですが、

何せ計算手順が複雑なのでまだベースプログラムすら出来上がってません・・・

（一応自分が過去にＨＳＰコンテストに出した中で一つkaratsuba法を使った計算プログラムがあるのですが、また新しくプログラムを打ち直すつもりです）

なので詳しいkaratsuba法の計算についてはまたの機会にしたいと思います

次回：技術解説「ババドン」その３　吸引アルゴリズム