bincombinations から n-combinations まで

R の e1071 ライブラリに bincombinations がある。

bincombinations(3) で，2^3 × 3 行列の結果が返される。

> bincombinations(3)
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    1
[3,]    0    1    0
[4,]    0    1    1
[5,]    1    0    0
[6,]    1    0    1
[7,]    1    1    0
[8,]    1    1    1

要するに，n 桁の 2 進数の各桁が 2^n 行で表される。ただ，これだと結果が行列で返されるので，不都合な場合もある。また，一般に n 進数の結果が必要になることもある。bincombinations を一般化するのは，このページの下の方に再掲載しておく。

まず，与えられた n 進数表記の次の数の表記を与える関数を書いてみる。

next.Ncombination = function(x, base=2) {
  n = length(x)
  carry = 0
  x[n] = x[n] + 1
  for (i in n:1) {
    x[i] = x[i] + carry
    carry = 0
    if (x[i] >= base) {
      x[i] = 0
      carry = 1
    }
  }
  return(x)
}

この関数を

n = 4
base = 2
x = integer(n)
for (i in seq.int(base ^ n)) {
  print(x)
  x = next.Ncombination(x, base)
}

のように指定すると全ての組み合わせが得られる。bincombinations(4) と同じであるが，それぞれの結果はベクトルである。

[1] 0 0 0 0
[1] 0 0 0 1
[1] 0 0 1 0
[1] 0 0 1 1
[1] 0 1 0 0
[1] 0 1 0 1
[1] 0 1 1 0
[1] 0 1 1 1
[1] 1 0 0 0
[1] 1 0 0 1
[1] 1 0 1 0
[1] 1 0 1 1
[1] 1 1 0 0
[1] 1 1 0 1
[1] 1 1 1 0
[1] 1 1 1 1

上の使用例の print 文の位置に必要な処理プログラムを書けばよい。

tricombinations に相当するプログラムは base = 3 とするだけでよい。

> n = 4
> base = 3
> x = integer(n)
> for (i in seq.int(base ^ n)) {
+   print(x)
+   x = next.Ncombination(x, base)
+ }
[1] 0 0 0 0
[1] 0 0 0 1
[1] 0 0 0 2
[1] 0 0 1 0
[1] 0 0 1 1
[1] 0 0 1 2
[1] 0 0 2 0
[1] 0 0 2 1
[1] 0 0 2 2
   中略
[1] 2 2 1 1
[1] 2 2 1 2
[1] 2 2 2 0
[1] 2 2 2 1
[1] 2 2 2 2

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

オリジナルの bincombinations のソースプログラム

bincombinations = function (p) {
    retval = matrix(0, nrow = 2^p, ncol = p)
    for (n in 1:p) {
        retval[, n] = rep(c(rep(0, (2^p/2^n)), rep(1, (2^p/2^n))),
            length = 2^p)
    }
    retval
}

これを書き換えて，3, 4, 5, ... 進数に対応させた tricombinations, quatrcombinations, ... quintcombinations など

tricombinations = function(p) {
    retval = matrix(0, nrow = 3^p, ncol = p)
    for (n in 1:p) {
        retval[, n] = rep(c(rep(0, (3^p/3^n)), rep(1, (3^p/3^n)), rep(2, (3^p/3^n))),
            length = 3^p)
    }
    retval
}

quatrcombinations = function(p) {
    retval = matrix(0, nrow = 4^p, ncol = p)
    for (n in 1:p) {
        retval[, n] = rep(c(rep(0, (4^p/4^n)), rep(1, (4^p/4^n)), rep(2, (4^p/4^n)), rep(3, (4^p/4^n))),
            length = 4^p)
    }
    retval
}

quintcombinations = function(p) {
    retval = matrix(0, nrow = 5^p, ncol = p)
    for (n in 1:p) {
        retval[, n] = rep(c(rep(0, (5^p/5^n)), rep(1, (5^p/5^n)), rep(2, (5^p/5^n)), rep(3, (5^p/5^n)), rep(4, (5^p/5^n))),
            length = 5^p)
    }
    retval
}

sextetcombinations = function(p) {
    retval = matrix(0, nrow = 6^p, ncol = p)
    for (n in 1:p) {
        retval[, n] = rep(c(rep(0, (6^p/6^n)), rep(1, (6^p/6^n)), rep(2, (6^p/6^n)), rep(3, (6^p/6^n)), rep(4, (6^p/6^n)), rep(5, (6^p/6^n))),
            length = 6^p)
    }
    retval
}

septetcombinations = function(p) {
    retval = matrix(0, nrow = 7^p, ncol = p)
    for (n in 1:p) {
        retval[, n] = rep(c(rep(0, (7^p/7^n)), rep(1, (7^p/7^n)), rep(2, (7^p/7^n)), rep(3, (7^p/7^n)), rep(4, (7^p/7^n)), rep(5, (7^p/7^n)), rep(6, (7^p/7^n))),
            length = 7^p)
    }
    retval
}

