消費増税ーーなんで，理性的な判断ができないか？

今回の消費増税にあたって，駆け込み購入が目立たないという報道が多数あったが，ここ数日，それを煽るような報道が多く見られるようになった。

たとえば，リクルートスーツ。定価 2万5千円 を駆け込み購入するとか。

従来の消費税なら税率 8%  で，税込み 2万7千円。消費税分は 2000 円。

増税後は税率 10% で，税込み 2万7千500円。消費税分は 2500 円。

差額は 500円。これは大きいか？

リクルートスーツ。就職活動の後も使うだろう。何年間？まあ，2年としても，2 年間の内，500 円くらいのラーメンを一回我慢するだけで 500 円はカバーできる。2 年間の内の一回のラーメン代も我慢できないというなら，それは仕方ない。普通の人ならできるだろう。

駆け込み購入は，トイレットペーパーから猫の餌までということだ。

前回の消費増税のとき，ふりかけ（ごはんにかけるふりかけですよ！！）を 3 年分くらい買い込んだ アファー な人がいたということを聞いて，思わず笑ってしまった。

個々人にとってはどうって事がない話のはずなのに，それが何十万，何百万，何千万人の話になると消費増税後の不況になる。

もし，10万人の消費者が駆け込み購入すると 10万 × 500円だから 5千万円になるのだ。増税後は 5千万の売上減になる（それ以前に 5千万売り上げたので，トントンなのだけど）。それが，増税不況だ。

お互い，不幸な話ではないですか。

もっと，冷静になった方がどれだけ幸せか。

 

多倍長精度による円周率 π の計算

「多倍長精度による計算」に，R による円周率 π の計算プログラムが掲載されている。R でも gmp パッケージを使えば多倍長精度による計算はできるので，このプログラムを書き換えればよい。

このプログラムは円周率を10の羃乗倍した固定小数点数として計算する。

（浮動小数点数で計算するプログラムは「多倍長精度計算パッケージ Rmpfr」に掲載した。）

末尾に Python と Julia のプログラムを追記した。

library(gmp)
calc.pi = function(accuracy=10000)    {
  FACTOR = as.bigz(10)^accuracy # FACTOR 倍した固定小数点数として計算
  pi = as.bigz(0)
  a = as.bigz(16*5*FACTOR)
  n = 1
  repeat {
    a = a %/% 25
    b = a %/% n
    if (b == 0) break
    pi = pi + b
    n = n + 2
    a = a %/% 25
    b = a %/% n
    if (b == 0) break
    pi = pi - b
    n = n + 2
  }
  a = as.bigz(4*239*FACTOR)
  n = 1
  repeat {
    a = a %/% 57121
    b = a %/% n
    if (b == 0) break
    pi = pi - b
    n = n + 2
    a = a %/% 57121
    b = a %/% n
    if (b == 0) break
    pi = pi + b
    n = n + 2
  }
  return(pi)
}

> calc.pi(50)
Big Integer ('bigz') :
[1] 314159265358979323846264338327950288419716939937510

Python でも gmpy モジュールを使えば多倍長精度による計算ができる。

移植に際して注意すべきことは，除算は / ではなく，整数除算 // を使うこと。また，π は浮動小数点数であるが，指定した桁の固定小数点数として扱う。具体的には小数点以下の桁数（acculacy で指定）の下に小数点を置いて，整数値 π x acculacy を計算する。

def calc_pi(acculacy=10000):
    FACTOR = 10**acculacy
    pi = 0
    n, a = 1, 16 * 5 * FACTOR
    while True:
        a //= 25
        b = a // n
        if b == 0:
            break
        pi += b
        n += 2
        a //= 25
        b = a // n
        if b == 0:
            break
        pi -= b
        n += 2
    n, a = 1, 4 * 239 * FACTOR
    while True:
        a //= 57121
        b = a // n
        if b == 0:
            break
        pi -= b
        n += 2
        a //= 57121
        b = a // n
        if b == 0:
            break
        pi += b
        n += 2
    return pi

>>> calc_pi(50)
314159265358979323846264338327950288419716939937510

実際の小数点は，言うまでもなく，先頭桁の 3 の次にある。

