LoginSignup
4
1

More than 3 years have passed since last update.

@bellbind

有限体のプログラミング中編: JavaScriptで実装する有限体

Last updated at Posted at 2020-03-28

前編内容

  1. 素数べき要素数の有限体のための基礎知識
  2. 素数べき要素数の有限体とその作り方
  3. 素数べき要素数の有限体の例
  4. 素数べき要素数の有限体のガロア群

中編内容

  1. JavaScriptで実装する素数べき要素数有限体
  2. 2べき要素数の有限体に特化した実装

後編内容

  1. 有限体の応用: リードソロモン符号とBCH符号

注: 記事内のJavaScriptコードについて

この記事内のコードで用いるJavaScriptの仕様は、Array.flatMap()のあるECMAScript2019仕様までを用います。JSファイルは、すべて拡張子は.jsのES moduleです。

非ECMAScriptな要素として、console.logとconsole.assertのみを用います。その他のWeb APIやnode.jsビルトイン等は一切使用していません。

この記事内のコードは、以下のgistに置いています。

多項式や有限体等の要素は、数値とArrayにマッピングし、それらのデータ処理をfunction式の中で記述します。体や多項式等は、このfunction実装をメンバーに持つObjectとして提供します。たとえば、有理数体Qでの加算は以下のようになります。

  • $c = \frac{1}{2} \times -\frac{3}{4}$ : const c = Q.mul([1, 2], [-3, 4]);

JavaScriptで素数べき要素数有限体

ここから、プログラム上で素数べき要素数の有限体GF($p^n$)の扱い方についての内容となります。

GF($p^n$)のプログラムは、

  • ユーティリティ
  • 素数要素数の有限体(Prime Field): PF(p)
  • 体Kを係数に持つ多項式: Polynomial(K)
  • 素数べき要素数の有限体 GF($p^n$)

の順に実装します。そのあとで、多項式表現としてビット列を用いるp=2に特化したGF($2^n$)実装について触れます。

この有限体実装を使う応用として、前編の加算表や乗算表のmarkdownテーブルを出すコードを、そして、後編でエラー訂正可能なリードソロモン符号の実装を行います。

ユーティリティ: field-utils.js


コード全体
field-utils.js
// Array utils
export const range = n => [...Array(n).keys()];
export const uniq = (a, eql = (e1, e2) => e1 === e2) =>
      a.filter((e1, i) => a.findIndex(e2 => eql(e1, e2)) === i);
export const transpose = a => range(a[0].length).map(k => a.map(e => e[k]));

// Field suppliment: generate isZero, sub, div, sum, prod
export const Field = K => {
  const {eql, zero, one, add, mul, neg, inv} = K;
  const sub = (a, b) => add(a, neg(b));
  const div = (a, b) => mul(a, inv(b));
  const sum = es => es.reduce((r, e) => add(r, e), zero());
  const prod = es => es.reduce((r, e) => mul(r, e), one());
  const isZero = e => eql(e, zero());
  return {sub, div, sum, prod, isZero, ...K};
};

// Example Field Q: Rational number
export const Q = () => {
  const gcd = (a, b) => a === 0 ? b : gcd(b % a, a);
  const normalize = ([e0, e1]) => {
    const f = gcd(e1, e0);
    return [e0 / f, e1 / f];
  };

  const eql = (a, b) => a[0] * b[1] === a[1] * b[0];
  const zero = () => [0, 1];
  const one = () => [1, 1];
  const neg = e => [-e[0], e[1]];
  const inv = e => normalize([e[1], e[0]]);
  const add = (a, b) => normalize([a[0] * b[1] + b[0] * a[1], a[1] * b[1]]);
  const mul = (a, b) => normalize([a[0] * b[0], a[1] * b[1]]);

  const times = (e, k) => normalize([e[0] * k, e[1]]);
  const pow = (e, k) => normalize([e[0] ** k, e[1] ** k]);
  const toStr = e => e[1] === 1 ? `${e[0]}` : `${e[0]}/${e[1]}`;
  const toNum = e => Math.round(e[0] / e[1]);
  const fromNum = n => [n, 1];

  return Field({
    eql, zero, one, neg, inv, add, mul, times, pow, toStr, toNum, fromNum});
};

// Field utilities
export const FieldUtils = K => {
  const {isZero, zero, one, neg, inv, sub, mul, div} = K;
  // linear equations solver on the field f
  const pivoting = (r, u, i) => {
    for (let v = u + 1; v < r.length; v++) {
      if (!isZero(r[v][i])) {
        [r[u], r[v]] = [r[v], r[u]];
        break;
      }
    }
  };

  const upperTriangle = r => {
    const n = r[0].length - 1, m = r.length;
    for (let i = 0, u = 0; i < n && u < m; i++) {
      if (isZero(r[u][i])) pivoting(r, u, i);
      if (isZero(r[u][i])) continue;
      for (let v = u + 1; v < m; v++) {
        const c = div(r[v][i], r[u][i]);
        for (let j = i; j < n + 1; j++) {
          r[v][j] = sub(r[v][j], mul(c, r[u][j]));
        }
      }
      u++;
    }
  };

  const rank = lineqs => {
    const r = lineqs.map(lineq => lineq.slice());
    upperTriangle(r);
    return r.filter(eq => !eq.slice(0, -1).every(isZero)).length;
  };

  const solve = lineqs => {
    const r = lineqs.map(lineq => lineq.slice());
    upperTriangle(r);
    const n = r[0].length - 1, m = r.length;
    for (let i = n - 1, u = m - 1; i >= 0 && u >= 0; i--) {
      while (u > 0 && isZero(r[u][i])) u--;
      if (isZero(r[u][i])) continue;
      r[u][n] = div(r[u][n], r[u][i]);
      r[u][i] = one();
      for (let v = u - 1; v >= 0; v--) {
        r[v][n] = sub(r[v][n], mul(r[v][i], r[u][n]));
        r[v][i] = zero();
      }
    }
    return r.slice(0, n).map(l => neg(l[n]));
  };
  return {rank, solve};
};


JavaScript配列用ユーティリティ関数群

field-utils.js[1]
export const range = n => [...Array(n).keys()];
export const uniq = (a, eql = (e1, e2) => e1 === e2) =>
      a.filter((e1, i) => a.findIndex(e2 => eql(e1, e2)) === i);
export const transpose = a => range(a[0].length).map(k => a.map(e => e[k]));
  • range(n): 0からn-1までの数の列挙配列を返す
  • uniq(a, eql): aのうちの重複要素を除いた配列を返す
  • transpose(m): 二階配列mの転置配列を返す

二階配列aの転置配列tというのは、a[i][j]の値がt[j][i]に入っている二階配列tのことです。
たとえば、transpose([[1,2,3], ["a","b","c"]])の結果は、[[1,"a"], [2,"b"], [3,"c"]]です。

体の機能補完

field-utils.js[2]
export const Field = K => {
  const {eql, zero, one, add, mul, neg, inv} = K;

  const sub = (a, b) => add(a, neg(b));
  const div = (a, b) => mul(a, inv(b));
  const sum = es => es.reduce((r, e) => add(r, e), zero());
  const prod = es => es.reduce((r, e) => mul(r, e), one());
  const isZero = e => eql(e, zero());

  return {sub, div, sum, prod, isZero, ...K};
};

Field(K)は、基本となる関数群eql(a,b)、zero()、one()、add(a,b)、mul(a,b)、neg(e)、inv(e)を備える体Kをもとに、以下の関数実装を補完します。

  • sub(a,b): 減算: $a - b$
  • div(a,b): 除算: $a / b$
  • sum(es): 体の要素の配列esの総和: $\sum_{i}{e_i}$
  • prod(es): 体の要素の配列esの総乗: $\prod_{i}{e_i}$
  • isZero(e): ゼロ判定: $e = 0$

これらの関数の実装に必要な、上述の関数群については、次の有限体Qで説明します。

ちなみに、この記事のプログラミング上での体は、多項式の係数など、積が可換であることを前提とした利用をします。たとえば、積が非可換になる4元数体(quaternion)は対象外です。

