リポジトリ / emssr / independent_common_compiler

# -*- coding: utf-8 -*-

"""lane64 の GPU 実行台（limb 版）。thread 一本が lane 一本になる。

仕様は docs/common-compiler/arch/20。式 (2.1) の解釈器と同じ意味を、CUDA の kernel の中で

lane ごとに走らせる。cuvcpu.py は語を bit 平面に持ち、命令一つが平面の演算の列になるので、

語が広いと launch の数で足が止まる。ここでは語を 64 bit の limb の並びで持ち、lane ごとの

thread が block を跨いで命令を辿る。分岐は thread ごとに好きに飛ぶ（mask は要らない）。

  並び   語 主。語 w の limb i、lane l は M[(w·NLIMB + i)·lanes + l] に置く。

         同じ語を読む lane の番地が連続するので、warp の memory access が一度で済む。

         lane 主に置くと lane ごとに words·NLIMB 語ずつ離れ、一命令が 32 回の転送になる

  語     NLIMB 本の limb。二の補数、little endian。既定 16 本 = 1024 bit

  landmark  R[0..7]。commonc の lane64 版は R1 R2 R3 だけを使う（back_lane64.cam の覆い）

  作業領域  lane ごとに words 語。読む語（watch）だけを host へ戻す

語の幅を越えた乗算は報告 43 で止まる。lane64 の語は ℤ なので（20 §0.3）、この道は

NLIMB·64 bit を越えた時点で core 意味論と離れる。

  run(prog, seeds, words, watch, need=None) → lane ごとの (報告, {語: 値})

"""

import numpy as np

INSN_FIELDS = 7

NREG_BYTES = 8 * 8          # shared memory へ置く register。NREG 本 × 8 byte

R_OK, R_KO, R_WIDE, R_BUDGET = 0, 42, 43, 44