なお，Python には，sympy があり，それを用いれば，任意の精度の π が簡単に得られる。

>>> import sympy
>>> sympy.pi.evalf(50)
3.1415926535897932384626433832795028841971693993751

50 桁までの指定では，sympy.pi と比べると，最後まで合っている。普通は，最後の2桁ぐらいに誤差が出る。

Julia の場合も裏で gmp ライブラリを呼んでいるので，ほぼ同じように書くことができる。
多倍長精度の整数は big(数) で定義し，整数除算は ÷，if (b == 0) break が b == 0 && break と書けるなど。

function calcpi(accuracy = 10000)
    factor = big(10) ^ accuracy
    pi = big(0)
    a = 16 * 5 * factor
    n = big(1)
    while true 
        a ÷= 25
        b = a ÷ n
        b == 0 && break
        pi += b
        n += 2
        a ÷= 25
        b = a ÷ n
        b == 0 && break
        pi -= b
        n += 2
    end
    a = 4 * 239 * factor
    n = big(1)
    while true 
        a ÷= 57121
        b = a ÷ n
        b == 0 && break
        pi -= b
        n += 2
        a ÷= 57121
        b = a ÷ n
        b == 0 && break
        pi += b
        n += 2
    end
    pi
end

長精度整数演算

Sum of three cubes for 42 finally solved -- using real life planetary computer

3 つの整数の 3 乗和が 42 になる解が見つかったとのこと。

なんで，42 なのかは興味がないので説明しない。

かなり大きな整数なので，R の普通の変数型 numeric(double) では扱えないので，gmp ライブラリを使う必要がある。

> library(gmp)
> x = as.bigz("-80538738812075974")
> y = as.bigz("80435758145817515")
> z = as.bigz("12602123297335631")
> ans = x ^ 3 + y ^ 3 + z ^ 3
> ans
Big Integer ('bigz') :
[1] 42

> x ^ 3
Big Integer ('bigz') :
[1] -522413599036979150280966144853653247149764362110424
> y ^ 3
Big Integer ('bigz') :
[1] 520412211582497361738652718463552780369306583065875
> z ^ 3
Big Integer ('bigz') :
[1] 2001387454481788542313426390100466780457779044591

Python の整数型は通常 64 ビット整数だが，必要に応じて精度が増やされる。ので，そのまま計算できる。

>>> x = -80538738812075974
>>> y = 80435758145817515
>>> z = 12602123297335631
>>> ans = x ** 3 + y ** 3 + z ** 3
>>> print("ans =", ans)
ans = 42

>>> x ** 3
-522413599036979150280966144853653247149764362110424
>>> y ** 3
520412211582497361738652718463552780369306583065875
>>> z ** 3
2001387454481788542313426390100466780457779044591

マシュマロを食べることができるだろうか？

本日(2019/09/08) の日本テレビ，世界の果てまで行ってQ で，100m の高さから 10 個くらいのマシュマロを一度に投げて，下で待ち受ける8人ほどが口に入れることができるだろうかという実験をやっていた。結果は失敗だったけど，どんな条件だと成功する確率がどれくらいになるのかシミュレーションしてみた。

library(MASS)

# 9 人を 40 センチ間隔の格子点に配置（単位はセンチメートル）
x = c(-40, -40, -40, 0, 0, 0, 40, 40, 40)
y = c(-40, 0, 40, -40, 0, 40, -40, 0, 40)
person = length(x)

sigma = 4000 # 二次元正規分布の分散
marshmallow = 10 # 一度に投げるマシュマロの数
trial = 1000 # 実験回数
result = integer(trial) # マシュマロが何人の口に入ったかの記録
for (loop in 1:trial) {
    # 二次元正規乱数の発生
    d = mvrnorm(n=marshmallow, mu=c(0, 0), Sigma=matrix(c(sigma, 0, 0, sigma), 2), empirical=TRUE)
    # plot(d, cex = 0.5)
    # points(x, y, pch=19, col=2, cex = 0.5)
    success = 0 # 何人の口に入ったか
    for (i in 1:marshmallow) { # それぞれのマシュマロの着地点について
        for (j in 1:person) { # 各人の口から半径 5 センチ以内なら口に入るだろうと判定
            if ((d[i, 1] - x[j])^2 + (d[i, 2] - y[j])^2 < 5^2) {
                success = success + 1
            }
        }
    }
    result[loop] = success # 結果の記録
}
table(result) # 結果集計
print(mean(result)) # 実験 1 回あたりの，マシュマロが口に入る回数