有限体Q実装

field-utils.js[3]
export const Q = () => {
  const gcd = (a, b) => a === 0 ? b : gcd(b % a, a);
  const normalize = ([e0, e1]) => {
    const f = gcd(e1, e0);
    return [e0 / f, e1 / f];
  };

  const eql = (a, b) => a[0] * b[1] === a[1] * b[0];
  const zero = () => [0, 1];
  const one = () => [1, 1];
  const neg = e => [-e[0], e[1]];
  const inv = e => normalize([e[1], e[0]]);
  const add = (a, b) => normalize([a[0] * b[1] + b[0] * a[1], a[1] * b[1]]);
  const mul = (a, b) => normalize([a[0] * b[0], a[1] * b[1]]);

  const times = (e, k) => normalize([e[0] * k, e[1]]);
  const pow = (e, k) => normalize([e[0] ** k, e[1] ** k]);
  const toStr = e => e[1] === 1 ? `${e[0]}` : `${e[0]}/${e[1]}`;
  const toNum = e => Math.round(e[0] / e[1]);
  const fromNum = n => [n, 1];

  return Field({
    eql, zero, one, neg, inv, add, mul, times, pow, toStr, toNum, fromNum});
};

体の実装例としての有理数体Qです。

  • 体の要素: ２要素の整数配列 [分子, 分母]