SRC = r"""

#define NLIMB %(nlimb)d

#define NREG 8

typedef unsigned long long u64;

typedef long long i64;

struct Big { u64 v[NLIMB]; };

__device__ __forceinline__ void big_zero(Big &r) { for (int i = 0; i < NLIMB; i++) r.v[i] = 0; }

__device__ __forceinline__ void big_imm(Big &r, i64 x) {

    r.v[0] = (u64)x; u64 f = (x < 0) ? ~0ULL : 0ULL;

    for (int i = 1; i < NLIMB; i++) r.v[i] = f;

}

__device__ __forceinline__ bool big_neg_p(const Big &a) { return (a.v[NLIMB - 1] >> 63) != 0; }

__device__ __forceinline__ void big_add(Big &r, const Big &a, const Big &b) {

    u64 c = 0;

    for (int i = 0; i < NLIMB; i++) {

        u64 s = a.v[i] + b.v[i]; u64 c1 = (s < a.v[i]) ? 1 : 0;

        u64 t = s + c; u64 c2 = (t < s) ? 1 : 0;

        r.v[i] = t; c = c1 | c2;

    }

}

__device__ __forceinline__ void big_not(Big &r, const Big &a) { for (int i = 0; i < NLIMB; i++) r.v[i] = ~a.v[i]; }

__device__ __forceinline__ void big_negate(Big &r, const Big &a) {

    Big n; big_not(n, a); Big one; big_imm(one, 1); big_add(r, n, one);

}

__device__ __forceinline__ void big_sub(Big &r, const Big &a, const Big &b) { Big nb; big_negate(nb, b); big_add(r, a, nb); }

__device__ __forceinline__ bool big_eq(const Big &a, const Big &b) {

    for (int i = 0; i < NLIMB; i++) if (a.v[i] != b.v[i]) return false; return true;

}

__device__ __forceinline__ bool big_lt(const Big &a, const Big &b) {

    bool na = big_neg_p(a), nb = big_neg_p(b);

    if (na != nb) return na;

    for (int i = NLIMB - 1; i >= 0; i--) { if (a.v[i] != b.v[i]) return a.v[i] < b.v[i]; }

    return false;

}

__device__ __forceinline__ void big_shl(Big &r, const Big &a, int s) {

    if (s <= 0) { r = a; return; }

    int w = s / 64, bts = s %% 64;

    Big t; big_zero(t);

    for (int i = NLIMB - 1; i >= w; i--) {

        u64 hi = a.v[i - w] << bts;

        u64 lo = (bts && i - w - 1 >= 0) ? (a.v[i - w - 1] >> (64 - bts)) : 0;

        t.v[i] = hi | lo;

    }

    r = t;

}

// |a|·|b| を 2·NLIMB limb で作り、上半分が零でなければ幅を越えたとみなす

__device__ bool big_mul(Big &r, const Big &a0, const Big &b0) {

    Big a, b; bool na = big_neg_p(a0), nb = big_neg_p(b0);

    if (na) big_negate(a, a0); else a = a0;

    if (nb) big_negate(b, b0); else b = b0;

    // 上位の零の limb は掛けない。小さい数どうしなら 1 回の乗算で済む

    int la = 0, lb = 0;

    for (int i = 0; i < NLIMB; i++) { if (a.v[i]) la = i + 1; if (b.v[i]) lb = i + 1; }

    u64 acc[2 * NLIMB];

    for (int i = 0; i < 2 * NLIMB; i++) acc[i] = 0;

    for (int i = 0; i < la; i++) {

        u64 carry = 0;

        for (int j = 0; j < lb; j++) {

            u64 lo = a.v[i] * b.v[j];

            u64 hi = __umul64hi(a.v[i], b.v[j]);

            u64 s = acc[i + j] + lo; u64 c1 = (s < lo) ? 1 : 0;

            u64 t = s + carry; u64 c2 = (t < s) ? 1 : 0;

            acc[i + j] = t;

            carry = hi + c1 + c2;

        }

        int k = i + lb;

        while (carry && k < 2 * NLIMB) { u64 s = acc[k] + carry; carry = (s < carry) ? 1 : 0; acc[k] = s; k++; }

    }

    bool wide = false;

    for (int i = NLIMB; i < 2 * NLIMB; i++) if (acc[i]) wide = true;

    if (acc[NLIMB - 1] >> 63) wide = true;

    for (int i = 0; i < NLIMB; i++) r.v[i] = acc[i];

    if (na != nb) { Big t = r; big_negate(r, t); }

    return wide;

}

// 語の番地 → 詰めた並びの位置。work が 0 なら詰めない

__device__ __forceinline__ i64 slot(i64 a, int need, int work, int words) {

    if (work == 0) return (a >= 0 && a < words) ? a : -1;

    if (a >= 0 && a < need) return a;

    if (a >= work && a - work + need < words) return a - work + need;

    return -1;

}

extern "C" __global__

void lane64_run(const i64 *prog, int ninsn, u64 *M, int words, int lanes, i64 budget,

                int need, int work, int *report, i64 *steps,

                const u64 *ROM, i64 romlo, i64 romn)

{

    int lane = blockIdx.x * blockDim.x + threadIdx.x;

    if (lane >= lanes) return;

    Big R[NREG];

    int pc = 0; i64 n = 0; int rep = %(r_budget)d;

    while (n < budget) {

        n++;

        if (pc < 0 || pc >= ninsn) { rep = %(r_end)d; break; }

        const i64 *ins = prog + (size_t)pc * %(fields)d;

        int op = (int)ins[0], a = (int)ins[1], b = (int)ins[2], c = (int)ins[3];

        int sb = (int)ins[4], sc = (int)ins[5]; i64 im = ins[6];

        if (op == 0) { rep = %(r_end)d; break; }                         // HALT

        else if (op == 1) { big_imm(R[a], im); pc++; }                  // MAT

        else if (op == 2) { Big x, y; big_shl(x, R[b], sb); big_shl(y, R[c], sc); big_add(R[a], x, y); pc++; }   // ADD

        else if (op == 3) { Big t; if (big_mul(t, R[b], R[c])) { rep = %(r_wide)d; break; } R[a] = t; pc++; }  // MUL

        else if (op == 4) { Big t; big_negate(t, R[b]); R[a] = t; pc++; }   // NEG

        else if (op == 5) { Big x, y; big_shl(x, R[b], sb); big_shl(y, R[c], sc); pc += big_eq(x, y) ? (int)im : 1; }   // EQ

        else if (op == 6) { Big x, y; big_shl(x, R[b], sb); big_shl(y, R[c], sc); pc += big_lt(x, y) ? (int)im : 1; }   // LT

        else if (op == 7) { pc += (int)im; }                            // JMP

        else if (op == 8) {                                             // LD

            i64 raw = (i64)R[b].v[0];

            if (raw >= romlo) {                                         // 共有帯。lane で割らない

                i64 o = raw - romlo; if (o >= romn) { rep = %(r_end)d; break; }

                const u64 *w = ROM + (size_t)o * NLIMB;

                for (int i = 0; i < NLIMB; i++) R[a].v[i] = w[i]; pc++;

            } else {

                i64 addr = slot(raw, need, work, words); if (addr < 0) { rep = %(r_end)d; break; }

                const u64 *w = M + ((size_t)addr * NLIMB) * (size_t)lanes + lane;

                for (int i = 0; i < NLIMB; i++) R[a].v[i] = w[(size_t)i * lanes]; pc++;

            }

        }

        else if (op == 9) {                                             // ST

            i64 raw = (i64)R[b].v[0];

            if (raw >= romlo) { rep = %(r_end)d; break; }               // 共有帯へは書けない

            i64 addr = slot(raw, need, work, words); if (addr < 0) { rep = %(r_end)d; break; }

            u64 *w = M + ((size_t)addr * NLIMB) * (size_t)lanes + lane;

            for (int i = 0; i < NLIMB; i++) w[(size_t)i * lanes] = R[c].v[i]; pc++;

        }

        else if (op == 10) { rep = %(r_ok)d; break; }                   // OK

        else if (op == 11) { rep = %(r_ko)d; break; }                   // KO

        else if (op == 12) { pc++; }                                    // PUT（GPU では捨てる）

        else { rep = %(r_end)d; break; }

    }

    report[lane] = rep; steps[lane] = n;

}

// ===== 語が 64 bit で足りるときの道 =====

//

// Big R[NREG] は添字が命令から来るので、compiler が register に置けず local memory へ落ちる。

// 一命令ごとに local への往復が入り、thread 一本の歩みが桁で遅くなる。

// 語を 64 bit に決め打てる場合は、register を shared memory の

//

//     R[k]  =  sreg[k * blockDim.x + threadIdx.x]

//

// へ置く。stride が blockDim.x なので bank が衝突せず、添字が動いても速い。

// 積が 64 bit を越えたら幅の報告で止める（意味は Big の道と同じ）。

extern "C" __global__

void lane64_run1(const i64 *prog, int ninsn, u64 *M, int words, int lanes, i64 budget,

                 int need, int work, int *report, i64 *steps,

                 const u64 *ROM, i64 romlo, i64 romn)

{

    extern __shared__ i64 sreg[];

    int lane = blockIdx.x * blockDim.x + threadIdx.x;

    if (lane >= lanes) return;

    i64 *R = sreg + threadIdx.x;

    const int rs = blockDim.x;

    int pc = 0; i64 n = 0; int rep = %(r_budget)d;

    while (n < budget) {

        n++;

        if (pc < 0 || pc >= ninsn) { rep = %(r_end)d; break; }

        const i64 *ins = prog + (size_t)pc * %(fields)d;

        int op = (int)ins[0], a = (int)ins[1], b = (int)ins[2], c = (int)ins[3];

        int sb = (int)ins[4], sc = (int)ins[5]; i64 im = ins[6];

        if (op == 0) { rep = %(r_end)d; break; }

        else if (op == 1) { R[a * rs] = im; pc++; }

        else if (op == 2) { R[a * rs] = (R[b * rs] << sb) + (R[c * rs] << sc); pc++; }

        else if (op == 3) {

            i64 x = R[b * rs], y = R[c * rs];

            u64 ux = (u64)(x < 0 ? -x : x), uy = (u64)(y < 0 ? -y : y);

            if (__umul64hi(ux, uy) || ((ux * uy) >> 63)) { rep = %(r_wide)d; break; }

            R[a * rs] = x * y; pc++;

        }

        else if (op == 4) { R[a * rs] = -R[b * rs]; pc++; }

        else if (op == 5) { pc += ((R[b * rs] << sb) == (R[c * rs] << sc)) ? (int)im : 1; }

        else if (op == 6) { pc += ((R[b * rs] << sb) < (R[c * rs] << sc)) ? (int)im : 1; }

        else if (op == 7) { pc += (int)im; }

        else if (op == 8) {

            i64 raw = R[b * rs];

            if (raw >= romlo) {                                         // 共有帯。lane で割らない

                i64 o = raw - romlo; if (o >= romn) { rep = %(r_end)d; break; }

                R[a * rs] = (i64)ROM[(size_t)o]; pc++;

            } else {

                i64 addr = slot(raw, need, work, words); if (addr < 0) { rep = %(r_end)d; break; }

                R[a * rs] = (i64)M[(size_t)addr * (size_t)lanes + lane]; pc++;

            }

        }

        else if (op == 9) {

            i64 raw = R[b * rs];

            if (raw >= romlo) { rep = %(r_end)d; break; }               // 共有帯へは書けない

            i64 addr = slot(raw, need, work, words); if (addr < 0) { rep = %(r_end)d; break; }

            M[(size_t)addr * (size_t)lanes + lane] = (u64)R[c * rs]; pc++;

        }

        else if (op == 10) { rep = %(r_ok)d; break; }

        else if (op == 11) { rep = %(r_ko)d; break; }

        else if (op == 12) { pc++; }

        else { rep = %(r_end)d; break; }

    }

    report[lane] = rep; steps[lane] = n;

}

// ===== 多数の lane が一つの仕事を分ける道 =====

//

// 段と段の間では、lane i が書いた値を lane j が読む。そのために

//

//   1. 共有帯へ書ける（既存の道は報告 45 で止める）

//   2. 書いた後に grid 全体で足並みを揃える

//

// の二つが要る。合図は「共有帯への書き込み」そのものにする。境界へ書く所が

// チェンバーの接合点なので、pc を外から指さなくて済む。全 lane が段ごとに

// ちょうど一度ずつ書く形の平文でだけ使える。

//

// 足並み揃えは世代を反転させる素朴な形。block が全部同時に載っていないと

// 止まるので、grid は device に載る数までに抑える（host 側で決める）。

__device__ __forceinline__ void gbar(int *cnt, int *gen, int nblocks)

{

    __syncthreads();

    if (threadIdx.x == 0) {

        int g = *((volatile int *)gen);

        __threadfence();

        if (atomicAdd(cnt, 1) == nblocks - 1) {

            *cnt = 0;

            __threadfence();

            atomicExch(gen, g + 1);

        } else {

            while (*((volatile int *)gen) == g) __threadfence();

        }

    }

    __syncthreads();

}

extern "C" __global__

void lane64_runb(const i64 *prog, int ninsn, u64 *M, int words, int lanes, i64 budget,

                 int need, int work, int *report, i64 *steps,

                 u64 *ROM, i64 romlo, i64 romn, int *cnt, int *gen, int nblocks)

{

    extern __shared__ i64 sregb[];

    int lane = blockIdx.x * blockDim.x + threadIdx.x;

    i64 *R = sregb + threadIdx.x;

    const int rs = blockDim.x;

    // 足並みを揃える以上、休んでいる thread を作れない。lanes は grid x block に

    // ちょうど合わせる（余りの lane は自分の席で同じ文を走らせる）。

    int pc = 0; i64 n = 0; int rep = %(r_budget)d;

    while (n < budget) {

        n++;

        if (pc < 0 || pc >= ninsn) { rep = %(r_end)d; break; }

        const i64 *ins = prog + (size_t)pc * %(fields)d;

        int op = (int)ins[0], a = (int)ins[1], b = (int)ins[2], c = (int)ins[3];

        int sb = (int)ins[4], sc = (int)ins[5]; i64 im = ins[6];

        if (op == 0) { rep = %(r_end)d; break; }

        else if (op == 1) { R[a * rs] = im; pc++; }

        else if (op == 2) { R[a * rs] = (R[b * rs] << sb) + (R[c * rs] << sc); pc++; }

        else if (op == 3) {

            i64 x = R[b * rs], y = R[c * rs];

            u64 ux = (u64)(x < 0 ? -x : x), uy = (u64)(y < 0 ? -y : y);

            if (__umul64hi(ux, uy) || ((ux * uy) >> 63)) { rep = %(r_wide)d; break; }

            R[a * rs] = x * y; pc++;

        }

        else if (op == 4) { R[a * rs] = -R[b * rs]; pc++; }

        else if (op == 5) { pc += ((R[b * rs] << sb) == (R[c * rs] << sc)) ? (int)im : 1; }

        else if (op == 6) { pc += ((R[b * rs] << sb) < (R[c * rs] << sc)) ? (int)im : 1; }

        else if (op == 7) { pc += (int)im; }

        else if (op == 8) {

            i64 raw = R[b * rs];

            if (raw >= romlo) {

                i64 o = raw - romlo; if (o >= romn) { rep = %(r_end)d; break; }

                R[a * rs] = (i64)ROM[(size_t)o]; pc++;

            } else {

                i64 addr = slot(raw, need, work, words); if (addr < 0) { rep = %(r_end)d; break; }

                R[a * rs] = (i64)M[(size_t)addr * (size_t)lanes + lane]; pc++;

            }

        }

        else if (op == 9) {

            i64 raw = R[b * rs];

            if (raw >= romlo) {                    // 共有帯へ書く。ここが接合点

                i64 o = raw - romlo; if (o >= romn) { rep = %(r_end)d; break; }

                ROM[(size_t)o] = (u64)R[c * rs];

                pc++;

                __threadfence();

                gbar(cnt, gen, nblocks);

                continue;

            }

            i64 addr = slot(raw, need, work, words); if (addr < 0) { rep = %(r_end)d; break; }

            M[(size_t)addr * (size_t)lanes + lane] = (u64)R[c * rs]; pc++;

        }

        else if (op == 10) { rep = %(r_ok)d; break; }

        else if (op == 11) { rep = %(r_ko)d; break; }

        else if (op == 12) { pc++; }

        else { rep = %(r_end)d; break; }

    }

    report[lane] = rep; steps[lane] = n;

}

// ===== warp を pc の小さい順にそろえて歩く道 =====

//

// thread ごとに好きに飛ぶと、warp の中で pc が割れた瞬間から命令が直列になり、

// 語の読みも散って memory に縛られる（電力 50 W で頭打ち）。

// 20 §6 命題 6.1: 合流は最も小さい pc から取れば足りる。warp の中で最小の pc を

// 一つ選び、そこに居る lane だけが一命令進む。分岐で割れても、後ろの lane が

// 先に追いつくので、loop の出口で自然に合流する。

extern "C" __global__

void lane64_run1_sync(const i64 *prog, int ninsn, u64 *M, int words, int lanes, i64 budget,

                      int need, int work, int *report, i64 *steps,

                      const u64 *ROM, i64 romlo, i64 romn, int stride)

{

    // stride は語の帯の中の lane の数（M2 の列数）。lanes はいま走らせる lane の数で、

    // 先頭から lanes 本だけが歩く。段を長さの順に並べておけば、短い段が終わった後の

    // 歩は先頭の生きている lane だけで済む

    extern __shared__ i64 sreg[];

    int lane = blockIdx.x * blockDim.x + threadIdx.x;

    bool alive = lane < lanes;

    i64 *R = sreg + threadIdx.x;

    const int rs = blockDim.x;

    int pc = 0; i64 n = 0; int rep = %(r_budget)d;

    bool done = !alive;

    while (true) {

        unsigned live = __ballot_sync(0xffffffffu, !done);

        if (live == 0) break;

        int mypc = done ? 0x7fffffff : pc;

        int pcmin = __reduce_min_sync(0xffffffffu, mypc);

        if (done || pc != pcmin) continue;

        n++;

        if (n > budget) { done = true; continue; }

        if (pc < 0 || pc >= ninsn) { rep = %(r_end)d; done = true; continue; }

        const i64 *ins = prog + (size_t)pc * %(fields)d;

        int op = (int)ins[0], a = (int)ins[1], b = (int)ins[2], c = (int)ins[3];

        int sb = (int)ins[4], sc = (int)ins[5]; i64 im = ins[6];

        if (op == 0) { rep = %(r_end)d; done = true; }

        else if (op == 1) { R[a * rs] = im; pc++; }

        else if (op == 2) { R[a * rs] = (R[b * rs] << sb) + (R[c * rs] << sc); pc++; }

        else if (op == 3) {

            i64 x = R[b * rs], y = R[c * rs];

            u64 ux = (u64)(x < 0 ? -x : x), uy = (u64)(y < 0 ? -y : y);

            if (__umul64hi(ux, uy) || ((ux * uy) >> 63)) { rep = %(r_wide)d; done = true; }

            else { R[a * rs] = x * y; pc++; }

        }

        else if (op == 4) { R[a * rs] = -R[b * rs]; pc++; }

        else if (op == 5) { pc += ((R[b * rs] << sb) == (R[c * rs] << sc)) ? (int)im : 1; }

        else if (op == 6) { pc += ((R[b * rs] << sb) < (R[c * rs] << sc)) ? (int)im : 1; }

        else if (op == 7) { pc += (int)im; }

        else if (op == 8) {

            i64 raw = R[b * rs];

            if (raw >= romlo) {

                i64 o = raw - romlo; if (o >= romn) { rep = %(r_end)d; done = true; }

                else { R[a * rs] = (i64)ROM[(size_t)o]; pc++; }

            } else {

                i64 addr = slot(raw, need, work, words); if (addr < 0) { rep = %(r_end)d; done = true; }

                else { R[a * rs] = (i64)M[(size_t)addr * (size_t)stride + lane]; pc++; }

            }

        }

        else if (op == 9) {

            i64 raw = R[b * rs];

            if (raw >= romlo) { rep = %(r_end)d; done = true; }

            else {

                i64 addr = slot(raw, need, work, words); if (addr < 0) { rep = %(r_end)d; done = true; }

                else { M[(size_t)addr * (size_t)stride + lane] = (u64)R[c * rs]; pc++; }

            }

        }

        else if (op == 10) { rep = %(r_ok)d; done = true; }

        else if (op == 11) { rep = %(r_ko)d; done = true; }

        else if (op == 12) { pc++; }

        else { rep = %(r_end)d; done = true; }

    }

    if (alive) { report[lane] = rep; steps[lane] = n; }

}

extern "C" __global__

void gather(const u64 *M, int words, int lanes, const int *watch, int nwatch, u64 *out)

{

    int lane = blockIdx.x * blockDim.x + threadIdx.x;

    if (lane >= lanes) return;

    u64 *o = out + (size_t)lane * nwatch * NLIMB;

    for (int k = 0; k < nwatch; k++) {

        const u64 *w = M + ((size_t)watch[k] * NLIMB) * (size_t)lanes + lane;

        for (int i = 0; i < NLIMB; i++) o[k * NLIMB + i] = w[(size_t)i * lanes];

    }

}

// 先頭 lanes 本だけを集める。stride は帯の中の lane の数

extern "C" __global__

void gather_s(const u64 *M, int stride, int lanes, const int *watch, int nwatch, u64 *out)

{

    int lane = blockIdx.x * blockDim.x + threadIdx.x;

    if (lane >= lanes) return;

    u64 *o = out + (size_t)lane * nwatch;

    for (int k = 0; k < nwatch; k++) o[k] = M[(size_t)watch[k] * stride + lane];

}

"""

