仮眠プログラマーのつぶやき

自分がプログラムやっていて、思いついたことをつぶやいていきます。

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

2012年12月

流体力学のプログラムを作りたい!その3

流体力学に関する記事はかなり久しぶりである。

約束通りCIP法の解説をしたいと思ったが、解説に使う労力があまりに多いため、また最近多忙なため、CIP法のソースの公開と、参考になったリンクを貼るのみとさせて頂く。

今までネット上に落ちているソースを元に独学で流体力学を学んできたが、ついにCIP法の難しさにより本を購入。

CIP法とJavaによるCGシミュレーション 単行本(ソフトカバー) – 2007/2/17
矢部 孝  (著), 尾形 陽一  (著), 滝沢 研二  (著)

といっても私でも理解できたのでかなりやさしいのでは。CIP法のソースが最初から最後まで乗っている解説本。
今回紹介する公開サンプルはそれのHSP移植かつ模倣にとどまる。

以下のスクリーンショットがHSP上で実行したものだ。
cip法流体.rar

cip有限音速
これは、音速が有限で、混合密度性流体である。
「cipcup格子法低音速、混合密度.hsp」の実行画面である。

画面の意味するものは、明るい部分が密度の高いもの、暗い部分が密度の低いもので、重力なし、境界は反対側につながっていて、初期条件は左半分が高密度で真下に向かう速度を持ち、右半分が低密度で真上に向かう速度を持つ。
流体を妨げる障害物はない。


そして以下のものは、本のソースを参考に自分で改変し、音速を無限、密度一定で非圧縮性流体にしたものだ。
当然CIP法なので、数値拡散がとても少なく、サンプル内のレンダリング動画を見てもらえれば分かる通り、160×120格子でも細かい渦が消えずに残っているのが分かる。
(もしCIP法でないと細かい渦が拡散してしまう)
これは障害物(壁)を3枚、左側に配置してある。
cip無限音速

壁の境界条件をどの段階で設定すればいいのかわからず、移流項計算、圧力項計算の前後に無駄に多くの壁境界条件設定をしているが、こんなにする意味は無い・・お陰で汚いソースである・・・・。



最後にお世話になりまくったサイト or PDFのご紹介。
特に参考になったリストは 
以下、その他参考になったサイト or PDF
http://www.4gamer.net/games/000/G000000/20070926041/
http://cfdsim.web.fc2.com/download_source.html
http://www.civil.hokudai.ac.jp/yasu/home.htm
http://www.exocortex.com/plugins/slipstreamvx
http://grantkot.com/Multiphase/Liquid.html
http://grantkot.com/MPM/Liquid.html
http://ja.wikipedia.org/wiki/CIP%E6%B3%95
http://news.mynavi.jp/articles/2007/10/09/cedec01/001.html
http://art-science.org/journal/v3n2/v3n2pp159/artsci-v3n2pp159.pdf
http://www.csp-consortium.org/action/csp_200910/nakane_200910.pdf

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=c」の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
w-
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初期設定」の3つを登録した。それぞれ離散フーリエ変換、逆フーリエ、バタフライ演算のインデックス作成、となっている。
「fft」 「ifft」ともに、とる引数は実数配列(double),虚数配列(double)で、変換結果がそれに返される。
配列の要素数は、2^beki個であることが条件、それ以上でも以下でも、1でもずれたらうまく計算されない。

実行すると、
実行画面の左には、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モジュールはもちろん乗算以外の用途にも使える。 

パズドラルート解析アプリの問題点と今後の抱負

PC用解析君はコチラ

新しいビットマップ イメージ




問題のアプリです
https://itunes.apple.com/jp/app/andrpzsolve/id580792881?l=ja&ls=1&mt=8
https://play.google.com/store/apps/details?id=net.onionsoft.hsptv.pdroute
アプリの名前柄、流行りに便乗して多くのユーザーにダウンロードしてもらったという形になったわけですが、
まぁ予想通りではあったが、計算速度があまりに遅すぎるというご指摘を多く頂きました。
その割にはちゃんと評価してくださっているユーザーもいて、親切だなとも思ったというのが感想です。

では早速、このアプリの問題点を全て挙げてみましょう
1、計算速度が遅すぎる
2、スクリーンショット取り込みができない
3、iphone、ipadだとメニュー画面に戻れない 
4、androidでは計算中に強制終了画面がでる 



1、計算速度が遅い
一番の原因は「全経路検索」をしているためです。
数学に「場合の数」という概念があります。
数学的に考えるなら、ドロップをあるルートで動かしたとき、一番コンボ数が多くなる経路を求めるには、すべての場合のルートを確認する必要があります。

ここでルートの最大距離を1マスとしたとき、全ての経路の場合の数は98通りとなります。
とりうる始点は30通り(横6マス×縦5マス)、そこからドロップ移動可能方向の上下左右4をかけ、30×4で120通りとなりますが、欄外にはみ出る経路の場合の数は考慮してはなりません。
例えば始点位置が一番左上の場合、そこから左、上へはもう移動できません。

同じように、ルートの最大距離が2マスとすれば、1マス移動の98通りのルートからさらに1マス上下左右へ移動可能なので4を掛けますが、先ほどと真逆の方向へは移動できません(戻る方向に行けばルート距離が1減るため)。さらに欄外へ行く経路を除外すれば236通りとなります。
以降ルートの場合の数は、 586通り、1452通り、3574通り、8764通り、21412通り、52220通り、127376通り、311224通り、760980通り・・と指数関数的に増えていくわけです。