実装した体の演算

  • eql(a, b): 同値関係 $a = b$
  • zero(): $0$
  • one(): $1$
  • `neg(e): 負数 $-e$
  • `inv(e): 逆数 $1/e$
  • add(a, b): 加算 $a + b$
  • mul(a, b): 乗算 $a * b$
  • `times(e, k): 整数倍 $ke$
  • pow(e, k): べき乗 $e^k$
  • toStr(e): $e$の文字列表現
  • toNum(e): $e$の数値表現
  • fromNum(n): 数値表現に対応するQの要素(2要素配列)

引数のa、b、eは体の要素の値で、kは整数値です。

timesとpowのkは、Qの整数[e0, 1]ではない点に注意です。このkは、addやmulで演算の要素の数です。たとえば、$e + e + e$を$3e$、$e * e * e$を$a^3$と表現したものです。

toStr、toNum、fromNumはJavaScriptの値との相互変換機能になります。

これにFieldで追加されるisZero(e)、sub(a, b)、div(a, b)、sum(es)、prod(es)を持ちます。

体のユーティリティFieldUtils(K)

field-utils.js[4]
export const FieldUtils = K => {
  const {isZero, zero, one, neg, inv, sub, mul, div} = K;
  // linear equations solver on the field f
  const pivoting = (r, u, i) => {
    for (let v = u + 1; v < r.length; v++) {
      if (!isZero(r[v][i])) {
        [r[u], r[v]] = [r[v], r[u]];
        break;
      }
    }
  };

  const upperTriangle = r => {
    const n = r[0].length - 1, m = r.length;
    for (let i = 0, u = 0; i < n && u < m; i++) {
      if (isZero(r[u][i])) pivoting(r, u, i);
      if (isZero(r[u][i])) continue;
      for (let v = u + 1; v < m; v++) {
        const c = div(r[v][i], r[u][i]);
        for (let j = i; j < n + 1; j++) {
          r[v][j] = sub(r[v][j], mul(c, r[u][j]));
        }
      }
      u++;
    }
  };

  const rank = lineqs => {
    const r = lineqs.map(lineq => lineq.slice());
    upperTriangle(r);
    return r.filter(eq => !eq.slice(0, -1).every(isZero)).length;
  };

  const solve = lineqs => {
    const r = lineqs.map(lineq => lineq.slice());
    upperTriangle(r);
    const n = r[0].length - 1, m = r.length;
    for (let i = n - 1, u = m - 1; i >= 0 && u >= 0; i--) {
      while (u > 0 && isZero(r[u][i])) u--;
      if (isZero(r[u][i])) continue;
      r[u][n] = div(r[u][n], r[u][i]);
      r[u][i] = one();
      for (let v = u - 1; v >= 0; v--) {
        r[v][n] = sub(r[v][n], mul(r[v][i], r[u][n]));
        r[v][i] = zero();
      }
    }
    return r.slice(0, n).map(l => neg(l[n]));
  };
  return {rank, solve};
};

体Kをパラメータに取るFieldUtils(K)は、指定した体Kの演算を用いて実装する汎用計算を提供します。
ここでは、後の応用で使用する連立方程式用の機能を提供します。

  • rank(lineqs): 体Kを係数とする連立方程式の係数行列のランク
    • lineqs: 1次方程式E[0]x0 + E[1]x1 + ... + E[N]xn + E[N+1] = 0を表す配列Eを要素とする連立方程式を表す配列
    • 戻り値: 係数行列のランク
  • solve(lineqs): 体Kを係数とする連立方程式を解くガウス除去ソルバ
    • lineqs: 1次方程式E[0]x0 + E[1]x1 + ... + E[N]xn + E[N+1] = 0を表す配列Eを要素とする連立方程式を表す配列
    • 戻り値: 解[x0, x1, ..., xn]の配列

rank(lineqs)での、連立方程式の係数行列とは、方程式a*x+b*y = fとc*x+d*y = gの、変数xとyの連立方程式では、その係数でつくる[[a, b], [c, d]]に相当する行列のことです。


ノート: 係数行列のランク、連立方程式ソルバー実装

係数行列のランクは、変数の実質の次元のことです。3つの方程式がx、y、zと３変数で構成されていても、係数行列のランクが2であれば、2変数で残りの1変数が表現できる係数構成になっていることを意味し、連立方程式として単一の3変数の解を求めることができなくなります。つまり、ランクと変数の数を比較することで、連立方程式が単一の解を持つかどうかの判定ができます。つまり、解ける連立方程式は、変数の数とランクが一致する、ということです。

solve(lineqs)は、連立方程式を満たす解を計算し、配列で出力します。
実装は、ガウス消去法です。
ただし、係数の数と方程式の数が一致しない、冗長な連立方程式でも解けるようにするために、行と列のためのループの変数を分離したため、対角行列特化の実装と比べて、若干複雑になっています。

ランクの計算は、ガウス消去法の前半の結果として得られる上三角行列で明らかになります。
この上三角行列の中の、すべて0でない行の数がランクとなります。



連立方程式ソルバ利用例
solve-example.js
import {Q, FieldUtils} from "./field-utils.js";

{
  const K = Q(), futils = FieldUtils(K);
  const lineqs = [
    [1, 2, -5], // x+2y-5=0
    [2, 3, -8], // 2x+3y-8=0
  ].map(eq => eq.map(K.fromNum));

  console.log("rank =", futils.rank(lineqs));
  console.log("[x, y] =", futils.solve(lineqs).map(K.toStr));
}

実行結果

rank = 2
[x, y] = [ '1', '2' ]


素数要素数の有限体PF(p): pf.js


コード全体
pf.js
import {range, uniq, Field} from "./field-utils.js";

// number utils
const egcd = (a, b) => {
  if (a === 0) return [b, 0, 1];    // b = 0*0 + b*1
  const r = b % a, q = (b - r) / a; // r = b - a*q
  const [g, s, t] = egcd(r, a);     // g = r*s + a*t;
  return [g, t - q * s, s];         // g = (b-a*q)*s + a*t = a*(t-q*s) + b*t
};
const primes = max => {
  const ps = [];
  let ns = range(max).slice(2);
  while (ns.length > 0) {
    const p = ns[0];
    ps.push(p);
    ns = ns.filter(n => n % p);
  }
  return ps;
};
const primeFactors = n => {
  const ps = [];
  for (const p of primes(n)) {
    if (n % p === 0) ps.push(p);
  }
  return ps;
};


// GF(p)
export const PF = (p, pe) => {
  const p1 = p - 1;
  const modp = n => (p + n % p) % p;
  const modp1 = n => (p1 + n % p1) % p1;

  const eql = (a, b) => modp(a) === modp(b);
  const zero = () => 0;
  const one = () => 1;
  const add = (a, b) => modp(a + b);
  const mul = (a, b) => modp(a * b);
  const neg = e => modp(-e);
  const inv = e => modp(egcd(e, p)[1]);

  const times = (e, k) => modp(e * k);
  const pow = (e, k) => {
    k = modp1(k);
    return k === 0 ? one() : k === 1 ? e : mul(e, pow(e, k - 1));
  };
  const toStr = e => `${e}`;
  const toNum = e => e;
  const fromNum = n => n;

  // find primitive element for GF(p)
  const isPrimitiveElement = (e, pfs = primeFactors(p1)) => {
    for (const pf of pfs) if (pow(e, p1 / pf) === 1) return false;
    return true;
  };
  const primitiveElement = p => {
    const pfs = primeFactors(p1);
    for (let e = 2; e < p1; e++) if (isPrimitiveElement(e, pfs)) return e;
    throw Error("never reached");
  };
  const a = pe ? pe : p === 2 ? 1 : p === 3 ? 2 : primitiveElement(p);
  const alpha = (k = 1) => pow(a, k);

  const size = () => p;
  const isPrimitive = () => isPrimitiveElement(a);
  const elems = () => range(p);

  return Field({
    p, eql, zero, alpha, one, add, mul, neg, inv, times, pow,
    toStr, toNum, fromNum, size, elems});
};


整数ユーティリティ関数

pf.js[1]
const egcd = (a, b) => {
  if (a === 0) return [b, 0, 1];    // b = 0*0 + b*1
  const r = b % a, q = (b - r) / a; // r = b - a*q
  const [g, s, t] = egcd(r, a);     // g = r*s + a*t;
  return [g, t - q * s, s];         // g = (b-a*q)*s + a*t = a*(t-q*s) + b*t
};
const primes = max => {
  const ps = [];
  let ns = range(max).slice(2);
  while (ns.length > 0) {
    const p = ns[0];
    ps.push(p);
    ns = ns.filter(n => n % p);
  }
  return ps;
};
const primeFactors = n => {
  const ps = [];
  for (const p of primes(n)) {
    if (n % p === 0) ps.push(p);
  }
  return ps;
};
  • egcd(a, b): a*s + b*t = g となる[g, s, t]の配列を返す拡張ユークリッド互除法
  • primes(max): max未満の素数の配列を返す
  • primeFactors(n): nの素因数の配列を返す

egcdはPF(p)での逆数inv(a)実装で使用します。
また、primeFactorsは、PF(p)での原始元を見つけるために使用します。

素数要素数の有限体PF(q)実装

pf.js[2]
export const PF = (p, pe) => {
  const p1 = p - 1;
  const modp = n => (p + n % p) % p;
  const modp1 = n => (p1 + n % p1) % p1;

  const eql = (a, b) => modp(a) === modp(b);
  const zero = () => 0;
  const one = () => 1;
  const add = (a, b) => modp(a + b);
  const mul = (a, b) => modp(a * b);
  const neg = e => modp(-e);
  const inv = e => modp(egcd(e, p)[1]);

  const times = (e, k) => modp(e * k);
  const pow = (e, k) => {
    k = modp1(k);
    return k === 0 ? one() : k === 1 ? e : mul(e, pow(e, k - 1));
  };
  const toStr = e => `${e}`;
  const toNum = e => e;
  const fromNum = n => n;

  // find primitive element for GF(p)
  const isPrimitiveElement = (e, pfs = primeFactors(p1)) => {
    for (const pf of pfs) if (pow(e, p1 / pf) === 1) return false;
    return true;
  };
  const primitiveElement = p => {
    const pfs = primeFactors(p1);
    for (let e = 2; e < p1; e++) if (isPrimitiveElement(e, pfs)) return e;
    throw Error("never reached");
  };
  const a = pe ? pe : p === 2 ? 1 : p === 3 ? 2 : primitiveElement(p);
  const alpha = (k = 1) => pow(a, k);

  const size = () => p;
  const isPrimitive = () => isPrimitiveElement(a);
  const elems = () => range(p);

  return Field({
    p, eql, zero, alpha, one, add, mul, neg, inv, times, pow,
    toStr, toNum, fromNum, size, elems});
};

PF(p)はパラメータとして素数値pを受け取り、モジュロ演算を用いた、素数要素数の有限体GF(p)を作ります。

  • PF(p)の要素: 非負整数値(JavaScriptのnumber)

PF(p)の演算実装(Qと同じ関数群を実装)

  • eql(a, b): a % q === b % q
  • zero(): 0
  • one(): 1
  • add(a, b): a + b % q
  • mul(a, b): a * b % q
  • neg(e): q - e
  • inv(e): [1, s, t] = egcd(e, q)のs % qがeの逆数
  • times(e, k): l * e % q
  • pow(e, k): 1にeをk回上記のmulでかける
  • toStr(e), toNum(e)、fromNum(n): 略

また、有限体として以下の関数を備えます:

  • alpha(k=1): PF(p)の原始元(のk乗)を返す
  • isPrimitive(): (パラメータpeで設定可能な)aが原始元かどうか
  • size(): 要素数を返す
  • elems(): 全要素の配列を返す

これらに、前述のField()によるisZero(a), sub(a, b)、div(a, b)、sum(es)、prod(es)が追加されたオブジェクトを返します。


ノート: 素数要素数の有限体PF(p)の原始元

PF(q)にも、べき乗していくことで、すべての数になる原始元(Primitive Element)が(複数)存在します。

ただし、有限体GF($p^n$)が原始元aを基準にして、その多項式を要素とするのと違い、要素は整数値を用いるため、数値の中から原始元aになるものを発見しなくてはいけません。

この発見方法は、2からp-1までの数の中から、$p-1$乗までべき乗することで初めて1になる数であるかをチェックすることです。

GF($p^n$)と同様、数eが原始元でなければ、eの指数表現$a^k$としたときの指数kはp-1と互いに素ではない数になります。そして、kとp-1の公約数をdとすると、どの公約数でも、$a^k$は(p-1)/d乗したときにも1になります。 k(p-1)/d mod p-1が0になり、$a^0=1$だからです。

ただし原始元が決まらなければこの公約数dもわかりません。
そこでdになりうる候補としてp-1の素因数を総当りし、どの(p-1)/d乗でも1にならなければ、それは原始元である、として実装することができます。

暗号技術のDH鍵交換などでは、PF(p)の原始元のことを、乗算の巡回群関係だけに注目して、(巡回群としての)生成元(generator)とよんでいます。同時にどの数も$g^0=1$になる、べき乗の指数p-1のことも(巡回群としての)位数(order)と呼びます。

また、暗号化技術では、原始元(生成元)は、解析を難しくするため、大きめ値のものを指定して用います。暗号技術の仕組みの理解するために、小さい値の原始元を用いて手計算で行っても、同様の結果が得られます。


体Kを係数とする多項式Polynomial(K): polynomial.js


コード全体
polynomial.js
import {range, uniq} from "./field-utils.js";

// Polynomial[F]
export const Polynomial = K => {
  const zeros = n => Array(n).fill(K.zero());

  const eql = (a, b) => {
    if (a.length <  b.length) [a, b] = [b, a];
    return a.every((c, k) => k < b.length ? K.eql(c, b[k]) : K.isZero(c));
  };
  const zero = () => [K.zero()];
  const one = () => [K.one()];
  const add = (a, b) => {
    if (a.length < b.length) [a, b] = [b, a];
    return a.map((c, k) => k < b.length ? K.add(c, b[k]) : c);
  };
  const neg = e => e.map(c => K.neg(c));
  const sub = (a, b) => add(a, neg(b));
  const sum = es => es.reduce((r, e) => add(r, e), zero());
  const scale = (e, c) => e.map(ci => K.mul(ci, c));
  const carry = (e, k) => [...zeros(k), ...e];
  const mul = (a, b) => sum(b.map((bk, k) => carry(scale(a, bk), k)));

  const order = e => {
    for (let i = e.length - 1; i > 0; i--) if (!K.isZero(e[i])) return i;
    return 0;
  };
  const mod = (a, b) => {
    const an = order(a), bn = order(b);
    console.assert(bn > 0, "poly.mod");
    if (an < bn) return a.slice(0, bn);
    const f = K.div(a[an], b[bn]);
    return mod(sub(a, carry(scale(b, f), an - bn)).slice(0, an), b);
  };
  const prod = es => es.reduce((r, e) => mul(r, e), one());
  const times = (e, k) => e.map(c => K.times(c, k));
  const diff = e => e.map((c, k) => K.times(c, k)).slice(1);

  const coef = (e, k) => e[k] || K.zero();
  const monomial = (c, k) => carry(scale(one(), c), k);
  const apply = (e, v) => K.sum(e.map((c, k) => K.mul(c, K.pow(v, k))));

  const toStr = (e, name = "x") => {
    const s1 = e.map((c, k) => {
      if (k === 0) {
        const coef = K.toStr(c);
        return coef.includes("+") ? `(${coef})` : coef;
      }
      if (K.isZero(c)) return "";
      const coef = K.eql(c, K.one()) ? "" : K.toStr(c);
      const cstr = coef.includes("+") ? `(${coef})` : coef;
      const pow = k === 1 ? "" : `^{${k}}`;
      return `${cstr}${name}${pow}`;
    }).filter(e => e);
    const s2 = s1.length > 1 && s1[0] === K.toStr(K.zero()) ? s1.slice(1) : s1;
    return s2.reverse().join("+");
  };
  const toArray = (e, len) =>
        e.length < len ? [...e, ...zeros(len - e.length)] : e.slice(0, len);
  const fromArray = cs => cs;
  const R = {
    K, eql, zero, one, add, scale, carry, neg, sub, sum, mul, mod, prod, times,
    order, diff, coef, monomial, apply, toStr, toArray, fromArray,
  };

  if (K.size && K.elems) {
    const toNum = e => e.reduceRight((r, c) => r * K.size() + K.toNum(c), 0);
    const fromNum = n => {
      const e = [];
      while (n > 0) {
        const r = n % K.size();
        e.push(K.fromNum(r));
        n = (n - r) / K.size();
      }
      return e;
    };
    const elems = k => {
      if (k < 0) return [[]];
      const base = elems(k - 1);
      return K.elems().flatMap(c => base.map(b => [...b, c]));
    };
    return {...R, toNum, fromNum, elems};
  } else return R;
};

// Polynomial utilities
export const PolynomialUtils = poly => {
  const {K, eql, one, add, sub, scale, carry, mul, coef, elems} = poly;
  // listing irreducible polynomials
  const kPolynoms = k => elems(k - 1).map(e => add(e, carry(one(), k)));
  const reducibles = n => {
    const rs = [];
    for (let i = 1, nh = n >> 1; i <= nh; i++) {
      const l = kPolynoms(i), r = kPolynoms(n - i);
      for (const lp of l) for (const rp of r) rs.push(mul(lp, rp));
    }
    return uniq(rs, eql);
  };
  const irreducibles = n => {
    const reds = reducibles(n);
    return kPolynoms(n).filter(f => !reds.find(r => eql(f, r)));
  };

  const findLinearRecurrence = s => {
    // Berlekamp-Massey algorithm
    let cx = one(), cl = 1, bx = one(), bl = 1, b = K.one(), m = 0;
    for (let i = 0; i < s.length; i++) {
      const d = K.sum(range(cl).map(k => K.mul(coef(cx, k), s[i - k])));
      m++;
      if (K.isZero(d)) continue;
      const tx = cx, tl = cl;
      cx = sub(cx, scale(carry(bx, m), K.div(d, b)));
      cl = Math.max(cl, bl + m);
      if (cl > tl) [bx, bl, b, m] = [tx, tl, d, 0];
    }
    return cx;
  };

  return {irreducibles, findLinearRecurrence};
};


多項式Polynomial(K)

多項式実装は、代数系としての演算以外にも、多項式固有の操作も提供するため、大きくなります。
いくつかに、分割しながら、説明していきます。

polynomial.js[1]
export const Polynomial = K => {
  const zeros = n => Array(n).fill(K.zero());

  const eql = (a, b) => {
    if (a.length <  b.length) [a, b] = [b, a];
    return a.every((c, k) => k < b.length ? K.eql(c, b[k]) : K.isZero(c));
  };
  const zero = () => [K.zero()];
  const one = () => [K.one()];
  const add = (a, b) => {
    if (a.length < b.length) [a, b] = [b, a];
    return a.map((c, k) => k < b.length ? K.add(c, b[k]) : c);
  };

多項式Polynomial(K)は、パラメータとして体Kを係数とする多項式です。

  • 要素としての多項式: 変数の指数を配列インデックスとする、係数値のJavaScript配列
    • 例 PF(5)係数の多項式 $x^{3}+2x$: [0,2,0,1]

多項式の配列要素数は任意長とし、eqlやaddでの実装のように、配列サイズを超えたインデックスの係数は(体Kでの)0とみなすように各演算を実装します。

polynomial.js[2]
  const neg = e => e.map(c => K.neg(c));
  const sub = (a, b) => add(a, neg(b));
  const sum = es => es.reduce((r, e) => add(r, e), zero());
  const scale = (e, c) => e.map(ci => K.mul(ci, c));
  const carry = (e, k) => [...zeros(k), ...e];
  const mul = (a, b) => sum(b.map((bk, k) => carry(scale(a, bk), k)));

  const order = e => {
    for (let i = e.length - 1; i > 0; i--) if (!K.isZero(e[i])) return i;
    return 0;
  };
  const mod = (a, b) => {
    const an = order(a), bn = order(b);
    console.assert(bn > 0, "poly.mod");
    if (an < bn) return a.slice(0, bn);
    const f = K.div(a[an], b[bn]);
    return mod(sub(a, carry(scale(b, f), an - bn)).slice(0, an), b);
  };
  const prod = es => es.reduce((r, e) => mul(r, e), one());
  const times = (e, k) => e.map(c => K.times(c, k));

次は乗算の実装ですが、まず多項式でのneg(e)、sub(a, b)、sum(es)、scale(e, c)、carry(e, k)、order(e)を実装し、それらを用いてmul(a, b)とmod(a, b)を実装します。

  • neg(e): 負多項式。e(x)の各係数を(Kでの)負数にした多項式
  • sub(a, b): 多項式の減算。a(x)からb(x)を係数ごとに引いた多項式
  • sum(es): 多項式の総和。各項ごとの係数での総和
  • scale(e, c): e(x)のc倍。e(x)の各係数へｃを(体K)で乗算した多項式。
  • carry(e, k): e(x)のk桁上げ。e(x)を$x^k$倍した多項式。e(x)配列の頭にk個の0を追加した配列を作る
  • order(e): e(x)の次数。e(x)配列の0でない最上位要素のインデックスが多項式の次数

多項式の乗算mul(a, b)の実装では、b(x)の各k次の項ごとに、$b_k * a(x) * x^k$し(carry(scale(a, bk), k))、それらの総和を取ります。

多項式のモジュロ算mod(a, b)は、b(x)の次数以下になるまで、a(x)の最上位項が0になるよう引き算で取り除くことを(再帰呼び出しで)繰り返すことで剰余を割り出します。
モジュロ算が配列としての桁を抑える実用効果として、結果の配列の要素数をbの要素数引く1に切り詰めています。

  • prod(es): esの多項式の総乗の実装
  • times(e, k): k個のe(x)の総和の実装

これらは応用時に使用するため、用意しました。

polynomial.js[3]
  const diff = e => e.map((c, k) => K.times(c, k)).slice(1);

  const coef = (e, k) => e[k] || K.zero();
  const monomial = (c, k) => carry(scale(one(), c), k);
  const apply = (e, v) => K.sum(e.map((c, k) => K.mul(c, K.pow(v, k))));

  const toStr = (e, name = "x") => {
    const s1 = e.map((c, k) => {
      if (k === 0) {
        const coef = K.toStr(c);
        return coef.includes("+") ? `(${coef})` : coef;
      }
      if (K.isZero(c)) return "";
      const coef = K.eql(c, K.one()) ? "" : K.toStr(c);
      const cstr = coef.includes("+") ? `(${coef})` : coef;
      const pow = k === 1 ? "" : `^{${k}}`;
      return `${cstr}${name}${pow}`;
    }).filter(e => e);
    const s2 = s1.length > 1 && s1[0] === K.toStr(K.zero()) ? s1.slice(1) : s1;
    return s2.reverse().join("+");
  };
  const toArray = (e, len) =>
        e.length < len ? [...e, ...zeros(len - e.length)] : e.slice(0, len);
  const fromArray = cs => cs;
  const R = {
    K, eql, zero, one, add, scale, carry, neg, sub, sum, mul, mod, prod, times,
    order, diff, coef, monomial, apply, toStr, toArray, fromArray,
  };

多項式固有の機能

  • diff(e): e(x)の微分
  • coef(e, k): e(x)のk次項の係数(Coefficient)
  • monimial(c, k): 単項式$cx^k$
  • apply(e, v): 体Kの値vをe(x)に代入した体Kの値e(x=v)

JavaScriptデータとの相互変換

  • toStr(e, name="x"): e(a)の文字列表現(LaTeX数式形式)
  • toArray(e, len): 係数配列化(len要素固定)
  • fromArray(cs): 係数配列から多項式値の復元

多項式文字列の変数名のデフォルト名は"x"にしています。

  • 例 $a^3+2a$: [0,2,0,1] => "a^{3}+2a"
polynomial.js[4]
  if (K.size && K.elems) {
    const toNum = e => e.reduceRight((r, c) => r * K.size() + K.toNum(c), 0);
    const fromNum = n => {
      const e = [];
      while (n > 0) {
        const r = n % K.size();
        e.push(K.fromNum(r));
        n = (n - r) / K.size();
      }
      return e;
    };
    const elems = k => {
      if (k < 0) return [[]];
      const base = elems(k - 1);
      return K.elems().flatMap(c => base.map(b => [...b, c]));
    };
    return {...R, toNum, fromNum, elems};
  } else return R;
};

Kが有限体の場合(K.size()とK.elems()が存在する)には、以下の関数群を追加します。

  • toNum(e): e(x)のJavaScript整数表現
  • fromNum(n): JavaScript整数から多項式を復元
  • elems(k): k次未満の多項式の列挙

係数が有限体の場合、係数の有限体Kは要素数n未満の整数値にマッピングでき、整数化した係数をn進数の各桁の数とすることで、有限体係数の多項式も単一の整数値にマッピングできます。

多項式ユーティリティPolynomialUtils(poly)

polynomial.js[1]
export const PolynomialUtils = poly => {
  const {K, eql, one, add, sub, scale, carry, mul, coef, elems} = poly;
  // listing irreducible polynomials
  const kPolynoms = k => elems(k - 1).map(e => add(e, carry(one(), k)));
  const reducibles = n => {
    const rs = [];
    for (let i = 1, nh = n >> 1; i <= nh; i++) {
      const l = kPolynoms(i), r = kPolynoms(n - i);
      for (const lp of l) for (const rp of r) rs.push(mul(lp, rp));
    }
    return uniq(rs, eql);
  };
  const irreducibles = n => {
    const reds = reducibles(n);
    return kPolynoms(n).filter(f => !reds.find(r => eql(f, r)));
  };

  const findLinearRecurrence = s => {
    // Berlekamp-Massey algorithm
    let cx = one(), cl = 1, bx = one(), bl = 1, b = K.one(), m = 0;
    for (let i = 0; i < s.length; i++) {
      const d = K.sum(range(cl).map(k => K.mul(coef(cx, k), s[i - k])));
      m++;
      if (K.isZero(d)) continue;
      const tx = cx, tl = cl;
      cx = sub(cx, scale(carry(bx, m), K.div(d, b)));
      cl = Math.max(cl, bl + m);
      if (cl > tl) [bx, bl, b, m] = [tx, tl, d, 0];
    }
    return cx;
  };

  return {irreducibles, findLinearRecurrence};
};

多項式ユーティリティPolynomialUtils(poly)は、パラメータpolyの多項式演算を利用した汎用機能を提供します。

  • irreducible(n): 有限体係数のn次既約多項式(irreducible polynomial)を全列挙した配列を返す関数
  • findLinearRecurrence(s): 数列s[k]から、任意階数の線形漸化式(linear recurrence equation)の係数を見つける関数

irreducibles(n)による有限体係数のn次既約多項式の全列挙は、全n次多項式から、n次の可約多項式(reducible polynomial)を取り除いたものです。n次の可約多項式は、n未満次の全多項式同士を掛け合わせることで導出します。

線形漸化式は、任意個の前項の線形和として構成する数列のことです。
たとえばフィボナッチ数の漸化式$s_n = s_{n-1}+s_{n-2}$は２つ前までの項の線形和であり、線形漸化式の一種となります。このフィボナッチ数の漸化式の係数列は、1, 1です。ちなみに、等比aの等比数列は$s_n = as_{n-1}$の線形漸化式であり、等差数列の等比はどのk項でも同じ値の$s_{k} - s_{k-1}$になるため、等比数列はすべて、$s_n = s_{n-1} + (s_{n-1}-s_{n-2}) = 2s_{n-1} - s_{n-2}$の線形漸化式になります。

ただし、findLinearRecurrence(s)では、数列sが線形漸化式によるものである場合、漸化式表現を$s_n + c_1s_{n-1} + c_2s_{n-2} + ... = 0$と表現したときの係数列[1, $c_1$, $c_2$, ...]を(多項式扱いで)返します。

たとえば、フィボナッチ数s = [1, 1, 2, 3, 5, 8]を渡した結果は、配列として[1, -1, -1]、つまり$s_n - s_{n-1} - s_{n-2} = 0$が得られます。
この結果は、sのうちの3つ目以降の項2,3,5,8すべてに対して成立する関係となっています。


参考: Berlekamp-Masseyアルゴリズム

このfindLinearRecurrenceの実装のアルゴリズムは、Berlekamp-Masseyアルゴリズムというものです。

このアルゴリズムは、本来は後述する応用のBCH符号でのエラー位置方程式を、有限体数列から算出するための手法でした。この有限体の演算を切り替えることで、そのまま有理数体Qなどの無限体でも適用可能なものとなり、漸化式を発見する手法として解釈されるものとなりました。

以下は、このアルゴリズムの内容です。

  • 注: wikipediaのコードでのLやその判定条件は直観的ではないので、より素直な多項式の長さを用いるよう、置き換えています

まず、変数cxは漸化式係数(の多項式c(x))であり、数列sをたどるループ内で更新していきます。変数clはこのcxの長さです。

ループ内のdは、s[i]に対し、漸化式cxを適用した結果です。

  • $d = \sum_{k=0}^{{cl}-1}{c(x)}_k s[i-k]$

dが0であれば、cxの漸化式がそのiの時点では成立していることを意味します。

もし、dが0でなければ、以下の式でcxを更新します。

  • $c(x) - \frac{d}{b}b(x)x^m$

変数bxは、cxを更新してclが伸びたときの、その直前のcxになります。blはbxの長さです。
bは、bxを更新したときのdです。

mは、cx更新時に使うbxの桁上げ値です。ループで増やし続け、bx更新時にリセットします。
よってi-mはbxを更新した時点のiであり、bとmとbxの関係は、常に以下になります。

  • $b = \sum_{k=0}^{{bl}-1}{b(x)}_k s[i-m-k]$

このため、cxを更新した時点のiで、この新しいcxを用い、再びdを計算すれば必ず0になります。

\begin{align}
\sum_{k=0}^{{cl}-1}{(c(x) - \frac{d}{b}b(x)x^m)}_{k} s[i-k] &= \sum_{k=0}^{{cl}-1}{c(x)}_k s[i-k] - \frac{d}{b}\sum_{k=0}^{{bl}-1}{b(x)}_k s[i-m-k] \\
&= d - \frac{d}{b}{b} \\
&= 0
\end{align}

よって、ループが数列sの最後まで到達したあとのcxによる漸化式は、そのs全体で成立するものとなっています。

  • $\sum_{k=0}^{{cl}-1}{c(x)}_k s[i-k] = 0$ (ただしcx.length <= i < s.length)


素数べき要素数の有限体GF(p^n): gf.js


コード全体
gf.js
import {range, uniq, transpose, Field, FieldUtils} from "./field-utils.js";
import {PF} from "./pf.js";
import {Polynomial} from "./polynomial.js";

// GF(p^n)
export const GF = (K, n, f, name="a") => {
  const pn = K.size() ** n, pn1 = pn - 1;
  const poly = Polynomial(K);
  const p2gf = pe => poly.toArray(pe, n);
  const modpn1 = k => (pn1 + k % pn1) % pn1;

  const {eql, add, times, neg, toNum} = poly;
  const zero = () => p2gf(poly.zero());
  const one = () => p2gf(poly.one());
  const a = p2gf(poly.carry(poly.one(), 1));
  const isZero = e => eql(zero(), e);

  const mul0 = (a, b) => poly.mod(poly.mul(a, b), f);
  const pow0 = (e, k) => k === 0 ? one() : mul0(e, pow0(e, k - 1));

  const powList = Object.freeze(
    range(pn1).map(k => Object.freeze(pow0(a, k))));
  const exponent = e => powList.findIndex(pe => eql(e, pe));
  const mul = (a, b) => isZero(a) || isZero(b) ? zero() :
        powList[modpn1(exponent(a) + exponent(b))];
  const pow = (e, k) => isZero(e) ? zero() : k == 1 ? e :
        powList[modpn1(exponent(e) * k)];
  const alpha = (k = 1) => pow(a, k);
  const inv = e => powList[modpn1(-exponent(e))];
  const pows = () => powList;

  const size = () => pn;
  const toStr = e => poly.toStr(e, name);
  const fromNum = num => p2gf(poly.fromNum(num));
  const elems = () => poly.elems(n - 1);
  const isPrimitive = () => uniq(powList, eql).length === powList.length;
  const fromSF = e => p2gf(poly.monomial(e, 0));
  const toSF = e => e[0];

  return Field({
    K, poly, n, f, eql, zero, one, alpha, isZero,
    add, times, neg, exponent, mul, pow, inv,
    pows, toStr, size, toNum, fromNum, elems, isPrimitive, fromSF, toSF,
  });
};

// GF Utilities: judge primitive polynomial, irreducible formula of a GF elem
export const GFUtils = gf => {
  const {poly, n, eql, pow, exponent, elems} = gf;
  const futils = FieldUtils(poly.K);
  const p = gf.K.size(), pn1 = gf.size() - 1;

  const cyclicOrder = e => {
    const k = exponent(e);
    for (let i = 1; i <= n; i++) {
      if (k * (p ** i) % pn1 === k) return i;
    }
    throw Error("never reached");
  };

  const findIrreducible = e => {
    const d = cyclicOrder(e);
    // gfeq: [1, e, e^2, ..., e^d] for c0 + c1e +...+ c(d-1)e^(d-1) + e^d = 0
    const gfeq = range(d + 1).map(k => poly.toArray(pow(e, k), n));
    // coefEqs: split gfeq with each coeddicients in e^0,e^1,...,e^d
    const coefEqs = transpose(gfeq);
    // modDim: [c0, c1, ...c(d-1)] satisfied above equation
    console.assert(futils.rank(coefEqs) === d, "findIrreducible");
    const cs = futils.solve(coefEqs);
    const modDim = poly.fromArray(cs);
    return poly.add(modDim, poly.carry(poly.one(), d));
  };

  const exponentGroups = () => {
    const used = new Set();
    const gs = [];
    for (const expo of range(pn1)) {
      if (used.has(expo)) continue;
      const g = uniq(range(n).map(k => (expo * (p ** k)) % pn1));
      g.forEach(expo => used.add(expo));
      gs.push(g);
    }
    return gs;
  };

  const generatorAutomorphism = () => elems().map(e => pow(e, p));

  return {findIrreducible, exponentGroups, generatorAutomorphism};
};


素数べき要素数の有限体GF(p, n, f)実装

gf.js[1]
export const GF = (K, n, f, name="a") => {
  const pn = K.size() ** n, pn1 = pn - 1;
  const poly = Polynomial(K);
  const p2gf = pe => poly.toArray(pe, n);
  const modpn1 = k => (pn1 + k % pn1) % pn1;

  const {eql, add, times, neg, toNum} = poly;
  const zero = () => p2gf(poly.zero());
  const one = () => p2gf(poly.one());
  const a = p2gf(poly.carry(poly.one(), 1));
  const isZero = e => eql(zero(), e);

  const mul0 = (a, b) => poly.mod(poly.mul(a, b), f);
  const pow0 = (e, k) => k === 0 ? one() : mul0(e, pow0(e, k - 1));

  const powList = Object.freeze(
    range(pn1).map(k => Object.freeze(pow0(a, k))));
  const exponent = e => powList.findIndex(pe => eql(e, pe));
  const mul = (a, b) => isZero(a) || isZero(b) ? zero() :
        powList[modpn1(exponent(a) + exponent(b))];
  const pow = (e, k) => isZero(e) ? zero() : k == 1 ? e :
        powList[modpn1(exponent(e) * k)];
  const alpha = (k = 1) => pow(a, k);
  const inv = e => powList[modpn1(-exponent(e))];
  const pows = () => powList;

  const size = () => pn;
  const toStr = e => poly.toStr(e, name);
  const fromNum = num => p2gf(poly.fromNum(num));
  const elems = () => poly.elems(n - 1);
  const isPrimitive = () => uniq(powList, eql).length === powList.length;
  const fromSF = e => p2gf(poly.monomial(e, 0));
  const toSF = e => e[0];

  return Field({
    K, poly, n, f, eql, zero, one, alpha, isZero,
    add, times, neg, exponent, mul, pow, inv,
    pows, toStr, size, toNum, fromNum, elems, isPrimitive, fromSF, toSF,
  });
};

GF($p^n$)の関数群を作る関数がGF(K, n, f)です。各パラメータは、

  • K: 有限体の要素の多項式の係数体
  • n: 2以上の整数値
  • f: 多項式aが原始元としての解となるn次既約多項式Polynomial(PF(p))
  • オプションname = "a": 有限体の値の多項式の変数名

となります。GF(PF(p),n,f)の内部では、Polynomial(PF(p))を使用します。

GF($p^n$)の要素は、Polynomial(PF(p))によるn-1次多項式となるJavaScript配列で表しています。
ここではGF($p^n$)の要素である確認がしやすいよう、p2gfにより、要素数をn個固定になるようにしています。

  • 例: GF($2^3$)の$a+1$: [1, 1, 0]

GF($p^n$)が提供する関数群は、有限体としてPF(p)の関数群と同じラインナップになります。

  • eql(a, b)、add(a, b)、neg(e)、times(e, k)、toStr(e)、toNum(e):Polynomial(PF(p))`の実装そのもの
  • zero()、one()、fromNum(n): Polynomialのzero()とone()、fromNum(n)をn要素配列化した実装
  • size(): $p^n$固定
  • elems(): Polynomialのelems(n - 1)`の全多項式値の配列n要素配列化した実装

乗算やべき乗の実装は、原始元の指数表を作り、指数の加算演算で行う実装にします。

  • mul0(a, b): Polynomialのmulとmodを使用した有限体要素用の乗算実装
  • pow0(e, k): mul0を使用したべき乗実装

このpow0を用い、指数値から多項式値をひける原始元の指数表powListを作ります。
そして、このpowListを用い、以下の関数群を実装します。

  • exponent(e): 0以外の要素の原始元の指数値を返す関数
  • mul(a, b): 指数値の和で実装した有限体の乗算
  • pow(e, k): 指数値の乗算で実装した有限体のべき乗
  • alpha(k=1): 原始元のべき乗$a^k$の多項式値
  • inv(e): 指数値の負数で実装した有限体の逆数
  • isPrimitime(): パラメータの既約多項式fが、原始元を解とする多項式かどうか

上記の関数群がaの指数表が成立することに基づいているため、isPrimitive()がfalseなら、このGF｀のインスタンスは機能しないものになります。
これらの関数群は、どれもmul0やpow0`を直接使う実装が可能です。ただし、GFの応用は、多項式aが原始元であることに依存しているものが普通であるため、有用ではありません。