R_END = 45

_MOD = {}

# ===== block 内の即値を畳む =====

#

# commonc が出す lane64 の命令列は、番地を one の倍化と加算で組み立てる。

# CPU 側の vcpu.translate は block ごとに値を記号で追い、即値どうしの算術を

# 一つの即値へ畳む（casm.lane64.vn）。GPU の解釈器は命令を一つずつ歩くので、

# 同じ畳み込みを先に済ませておかないと、番地を綴る鎖がそのまま実行費用になる。

#

#   R1 = 1 ; R1 = R1+R1 ; ... ; R2 = R1+R3   →   R1 = 65554 の一命令

#

# block の入口で値が分かっている register は無い、として始める。分岐は畳まない。

MAT, ADD, MUL, NEG, EQ, LT, JMP, LD, ST = 1, 2, 3, 4, 5, 6, 7, 8, 9

_BR = (EQ, LT, JMP)

_TERM = (0, 10, 11)

def leaders(prog):

    """block の先頭になる命令の番号"""

    out = set([0])

    for i, ins in enumerate(prog):

        op, im = ins[0], ins[6]

        if op in _BR:

            out.add(i + int(im))

            out.add(i + 1)

        elif op in _TERM:

            out.add(i + 1)

    return sorted(a for a in out if 0 <= a < len(prog))