あとは，ご自由に定義してください

insert, pop, count, index, concat, remove を定義する

R には，ベクトルに要素を挿入する append があるが，関連する他の関数がないので作ってみた。

 

1. 要素の挿入　index(x, a, pos=length(x))　ベクトル x の pos 位置に a を挿入する。先頭に加える場合は pos=0 とする。

insert = function(x, a, pos=length(x)) {
  assign(deparse(substitute(x)), append(x, a, after=pos), envir = .GlobalEnv)
}

> x = c('w', 'q', 'p', 'w')
> insert(x, 'first', 0)
> insert(x, 'last')
> x
[1] "first" "w"     "q"     "p"     "w"     "last"

2. 要素の取り出し　pop(x, pos=length(x))　ベクトル x の pos 位置にある要素を取除く。取り除かれた値を返す。

pop = function(x, pos=length(x)) {
  assign(deparse(substitute(x)), x[-pos], envir = .GlobalEnv)
  return(x[pos])
}

> a = c('1', '2', 'a', 'x')
> pop(a, 1)
[1] "1"
> a
[1] "2" "a" "x"
> pop(a)
[1] "x"
> a
[1] "2" "a"

3. 要素の個数　count(x, a)　ベクトル x の中に a がいくつあるかを返す。

count = function(x, a) {
  sum(x %in% a)
}

> x = c('w', 'q', 'p', 'w')
> count(x, "w")
[1] 2

4. 要素の位置　index(x, a)　ベクトル x において， a がある位置を返す。なければ 0，複数ある場合は最初のものの位置を返す。

index = function(x, a) {
  ans = which(x == a)
  if (length(ans) >= 1) {
    return(ans[1])
  } else {
    return(0)
  }
}

> a = c('1', '2', 'a', 'x')
> index(a, '3')
[1] 0
> index(a, '1')
[1] 1
> index(a, 'a')
[1] 3

5. 連結　指定した複数のベクトルを連結したものを返す。

concat = function(a, b, ...) {
  others = list(...)
  ret = c(a, b)
  len = length(others)
  if (len != 0) {
    for (i in seq_along(others)) {
      ret = c(ret, others[[i]])
    }
  }
  return(ret)
}

> alp = c("a", "b", "c")
> num = c("1", "2", "3")
> sym = c("@", "*", "/")
> foo = c("x", "y", "z", "w")
> concat(alp, num)
[1] "a" "b" "c" "1" "2" "3"
> concat(alp, num, sym, foo)
 [1] "a" "b" "c" "1" "2" "3" "@" "*" "/" "x" "y"
[12] "z" "w"

6. 消去　remove(x, a)　ベクトル x から，要素 a を取り除く。

remove = function(x, a) {
  ans = which(x == a)
  if (length(ans) >= 1) {
    assign(deparse(substitute(x)), x[-ans], envir = .GlobalEnv)
  }
}

> x = c('a', 'b', 'c')
> remove(x, "w")
> x
[1] "a" "b" "c"
> remove(x, "a")
> x
[1] "b" "c"

 

連続変数をカテゴリー化する

連続変数をカテゴリー化する関数としては findInterval 関数がある。

set.seed(123)
x = rnorm(100000)
v = c(-Inf, -1.2, -0.8, 0, 0.5, 1.4, 2.2, Inf)

system.time({
a = findInterval(x, v)
table(a)
})

もし，以下のような関数を作ってみると，findInterval 関数より，15倍遅い。

system.time({
b = sapply(x, function(y) sum(y > v))
table(b)
})

四捨五入関数は round ではない（かも）

三省堂 大辞林によれば，四捨五入とは，「必要とする位の次の位の数字が四以下ならばこれを切り捨て、五以上ならば切り上げて、上の位に一を加える方法」とある。

R では，round 関数がこれを行う。

1.5 を，小数第1位で四捨五入を行うと 2 である。
> round(1.5)
[1] 2

しかし 0.5 を，小数第1位で四捨五入を行うと 0 である。

 > round(0.5)
[1] 0

よって，round は四捨五入を行う関数ではない。

round は，偶数への丸め（最近接丸め，JIS 丸め，ISO 丸め，銀行丸めなどとも呼ばれる）を行う関数である。これは，「端数が 0.5 より小さいなら切り捨て、端数が 0.5 より大きいならは切り上げ、端数がちょうど 0.5 なら切り捨てと切り上げのうち結果が偶数となる方へ丸める。（JIS Z 8401 の規則A）」として定められており、規則B（四捨五入）より「望ましい」とされている。