PolynomialUtilsのirreducibles()で得た既約多項式が、原始元の多項式になるかどうかを調べるためにも使えます。

最後に、要素の多項式の係数である部分体(sub-field)との間の変換機能をもたせています。

  • fromSF(e): 部分体Kの要素eを定数項とするGFでの多項式要素に変換する関数
  • toSF(e): GFの多項式要素eの定数項を返す関数

指数べき要素数の有限体ユーティリティ

gf.js[2]
export const GFUtils = gf => {
  const {poly, n, eql, pow, exponent, elems} = gf;
  const futils = FieldUtils(poly.K);
  const p = gf.K.size(), pn1 = gf.size() - 1;

  const cyclicOrder = e => {
    const k = exponent(e);
    for (let i = 1; i <= n; i++) {
      if (k * (p ** i) % pn1 === k) return i;
    }
    throw Error("never reached");
  };

  const findIrreducible = e => {
    const d = cyclicOrder(e);
    // gfeq: [1, e, e^2, ..., e^d] for c0 + c1e +...+ c(d-1)e^(d-1) + e^d = 0
    const gfeq = range(d + 1).map(k => poly.toArray(pow(e, k), n));
    // coefEqs: split gfeq with each coeddicients in e^0,e^1,...,e^d
    const coefEqs = transpose(gfeq);
    // modDim: [c0, c1, ...c(d-1)] satisfied above equation
    console.assert(futils.rank(coefEqs) === d, "findIrreducible");
    const cs = futils.solve(coefEqs);
    const modDim = poly.fromArray(cs);
    return poly.add(modDim, poly.carry(poly.one(), d));
  };

  const exponentGroups = () => {
    const used = new Set();
    const gs = [];
    for (const expo of range(pn1)) {
      if (used.has(expo)) continue;
      const g = uniq(range(n).map(k => (expo * (p ** k)) % pn1));
      g.forEach(expo => used.add(expo));
      gs.push(g);
    }
    return gs;
  };

  const generatorAutomorphism = () => elems().map(e => pow(e, p));

  return {findIrreducible, exponentGroups, generatorAutomorphism};
};