def _reads(ins):

    op = ins[0]

    if op in (ADD, EQ, LT, MUL):

        return (ins[2], ins[3])

    if op in (NEG, LD):

        return (ins[2],)

    if op == ST:

        return (ins[2], ins[3])

    if op == 12:

        return (ins[2],)

    return ()

def _writes(ins):

    op = ins[0]

    return ins[1] if op in (MAT, ADD, MUL, NEG, LD) else None

def _walk(prog, lo, hi, live, dead, rewrite):

    """block をひと通り追う。live には「自分の Store より前に読んだ定数番地」を集め、

    rewrite が真なら即値の算術と語への往復を畳み、dead の番地への Store を落とす。"""

    val, mem = {}, {}

    for i in range(lo, hi):

        ins = prog[i]

        if ins is None:

            continue

        op = ins[0]

        if op == MAT:

            val[ins[1]] = int(ins[6])

        elif op in (ADD, MUL, NEG):

            b, c, sb, sc = ins[2], ins[3], ins[4], ins[5]

            vb = val.get(b)

            vc = val.get(c) if op != NEG else 0

            if vb is not None and vc is not None:

                v = (vb << sb) + (vc << sc) if op == ADD else (vb * vc if op == MUL else -vb)

                if rewrite:

                    prog[i] = [MAT, ins[1], 0, 0, 0, 0, v]

                val[ins[1]] = v

            else:

                val.pop(ins[1], None)

        elif op == LD:

            a = val.get(ins[2])

            if a is not None and a in mem:

                if rewrite:

                    prog[i] = [MAT, ins[1], 0, 0, 0, 0, mem[a]]

                val[ins[1]] = mem[a]

            else:

                if a is not None:

                    live.add(a)

                val.pop(ins[1], None)

        elif op == ST:

            a, v = val.get(ins[2]), val.get(ins[3])

            if a is None:

                mem.clear()

            elif v is None:

                mem.pop(a, None)

            else:

                mem[a] = v

            # 同じ block の Load は先に繋いでから落とす。落とすのは block を跨ぐ分だけ

            if rewrite and a is not None and a in dead:

                prog[i] = None

        else:

            w = _writes(ins)

            if w is not None:

                val.pop(w, None)

