天泣記

2021-02-28 (Sun)

#1 OCaml 4.12.0 と opam

OCaml 4.12.0 が出ていたので、 いつものように opam で install しようとしたが、 opam switch list-available としても、 4.11.2+flambda みたいな flambda を有効にして install する選択肢が出てこない

調べてみるとリリースノートに説明があった

OCaml 4.12.0

flambda (だけ) を有効にするには以下のようにすればいいようだ

opam switch create 4.12.0+flambda --package=ocaml-variants.4.12.0+options,ocaml-option-flambda

まぁ、組み合わせが自由になっていいことだ、というところかな

2021-01-18 (Mon)

#1 PEPM 2021 で発表した

Akira Tanaka. Coq to C Translation with Partial Evaluation. Proceedings of the 2021 ACM SIGPLAN Workshop on Partial Evaluation and Program Manipulation (PEPM '21), 2021-01-18, Virtual. (paper) (slides) (DOI)

2020-12-30 (Wed)

#1 線形合同法は乱数としてどう悪いのか

線形合同法 (Linear Congruential Generator, LCG) は擬似乱数のアルゴリズムだが、 下位ビットがあまり乱数になってないとか、 良くないので使わない方がいいといわれている

なぜ良くないのか、ということを真面目に考えたことがなかったのだが、 ふと考えてみた

線形合同法は以下によって生成される数列である

x{0} = seed
x{n+1} = (a * x{n} + c) % m

wikipedia の線形合同法の英語版 によれば、 m は modulus, a は multiplier, c は increment という名前である

さて、x{n+1} とか x{n} と書くのは長いので、 y = x{n+1}, x = x{n} とおくことにする

y = (a * x + c) % m

ここで、よくある実装では m は 2^32 など、2のべき乗が選ばれる (wikipedia の線形合同法の英語版 にはそういう例がたくさん載っている)

y = (a * x + c) % (2^n)

この場合、以下が成り立つ

y % 2 = (a * (x % 2) + c) % (2^n)

つまり、y の最下位ビット y % 2 は x の直前の状態 x の最下位ビット x % 2 だけに依存して決まる

ということは、状態が2種類しかないので、 最下位ビットだけを取り出した数列は (定常状態としては) 0 -> 1 -> 0 -> 1 -> ... と 0 と 1 が交互に現れるか、 0 -> 0 -> ... と 0 がずっと続くか、 1 -> 1 -> ... と 1 がずっと続くか、 という 3種類に限られる

つまり、最下位ビットが、最大周期2 の繰返しになるということで、これはよろしくない

2^n が悪いなら、3^n にすればいいのか、というと、

y = (a * x + c) % (3^n)
y % 3 = (a * (x % 3) + c) % (3^n)

となるので、3で割った余り、つまり、3進数における最下位桁が、直前の状態の最下位桁だけに依存して決まる、 となる

そうすると、状態は3種類になるので、最大周期3 の繰返しになるので、これもよろしくない

では、(2^n)*(3^n) とか素因数をふたつ以上混ぜればいいのかというと、

y = (a * x + c) % (2^n)(3^n)
y % 2 = (a * (x % 2) + c) % (2^n)(3^n)
y % 3 = (a * (x % 3) + c) % (2^n)(3^n)

2進数として解釈したときの最下位桁が最大周期2の繰返しになり、かつ、 3進数として解釈したときの最下位桁が最大周期3の繰返しになる

つまりぜんぜん解決になっていない

一般に、y = (a * x + c) % m のとき、m を割りきれる d があると、 y % d = (a * (x % d) + c) % d が成り立つ

これは、Coq/SSReflect では以下のように証明できる (d %| m は d が m を割り切れる、という意味で、あと、%% は剰余で C の % に対応する)

From mathcomp Require Import all_ssreflect.

Goal forall a c m d x y, d %| m ->
  y = (a * x + c) %% m ->
  y %% d = (a * (x %% d) + c) %% d.
Proof.
  move=> a c m d x y Hdvdn ->.
  rewrite modn_dvdm; last by [].
  apply /eqP.
  by rewrite eqn_modDr modnMmr.
Qed.

ということで、modulus が 2^32 なら、2進数の最下位桁は周期が最大2になるし、 下位2桁は周期が最大4になるし... 下位i桁は周期が最大 2^i になるし... となる

ということで、線形合同法の modulus は、 modulus を割り切る小さな数がないもの、つまり素数を選ぶことが望ましい

Wikipedia の線形合同法の日本語版の項には 「Stephen K. Park と Keith W. Miller が、彼らのサーベイ中で「最低基準」として示したもの」として、

x{i+1} = (48271 * x{i}) % (2**31 - 1)

が示されており、この 2**31 - 1 = 2147483647 は素数である

wikipedia の線形合同法の英語版 には、 いろいろな線形合同法の実装で使われているパラメータがのっているが、 modulus はだいたいが 2のべき乗で、それ以外は 2**31 - 1 = 2147483647 が 3個、 2**3 * 7**5 = 134456 が 1個 載っている

ところで、線形合同法における乱数の内部状態は直前に生成した乱数そのものなわけであるが、 その状態は、32ビット整数とか、Cの整数型として普通に利用可能なものに保存されるだろう

また、乱数は、modulus 未満の整数である

とすると、(modulus が 2**n でない場合) 状態として表現可能な数の中には、乱数列に含まれない値が存在する。 そういう値を種に設定したりすると何が起こるのだろうか、と考えると、 まぁ、m は m を割りきるから、y = (a * (x % m) + c) % m であり、modulus 以上の値を指定することは、 mod m の値を指定することに対応するのだろう

あと気になるのは、modulus 未満の値を種にしても、周期が modulus でない場合はありうるという点である

ということで、

x{i+1} = (48271 * x{i}) % (2**31 - 1)

について、状態のなかにどういう列が入っているか (力ずくで) 調べてみた

% cat tst.c
#include <stdio.h>
#include <stdint.h>

uint32_t bitary[1 << (32-5)];
#define SETBIT(i) (bitary[(i) >> 5] |= 1UL << ((i) & 0x1f))
#define TESTBIT(i) ((bitary[(i) >> 5] & (1UL << ((i) & 0x1f))) != 0)

#define A 48271UL
#define C 0
#define MODULUS ((1UL << 31) - 1)

int main(int argc, char *argv)
{
  uint64_t i;
  for (i = 0; i < MODULUS; i++) {
    if (!TESTBIT(i)) {
      uint64_t j = i;
      uint64_t n = 0;
      printf("start=%lu ", i);
      fflush(stdout);
      while (!TESTBIT(j)) {
        SETBIT(j);
        n++;
        j = (j * A + C) % MODULUS;
      }
      printf("length=%lu end=%lu\n", n, j);
      fflush(stdout);
    }
  }
}
% gcc -O3 tst.c
% ./a.out
start=0 length=1 end=0
start=1 length=2147483646 end=1

つまり、状態が 0 というのが不動点で、これはずっと 0 で変わらないので乱数にならず、 それ以外のすべての状態がひとつの乱数列をなしていることがわかる (まぁ、c = 0 なら 0 が不動点になるのは仕方ない)

また、Wikipedia に random0 として載っている

#define A 8121UL
#define C 28411UL
#define MODULUS 134456UL

というのを試すと、

start=0 length=134456 end=0

となった これは 0..134455 の範囲がひとつの周期的な乱数列になっている

もちろん、134456 = 2**3 * 7**5 なので、 この乱数は 2進数や7進数での下位桁が短い周期になっていてよろしくはない

では、modulus が素数であり、かつ、c が 0 でない例があるかと思って wikipedia の線形合同法の英語版 を 探すと、RtlUniform from Native API というのがあげられていた

#define A 2147483629UL
#define C 2147483587UL
#define MODULUS ((1UL << 31) - 1)

これを実行すると、以下のようになった

start=0 length=715827882 end=0
start=1 length=715827882 end=1
start=5 length=715827882 end=5
start=1243280003 length=1 end=1243280003

つまり、0..(MODULUS-1) という中に、4つの周期数列が含まれていて、 3個は周期 715827882 で、1個は周期 1 というわけである

まぁ、好んで使うものではない気がする

では、c = 0 だと種として 0 を選ぶと乱数でなくなってしまうので c を 0 以外に選び、 下位桁の周期性を避けるために m を素数として選び、 さらに周期が m になる、というのはあり得るだろうか、 と考えたが、これはうまくいかないようだ

wikipediaには、 c != 0 の場合で周期を m になる場合について、 m が square-free integer の場合は a = 1 にならざるを得ない、と書いてある。

square-free integer というのは素因数分解したときに、同じ素数は高々1つしか現れないというやつで、 素数は square-free integer である。

ということで、a = 1 にならざるを得ないということで、y = (x + c) % m になるわけだが、 これは乱数としてはまったくよろしくない

2020-11-24 (Tue)

#1 学術出版の難しさ

ふと興味があって、以下の本をざっと眺めた

学術出版の技術変遷論考 : 活版からDTPまで, 中西 秀彦, 印刷学会出版部, 2011

これを読むと、学術出版の難しさは、えんえんと続く校正の多さ、 多種の文字 (西夏文字とかもう使われていない文字を含む)、 数式あたりにあるようだ

校正への対応は写植よりも活版印刷のほうが楽、というのは興味深かった

2020-10-30 (Fri)

#1 italic と roman の使い分け

論文に式を書いていて、italic と roman の使い分けが 気になったのでちょっと調べた。

SIパンフレットには、「量の記号は斜体で書く」「単位記号は直立体で書く」とある。 JIS Z 8000-1:2014「量および単位 - 第1部: 一般」にも「量を表す記号は, 斜体文字で記載する」という記述がある。

しかし、ソフトウェアの話だと、 ベンチマークで測定する時間あたりは SI の記述に従えばいいが、 適用しがたいものも多い。 関数 eval はプログラム p と環境 e を受け取って値 v を返します、 というときに、プログラム、環境、値は量なのか、またその単位は何か、 などといわれても困ってしまう。

JIS Z 8201-1981「数学記号」には 「数値が一般的に定められている定数の記号は, 原則として立体とする」 「演算の記号は, 原則として立体とする」 「変数の記号は, 斜体とする」 とある。 ここで、定数の例として π, 演算の記号の例として log とかがあげられている。

ソフトウェアを考えると、 定数と変数の区別はつけやすい気がする。 ただ、高階関数を扱うと、演算 (関数) が変数に束縛されるので 演算と変数を区別しろといわれるとそれは困る。

Wikipedia:イタリック体には、 「変数(内容の変化する関数、作用素、物理定数などを含む)、すなわち内容が変化するものを表す記号は原則的にイタリック体で表記される」 「これに対して、2文字以上の標準的な関数名(log, sin, exp など)、演算記号(lim や d など)、数学定数(円周率 π, 虚数単位 i など)、物理単位など、内容が変化しない記号は立体で表記することが国際標準化機構 (ISO)、日本工業規格 (JIS)、日本物理学会などによって定められている」 とあるが、いろいろと例外も書いてある。

まぁ、定数と変数で区別するのが妥当かなぁ。

2020-09-30 (Wed)

#1 coq pull request

GitHub coq/coq/pull/13111: Small document fixes.



田中哲