GFUtils(gf)は、指数べき要素数の有限体GFを用い、以下の関数を実装したものです。

  • findIrreducible(e): 要素eを解に持つ既約多項式を返す関数
  • expoentGroups(): 同一既約多項式を満たすGFの元のグルーピング(原始元の指数の配列の配列)
  • generatorAutonorphism(): GF.elems()への(自己同型群の生成元となる)自己同型変換

findIrreducible(e)は、eを解に持つ既約多項式を見つける関数です。
内部で既約多項式の係数を求める(冗長な)連立方程式を解くため、前述のFieldUtilsを使用します。

exponentGroups()は、たとえばGF(PF(2), 2)であれば、[[0], [1, 2]]を返します。
この結果は、$a^0$で一つの1次既約多項式の解、$a^1$と$a^2$とで同じ2次既約多項式の解であることを示したものです。


findIrreduciblesとexponentGroupsの利用コード例
find-irreducible-example.js
import {PF} from "./pf.js";
import {GF, GFUtils} from "./gf.js";

const p = 2, n = 2, f0 = [1, 1, 1];
const gf = GF(PF(p), n, f0), gfutils = GFUtils(gf);

const eg = gfutils.exponentGroups()[1]; //=> [1, 2]
const f = gfutils.findIrreducible(gf.alpha(eg[0]));
console.log(gf.poly.toStr(f)); //=> "x^{2}+x+1"