def fold(prog, wires=None):

    """block の中で値を記号で追い、即値どうしの算術と語への往復を畳む。意味は変わらない。

    lane64 の一箱は MAT・ST・LD の往復に開くので、register だけを追っても畳めない。

    定数番地への Store と、その番地からの Load を繋いで初めて鎖が一つの即値になる。

    """

    prog = [list(x) for x in prog]

    heads = leaders(prog) + [len(prog)]

    # 一巡目。block を跨いで読まれる定数番地を集める

    live = set()

    for k in range(len(heads) - 1):

        _walk(prog, heads[k], heads[k + 1], live, (), False)

    # commonc が出すワイヤの番地は定数だけで、本文の番地は WORK 以上にある。

    # だから [64, need) のうち誰も跨いで読まない番地への Store は落としてよい

    dead = set()

    if wires:

        lo, hi = wires

        dead = set(a for a in range(lo, hi) if a not in live)

    # 二巡目。畳んで、死んだ Store を落とす

    for k in range(len(heads) - 1):

        lo, hi = heads[k], heads[k + 1]

        _walk(prog, lo, hi, set(), dead, True)

    # block の中だけの死んだ書き込みを落とす

    for k in range(len(heads) - 1):

        lo, hi = heads[k], heads[k + 1]

        last, drop = {}, set()

        for i in range(lo, hi):

            ins = prog[i]

            if ins is None:

                continue

            for r in _reads(ins):

                last.pop(r, None)

            w = _writes(ins)

            if w is not None:

                if ins[0] == MAT and w in last:

                    drop.add(last[w])

                last[w] = i if ins[0] == MAT else None

                if last[w] is None:

                    last.pop(w)

        for i in drop:

            prog[i] = None

    # 落とした分だけ番号が詰まるので、分岐の変位を貼り直す

    keep = [i for i, x in enumerate(prog) if x is not None]

    newpos, n = {}, 0

    for i in range(len(prog)):

        newpos[i] = n

        if prog[i] is not None:

            n += 1

    newpos[len(prog)] = n

    out = []

    for i in keep:

        ins = list(prog[i])

        if ins[0] in _BR:

            ins[6] = newpos[min(i + int(ins[6]), len(prog))] - newpos[i]

        out.append(tuple(ins))

    return out