インターネット上の多くのページでは，プログラム言語（処理系）によって違いがあるとかの記述はあるが，実際に四捨五入を行うにはどのような関数を定義すればよいかはほとんど書かれていない。

そこで，R においては以下の関数を定義する。

round2 = function(x, d=0) {
    p = 10^d
    return((x * p * 2 + 1) %/% 2 / p)
}

この関数はまず，四捨五入を行う位置 d により，対象とする数 x に p = 10^d をかける。たとえば x = 0.45, d = 1 ならば x = 0.45 * 10^1 = 4.5 となる。ここで小数点以下を四捨五入すれば 5 となる。同じことをするのに x = (x * 2 + 1) %/% 2 としている。x = (4.5 * 2 + 1) = 10 を 2 で割って整数部を取る。つまり x = 5 となるので四捨五入していることになる。なお，(x * p * 2 + 1) %/% 2  は as.integer(x * p + 0.5) と同じである。x * p の整数部 と 0.5 は共にコンピュータ上で正確に表現できる数なので，誤差はない。最後に得られた数を p で割る。つまり x = 5 / 10 = 5 が答えである。

> round2(0.45, 1)
[1] 0.5

これに対して，元からある round は以下のような結果になる。
> round(0.45, 1)
[1] 0.4

解釈は，「0.45 という数値がコンピュータ上で正確に表されているなら，小数点第2位を丸めるときに，切り上げて 0.5 とする場合も，切り下げて 0.4 とする場合もその差の絶対値は 0.05 である。このように丸めの候補値が 2 つあるときは，結果の末尾が偶数のほうを採用する。つまり，この場合だと 0.4 になる」ということである。

注意すべきは「0.45 という数値がコンピュータ上で正確に表されているなら」というところ。0.45 は 10 進数では循環小数ではないが 2 進数では循環小数で内部的には誤差を含む数値であることもある。

0.45 ≒ 0x1.ccccccccccccdp-2 = (1+12/16^1 + 12/16^2 + … + 12/16^12 + 13/16^13) *2 ^-2

> round(0.8500000000000000, 1)
[1] 0.8

となるのは正しい結果であるが，

> round(0.8500000000000001, 1)
[1] 0.8

は 0.8500000000000001 が内部的に正しく表現されているなら，本来は 0.9 にならないといけない。

> round(0.8500000000000002, 1)
[1] 0.9

は正しい結果である。

つまり，コンピュータ内部では，0.8500000000000001 と 0.8500000000000000 は異なる数であるが，

> 0.8500000000000001 == 0.8500000000000000
[1] FALSE

その差は末尾 1 ビット

> 0.8500000000000001 - 0.8500000000000000 < EPSILON
[1] TRUE

> sprintf("%a %a %a %a %a", EPSILON, 0.8499999999999999, 0.8500000000000000, 0.8500000000000001, 0.8500000000000002)

EPSILON    =　0x1p-52
0.8499999999999999　=　0x1.b333333333332p-1
0.8500000000000000　=　0x1.b333333333333p-1
0.8500000000000001　=　0x1.b333333333334p-1
0.8500000000000002　=　0x1.b333333333335p-1"

いずれの場合も，四捨五入関数 round2 は正しく機能している。

> round(0.8499999999999999, 1); round2(0.8499999999999999, 1)
[1] 0.8
[1] 0.8
> round(0.85, 1); round2(0.85, 1)
[1] 0.8
[1] 0.9
> round(0.8500000000000001, 1); round2(0.8500000000000001, 1)
[1] 0.8
[1] 0.9
> round(0.8500000000000002, 1); round2(0.8500000000000002, 1)
[1] 0.9
[1] 0.9

なかなか頷きがたい

数学的に考える「じゃんけん」で有利な手は何か

> 725人から集めたのべ11567回のじゃんけんデータからみたじゃんけんの確率はどうだったのか？　　集計結果は以下である。グーが4054回、チョキが3664回、パーが3849回。したがって、「グーを出す人が多いことから一般にじゃんけんではパーが有利」だと導かれた。

それぞれ，35.0%, 31.7%, 33.3% である。一番弱いチョキと比べてグーだと3%しか有利ではない。

確かに，のべ11567回の結果からは p 値 = 5.166e-05 などと言うことになるんだが，これはいかに何でもサンプルサイズが大きすぎる。それぞれの手の割合が同じでも，サンプルサイズを 1/4 にしたら，p 値 = 0.0826 で有意ではない。

サンプルサイズが 2892 でも相当大きい。ねえ！世論調査や視聴率調査のときのサンプルサイズは 1000 以下なんだから。

サンプルサイズが 1 万を超えると，母比率のわずかな違いを検出することになる。

そんな調査結果も，

> 心理学的には、「人間は警戒心をもつと拳を握る傾向がある」という説明のほか、「チョキはグーやパーと比べて出しがたい手である」という説明もある。