GFの応用: GF(p^n)のマークダウン形式の演算表等を作る

以下のコードは、Qiitaに埋め込めるMathJax&Markdown table形式の有限体の加算と乗算の演算表を作る、gf.js実装の応用例です。

calc-table-example.js
import {PF} from "./pf.js";
import {GF} from "./gf.js";

const calcTable = (elems, op) => elems.map(e1 => elems.map(e2 => op(e1, e2)));

const mdTable = (heads, elements) => {
  const maxes = heads.map(
    (h, i) => Math.max(h.length, ...elements.map(l => l[i].length)));
  const top = `| ${heads.map((h, i) => h.padEnd(maxes[i])).join(" | ")} |`;
  const guide = `|-${heads.map((_, i) => "-".repeat(maxes[i])).join("-|-")}-|`;
  const lines = elements.map(
    l => `| ${l.map((e, i) => e.padEnd(maxes[i])).join(" | ")} |`);
  return [top, guide, ...lines].join("\n");
};

const mdCalcTable = (f, table, mark) => {
  const elems = f.elems().map(e => f.toStr(e));
  const len = Math.max(...elems.map(e => e.length));
  const strList = elems.map(e => `$${e.padEnd(len)}$`);
  const strTable = table.map(
    l => l.map(e => `$${f.toStr(e).padEnd(len)}$`));
  const heads = [mark, ...strList];
  const elements = strList.map((e, i) => [`**${e}**`, ...strTable[i]]);
  return mdTable(heads, elements);
};