ルート距離が12マスの場合、全ての経路のとりうる場合の数は1860172通りあります。
ハード的にはだいたいここらまでが限界です。
なぜならここらでスマホの使用メモリが200Mbを超えてしまうからです。
ルート1通りあたりに使うバイト数は、95byte。
内訳は、30個のドロップ(6種類)配列の格納で30byte。
同色ドロップが横2つ続いているときのみ0を表すいわゆる横差分値を格納する変数で30byte、同じく縦でも同じのを作り30byte。
これに加えルートの移動ログを格納する変数で4byte(1マス移動につき2bit)
さらに、ルートの先端部分つまり現在位置を格納する変数で1byte

スマホのメモリ容量は機種によりけりだが、スワップが発生すると劇的に遅くなるのでこれ以上のメモリ使用は難と考えられます。
また計算時間も非常識なものとなります。
1860172通りのルート全てで、コンボが発生するか、発生する場合は何色のドロップが何個何連鎖コンボが発生するかを計算する必要があるので、2.4GhzのノートパソコンのCPUでも3分近くかかります。
スマホのCPUは約1Ghz前後であるしキャッシュもあるのか無いのかよくわからないので、最大12マスの設定で全経路検索をすれば10分近くかかることとなります。

まぁただでさえ非常識な時間なのでこれ以上は不可能でしょう。
全経路検索にしなければ計算時間は劇的に早くなるわけですが、それではコンボの精度が悪化してしまうという欠点があります。
それでも計算速度が早くなる利点は非常に大きいですが、それは以下に示す問題点により、あえて精度をとったという経緯がありました。


2、スクリーンショット取り込みができない
これは、ただただ私の技術不足のため、という一言につきます。

アプリを開発するにあたり、取る方法は2通りあります。
ひとつはHSPdishを使い開発すること、もうひとつはxcodeでobject c、eclipseでjavaを使い開発すること。
後者ではスクリーンショット読み込み機能が実装できます。
私にはHSPで作るしか能力がなく、残念ながらスクリーンショット読み込み機能が実装できませんでした。

また私の機種にはそもそもスクリーンショットを取る機能がなく、そのような機種でもルート解析ができるようにしたいという思いがあったのも事実ですが、やはりスクショ機能はあったに越したことはありません

つまり、ドロップ配列を手動で入力するしかなく、ここでユーザーに多くの労力を強いることが容易に想像できました。
時間と労力をかけて入力しても、計算が一瞬で終わってしまいその割にはたいしたコンボ精度が得られなければ、意味がありません。

時間と労力をかける価値がある状況というのは通常の状況ではあまりありえなく、詰みそうな場合やこの1ターンで戦局が大きく変わる場合など、ここ一番という場面に限られます。
そういうときは、一回だけすごい時間をかけても最大コンボルートを導き出せればいいのですから、計算時間より精度を優先したシステムの方を選択し、結局このような非常識な計算時間になりました。


3、iphone、ipadだとメニュー画面に戻れない
※最新情報で、「ホームボタン二度押し、下の段に出てくるアプリアイコンを消してあげると再度アプリ起動した時に解析できる」ようになっているそうです。

これは評価コメントを見てはじめて気がついたのですが、iphone、ipadの場合、ホームボタンを押したときの動作として、アプリが終了するのではなく、アプリが待機状態になるという性質を理解していないために生まれてしまった不具合でした。申し訳ありません


4、androidでは計算中に強制終了画面がでることがある
これは私が、計算プログラムにwaitを入れるのを忘れたため、タスクがビジー状態になっているためです。
なんという初歩的なミス・・
対策方法は、エラーメッセージが出たら「待機」ボタンを押して待つ、です。いつか結果画面が表示されるでしょう

毎年のHSPプログラムコンテストの作品でもそうですが、私はwaitを 入れるのを忘れる癖があるようです
いい加減直したほうがいいですね
申し訳ありません



なにはともあれ、リリースして3週間くらい経とうとしていますがとても多くのダウンロードと評価をありがとうございます。そして沢山の不具合に関しては、申し訳ございません。

HSPdish製のアプリということで、アプリランキングでそこそこいい所までいけたというのは、HSPの評価を上げるまぁまぁいい話題になるのかなと思う反面、計算が遅かったりと評判が悪いのを見ると逆にHSPの評価を下げてしまう結果になるのではと少し恐れてもいます・・・

まぁアプリの名前は少し誤りがありますからね
もっと誠実につけるとしたら「移動距離が12マス以内で、コンボ数がそれ以上多くなるようなルートがないことを証明するアプリ」とでもなるのでしょうか。そのルートには移動距離制限がつくため、人間がルートを考えた場合13マス以上でもありうるので、人間が考えたほうがコンボ数が多くなる可能性もあります。

だから役に立つ場面なんてのはほとんどないかもしれません。



そういえばHSPdishが出た当初どこかで、その便利さを称えるコメントかなんかで、「これのお陰でしょうもないアプリがたくさん世に放たれるのか」といったのを見た記憶があるのですが、自分でまさにその通りにしてしまいました。



そして今後の抱負ですが・・
そろそろツールじゃなくてゲームを作りたいな・・

もちろんプラグイン開発も同時進行で
むしろプラグインを積極的に使ったゲームとか
プロフィール

toropippi

記事検索
アクセスカウンター

    QRコード
    QRコード
    • ライブドアブログ