まあ，どうとでも解釈のしようはあるとは言うものの，『心理学』がからむと，眉唾感が満載である。

また，

> 4桁の数字は全部で10000個あるので、4桁の数字が全部異なる確率は50.4％しかない。そして、重複のある4桁の数字が当せん番号になる確率は49.6％もあるのだ。

「ナンバーズ4」について書いている部分だが，0.8 % の差が実質的にどの程度の意味があると？どうせなら 49.6% の方を買おうと？みんながそちらを買うようになると，裏をかいて，50.4% の方を買う方が当たったときの当選金は高くなる。そのまた裏を考えて 49.6% の方を買う人が多くなるので... 以下，永遠に続く（笑）

 

Mac でのファイル名（UTF-8-MAC）

こっちにも控えておこう

list.files() で読み込んだファイル名に濁音，半濁音が含まれるときの対処法

> y = list.files()
> x = c('織田信長', 'パンダ')
> x %in% y
[1]  TRUE FALSE
> x = iconv(x, "UTF-8", "UTF-8-MAC")
> x %in% y
[1] TRUE TRUE

 

二次方程式の解を求めるプログラムのレビュー

二次方程式の解を求めるプログラム例はネット上でもよく見られる。
本来（？）お手本を示すべきものであるが，問題を含むものが少なくない（多いということだ）。

結論を先に述べておく。

(0) 全部で 40 サイトを対象とした
(1) float を使っているのは 8 サイト。そのうち C によるものは 6 サイト。C は 12 サイトだが，そのうち float を使っているのは 6 サイト。C はよく使われているが，不適切なサイトも多い。
(2) 桁落ちを考慮しているサイトは 8 サイト（20%）。
(3) 虚数型を使用しているのは 3 サイトに過ぎない。
(4) 虚数解の場合を対象としていないのは 17 サイト。実部虚部を表示するのは 19 サイト。

＝＝以下は各サイトのレビュー＝＝

●● https://webkaru.net/clang/quadratic-formula/

 は C 言語によるものだ。

いくつか問題点がある。
(1) 使用しているのが double ではなく float --- float を使うのに何の意味があるのか。float など使う癖がつくと，後々困る。
(2) 実数解を求めるとき，2 つの解を解の公式だけで求めている。
    kai1 = ( -b + sqrt(discriminant) ) / (2*a);
    kai2 = ( -b - sqrt(discriminant) ) / (2*a);
これは，次のような方程式を正しく解けない。
2次方程式の定数を入力してください。
a = 1
b = 1.0000000001
c = 0.0000000001
2次方程式の解: x = 0.00, -1.00
答えは printf で書かれているので少数以下を 2 桁に指定しているのでこうなる。また，上述のように float を使っているので有効桁数は 6 桁ほどしかない。そこで，float は全て double にする。係数を端末から読み込むのに scanf("%f", &a) などと，お勧めできない書法で書かれているが，最小限の手当で，scanf("%lf", &a) にする。さらに，出力書式もより広範囲の数値を表現できるよう %g を使う。
正しい答えはいうまでもなく，x = -1, -1e-10 である。ここまでの改変後のプログラムでは
2次方程式の定数を入力してください。
a = 1
b = 1.0000000000000001
c = 0.0000000000000001
2次方程式の解: x = -1.11022e-16, -1
が得られる。絶対値の小さい方の解の誤差がひどい。原因は先に述べたように，2つの解を解の公式をそのまま使って求めているからである。大きい方の解をまず求めて，絶対値の小さい方の解は解と係数の関係から求めるようにしなければならない。
    if (b < 0) {
      kai1 = ( -b + sqrt(discriminant) ) / (2*a);
    }
    else {
      kai1 = ( -b - sqrt(discriminant) ) / (2*a);
    }
    kai2 = (c/a)/kai1;
これでやっと，正しい答えが得られる。
2次方程式の定数を入力してください。
a = 1
b = 1.0000000000000001
c = 0.0000000000000001
2次方程式の解: x = -1, -1e-16

このプログラムではちゃんと虚数解も求めるように作られている。

2次方程式の定数を入力してください。
a = 1    
b = 1
c = 1
2次方程式の解: -0.5+0.866025i, -0.5-0.866025i
前にも書いたが，解は printf で書かれている。実部と虚部を文字列としてつなぎ合わせているだけである。
printf("2次方程式の解: %.2f+%.2fi, %.2f-%.2fi\n", real, imag, real, imag);
普通，二次方程式の解は別の計算のために使われるのだから，文字列でコンソールに表示するだけでは意義は半減する。C は複素数も扱えるのだから，複素数にすればよいのだ。もっとも，その場合に結果を出力するときは，creal, cimag で実部と虚部を取り出して出力しなければならないが。
#include
:
  double _Complex x1, x2;