マシュマロを10個投げる。下で 9 人が待ち受ける。マシュマロの着地点が，待ち受ける人口の半径5センチ以内なら，口に入れることができるとの条件で計算してみた。

ちなみに，マシュマロの着地点と 9 人の位置関係は，下図のような感じ。ここでは，マシュマロを 100 個投げたときの状況を示す。白抜きの丸がマシュマロ。赤丸が待ち受けている人。これだと，けっこう入りそうな気がする。

結果は，ほぼ 20% の成功確率。

マシュマロの数を倍の 20 個にすると当然（？？）倍くらい（40%）の成功確率となる。

判定条件を各人の口から 3 センチ以内とすると，成功確率は 6% ぐらいになる。

まあ，大まかでもいいけどシミュレーションしておかないと，「成功しました！！」とはなかなか行かないものだということだ。

Python って奴は，なんて自由な奴なんだ！

プログラミングの初歩の段階で，「1 以上，10 以下の場合は 'OK'，そうでない場合は 'NO' を書くプログラムを作りなさい」などというのがあり，初心者の中には，

if (1 < x < 10) print('OK')

などと書きたくなる人がいる。大抵の言語はこれはだめで，

if (1 < x && x < 10) print('OK')

としないといけない。

 

しかし，である。Python はこれを許す。なんですと!! ちっとも知らなかった。

def if_test(x):
    if 1 <= x <= 10:
        print('OK')
    else:
        print('NG')

if_test(0) # NG
if_test(5) # OK
if_test(10) # OK

ビットカウント（ビットが 1 である個数）

「ビットカウントする高速アルゴリズムをPythonで実装しながら詳しく解説してみる」というのを見つけた。

Python の整数は 64 ビット整数，R は 32 ビット整数であることを考慮し，以下のように翻訳した。

bitcount = function(n) {
  # nの立っているビット数をカウントする関数(nは64bit整数)
  # 2bitごとの組に分け、立っているビット数を2bitで表現する
  n = n - bitwAnd(bitwShiftR(n, 1), 0x55555555)
  # 4bit整数に 上位2bit + 下位2bit を計算した値を入れる
  n = bitwAnd(n, 0x33333333) + bitwAnd(bitwShiftR(n, 2), 0x33333333)
  n = bitwAnd((n + bitwShiftR(n, 4)), 0x0f0f0f0f) # 8bitごと
  n = n + bitwShiftR(n, 8) # 16bitごと
  n = n + bitwShiftR(n, 16) # 32bitごと = 全部の合計
  return(bitwAnd(n, 0x0000007f))
}

正しく移植できたか，確認してみる。

> bitcount(7)
[1] 3
> bitcount(127)
[1] 7
> bitcount(2^31-1)
[1] 31
> bitcount(-1)
[1] 32

R の intToBits() は，整数の二進表記を得るものである。これを用いてビット 1 を数えるのは自然である。 intToBits() のクラスは "raw" なので，as.integer() で整数化してから合計する。

bitcount2 = function(n) {
  return(sum(as.integer(intToBits(n))))
}

同じ結果になることが確認できる。

> bitcount2(7)
[1] 3
> bitcount2(127)
[1] 7
> bitcount2(2^31-1)
[1] 31
> bitcount2(-1)
[1] 32

では，両者の速度比較をしてみる。intToBits() がベクトルに対応していないので，心ならずも for ループで検証する。

> mx = 1000000000
> n = sample(mx, 1000000) # 100 万個の整数
> a = integer(mx)
> b = integer(mx)
> system.time({for (i in n) a = bitcount(i)})
   ユーザ   システム       経過  
     6.915      0.198      7.118
> system.time({for (i in n) b = bitcount2(i)})
   ユーザ   システム       経過  
     1.620      0.148      1.769
> all(a == b) # どちらの結果も同じか？
[1] TRUE