const outputCalcTables = (gf) => {
  const elems = gf.elems();
  console.log(mdCalcTable(gf, calcTable(elems, gf.add), "+"));
  console.log();
  console.log(mdCalcTable(gf, calcTable(elems, gf.mul), "*"));
  console.log();
};

// example
{
  const p = 2, n = 3, f = [1, 1, 0, 1];
  outputCalcTables(GF(PF(p), n, f));
}

実行結果

| +               | $0        $ | $1        $ | $a        $ | $a+1      $ | $a^{2}    $ | $a^{2}+1  $ | $a^{2}+a  $ | $a^{2}+a+1$ |
|-----------------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|
| **$0        $** | $0        $ | $1        $ | $a        $ | $a+1      $ | $a^{2}    $ | $a^{2}+1  $ | $a^{2}+a  $ | $a^{2}+a+1$ |
| **$1        $** | $1        $ | $0        $ | $a+1      $ | $a        $ | $a^{2}+1  $ | $a^{2}    $ | $a^{2}+a+1$ | $a^{2}+a  $ |
| **$a        $** | $a        $ | $a+1      $ | $0        $ | $1        $ | $a^{2}+a  $ | $a^{2}+a+1$ | $a^{2}    $ | $a^{2}+1  $ |
| **$a+1      $** | $a+1      $ | $a        $ | $1        $ | $0        $ | $a^{2}+a+1$ | $a^{2}+a  $ | $a^{2}+1  $ | $a^{2}    $ |
| **$a^{2}    $** | $a^{2}    $ | $a^{2}+1  $ | $a^{2}+a  $ | $a^{2}+a+1$ | $0        $ | $1        $ | $a        $ | $a+1      $ |
| **$a^{2}+1  $** | $a^{2}+1  $ | $a^{2}    $ | $a^{2}+a+1$ | $a^{2}+a  $ | $1        $ | $0        $ | $a+1      $ | $a        $ |
| **$a^{2}+a  $** | $a^{2}+a  $ | $a^{2}+a+1$ | $a^{2}    $ | $a^{2}+1  $ | $a        $ | $a+1      $ | $0        $ | $1        $ |
| **$a^{2}+a+1$** | $a^{2}+a+1$ | $a^{2}+a  $ | $a^{2}+1  $ | $a^{2}    $ | $a+1      $ | $a        $ | $1        $ | $0        $ |