:
    x1 = real+I*imag;
    x2 = conj(x1);

●● https://na-inet.jp/research/equation_cpp.pdf

は，ここに挙げたものの中では，一番優れていると思う。
C++ で，高校生向けの授業用資料であるが，「「桁落ち」現象と呼びます。これを解決する手段はいくつかありますが，ここでは最後の節で，倍精度浮動小数点数よりずっと多い桁数で計算する「多倍長浮動小数点数」というものを使って強引に解決していきます」と勇ましい（^_^;）。そこまで行かなくても，解の公式と解と係数の関係を使えばよいとは思う。
また，虚数解にも対処している。解を虚数型変数で求めている。
complex<double> a, b, c;
complex<double> d, sol[2];
cout << "Solutions: " <
虚数解も，
Solutions: (1,1.41421), (1,-1.41421)
のように自然に表示される。

●● https://python-programming-iot.com/2018/08/26/61/

は，Python によるもの。sqrt(x) を x**0.5 とするなど，のけぞるが，Python では，x が負の場合には（自動的に）虚数になる。つまり，つねに虚数型を使って解の公式により解を求める（これはユニーク）。x が正の場合には x**0.5 は実数であるので，実数型の解が求まるが，桁落ちを考えていないので小さい方の解の精度が低くなることがある。

●● https://books.google.co.jp/books?id=bpYhLwff3zgC&pg=PA58&lpg=PA58&dq=%E4%BA%8C%E6%AC%A1%E6%96%B9%E7%A8%8B%E5%BC%8F+%E3%83%97%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%A0&source=bl&ots=8B4TIiQMiT&sig=ACfU3U3i-YLWrRxyy1T0PtBi6wZnosrGqA&hl=ja&sa=X&ved=2ahUKEwjTpdiRtuPjAhWtUN4KHT5nC3w4PBDoATASegQIBxAB
書名: Javaで学ぶ数値解析，著者: 和光システム研究所

は，Java によるもの。double を使う，解の公式と解と係数の関係を使う。虚数解は実部と虚部を表示する。

●● https://blog.masuqat.net/2014/04/28/cancellation-of-significant-digits-of-quadratic-formula-with-csharp/

は，Java によるもの。float を使う。解の公式のみを使うのでは桁落ちがあるので，解と係数の関係を使うべしとの指摘をしている。虚数解の場合には，実部と虚部を表示する。

●● http://aoki2.si.gunma-u.ac.jp/Hanasi/Algo/letsc/1-body.html

は，double を使う，解の公式と解と係数の関係を使う。虚数解の場合は「実数解はありません」と表示する。

●● http://www.ritsumei.ac.jp/se/~osaka/rejime/suuti/equation.pdf

は，桁落ちの起きる原因（float を使う，解の公式をそのまま使う）を説明し，それを避ける方法（float のままであったとしても，解の公式と解と係数の関係を使う）を示している。

●● http://www.crl.nitech.ac.jp/~ida/education/MaterialsDesign/1453.html

は，プログラムの提示はなく，アルゴリズムを示すもの。解の公式のみを使うのではなく，解と係数の公式を使うべきであると述べている。

●● http://www.ne.phen.mie-u.ac.jp/misc/equation.pdf

では，b^2 - 4ac が「正確に0であるかどうか」に注意するようにアドバイスをしている。
> d = b^2 - 4acが0かどうかをを正しく判断するためには，fabs(d/(b*b+fabs(4.0*a*c)))

しかし，2つの解を共に解の公式をそのままあてはめたり，虚数解を扱うとき複素数型を使わないという欠点は最初のプログラムと同じである。変数としては double を使い，出力にも %lg を使うことを勧めているのは当たり前のことではあるが，よい。

●● https://takuyab.com/archives/2613

は，「プログラムに解かせる為にはおなじみ“解の公式”を使います」と堂々と書いてある，float を使っている，虚数解は出力として形式的表示。

●● https://qiita.com/watarumohawk/items/c6fa8dfc9e2f2a23a53b

は，変数の型は double であるが，2つの解を共に解の公式をそのままあてはめたり，虚数解を扱うとき複素数型を使わないという欠点は最初のプログラムと同じである。

●● https://seaotter.cite.tohoku.ac.jp/coda/clang/c-3-branching.html

は，float を使い，2つの解を共に解の公式をそのままあてはめている。虚数解は考慮されていない（「実数解はない」と表示するのみ）。

●● https://algorithm.joho.info/programming/c-language/equation-solution-c/

は，関数形式のプログラムであるが，変数の型は double であるが，2つの解を共に解の公式をそのままあてはめたり，虚数解を扱うとき複素数型を使わない。

●● https://keisan.casio.jp/exec/system/1161228770

は，さすが CASIO だけあって，a = 1，b = 1.0000000000000001，c = 0.0000000000000001 のような場合であっても，誤差のない計算結果を返します。