def compact(prog, wires, work=0):

    """ワイヤの番地を詰める。返りは (prog, 新しい need)。

    commonc の一箱はワイヤ一本を語 64 + w に取るので、image が大きいと lane ごとの

    ワイヤの帯が数万語になる。畳んだ後に生き残る番地はその一部なので、番地として

    使われている即値だけを 64 から順に付け直す。

    即値が番地か値かは block の中の使われ方で決める。MAT で置いた register が、

    上書きされるまでに LD / ST の番地としてしか使われなければ番地。値として

    使われた即値は触らない。触らない即値がワイヤの帯に残れば、need は

    その最大までしか下がらない（返りの need がそれを示す）。

    """

    lo, hi = wires

    # 二値の八箱の退避語は定数番地の最大の上に置かれ、ワイヤの帯と本文の間の穴に落ちる

    # ことがある（16 §5.3）。穴は lane に無いので、そこも一緒に詰める

    hi = max(hi, work) if work else hi

    prog = [list(x) for x in prog]

    heads = leaders(prog) + [len(prog)]

    addr_mats, keep_hi = [], lo

    for k in range(len(heads) - 1):

        a, b = heads[k], heads[k + 1]

        for i in range(a, b):

            ins = prog[i]

            if ins[0] != MAT or not (lo <= ins[6] < hi):

                continue

            reg, ok, closed = ins[1], True, False

            for j in range(i + 1, b):

                jn = prog[j]

                op = jn[0]

                if op in (LD, ST) and jn[2] == reg and not (op == ST and jn[3] == reg):

                    pass

                elif reg in _reads(jn):

                    ok = False

                    break

                if _writes(jn) == reg:

                    closed = True

                    break

            # block の外へ出る register は次の block でも番地として使われる（commonc は

            # ワイヤの番地を定数でしか作らない）。値として読まれた即値だけを残す

            if ok:

                addr_mats.append(i)

            else:

                keep_hi = max(keep_hi, ins[6] + 1)

    used = sorted(set(prog[i][6] for i in addr_mats))

    base = keep_hi

    remap = {a: base + k for k, a in enumerate(used)}

    for i in addr_mats:

        prog[i][6] = remap[prog[i][6]]

    return [tuple(x) for x in prog], base + len(used)

def _module(cp, nlimb):

    if nlimb not in _MOD:

        src = SRC % {"nlimb": nlimb, "fields": INSN_FIELDS, "r_ok": R_OK, "r_ko": R_KO,

                     "r_wide": R_WIDE, "r_budget": R_BUDGET, "r_end": R_END}

        _MOD[nlimb] = cp.RawModule(code=src, options=("-std=c++11",))

    return _MOD[nlimb]

def to_limbs(v, nlimb):

    """整数を limb の並びへ（二の補数）"""

    mask = (1 << (64 * nlimb)) - 1

    x = v & mask

    return [(x >> (64 * i)) & 0xFFFFFFFFFFFFFFFF for i in range(nlimb)]

def from_limbs(limbs, nlimb):

    x = 0

    for i, l in enumerate(limbs):

        x |= int(l) << (64 * i)

    if x >> (64 * nlimb - 1):

        x -= 1 << (64 * nlimb)

    return x

_NOROM_BUF = {}

def _NOROM(cp):

    """共有帯を使わないときに渡す一語の置き場。kernel は romlo = 2^62 で触らない"""

    if "b" not in _NOROM_BUF:

        _NOROM_BUF["b"] = cp.zeros((1,), dtype=cp.uint64)

    return _NOROM_BUF["b"]

def pack(addr, need, work):

    """語番号 → 詰めた並びの位置。kernel の slot() と同じ式（式は device 側 §slot）。"""

    if work == 0:

        return addr

    if 0 <= addr < need:

        return addr

    return addr - work + need

