Skip to content

logicia32/z80-parallel-pi

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

2 Commits
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

z80-parallel-pi — FPU の無い 8 ビット Z80 で円周率を、並列に求めてみる

1970 年代の 8 ビット CPU Z80 に、浮動小数点を一切使わず(FPU が無いので) 固定小数点だけで円周率 π を計算させ、さらに Z80 を 1〜4 台つないで 「分担すると速くなるのか/ならないのか」を実測するための実験ノートです。 大学の数値計算・並列計算の授業で、手を動かして確かめる教材になればと思って まとめました。再現手順はこの README に全部書いてあります。

題材は並列計算の教科書でいちばん最初に出てくる定番です。

$$\int_{0}^{1}\frac{4}{1+x^{2}},dx=\pi$$

右辺がきっかり π になるので、左辺を数値積分すれば円周率が出ます。これを 3 つの数値解法と 2 つの並列方式で実際に走らせ、グラフにします。


何が分かる教材か(先に結論)

観点 見えること
台形公式 区間を分けるだけなので素直に並列化でき、台数で速くなる
シンプソン公式 同じ分割数でも精度が桁違いに高い(誤差 ~h⁴ 対 ~h²)
ニュートン法 反復が逐次の鎖なので並列にしても速くならない(解析的に自明=グラフの線は実測でなく解析値)
①分散メモリ 送って持ち寄る型。速くなるが 2倍そこそこで頭打ち、しかも走るたびに振れる(rank0 集約が直列+ホスト依存)
②共有メモリ 取り合う型。4 台で 3.73 倍、しかも決定論的に毎回同じ。ロックを忘れると答えは 0 ではなく過小集計に静かに狂う(P=4 で π≈0.94)
共通 同じカーネル pi_kernel.h(単一ソース)。固定小数点の整数和は分割の仕方・方式・台数に依らずビット単位で同じ答え π=3.14148

「並列にしても速くならない計算がある。それは ~ だから」を、ニュートン法を 題材に実機相当のクロック数で示すのが、この教材のいちばんの狙いです。


ディレクトリ構成