●● http://pc-physics.com/secondequation1.html

は，「特にコンピューターでは虚数が扱えないので、虚数解になる場合注意が必要です」という，不穏当な注意書きがある。変数の型は double であるが，2つの解を共に解の公式をそのままあてはめたり，虚数解を扱うとき複素数型を使わない。

●● https://www.kkaneko.jp/cc/program/equation.html

は，解の公式と解と係数の公式を使うプログラムを示しています。しかし，虚数解の場合は，
    printf("x1 = 虚数解\n");
    printf("x2 = 虚数解\n");
という，ビックリするような出力をする。

●● http://roko.eng.hokudai.ac.jp/studentadm/chiba_data/pro/jp1.pdf

は，C，C++，FORTRAN の場合を示している。しかし，基本的に float を使う，解の公式のみを使う，虚数解の場合は「実数解はない」と表示するのみ。

●● http://kentiku-kouzou.jp/keisan-2.html

は，前述の CASIO の計算サイトと同じく，計算サービスのみでソースプログラムは提示していない。しかし，
a = 1，b = 1.0000000000000001，c = 0.0000000000000001 の場合，-1.1102230246251565e-16，-0.9999999999999999 を返すので，2つの解とも解の公式で求めているのだろう。「建築学生が学ぶ構造力学」というページなので，ちょっと不安になる。

●● https://okayan08.hatenadiary.org/entry/20081219/1229677932

でギョッとさせられるのは，dの符号変換を d *=-1 とやっているような所（^_^;）この人は -b というような書き方ができないと思っているらしい。いつも -1*b+d と書いている。単項演算の - の存在を知らない（こういう人は案外多い）。
解の公式のみを使う，虚数解は実部と虚部を表示する。

●● http://www.u.tsukuba.ac.jp/~watanabe.norio.fw/H24-2FSubject.pdf

は，C によるもの。関数型で，double を使い，解の公式をそのまま使う，虚数解は実部と虚部を返すという仕様。「解答例」と称するにはちょっと不十分。

●● http://www.hino.meisei-u.ac.jp/is/iga/hit-u/No9/C3.pdf

は，float を使う，解の公式をそのまま使う，虚数解の場合は「実数解はありません」とのみ表示。

●● http://masudahp.web.fc2.com/jfbas/kiso/mon44.html
●● http://masudahp.web.fc2.com/n88basic/kiso/mon33.html

は，珍しく BASIC。解の公式をそのまま使う，虚数解の場合は「解なし」とのみ表示。

●● https://riptutorial.com/ja/fortran/example/3115/二次方程式

は，珍しく FORTRAN。real(float) を使う，解の公式そのままを使う，虚数解の場合は "No real roots." とのみ表示。

●● http://nalab.mind.meiji.ac.jp/~mk/labo/studying-C/Programing-in-C.pdf

は，C によるもの。double を使うが，解の公式そのままを使う，虚数解は実部と虚部を表示する。

●● https://wwws.kobe-c.ac.jp/deguchi/c/prob/eq.html

は，C によるもの。float を使う，解の公式そのままを使う，虚数解は「解なし」と表示する。

●● https://www.is.kochi-u.ac.jp/kyoko/edu/ip2010/python.html

は Python。必然的に double を使う。しかし，解の公式そのままを使う，虚数解は「実数解はありません」と表示する。

●● http://blog.livedoor.jp/add20/archives/897767.html

は，C, Perl, JavaScript。基本的にアルゴリズムが同じなら，同じような結果が表示される。解の公式そのままを使う，虚数解は「実数解は存在しない」と表示する。

●● http://macroscope.world.coocan.jp/ja/edu/computer/sabun/nizi.html

は，BASIC, Pascal, FORTRAN。基本的にアルゴリズムが同じなら，同じような結果が表示される。改定版として，double を使い，解の公式と解と係数の公式を使う方法が書かれている。虚数解の場合には実部，虚部を書くだけで，虚数型は使っていない。

●● http://nlp.dse.ibaraki.ac.jp/~shinnou/zemi2009/php/php-saito-1106.pdf

は，珍しい PHP でのプログラム。しかし，解の公式そのままを使う，虚数解は「実数解を持たない」と表示する。

●● https://excelmath.atelierkobato.com/quadratic-macro/

は，VBA によるもの。解の公式のみを使う。虚数解の場合は実部と虚部を表示する。

●● http://stroll.hatenablog.com/entry/2015/08/16/212016

は，Java によるもの。ここでも，-1 * b + Math.sqrt(d) のように，単項演算子 "-" を知らない表記をしているので，あきれる。いろいろゴタゴタ飾り立て，119 行のプログラムになっている。double を使っているものの，解の公式のみを使う。虚数解の場合は実部と虚部を表示する。