| *               | $0        $ | $1        $ | $a        $ | $a+1      $ | $a^{2}    $ | $a^{2}+1  $ | $a^{2}+a  $ | $a^{2}+a+1$ |
|-----------------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|-------------|
| **$0        $** | $0        $ | $0        $ | $0        $ | $0        $ | $0        $ | $0        $ | $0        $ | $0        $ |
| **$1        $** | $0        $ | $1        $ | $a        $ | $a+1      $ | $a^{2}    $ | $a^{2}+1  $ | $a^{2}+a  $ | $a^{2}+a+1$ |
| **$a        $** | $0        $ | $a        $ | $a^{2}    $ | $a^{2}+a  $ | $a+1      $ | $1        $ | $a^{2}+a+1$ | $a^{2}+1  $ |
| **$a+1      $** | $0        $ | $a+1      $ | $a^{2}+a  $ | $a^{2}+1  $ | $a^{2}+a+1$ | $a^{2}    $ | $1        $ | $a        $ |
| **$a^{2}    $** | $0        $ | $a^{2}    $ | $a+1      $ | $a^{2}+a+1$ | $a^{2}+a  $ | $a        $ | $a^{2}+1  $ | $1        $ |
| **$a^{2}+1  $** | $0        $ | $a^{2}+1  $ | $1        $ | $a^{2}    $ | $a        $ | $a^{2}+a+1$ | $a+1      $ | $a^{2}+a  $ |
| **$a^{2}+a  $** | $0        $ | $a^{2}+a  $ | $a^{2}+a+1$ | $1        $ | $a^{2}+1  $ | $a+1      $ | $a        $ | $a^{2}    $ |
| **$a^{2}+a+1$** | $0        $ | $a^{2}+a+1$ | $a^{2}+1  $ | $a        $ | $1        $ | $a^{2}+a  $ | $a^{2}    $ | $a+1      $ |

前半の表もこの実装で生成したものを使用しています。

2べき要素数の有限体に特化したプログラム実装

2要素有限体PF(2)の多項式値として(数値配列ではなく)ビット列を割り当てることで、多項式の演算(addやcarryなど)のために、CPUレベルで実装されるビット演算が適用できるようになります。

以下のgf2n.jsは、ビット列で実装するPF(2)特化のPolynomialを実装したPF2Polynomail、および、このPF2Polynomialを使用して実装したGF2n(p, f)のコードとなります。


コード: gf2n.js
gf2n.js
import {Field, range, uniq} from "./field-utils.js";
import {PF} from "./pf.js";

const msb32 = n => 31 - Math.clz32(n);
const popcnt32 = n => {
  n = (n & 0x55555555) + ((n >>> 1) & 0x55555555);
  n = (n & 0x33333333) + ((n >>> 2) & 0x33333333);
  n = (n & 0x0f0f0f0f) + ((n >>> 4) & 0x0f0f0f0f);
  n = (n & 0x00ff00ff) + ((n >>> 8) & 0x00ff00ff);
  return (n & 0x0000ffff) + ((n >>> 16) & 0x0000ffff);
};

//GF(2) coefficient Polynomial as bits
export const PF2Polynomial = () => {
  const K = PF(2);
  const eql = (a, b) => a === b;
  const zero = () => 0;
  const one = () => 1;
  const add = (a, b) => a ^ b;
  const neg = e => e;
  const sub = (a, b) => add(a, b);
  const sum = es => es.reduce((r, e) => add(r, e), 0);
  const scale = (e, c) => (c & 1) ? e : 0;
  const times = (e, k) => (k & 1) ? e : 0;
  const carry = (e, k) => e << k;
  const mul = (a, b) => {
    let r = 0;
    for (; b > 0; b >>>= 1, a <<= 1) {
      if (b & 1) r ^= a;
    }
    return r;
  };
  const order = e => Math.max(msb32(e), 0);
  const mod = (a, b) => {
    const mb = msb32(b);
    for (let i = msb32(a); i >= mb; i--) {
      if (a & (1 << i)) a ^= b << (i - mb);
    }
    return a;
  };
  const prod = es => es.reduce((r, e) => mul(r, e), one());
  const diff = e => (e & 0xaaaaaaaa) >>> 1;
  const coef = (e, k) => e & (1 << k);
  const monomial = (c, k) => carry(scale(one(), c), k);
  const apply = (e, v) => (v & 1) ? popcnt32(e) & 1 : 0;
  const toStr = (e, name = "x") => {
    const s1 = [...e.toString(2).split("")].reverse().map((v, i) => {
      if (i === 0) return v;
      if (v === "0") return "";
      const pow = i === 1 ? "" : `^{${i}}`;
      return `${name}${pow}`;
    }).filter(e => e);
    const s2 = s1.length > 1 && s1[0] === "0" ? s1.slice(1) : s1;
    return s2.reverse().join("+");
  };
  const toArray = (e, len) => {
    const r = Array(len).fill(0);
    for (let i = 0; i < len; i++) r[i] = (e >>> i) & 1;
    return r;
  };
  const fromArray = bs => {
    let e = 0;
    for (let i = 0; i < bs.length; i++) e |= (bs[i] & 1) << i;
    return e;
  };
  const toNum = e => e;
  const fromNum = e => e;
  const elems = k => range(2 ** (k + 1));
  return {
    K, eql, zero, one, add, scale, neg, sub, sum, times, carry, mul, mod, prod,
    order, diff, coef, monomial, apply,
    toStr, toArray, fromArray, toNum, fromNum, elems,
  };
};

// GF(2^n) as bits
export const GF2n = (n, f, name="a") => {
  const poly = PF2Polynomial();
  const {K, eql, zero, one, add, times, neg, toNum, fromNum} = poly;
  const pn = 2 ** n, pn1 = pn - 1;
  const modpn1 = n => (pn1 + n % pn1) % pn1;

  const isZero = e => e === 0;
  const mul0 = (a, b) => poly.mod(poly.mul(a, b), f);
  const pow0 = (e, k) => k === 0 ? one() : mul0(e, pow0(e, k - 1));

  const powList = Object.freeze(range(pn1).map(k => pow0(2, k)));
  const expoList = range(pn).map(k => k === 0 ? NaN : powList.indexOf(k));
  const exponent = e => expoList[e];
  const mul = (a, b) => isZero(a) || isZero(b) ? zero() :
        powList[modpn1(exponent(a) + exponent(b))];
  const pow = (e, k) => isZero(e) ? zero() : k === 1 ? e :
        powList[modpn1(exponent(e) * k)];
  const alpha = (k = 1) => pow(2, k);
  const inv = e => powList[modpn1(-exponent(e))];
  const pows = () => powList;

  const size = () => pn;
  const elems = () => poly.elems(n - 1);
  const isPrimitive = () => uniq(powList, eql).length === powList.length;
  const toStr = e => poly.toStr(e, name);
  const fromSF = e => e % 2;
  const toSF = e => e % 2;

  return Field({
    K, poly, n, f, eql, zero, one, alpha, add, times, neg,
    exponent, mul, pow, inv,
    pows, toStr, size, toNum, fromNum, elems, isPrimitive, fromSF, toSF,
  });
};


詳細は省きますが、PF2Polynomialが扱うPF(2)係数の多項式はビット列として扱う整数値です。同様にGF2n(n, f)の要素もビット列の整数値です。既約多項式fもビット列の整数値で表現します。

  • 要素 $a^2+a$ : 0b110
  • 既約多項式 $x^3 + x + 1$: 0b1011

多項式の加減算addとaubはどちらもXOR演算(a ^ b)で実装できます。
桁上げcarry(e, k)は、左シフト演算(e << k)になります。

多項式の係数は0か1のみであり、加算と減算の区別がないことなどから、乗算mulやモジュロ残mod等で、ビット列であることに特化させた実装ができます。

Polynomialと同じ関数群を実装しているため、PF2Polynomialは、polynomial.jsで実装したPolynomialUtilsで利用可能となっています。
同様に、GF2nも、gf.jsのGFUtilsで利用可能となっています。

4
1
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
1