z80-parallel-pi/
├── pi_kernel.h               ★ 単一の台形カーネル(下記3つが #include)
├── methods/                  単一ノードの 3 解法(まずここ)
│   ├── trapezoid.c             台形(pi_kernel.h を使う)
│   ├── simpson.c               シンプソン公式(同じく Q14・偶数N)
│   ├── newton.c                ニュートン法 = ヘロンの平方根(逐次の鎖)
│   ├── console.c / console.h   ucsim_z80 への文字出力(I/O ポート 0x80)
│   └── Makefile
├── parallel/
│   ├── distributed-memory/   ① メッセージパッシング型
│   │   ├── node_pi.c           SPMD ノード本体(pi_kernel.h を使う)
│   │   ├── console_mm.c        メモリマップ simif バックエンド
│   │   ├── router.py           ホスト側ルータ(= ネットワークそのもの)
│   │   └── Makefile
│   └── shared-memory/        ② 共有メモリ型(自作マルチコア)
│       ├── z80.h               1 ヘッダ Z80 コア(第三者・zlib、NOTICE 参照)
│       ├── z80impl.c           z80.h を実体化する1行(原作・MIT)
│       ├── machine.h mc.c cpu.c バス調停器 + ピン配線 + IM2 割り込み
│       ├── node_pi2.c          Z80 側ワークロード(pi_kernel.h・ロック有/無)
│       ├── t63.c               N コアを回すホスト harness(pi 集約)
│       ├── crt0_irq.s isrdemo.c 割り込みデモの Z80 側(IM2・手書きISR)
│       ├── irq_host.c          割り込みデモのホスト("NIC"役・検証)
│       └── Makefile
├── graphs/                   速度向上・精度のグラフ生成(matplotlib)
├── writeups/                 解説(読み物・01〜04)
├── LICENSE                   本リポジトリのオリジナル部分(MIT)
└── NOTICE                    第三者コードの帰属(z80.h のみ・zlib)

pi_kernel.hただ 1 ファイルgrep -rn pi_kernel.h で、単一 ノード版・①・② が同じカーネルを #include しているだけ=「同じ 計算、束ね方だけ違う」を自分の目で確認できます。


必要なもの

  • SDCC(Z80 向け C コンパイラ)と ucsim_z80(同梱のシミュレータ) 多くのディストリで SDCC パッケージに ucsim_z80 が同梱されます。
  • makepython3、グラフを描くなら matplotlib

確認:

sdcc --version          # z80 ターゲットが含まれること
ucsim_z80 -t Z80 -h     # 起動すること
python3 -c "import matplotlib"   # グラフを描く場合のみ

処理系を PATH 以外に置いている場合は、各 makeSDCC=/path/to/sdcc UCSIM=/path/to/ucsim_z80 を渡せます。

本リポジトリの絶対クロック数SDCC 4.5.24 + ucsim_z80 0.9.9 での実測値です。別の版でもビルド・実行できますが、クロックは数 % ずれます。π=3.14148・「①は走るたびに振れる/②は決定論的」・ 並列で頭打ちする/しない、という"形と結論"は版に依りません—— 手元で確かめてほしいのは、個々の数字ではなくそこです。


1. まず単一ノードで 3 解法(methods/)

cd methods
make run                       # 3 つを順に実行

期待される出力(抜粋):

trapezoid: pi~=3.1416 (Q14, N=1024)
simpson:   pi~=3.1416 (Q14, N=1024)
newton: sqrt(2) by x_{n+1}=(x_n + S/x_n)/2  (Q14)
  iter 0: x=1.00000
  iter 1: x=1.50000
  iter 2: x=1.41662
  iter 3: x=1.41418
  ...

分割数 N を変えて精度の出方を見られます(台形とシンプソンの差が出ます):

make N=4   run-trapezoid run-simpson     # 台形=3.1311 / シンプソン=3.1415
make N=64  run-trapezoid run-simpson     # 台形=3.1414 / シンプソン=3.1415

シンプソンは N=4 の時点で台形の N=256 相当の精度に達します。あとは Q14(= 1/16384 ≒ 6×10⁻⁵)の量子化の床に当たって、それ以上は詰まりません。 「もっと細かく刻む」より「いい公式を選ぶ」ほうが効く、という体験です。

ニュートン法は x_{n+1} を出すのに x_n が要ります。前の値が無いと次が 計算できない=逐次の鎖。これが後で「並列にしても速くならない」話に効きます。


2-①. 分散メモリ:送って持ち寄る(parallel/distributed-memory)

各 Z80 は独立した ucsim_z80 プロセス。ノード間に共有メモリは無く、やり取りは 全部「フレーム」。router.py がその配線(= ネットワーク)そのものです。

cd parallel/distributed-memory
make run                       # P=1,2,3,4 を走らせ表を出す(数回どうぞ)

実測例(クリティカルパス = 一番遅いノードのクロック数基準)。P=1 は 決定論的に 1324206。P≥2 は rank0 の集約待ちがホストのスケジューリング 次第で大きく振れます。下はある 1 回の観測で、speedup 列はその行の 1324206 ÷ crit-path(このサンプル自身の値)。カッコの帯は別の試行も 含めて観測した散らばりの目安で、厳密な上下限ではなく走るたび・環境ごと に変わります(グラフ speedup.png は帯の中央値を線で描きます):

 P | per-rank Z80 clks            | crit-path |   pi~= | speedup (=1324206/crit) | 観測帯の目安
 1 | 1324206                      |   1324206 | 3.14148|   1.00x                 |   —
 2 | 752287 608568                |    752287 | 3.14148|   1.76x                 |   1.6〜1.8
 3 | 539105 413728 423173         |    539105 | 3.14148|   2.46x                 |   1.7〜2.5
 4 | 681382 321060 319924 312529  |    681382 | 3.14148|   1.94x                 |   1.8〜2.9

この 1 回では P=4 が P=3 より遅い(2.46x → 1.94x と逆転している)こと に注目してください。台数を増やしたのに遅くなる——これこそ「①は走るたび に振れて読めない」の生の姿で、別の回では P=4 が 2.9x まで伸びることも あります。非 rank0 の計算クロックは決定論的(毎回同じ)。動くのは集約 担当 rank0 だけで、台数を増やしても縮まない(アムダール)うえに走る たびに振れる——「速いが読めない」のが①の正体です。②(次節)の決定論 との対比が芯。

2-②. 共有メモリ:取り合う(parallel/shared-memory)

検証済みの 1 ヘッダ Z80 コアを N 個インスタンス化し、バス調停器・ハードウェア セマフォ・バリア・IM2 割り込みを自前で足した「1 つの SMP マシン」。Z80 側の ワークロードは第1回と同じカーネル pi_kernel.h(束ね方が違うだけの 別 binary)。集約は共有 RAM 窓 0xC000 に置いた 32bit 値の RMWです。

cd parallel/shared-memory
make run    # ロック版(P=1..4一致) / 無ロック版(過小集計) / 割り込みデモ

実測例(自作シミュレータは決定論的=毎回この値):

=== (2) shared-RAM reduction, locked (HW semaphore) : P=1..4 ===
  P=1  shared acc=3294113 raw=51470 pi~=3.14148  MATCH  (crit 1205949)
  P=4  shared acc=3294113 raw=51470 pi~=3.14148  MATCH  (crit  323096, 3.73x)
       lock spins: P=2:9  P=3:23  P=4:46
=== (2) NO lock : the classic lost update (measured, not 0) ===
  P=2 acc=1349508 pi~=1.28699   P=3 acc=1106018 pi~=1.05475
  P=4 acc= 988371 pi~=0.94257   (正解は 3.14148)  ... MISMATCH
=== interrupt-driven receive (IM2, hand-written ISR) ===
  bytes by interrupt=4  interrupts taken=4  checksum 0xAA  -> PASS

ロック有りは 4 台で約 3.73 倍、しかも何回走らせても 1 clk も動きません (①と対照的)。ロックを外すと共有 RAM の read-modify-write が競合して 答えは 0 ではなく過小集計に静かに狂う(P=4 で π≈0.94。決定論的に毎回 同じ値)。「ロストアップデートは 0 になる」ではなく「気づきにくい中途半端 な値になる」——授業で習うやつを実数で。割り込みデモは NIC が /INT を 上げ、IM2+手書き ISR で受ける本物で、main は NIC を一切ポーリングせず、 受理割り込み回数とチェックサムを突き合わせて PASS/FAIL を出します。


3. グラフ(graphs/)

cd graphs
make                           # speedup.{png,svg} / accuracy.{png,svg}
# 別の python を使う場合: make PYTHON=/path/to/python3
  • speedup.png … ②が線形に近く決定論的(実測・毎回同じ)、①は代表値 +観測した散らばり帯(host 依存で走るたび変わる。帯は目安で厳密境界 ではない)、ニュートンは 1.00 の水平線(解析値:鎖は分割不能なので 測るまでもない。実測の②との対比用)。
  • accuracy.png … 解法ごとの誤差。シンプソンが早々に Q14 の床に達する。

graphs/plot.py の数値は上記 harness の実測(SDCC 4.5.24 + ucsim_z80 0.9.9)。②は決定論的なので固定値。①は host 依存で振れるため代表値+ 観測した散らばり帯(試行・環境依存で、厳密な上下限ではない)。 ニュートンの線だけは解析値(並列実装は無い)と docstring に明記。 SDCC の版が変わるとクロックは数 % 動くが「形」は不変。


なぜニュートン法は並列にしても速くならないのか

台形・シンプソンは「区間(ペア)ごとに独立」なので、範囲をコアに配るだけで 分割できます。足し算は結合的なので、どう分けてもビット単位で同じ答え。

ニュートン法は違います。

$$x_{n+1}=\frac{1}{2}\left(x_{n}+\frac{S}{x_{n}}\right)$$

x_{n+1} を計算するには x_n が要り、x_{n+2} には x_{n+1} が要る。 反復は 1 本の鎖(逐次の依存関係)で、途中の輪を別のコアに渡せません。 1 回ぶんの評価(掛け算・割り算 1 個ずつ)を分割しても割に合わない。 アムダールの法則でいう「直列部分の割合 ≈ 1」で、4 台ぶら下げても実際に 働けるのは 1 台・残りはただ待つだけ=速度向上は 1 倍に張り付きます (speedup.png の灰色の線。これは分割不能ゆえ自明なので解析値として 描いてあり、台形/シンプソンの実測曲線との対比用です)。

「並列化が効くかどうかは、計算の中身ではなく依存関係の形で決まる」 —— この教材でいちばん持ち帰ってほしいのはここです。


ライセンスと謝辞

  • 本リポジトリのオリジナル部分(pi_kernel.h、バス調停器・マルチコア harness、IM2 割り込み層・割り込みデモ、ルータ、ビルドスクリプト、 グラフ、z80impl.c の 1 行ラッパー)は MITLICENSE)。
  • parallel/shared-memory/z80.h のみが第三者の 1 ヘッダ Z80 エミュレータで、zlib/libpng ライセンス(Andre Weissflog)。帰属は NOTICE を参照してください。素晴らしいコアを公開してくださっている ことに感謝します。

詳しい背景と読み物は writeups/ にあります(第1〜4回)。

About

FPU の無い 8 ビット Z80 で固定小数点の円周率を、1〜4 台で並列に求めて実測する教材。台形/シンプソン/ニュートン+分散メモリ・共有メモリの2方式。

Resources

License

Stars

Watchers

Forks

Releases

No releases published

Packages

 
 
 

Contributors