●● http://www.kodama-lab.com/seminar/javalang/summer/ex001.html

は，Java によるもの。解の公式そのままを使う，虚数解は「None」と表示する。

●● https://www.kitasato-u.ac.jp/sci/resea/buturi/hisenkei/sogo/poly4.pdf

は，Python によるもの。
import sys, math, cmath
I=cmath.sqrt(-1)
omega1=(-1+cmath.sqrt(-3))/2
などの記述がある。I はこのように定義せずとも 1j でよい。
s1*s1*s1 などが頻繁に出てくるが s1**3 などの方が良いだろう。
3次方程式で，間接的ではあるが虚数型を使って解の公式のみを使っている。

●● http://newtral.blog.jp/archives/1010749637.html

は，PHP によるもの。解の公式のみを使う。虚数解の場合は実部と虚部を表示する。

●● https://codeday.me/jp/qa/20190425/696210.html

は珍しい Prolog によるもの。解の公式のみを使う。虚数解の場合は「空リスト」を返す。

●● http://vipprog.net/wiki/解答例/二次方程式の解/Excel%20VBA.html

は，VBA(Visual BASIC) によるもの。Double を使う。解の公式のみを使う。虚数解の場合は実部と虚部を表示する。

●● https://www.e-bridge.jp/eb/tcontents/plain_basic/ChapC21.html

は，Plain Basic によるもの（行番号がある）。あれは double だったっけ？解の公式のみを使う。虚数解の場合は「根がありません」と表示する。Plain Basic といい，「根」というだけあって，古い。

●● http://fullbasic.web.fc2.com/510_two_degree.html

は，full Basic によるもの（行番号がある）。穴埋め問題。「センター試験数学」なのでやむを得ないが，解の公式のみを使う。虚数解の場合は「解なし」と表示する。

●● https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=78&ved=2ahUKEwjTpdiRtuPjAhWtUN4KHT5nC3w4PBAWMBF6BAgEEAI&url=http%3A%2F%2Fwww.ocw.titech.ac.jp%2Findex.php%3Fmodule%3DGeneral%26action%3DDownLoad%26file%3D201017211-90-0-50.pdf%26type%3Dcal%26JWC%3D201017211&usg=AOvVaw2tSEKr0ecdF7qKxKG59MlB

は，C プログラムによるもの。double を使い，解の公式のみを使う。虚数解は実部と虚部を表示する。

●● https://teenaka.at.webry.info/201701/article_5.html

は，Python によるもの。解の公式のみを使う。虚数解は実部と虚部を表示する。

●● http://www.akita-pu.ac.jp/system/elect/ins/kusakari/japanese/teaching/Old/Programming/2005/note/5/Slide03.html

は，Java によるもの。解の公式のみを使う。虚数解の場合は「実数解はない」と表示する。

●● http://www.asahi-net.or.jp/~YY8A-IMI/20040913/x68/c_1.htm

は，C によるもの。解の公式のみを使う。虚数解は実部と虚部を表示する。

●● http://www.avis.ne.jp/~ammonite/javatest1.htm

は，Java によるもの。解の公式のみを使う。虚数解の場合は「回答なし」。

 

紙は何回折れるか？

かなり有名な事実で，「ほぼ8回」というのが正解。

厚さ 0.1 mm，一辺 100 m の紙を，偶数回ごとに正方形になるように 8 回折れば，厚さは 0.1 / 1000 * 2^8  = 0.0256 m で，正方形の一辺の長さは 6.25 m。まだ十分に折れそうではあるが，10 回折ると，厚さは 0.1024 m，一辺の長さ3.125 m だが，折られた紙のそれまでの折り線付近が折りにくくなる。フットボール競技場大の紙を折るときは 11 回折れたという動画がある。

柔らかい紙を一方向に折っていくともう少し多く折れる。1200メートルのトイレットペーパーを半分ずつに折るのは 12回できたという写真がある（ワンカットの動画でないので正しさの証明はないが）。このとき，最終のトイレットペーパーの長さは 1200*0.5^12 =  0.2929688 で，まあ写真にあるように肘の長さ。

以上のように折るのは難しいが，半分ずつに切って，それを積み重ねるというのは思考実験としてはもう少し可能性が高い。

上のように，厚さ 0.1 mm，一辺 100 m の紙を，偶数回ごとに正方形になるように 8 回切って積み重ねる。もちろん，カッターなどを使って一回で切る。あるいはそれをシミュレートするように分割して切る。こうすれば 25 回切ると 0.1 / 1000 * 2^25 = 3355.443 m でほぼ富士山と同じ高さ。しかし，このときの一枚の紙切れは 1.22 cm × 2.44 cm だ。とても積み重ねることはできない。25 回目の紙を切る長さは延べ 409600 m。一枚一枚ハサミで切っていたら大変。全部で 1228600 m = 1228.6 km。札幌〜福岡の直線距離が約 1400 km。

