不定方程式(ディオファントス方程式)

YouTubeのおすすめ動画でたまたま「整数問題」として以下の問題が出てきた。

a^2 + b^2 + c^2 = 292を満たす自然数(a, b, c)を求めよ

おとなになったら、ノートにガリガリ書いて解いたりしないよ。
という訳で、Python(Sympy)で考えると、
これは、不定方程式なので、sympy.solveでは解けないが、
ディオファントス方程式なので、ソルバを使えば解ける。

import sympy as sym
from sympy.abc import a, b, c, n
from sympy.solvers.diophantine import diophantine


f = sym.Eq(a**2 + b**2 + c**2, n)
diophantine(f.subs(n, 292))
{(0, 6, 16), (2, 12, 12)}

この組み合わせになるので、9通りの答えになる。

カテゴリー: 未分類 | コメントをどうぞ

悪夢

頭が禿げ上がる夢をみて、非常に気分が悪い。しかも、なんかすごい汚いハゲ散らかし方をしているという。夢というのは大きく「記憶のエンコーディング(いわゆる夢)」と浅い睡眠時の「漠然とした考え(多くが夢と思っている方)」があるが、後者は思い煩いが大きく影響するので、週末の一件に関して(自分は直接関係していないけど)何かしらストレスを感じているのかなあ。

まあ、僕があれこれ憂慮しても仕方がないんだけど。

カテゴリー: 未分類 | コメントをどうぞ

一人10万は少ないし非課税にするのは愚か

経済支援よりコロナ対策として給付すべきで、気兼ねなく在宅出来る様、という意味を考えないといけない。十万円やるから仕事に行くなと言われて、それで満足する人間なんてごくしょうすうだろう。

また、非課税にするのも馬鹿げている。一律給付で迅速に行い、かつ税システムに基づく再分配によって、所得水準に応じた合理化が出来る。

何もかも、後手後手な上、悪手で、明らかに財務省の考えが透けて見え、日本の経済対策で最も重要なのは財務省解体だろって思う。

カテゴリー: 未分類 | コメントをどうぞ

シェルピンスキーのギャスケット

家のデスクトップで久しぶりにRStudioを立ち上げたら、
ワークスペースに面白いものが残っていたので改めて。
(RのコードはどこかのWebサイトを参考にしたもの)

gasket <- function(n=100000) {
 t <- proc.time()
 m <- matrix(c(0.5, 0, 0, 0.5) , 2)
 lst <- matrix(0, n, 2) # matrixを用意
 lst[1,] <- c(100, 0)   # 初期値

 prb <- sample(1:3, n,rep=TRUE)

 for(i in 2:n){
   lst[i,] <- switch(prb[i], # switchで振り分け
                     m %*% lst[i-1,],
                     m %*% lst[i-1,] + c(100,0),
                     m %*% lst[i-1,] + c(50, 30))
 }
 plot(lst, pch=".", axes=F, asp=sqrt(3)*1.5)
 timeval <- proc.time() - t
 print(timeval)
}

## 10年以上前のPC(多分Duronだった気がする) 
# > gasket()
#  user  system elapsed
#  0.94    1.07    2.29
## 今年買ったノートPC(Core i7 10th)
# > gasket()
#   user  system elapsed 
#   0.27    0.50    0.79  
# %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import time
import statistics


proc_time_l = []
for _ in range(100):
#856 ms ± 8.06 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
    s = time.perf_counter()
    n=100000
    m = np.array(
       [[0.5, 0],
       [0, 0.5]]
    )
    lst = np.zeros((n, 2))
    lst[0] = [100, 0]
    prb = np.random.randint(3, size=n)

    for i in range(1, n):
       if prb[i] == 0:
           lst[i] = np.dot(m, lst[i-1])
       elif prb[i] == 1:
           lst[i] = np.dot(m, lst[i-1]) + [100, 0]
       else:
           lst[i] = np.dot(m, lst[i-1]) + [50, 30]

    fig = plt.figure(figsize=(13, 10))
    plt.plot(lst[:, 0], lst[:, 1], 'k.')
    plt.axis('off')
    proc_time_l.append(time.perf_counter()-s)
print(statistics.mean(proc_time_l))
print(statistics.stdev(proc_time_l))
0.49452882845999396
0.014653359037954079

ちょっとだけ改良。

%%timeit
import numpy as np
import matplotlib.pyplot as plt


n=100000
m = np.array(
   [[0.5, 0],
   [0, 0.5]]
)
out = np.empty((n, 2))
a = np.array([[0, 0], [100, 0], [50, 30]])
r = np.random.randint(3, size=n-1)
out[0] = [100, 0]
out[1:] = a[r]

for i in range(1, n):
    out[i] = out[i] + m.dot(out[i-1])

fig = plt.figure(figsize=(13, 10))
plt.plot(out[:, 0], out[:, 1], 'k.')
plt.axis('off')
1 loop, best of 3: 240 ms per loop

行列積の(線形じゃ無い)部分が再帰構造的処理になっているので、ループは無くせない(と思う)。
単純なループのみが残っているので、NumbaやCythonが効果的だろう。

import numpy as np
import matplotlib.pyplot as plt
from numba import jit


@jit(nopython=True, parallel=True, fastmath=True)
def f(out, n):
    m = np.array([[0.5, 0], [0, 0.5]])
    for i in range(1, n):
        out[i] = out[i] + m.dot(out[i-1])
%%timeit
n=100000
out = np.empty((n, 2))
a = np.array([[0, 0], [100, 0], [50, 30]])
r = np.random.randint(3, size=n-1)
out[0] = [100, 0]
out[1:] = a[r]
f(out, n)
fig = plt.figure(figsize=(13, 10))
plt.plot(out[:, 0], out[:, 1], 'k.')
plt.axis('off')
1 loop, best of 3: 139 ms per loop
%%cython
import cython
import numpy as np
import matplotlib.pyplot as plt
cimport numpy as np