def run(prog, seeds, words, watch, nlimb=16, budget=1 << 26, chunk=64, need=0, work=0,

        base=None):

    """image（decode 済み）を seeds の lane 分だけ走らせる。

    seeds は lane ごとの {語: 値}、watch は host へ戻す語の並び。

    返りは lane ごとの (報告, {語: 値})。報告は 0 = tu、42 = te、43 = 幅を越えた、44 = 歩数の上限。

    base は全 lane へ同じ値を置く語（{語: 値}）。模型の表のように lane ごとに同じ物を

    渡すとき、seeds へ毎回入れると host の側だけで lane × 語 の手間になる。ここでは

    語ごとに一度だけ device の帯へ書くので、手間は語の数で済む。

    need と work を渡すと、ワイヤ（0 .. need）と本文の作業領域（work 以上）の間の穴を詰める。

    commonc の lane64 版は本文を WORK = 2^16 から置くので、詰めないと lane ごとに

    65536 語ぶんの空きを確保することになる（20 §1.2 の lane の作業領域は交わらない）。

        words = need + (本文が使う語の数)

    穴を詰めると lane あたりの領域が減り、同じ memory でより多くの lane を並べられる。

    """

    import cupy as cp

    mod = _module(cp, nlimb)

    k_run = mod.get_function("lane64_run1" if nlimb == 1 else "lane64_run")

    k_gather = mod.get_function("gather")

    prog = fold(prog)

    table = np.array([list(p) for p in prog], dtype=np.int64).reshape(-1)

    d_prog = cp.asarray(table)

    d_watch = cp.asarray(np.array([pack(w, need, work) for w in watch], dtype=np.int32))

    out = []

    for start in range(0, len(seeds), chunk):

        part = seeds[start:start + chunk]

        lanes = len(part)

        M = cp.zeros((lanes * words * nlimb,), dtype=cp.uint64)

        if base:

            for addr, v in base.items():

                a = pack(addr, need, work)

                v = int(v)

                if 0 <= v < (1 << 63):

                    M[(a * nlimb) * lanes:(a * nlimb) * lanes + lanes] = np.uint64(v)

                else:

                    for i, limb in enumerate(to_limbs(v, nlimb)):

                        if limb:

                            o = (a * nlimb + i) * lanes

                            M[o:o + lanes] = np.uint64(limb)

        idx, vals = [], []

        for l, seed in enumerate(part):

            for addr, v in seed.items():

                a = pack(addr, need, work)

                v = int(v)

                if 0 <= v < (1 << 63):

                    idx.append((a * nlimb) * lanes + l)

                    vals.append(v)

                else:

                    for i, limb in enumerate(to_limbs(v, nlimb)):

                        if limb:

                            idx.append((a * nlimb + i) * lanes + l)

                            vals.append(limb)

        if idx:

            M[cp.asarray(np.array(idx, dtype=np.int64))] = cp.asarray(np.array(vals, dtype=np.uint64))

        report = cp.zeros((lanes,), dtype=cp.int32)

        steps = cp.zeros((lanes,), dtype=cp.int64)

        threads = 128

        grid = (lanes + threads - 1) // threads

        sh = threads * NREG_BYTES if nlimb == 1 else 0

        k_run((grid,), (threads,), (d_prog, np.int32(len(prog)), M, np.int32(words), np.int32(lanes),

                                    np.int64(budget), np.int32(need), np.int32(work), report, steps,

                                    _NOROM(cp), np.int64(1 << 62), np.int64(0)),

              shared_mem=sh)

        got = cp.zeros((lanes * len(watch) * nlimb,), dtype=cp.uint64)

        if watch:

            k_gather((grid,), (threads,), (M, np.int32(words), np.int32(lanes), d_watch, np.int32(len(watch)), got))

        rep = report.get()

        raw = got.get()

        for l in range(lanes):

            vals_l = {}

            for k, w in enumerate(watch):

                o = (l * len(watch) + k) * nlimb

                vals_l[w] = from_limbs(raw[o:o + nlimb], nlimb)

            out.append((int(rep[l]), vals_l))

        del M, got

    return out

def run_shared(prog, base, lanes, words, watch, need, work, nlimb=4, budget=1 << 26,

               chunk=2048, romlo=None, romn=0):

    """lane が共有する種 base に、lane ごとの少しの語 lanes[k] を重ねて走らせる。

    語の並びは詰める: ワイヤ [0, need) をそのまま、作業領域 [work, ...) を need の後ろに置く。

    words は詰めた後の語数（need + 作業領域の幅）。平文の側の語は全部 work 以上にあるので、

    間の番地は使われない。返りは lane ごとの (報告, {語: 値})。

    romlo を渡すと、番地 [romlo, romlo + romn) は**全 lane で一枚だけ**持つ。

    読むだけの表（値段の表など）をここへ置くと、lane あたりの重さが桁で落ちる。

        lane あたりの語数 = need + (romlo - work)

        共有帯の語数     = romn

    この帯へ Store すると報告 45 で止まる。書く語を入れてはいけない。

    """

    import cupy as cp

    mod = _module(cp, nlimb)

    k_run = mod.get_function("lane64_run1" if nlimb == 1 else "lane64_run")

    k_gather = mod.get_function("gather")

    prog = fold(prog, wires=(64, need))

    table = np.array([list(p) for p in prog], dtype=np.int64).reshape(-1)

    d_prog = cp.asarray(table)

    hi = romlo if romlo is not None else (1 << 62)

    def pos(a):

        return a if a < need else a - work + need

    d_watch = cp.asarray(np.array([pos(w) for w in watch], dtype=np.int32))

    # 共有の種を一枚の並びにする。共有帯の分は別の並びへ分ける

    row = np.zeros((words * nlimb,), dtype=np.uint64)

    rom = np.zeros((max(1, romn) * nlimb,), dtype=np.uint64)

    for addr, v in base.items():

        v = int(v)

        p0 = (addr - hi) * nlimb if addr >= hi else pos(addr) * nlimb

        dst = rom if addr >= hi else row

        if 0 <= v < (1 << 63):

            dst[p0] = v

        else:

            for i, limb in enumerate(to_limbs(v, nlimb)):

                dst[p0 + i] = limb

    d_row = cp.asarray(row)

    d_rom = cp.asarray(rom)

    out = []

    for start in range(0, len(lanes), chunk):

        part = lanes[start:start + chunk]

        n = len(part)

        M = cp.repeat(d_row, n)

        idx, vals = [], []

        for l, seed in enumerate(part):

            for addr, v in seed.items():

                idx.append((pos(addr) * nlimb) * n + l)

                vals.append(int(v))

        if idx:

            M[cp.asarray(np.array(idx, dtype=np.int64))] = cp.asarray(np.array(vals, dtype=np.uint64))

        report = cp.zeros((n,), dtype=cp.int32)

        steps = cp.zeros((n,), dtype=cp.int64)

        threads = 128

        grid = (n + threads - 1) // threads

        sh = threads * NREG_BYTES if nlimb == 1 else 0

        k_run((grid,), (threads,), (d_prog, np.int32(len(prog)), M, np.int32(words), np.int32(n),

                                    np.int64(budget), np.int32(need), np.int32(work), report, steps,

                                    d_rom, np.int64(hi), np.int64(romn)),

              shared_mem=sh)

        got = cp.zeros((n * len(watch) * nlimb,), dtype=cp.uint64)

        k_gather((grid,), (threads,), (M, np.int32(words), np.int32(n), d_watch, np.int32(len(watch)), got))

        rep = report.get()

        raw = got.get().reshape(n, len(watch), nlimb)

        for l in range(n):

            out.append((int(rep[l]), dict((w, from_limbs(raw[l, k], nlimb)) for k, w in enumerate(watch))))

        del M, got

    return out