結論は，「intToBits() を使う方が 4 倍ほど速い」ということだが，100 万個の整数を処理して 2 〜 7 秒なので，1 個や 100 個なら，どちらの方法でも，ゴミのような差しかない。

======================================

ついでに，Python でも代替関数について調べてみた。

Python には二進表記をするための bin() がある。これは文字列として二進数を返すので，文字列の中に '1' がいくつあるか数える。

大元のビットカウントのプログラム

def bitcount(x):
    '''xの立っているビット数をカウントする関数
    (xは64bit整数)'''
    # 2bitごとの組に分け、立っているビット数を2bitで表現する
    x = x - ((x >> 1) & 0x5555555555555555)
    # 4bit整数に 上位2bit + 下位2bit を計算した値を入れる
    x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)
    x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f # 8bitごと
    x = x + (x >> 8) # 16bitごと
    x = x + (x >> 16) # 32bitごと
    x = x + (x >> 32) # 64bitごと = 全部の合計
    return x & 0x0000007f

対する，bin() を使う関数。短い。

def bitcount2(x):
    return sum([i == '1' for i in bin(x)])

100 万個の整数を対象とする。

import scipy as sp
a = sp.random.randint(-10, 100000000, 1000000)


def func1(): # 4.383 sec.
    x = [bitcount(i) for i in a]

def func2(): # 3.416 sec.
    y = [bitcount2(i) for i in a]

やはり，bin() を使う方が速い。
sp.all(x == y) は True である。

utils:::combn の代替関数 next.combn

utils:::combn は組み合わせを配列で返すが，場合によっては一つずつ返してくれる方がうれしいこともある。

そこで，以下に示す next.combn を書いてみた。

next_combn = function(n, r, a) {
  t = r
  while (t >= 1 && a[t] == n - r + t) {
    t = t - 1
  }
  if (t ==  0) {
    return(FALSE)
  }
  a[t] = a[t] + 1
  for (i in t:r) {
    a[i] = a[t] + i - t
  }
  return(a)
}

使い方は，n, r を指定して a を 1,2, ..., n のベクトルに定義して

n = 5
r = 3
a = 1:n

以下のように引用する。print(a[1:r]) の所に，本来意図している処理プログラムを書く。

while (TRUE) {
  print(a[1:r])
  a = next_combn(n, r, a)
  if (length(a) == 1) {
    break
  }
}

 

e1071:::permutations の代替関数 next.permutation

e1071 ライブラリには permutations 関数がある。

permutations は処理が速いが，配列で返すために大きなメモリーを消費する。

時間はかかってもよいので，メモリー消費のない関数が必要なこともある。そのようなときのために以下の関数を書いた。

next.permutation = function(x) {
  n = length(x)
  pos1 = 0
  for (i in (n - 1):1) {
    if (x[i] < x[i + 1]) {
      pos1 = i
      break
    }
  }
  if (pos1 == 0) {
    return(0)
  }
  for (i in n:(pos1 + 1)) {
    if (x[pos1] < x[i]) {
      pos2 = i
      break
    }
  }
  temp = x[pos1]
  x[pos1] = x[pos2]
  x[pos2] = temp
  len = (n - pos1) %/% 2
  if (len >= 1) {
    for (i in 1:len) {
      temp = x[pos1 + i]
      x[pos1 + i] = x[n - i + 1]
      x[n - i + 1] = temp
    }
  }
  return(x)
}

使用法は以下のように，最初に x に 1, 2, 3, ..., n をセットして，「x を使った処理」の次に必要な処理をプログラムする。

単に print(x) とするならば，可能な順列を全て列挙してくれる。

x = 1:4
while (TRUE) {
  # x を使った処理
  print(x)
  x = next.permutation(x)
  if (length(x) == 1) {
    break
  }
}

> x = 1:4
> while (TRUE) {
+   # x を使った処理
+   print(x)
+   x = next.permutation(x)
+   if (length(x) == 1) {
+     break
+   }
+ }
[1] 1 2 3 4
[1] 1 2 4 3
[1] 1 3 2 4
[1] 1 3 4 2
　　　:
[1] 4 2 1 3
[1] 4 2 3 1
[1] 4 3 1 2
[1] 4 3 2 1