@cython.ccall
@cython.locals(out=np.ndarray, m=np.ndarray, i=cython.uint, n=cython.uint)
def f(out, n):
    m = np.array([[0.5, 0], [0, 0.5]])
    for i in range(1, n):
        out[i] = out[i] + m.dot(out[i-1])

@cython.ccall
@cython.locals(out=np.ndarray, a=np.ndarray, r=np.ndarray, n=cython.uint)
def main():
    n=100000
    out = np.empty((n, 2))
    a = np.array([[0, 0], [100, 0], [50, 30]])
    r = np.random.randint(3, size=n-1)
    out[0] = [100, 0]
    out[1:] = a[r]
    f(out, n)
    fig = plt.figure(figsize=(13, 10))
    plt.plot(out[:, 0], out[:, 1], 'k.')
    plt.axis('off')
%%timeit
main()
1 loop, best of 3: 203 ms per loop

ほぼCの様な書き方にしたらもっと速くなるが、Pythonのメリットを失うので、スケーラビリティを考えれば、ループ部分の高速化だけであればNumbaがリーズナブルチョイス。

カテゴリー: 未分類 | コメントをどうぞ

ビットコインはまだ下げるか

当初の予想通り、5,000前後でサポートされてからは行ったり来たりだが、再度下げてもそこを割らないだろうと考えている。基本的には、ダウとの相関を意識し、ダウがここから大きく崩れる様な事にならない限りは、ビットコインも5,500から上だろうと。

4月まではまだまだ不安定とみるが、そろそろ仕込み時で、次の押し目はドルやダウ、日経の買い場だと思っている。

カテゴリー: 未分類 | コメントをどうぞ

今は株のバーゲンセール……と考えるのは早計では

僕は楽観論者なので、戒めも兼ねて悲観論を考える。

多額のマージンコールでドルクランチに陥り、金利が上がっている。今は信用収縮の話に過ぎないが、このままでは借り入れコストの上昇から企業の資金繰り悪化が懸念される。

既にゼロ金利以下の国々では、金融政策の実効性は限りなく低く、(経済学の定石通りに)財政政策しか無い。問題はその金額で、効果的なな額を投入出来無ければ、ただのサンクコストになる。米国の場合、およそ2兆ドルが必要になると推計されているし、日本の場合も50-100兆円は必要になるだろう。たかだか10兆20兆では、お金を無駄にするだけだ。協調的にやれれば、シナジーが得られれば必要な予算は減らせるだろうが。

仮に適切に対処されるならば、今は間違いなく底値で、バーゲンセール会場かも知れないが、そうでない場合、二番底をつけにいく事になる。つまり、今よりももっと下がる。ダウは16,000ドルまで、日経は8,800円がみえてくる事になる。

カテゴリー: 未分類 | コメントをどうぞ

12兆ドルのマージンコールとか笑うしか無い

12兆ドルって凄いなあ、としか。リーマン・ショック時の2倍の規模で、通常時の外為より遥かに大きい規模でドル需要が生じて、完全に需給が逼迫している。それをもってどれくらいの圧力が生まれるものなのか分からないけど。リーマン・ショックの時、ドル需要で10円位円安に振れたと考えると、その倍位振れてもおかしくないのかなと。まあ、元々120円位までは振れる事もあると考えているので、101円から20円振れて121円だと、割とこれまでの考えとも整合性があるので、4月までにその水準まで動く事は十分ありえるかなって思っている。

カテゴリー: 未分類 | コメントをどうぞ

Fusion 360使いこなせる様になりたい

多分、使いこなせている人なら数分から数十分で作る事ができるのかなっていう造形物を、僕は独学で必要に迫られながらやっているだけなので、数十分、下手をしたら数時間掛けて作っている気がする。3Dプリンタと合わせて、色々と作ってみたいものはあるし、そのイメージは持っているんだけど、それをアウトプットする能力に欠けるのは歯がゆい。

ちょっと本屋に行って、良さそうな本でも買って真剣に学びたいなあ。暇になれば。

カテゴリー: 未分類 | コメントをどうぞ

右足が痛い

内側とうじょう骨の辺りが痛い。歩くのも結構辛いレベル。23日前からだが、何が原因やら……。

カテゴリー: 未分類 | コメントをどうぞ

ダウ平均とビットコイン

既に買い場だと考えるが、仮に下げるとすると、ダウは16,000ドルまで、ビットコインは3,500-3,800ドルまで下げる計算になる。仮にその時、日経平均はどうなっているか。15,000円をそうそう割らないと考えているが、これを割ると過去の安値8,000円アラウンドがみえてくる。過去のケースでもそうだった様に、いきなりそこまで下げるとは思わないが。そろそろ持ち直し、或いはズルズルとくすぶり続け、コロナ収束の兆しがみえないと判断された所でドカン、或いは幸運にも収束しリスクオンと考えている。株よりもドル円の方が面白そうで、既に日経平均は日銀の損失ラインを超えている蓋然性が高く、この水準に留まれば意識されるだろう。日銀が損失を抱えようと大した問題では無いが、それが気にされる事に意味があって、そうなれば円安動意が発生し得る。株安、円安の動意となれば、円高以上にコントロールし難いのではないか。とはいっても、減少傾向にあるとは云え、日本は未だ世界一の純債権国なので、コントロール不能な円安に陥る蓋然性は低いが。せいぜい120円オーバーで止まると思っている。

カテゴリー: 未分類 | コメントをどうぞ