class Lanes(object):

    """lane の作業領域と共有帯を device に置いたまま、image を何度も流す。

    host は語を一つずつ触らない。種は (lane, 語, 値) の並びを一度に散らし、読みは

    (lane × 語) の並びを一度に集める。どちらも cupy の一手で済むので、lane の数にも

    語の数にも Python の反復が付かない。

        L = Lanes(prog, need, work, romlo, romn, nlanes, base={語: 値})

        L.put(lanes, addrs, vals)          種を散らす（三つとも同じ長さの並び）

        L.rom_put(offset, values)          共有帯の一部を書く

        rep, W = L.run(watch)              rep[lane] が報告、W[lane, k] が watch[k] の値（device）

        L.M2                               語 × lane の二次元の眺め（device）

        L.reset()                          作業領域を base の並びへ戻す

    語の並びは run_shared と同じ: ワイヤ [0, need) をそのまま、作業領域 [work, romlo) を

    need の後ろに詰める。M は語主で M2[pos(語), lane]。

    共有帯 rom は cupy の並び一枚。別の Lanes の M をそのまま rom に渡せる（rom=）ので、

    lane の畳みのように「全 lane の語を読む」kernel は、その M を共有帯として読めば書ける。

    """

    def __init__(self, prog, need, work, romlo, romn, nlanes, base=None, nlimb=1, rom=None, sync=True):

        """sync が真なら warp を pc の小さい順にそろえて歩く（lane64_run1_sync）。"""

        import cupy as cp

        self.cp = cp

        if nlimb != 1:

            raise ValueError("Lanes は 64 bit の語だけを扱う（nlimb=1）")

        self.need, self.work, self.romlo = need, work, romlo

        self.words = need + (romlo - work)

        self.n = nlanes

        mod = _module(cp, 1)

        self.sync = sync

        self.k_run = mod.get_function("lane64_run1_sync" if sync else "lane64_run1")

        self.k_gather = mod.get_function("gather_s")

        prog = fold(prog, wires=(64, need))

        prog, need = compact(prog, (64, need), work=work)

        self.need = need

        self.words = need + (romlo - work)

        self.nprog = len(prog)

        self.d_prog = cp.asarray(np.array([list(p) for p in prog], dtype=np.int64).reshape(-1))

        if rom is not None:

            self.d_rom = rom

            self.romn = int(rom.size)

        else:

            self.romn = romn

            self.d_rom = cp.zeros((max(1, romn),), dtype=cp.uint64)

        row = np.zeros((self.words,), dtype=np.uint64)

        for addr, v in (base or {}).items():

            row[self.pos(addr)] = int(v)

        self.d_row = cp.asarray(row)

        self.M = cp.zeros((self.words * nlanes,), dtype=cp.uint64)

        self.reset()

    def pos(self, a):

        return a if a < self.need else a - self.work + self.need

    @property

    def M2(self):

        return self.M.reshape(self.words, self.n)

    def reset(self):

        self.M2[:, :] = self.d_row.reshape(-1, 1)

    def put(self, lanes, addrs, vals):

        """(lane, 語, 値) の並びを一度に散らす。addrs は語番号（pos は中で当てる）。"""

        cp = self.cp

        lanes = cp.asarray(np.asarray(lanes, dtype=np.int64))

        addrs = np.asarray(addrs, dtype=np.int64)

        addrs = np.where(addrs < self.need, addrs, addrs - self.work + self.need)

        idx = cp.asarray(addrs) * self.n + lanes

        self.M[idx] = cp.asarray(np.asarray(vals, dtype=np.uint64))

    def put_block(self, addr, block):

        """語 addr から続く語の帯へ、(語数 × lane) の並びを一度に置く。block は device でも host でもよい。"""

        cp = self.cp

        b = cp.asarray(block, dtype=cp.uint64)

        self.M2[self.pos(addr):self.pos(addr) + b.shape[0], :] = b

    def get_block(self, addr, n):

        """語 addr から n 語の帯を (n × lane) で返す（device）。"""

        return self.M2[self.pos(addr):self.pos(addr) + n, :]

    def rom_put(self, offset, values):

        self.d_rom[offset:offset + len(values)] = self.cp.asarray(np.asarray(values, dtype=np.uint64))

    def run(self, watch=(), budget=1 << 36, active=None):

        """先頭 active 本の lane を走らせる（省くと全部）。返りは (報告, watch の値) で、どちらも active 本ぶん。"""

        cp = self.cp

        n = self.n if active is None else min(active, self.n)

        report = cp.zeros((n,), dtype=cp.int32)

        steps = cp.zeros((n,), dtype=cp.int64)

        threads = 128

        grid = (n + threads - 1) // threads

        if self.sync:

            args = (self.d_prog, np.int32(self.nprog), self.M, np.int32(self.words), np.int32(n),

                    np.int64(budget), np.int32(self.need), np.int32(self.work), report, steps,

                    self.d_rom, np.int64(self.romlo), np.int64(self.romn), np.int32(self.n))

        else:

            if n != self.n:

                raise ValueError("先頭だけを走らせるには sync の道が要る")

            args = (self.d_prog, np.int32(self.nprog), self.M, np.int32(self.words), np.int32(n),

                    np.int64(budget), np.int32(self.need), np.int32(self.work), report, steps,

                    self.d_rom, np.int64(self.romlo), np.int64(self.romn))

        self.k_run((grid,), (threads,), args, shared_mem=threads * NREG_BYTES)

        W = None

        if watch:

            d_watch = cp.asarray(np.array([self.pos(w) for w in watch], dtype=np.int32))

            got = cp.zeros((n * len(watch),), dtype=cp.uint64)

            self.k_gather((grid,), (threads,), (self.M, np.int32(self.n), np.int32(n), d_watch,

                                                np.int32(len(watch)), got))

            W = got.reshape(n, len(watch)).astype(cp.int64)

        return report.get(), W