高さはものすごい勢いで大きくなり，面積は逆にものすごい勢いで小さくなる。

一休さんかだれかが，殿様が褒美をやろうというのに答えて，「今日は碁盤の1マスに1粒のお米，明日は2倍の2粒，その次にはさらに2倍の4粒というようにして最後の一マスまでお願いします」というのも有名な話。一般的な碁盤のマス目は 18*18 = 324。324 番目のマス目には 2^(324-1) =1.70879e+97 個の米粒が載る（わけない）。碁盤全部では sum(2^(0:323)) = 3.417579e+97 粒。

これを逆方向に使ったのが，f(x) = 0 の解を求める二分法というアルゴリズム。-1000 から 1000 までの区間を半分にし，解のある方の区間を更に半分にしていく。30 回繰り返すと解の存在する区間幅は 2000*0.5^30 = 1.862645e-06 になる。

もし 2000 ページの辞書があって，調べたい単語が載っているページを探すとき，ほぼ半分ぐらいの所を開き，前のほうにあるか後ろのほうにあるかだけの情報で単語のある方の更に半分のページを開き...ということを繰り返すと 10 回目には 2ページが残る。 000*0.5^10 = 1.953125。11 回目で単語を含むページに達する。

この例は，単語がまさに辞書順に載っているからできることである。2000人の名前が名簿がコンピュータに記憶されているとき，でたらめに並んでいると，ある人を探すとき，最短では 1 回，最長では 2000 回，平均でも 1000 回の探索が必要だが，辞書順に並べ替えておくと最長でも 11 回で探し出せる。

 N.PRIME 個までの素数の配列 prime がある。整数 n が何番目の素数であるか探索するプログラムを書く。ただし，prime に含まれない数が与えられた場合は not found と表示する。

N.PRIME = 2000

素数配列を作成する

library(gmp)
prime = integer(N.PRIME)
a = 1
for (i in 1:N.PRIME) {
  prime[i] = a = as.integer(nextprime(a))
}

head(prime)
tail(prime)

線形探索（端から順番に探していく）

find1 = function(n) {
  found = 0
  for (i in 1:N.PRIME) {
    if (n == prime[i]) {
      found = i
      break
    }
  }
  if (found == 0) {
    print("not found")
  } else {
    print(sprintf("%d = prime[%d]", n, found))
  }
}

find1(1)
find1(2)
find1(7919)
find1(7909)
find1(17389)
find1(20000)

バイナリーサーチ（二分探索）

find2 = function(n) {
  found = 0
  begin = 1
  end = N.PRIME
  n.begin = prime[begin]
  n.end = prime[end]
  while (n.begin <= n && n.end >= n) {
    if (n.end == n) {
      found = end
      break
    }
    if (end == begin + 1 && prime[end] != n && prime[begin] != n) {
      break
    }
    middle = (begin + end) %/% 2
    n.middle = prime[middle]
    if (n.middle == n) {
      found = middle
      break
    } else if (n.middle > n) {
      end = middle
    } else {
      begin = middle
    }
  }
  if (found == 0) {
    print("not found")
  } else {
    print(sprintf("%d = prime[%d]", n, found))
  }
}

find2(1)
find2(2)
find2(7919)
find2(7909)
find2(17389)
find2(20000)

R の機能を最大限利用するプログラム

find3 = function(n) {
  pos = prime == n
  if (sum(pos) == 0) {
    print("not found")
  } else {
    found = (1:N.PRIME)[pos]
    print(sprintf("%d = prime[%d]", n, found))
  }
}

find3(1)
find3(2)
find3(7919)
find3(7909)
find3(17389)
find3(20000)

2直線の交点座標を求めるプログラム（Python 版）

何故か相も変わらず「2直線の交点座標を求めるプログラム」がよく参照されている

そこにあるのは R プログラムなので，Python で書いたものをおまけに載せておこう

import numpy as np

def prog(x1, y1, x2, y2, x3, y3, x4, y4):
    if x1 == x2 and x3 == x4:
        x = y = np.nan
    elif x1 == x2:
        x = x1
        y = (y4 - y3) / (x4 - x3) * (x1 - x3) + y3
    elif x3 == x4:
        x = x3
        y = (y2 - y1) / (x2 - x1) * (x3 - x1) + y1
    else:
        a1 = (y2-y1)/(x2-x1)
        a3 = (y4-y3)/(x4-x3)
        if a1 == a3:
            x = y = np.nan
        else:
            x = (a1*x1-y1-a3*x3+y3)/(a1-a3)
            y = (y2-y1)/(x2-x1)*(x-x1)+y1
    return (x, y)

なんでこんなものが検索されるのやら？