　　　　On a Thread of the Web: 計算

ラベル計算の投稿を表示しています。すべての投稿を表示

2019年11月4日月曜日

反対称テンソルと対称和

ある量の対称和が必要な場合がある。これを，レヴィ=チヴィタ記号$\varepsilon_{ijk}$で表現する例をいくつか考えてみた。$\varepsilon_{ijk}^2 = (1-\delta_{ij}) (1-\delta_{jk}) (1-\delta_{ki})$が成り立つような気がする。
\begin{equation}
\begin{aligned}
A_1 B_2 C_3 + A_2 B_3 C_1 + A_3 B_1 C_2 &= \dfrac{1}{2} \sum_{ijk} (\varepsilon_{ijk} + \varepsilon_{ijk}^2) A_i B_j C_k\\
A_1 B_2 B_3 + A_2 B_3 B_1 + A_3 B_1 B_2 &= \dfrac{1}{2} \sum_{ijk} \varepsilon_{ijk}^2 A_i B_j B_k\\
A_1 + A_2 + A_3 &= \dfrac{1}{2} \sum_{ijk}\varepsilon_{ijk}^2 A_i
\end{aligned}
\end{equation}

2019年9月19日木曜日

平方ピラミッド問題

京都の烏丸御池にあるNHK京都支局の１Ｆには８Ｋビジョンが設置されていて，パブリックビューイングができるが，その横には世界一大きな（？）ドーモ君が鎮座している。ちょうどお月見の季節なので，１＋４＋９個の月見団子が設置されていた。

このようにピラミッド状にお団子を$\ m \ $段積んでいくとその総数$\ \sum_{k=1}^m k^2\ $が平方数$\ n ^2 \ $になることがある。それは，$\ m=24，n=70 \ $に限られることが証明されている。

石井夕紀子さんによる初等的な証明の解説がある。

2019年9月4日水曜日

ネイピア数と素数

（天才数学少女？からの続き）

ネイピア数を10進法で表現し，小数点以下650桁までとったものから2.7を減じたものを整数とみなした数が素数になるという話について。なぜ，最初の２桁を除いたものなのだろうか。その他ではどうなるのだろうか。そこで，上の桁から順に取り除いて整数化したもののうちはじめて素数が現れる小数点以下桁数を探してみた。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
N[E, 20]
2.7182818284590452354

g1[k_] := PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 2*10^(k + 1))]]
g2[k_] := PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 27*10^k)]]
g3[k_] := PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 271*10^(k - 1))]]
g4[k_] := PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 2718*10^(k - 2))]]
g5[k_] := PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 27182*10^(k - 3))]]
g6[k_] := PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 271828*10^(k - 4))]]
g7[k_] := PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 2718281*10^(k - 5))]]
g8[k_] := PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 27182818*10^(k - 6))]]
g9[k_] :=
PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 271828182*10^(k - 7))]]
g10[k_] :=
PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 2718281828*10^(k - 8))]]
g11[k_] :=
PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 27182818284*10^(k - 9))]]
g12[k_] :=
PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 271828182845*10^(k - 10))]]
g13[k_] :=
PrimeQ[Floor[(N[E, k + 3]*10^(k + 1) - 2718281828459*10^(k - 11))]]

{g1[1], g2[649], g3[26], g4[3], g5[25], g6[34], g7[18], g8[7], g9[273], g10[44], g11[10], g12[16], g13[16]}

{True, True, True, True, True, True, True, True, True, True, True, True, True}
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

最後の数字はgi[k]のテーブルを目の子で作って探した結果だけを書いたものであり，２桁の2.7を差し引いた場合が，他の場合に比べて際立っていることがわかる。が，それほど意味はない。

【付録】なお，素数定理とネイピア数には関係があるらしい。
f[n_] := Apply[LCM, Range[n]]^(1.0/n)
a = Table[f[i], {i, 1, 10000}]; b=Table[f[i],{i,1,100000}];
ListPlot[a]
ListPlot[b]

図　nまでの数の最小公倍数のn乗根の収束性

2019年8月22日木曜日

multiplicative persistence

しばらく前，TwitterでUniversity of WashingtonのKeiko Toriiさんの12歳の娘さん（gifted class）の算数の自由課題が話題になっていた。任意の自然数の各桁の積を計算して新たな自然数を求める。答えの自然数に対してさらにこの計算を反復しする。これを繰り返すと最後には１桁の自然数が得られる。１桁の自然数に到達するまでの反復回数がなるべく大きくなるもとの自然数を探せというのが賞品付きの課題だ。8th grade にも関らず，pythonでプログラムを組んで反復数11の数を見つけたとのこと。

飛び級小学生（＝中学生相当）には負けてられませんと，Brute Forceで計算してみた。残念ながら時間ばかりかかって反復数が10にしかとどかなかず老人は小学生に完敗だ。調べてみると，Persistence of a NumberというWikipediaの項目が見つかった。日本語版はない。自然数に一定のアルゴリズムでの操作を行って不変になるものを探すという問題群らしい。OEIS（On Line Encyclopedia of Integer Sequences ）A003001 にもある。

さらに調べると，Rubyで効率的な探索をしている人がいたので，この方法を Juliaにあてはめてみた。考え方としては，２と３と７のベキからなる数の集合から，反復数が最も大きいものを探し，その数を構成する因子の２と３と７の数字が並んだ数が，反復数＋１の数を与えるというものだ。反復回数11回の数が２つ見つかった。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function pcount(n::Int128)
local m::Int128
println(n)
c = 0
while n >= 10
m = 1
while n>0 && m>0
m *= n%10
n ÷= 10
end
n = m
c += 1
println(n)
end
end

function count(n::Int128)
local m::Int128
c = 0
while n >= 10
m = 1
while n>0 && m>0
m *= n%10
n ÷= 10
end
n = m
c += 1
end
return c
end

function search(n2,n3,n7)
local n,m3,m7::Int128
d = (0,0,0,0,0)
max = 0
for k in 1:n7
m7 = 7^k
for j in 1:n3
m3 = 3^j
for i in 1:n2
n = 2^i*m3*m7
c = count(n)
if c > max
d = (i,j,k,c,n)
max = c
end
end
end
end
println(d)
return d
end

function next(p)
s="23789"
m=[0,0,0,0,0]
m[1]=p[1]%3
m[2]=p[2]%2
m[3]=p[3]
m[4]=p[1]÷3
m[5]=p[2]÷2
c=""
for i in 1:5
for j in 1:m[i]
c=c*s[i]
end
end
return parse(Int128,c)
end

a=search(20,10,10)
pcount(next(a))
b=search(10,20,10)
pcount(next(b))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(19, 4, 6, 10, 4996238671872)
277777788888899
4996238671872
438939648
4478976
338688
27648
2688
768
336
54
20
0
(4, 20, 5, 10, 937638166841712)
27777789999999999
937638166841712
438939648
4478976
338688
27648
2688
768
336
54
20
0

2019年4月23日火曜日

5! * 6! = 10!

Fermer's Libraryの先日のツイートを受けて，tsujimotterさんが階乗数の間の関係式というブログ記事を書いている。そこでは，rubyのコードで n=200以下で n!=a! b! となる自明でない，(n, a, b) の組を計算している。自明な組とは， (n, 0, n)， (n, 1, n)，(n!, n!-1, n) などのことである。これを juliaで計算してみた。

rubyのコードがわかりにくいと思った原因は，階乗の計算のオーバーフローを避けるために，n! を素因数分解した素数の積で表現していたからだったか。以下のjuliaプログラムではBigIntを使ってこの課題をスキップしてしまった。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f=ones(BigInt,1000)

function fct(n,f)
for i in 2:n
f[i]=BigInt(i)*f[i-1]
end
end

function bin(n,f)
for k in 2:n
for i in 2:n-1, j in i:n
if(f[k]==f[i]*f[j])
println(i,"! * ",j,"! = ",k,"!")
end
end
end
end

n=1000
fct(n,f)
@time bin(n,f)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

3! * 5! = 6!
6! * 7! = 10!
4! * 23! = 24!
5! * 119! = 120!
6! * 719! = 720!
2017.072241 seconds (1.50 G allocations: 476.873 GiB, 2.13% gc time)

2019年3月17日日曜日

素数日

ほぼ毎日（多少前後して編集しているけれど）ブログを書き始めて，今日で100記事になった。いつまで続くのだろう？何のために続けているのだろう？惰性なのだろうか？

西暦（グレゴリオ暦）のyyyy年mm月dd日を8桁の整数 yyyymmdd で表したものが素数であるものを素数日と言い習わしているようだ。あちらこちらで目にします。で，Juliaで素数日を求める関数を書いてみた。あらかじめ，Primesパッケージを読み出すので，自分の仕事はほとんどない。

using Pkg
Pkg.add("Primes")
using Primes
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function primeday(y)
u=28
if(y%4==0 && y%100!=0 || y%400==0)
u=29
end
# println(u)
mx=(31,u,31,30,31,30,31,31,30,31,30,31)
for m in 1:12, d in 1:mx[m]
if(isprime(y*10000+m*100+d))
println(y*10000+m*100+d)
end
end
end

primeday(2019)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20190221
20190227
20190301
20190319
20190323
20190421
20190523
20190529
20190601
20190613
20190719
20190811
20190823
20190913
20191009
20191027
20191109
20191117
20191231

2019年3月15日金曜日

円周率の分数近似

円周率の日の余震が続いている。3月14日，Googleは円周率の計算の世界記録を更新したことを発表したとの報道があった。31.4兆桁（31,415,926,535,897桁）までの計算を111.8日（start to end で121.1日）かけて計算したそうだ。このプロジェクトを主導したのは，筑波大学出身の技術者 Emma Haruka Iwao （岩尾・エマ・はるか）さん。今年の1月には終わっていたが，このよき日に発表された。つぎに記録に挑戦するグループは，314兆桁を越えないといけないのか。

さて，だいぶ昔（2011.6.30）のWolfram Blogで，Jon McLooneが，"All Rational Approximations of Pi Are Useless"という記事を書いている。

どういうことかというと・・・円周率のよくある近似値22/7を例にとってみる。22/7 = 3.1426... は小数点以下2桁まで正しいので，小数点も含め，3 . 1 4 の4文字が表現されたことになる。このために必要な数は 2 2 / 7 の4文字である。これじゃ文字数が節約できていない。どちらを覚えてもいっしょじゃないか，ということらしい。

もっとやってみると，355/113 = 3.141592... これは小数点以下6桁すなわち8文字まで正しい円周率に対応している。一方，355/113は7文字，さきほどよりちょっとだけましだ。

McLooneの記事にあるMathematicaの関数を用いて，小数点以下50桁のπとその分数表示を求めると，

\begin{equation}
\frac{26151465932107044561886949}{8324270144388272579650158} = \\
　\\
3.14159265358979323846264338327950288419716939937510...
\end{equation}
ということで共に52文字必要としている。でこぼこはあるが，おおむねこの比が1で推移することから，分数表示意味ないじゃんということになったようだ。

必要なMathematica関数は以下のようになっている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
In[1]:= characterCount[expr_] := StringLength[
StringReplace[ToString[Unevaluated[expr], InputForm],
" " -> ""]]; SetAttribute[characterCount, HoldAll];

In[2]:= matchingCharacters[expr_, target_] :=

Ceiling[-Log10[Abs[expr - target]]] + 1

In[3]:= $MaxExtraPrecision = Infinity;

In[4]:= candidates = Table[Rationalize[Pi, 10^(-n)], {n, 1001}];

In[5]:= a = Table[candidates[[k]], {k, 1, 10}]

Out[5]= {22/7, 22/7, 201/64, 333/106, 355/113, 355/113,
75948/24175, 100798/32085, 103993/33102, 312689/99532}

In[6]:= b =
Table[characterCount[candidates[[k]]], {k, 1, 1001, 50}]

Out[6]= {4, 54, 103, 154, 203, 253, 303, 353, 404, 455, 504,
553, 603, 652, 704, 753, 803, 853, 904, 953, 1003}

In[7]:= c =
Table[matchingCharacters[candidates[[k]], Pi], {k, 1, 1001, 50}]

Out[7]= {4, 53, 103, 153, 204, 253, 303, 353, 404, 453, 503,
553, 604, 653, 704, 753, 803, 853, 903, 953, 1003}
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

2019年3月6日水曜日

ネイピア数の計算

（円周率の計算からの続き）

ネイピア数は自然対数の底 $e$ のこと。欧米では Euler's Number とよばれることが多いとあるが，Eulerが関わる数はたくさんあってややこしいので，ここではネイピア数とする。

円周率の計算方法というか表現方法はたくさんあって，その中には与えられた級数の少数の項で良い精度の近似を与えるものがあった。ところで，ネイピア数の表式をざっと調べても，なんだかぱっとしないのである。
\begin{equation}
\begin{aligned}
e^x &= 1 + x + \frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\cdots\\
e &= 1 + 1 + \frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\cdots
\end{aligned}
\end{equation}
を多少変形したものか，連分数で表現したものはあるけれど，ラマヌジャンの公式のようにはっとするものが見当たらない。単に調査不足だけかもしれないし，そもそも指数関数を微分や積分しても元に戻るだけなのが原因かもしれないが，よくわからないのだった。

５項くらいで精度が70桁程度になる公式はないものなのだろうか？Improving the Convergence of Newton's Series Approximation for $e$ とかあったけど，目的のものではない。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f=ones(BigFloat,1000)

function fact(f,n)
# calculation of factorial 1!..n!
for i in 2:n
f[i]=BigFloat(i)*f[i-1]
end
# println(f[n])
end

function cf(x,a,b,n)
# continued fraction
for k = n:-1:1
x = a[k] + b[k]/x
end
return x
end

function napier1(n)
# by series expansion
o=BigFloat(1)
sum=o
for k in 1:n
sum=sum+o/f[k]
end
println(sum-exp(o))
end

function napier2(n)
# by series expansion
o=BigFloat(1)
sum=o
for k in 1:n
sum=sum+o/((BigFloat(4)^k)*f[k])
end
println(sum*sum*sum*sum-exp(o))
end

function napier3(n)
# by continued fraction
o=BigFloat(1)
z=BigFloat(2)
a=ones(BigFloat,1000)
for k in 1:1000
a[k]=z*(z*k-z-o)
end
a[1]=o
a[2]=o
b=ones(BigFloat,1000)
b[1]=z
println(cf(z,a,b,n)-exp(o))
end

function napier4(n)
# sqrt(e) by series expansion
o=BigFloat(1)
z=BigFloat(2)
sum=BigFloat(0)
for k in 0:n
sum=sum+BigFloat(4*k+3)/(z^(2*k+1)*f[2*k+1])
end
println(sum*sum-exp(o))
end

fact(f,160)
@time napier1(57)
@time napier2(41)
@time napier3(26)
@time napier4(24)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
0.018460 seconds (30.73 k allocations: 1.488 MiB)
-2.763573937630222280123632596096127862757120116619610043207585110453949377011977e-76
0.029598 seconds (49.70 k allocations: 2.358 MiB)
-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
0.054272 seconds (100.99 k allocations: 5.090 MiB)
1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-76
0.032685 seconds (72.55 k allocations: 3.431 MiB)

（Juliaの数学的定数と精度へ続く）

2019年3月5日火曜日

円周率の計算

（電子計算機と人間からの続き）

現代科学シリーズの「電子計算機と人間」では，円周率の計算のためのFORTRANプログラムが書かれていたが，高校１年生の自分にはなぜそれで円周率が求まるのかがさっぱり分からなかった。50年ぶりに，図書館でその本を手にとると，たしかに分かりやすくはなかったが，アークタンジェントのマクローリン展開に１を代入した単純な公式なのであった。
\begin{equation}
\begin{aligned}
\arctan(x) & = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \cdots \\
\frac{\pi}{4} & = 1 - \frac{1}{3!} + \frac{1}{5!} - \frac{1}{7!} + \cdots
\end{aligned}
\end{equation}
この公式は収束がとても悪い見本のような式なので，1000項でようやく3.14が0.1%精度で求まる始末である。

そこで，いくつかの代表的な公式をJuliaプログラムに落としてみた。BigFloatの精度がeps(BigFloat) = $10^{-77}$程度なので，有効数字は77桁弱である。そのBigFloatの使い方が，いまひとつわかっておらず適当に入れているが，たぶん冗長だと思われる。

円周率の計算で，有効数字77桁程度の精度を出すには，昔よくみかけたマチンの公式で 54項，なんだかすごいのだけれど今回初めて使ったラマヌジャンの公式で９項，スマートなガウス＝ルジャンドルの方法で５項が必要であった。もっと楽しそうなJuliaでのπの話はこちら。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f=ones(BigFloat,1000)

function fact(f,n)
# calculation of factorial 1!..n!
for i in 2:n
f[i]=BigFloat(i)*f[i-1]
end
# println(f[n])
end

function pie1(n)
# Gregory-Leibniz Formula
# pi=4*arctan(1)
sum=BigFloat(0)
for i in 1:n
sum=sum+(-1)^(i-1)/BigFloat(2*i-1)
end
println(4*sum-pi)
end

function pie2(n)
# Machin's Formula
# pi = 16*arctan(1/5)-4*arctan(1/239)
sum=BigFloat(0)
sun=BigFloat(0)
for i in 1:n
sum=sum+(-1)^(i-1)/((2*i-1)*BigFloat(5)^(2*i-1))
sun=sun+(-1)^(i-1)/((2*i-1)*BigFloat(239)^(2*i-1))
end
println(16*sum-4*sun-pi)
end

function pie3(n)
# Takano-Kikuo's Formula
# pi = 48*arctan(1/49)+128*arctan(1/57)-30*arctan(1/239)+48*arctan(1/110443)
sum1=BigFloat(0)
sum2=BigFloat(0)
sum3=BigFloat(0)
sum4=BigFloat(0)
for i in 1:n
sum1=sum1+(-1)^(i-1)/((2*i-1)*BigFloat(49)^(2*i-1))
sum2=sum2+(-1)^(i-1)/((2*i-1)*BigFloat(57)^(2*i-1))
sum3=sum3+(-1)^(i-1)/((2*i-1)*BigFloat(239)^(2*i-1))
sum4=sum4+(-1)^(i-1)/((2*i-1)*BigFloat(110443)^(2*i-1))
end
println(48*sum1+128*sum2-20*sum3+48*sum4-pi)
end

function pie4(n)
# Ramanujan Formula
sum=BigFloat(1103)
for k in 1:n
sum=sum+f[4*k]/(f[k]^4*BigFloat(396)^BigFloat(4*k))*BigFloat(1103+26390*k)
end
println(BigFloat(9801)/(sqrt(BigFloat(8))*sum)-pi)
end

function pie5(n)
# Chudnovsky formula
sum=BigFloat(13591409)
for k in 1:n
sum=sum+f[6*k]/(f[3*k]*f[k]^3*BigFloat(-640320)^BigFloat(3*k))*BigFloat(13591409+545140134*k)
end
println(BigFloat(640320)^1.5/(BigFloat(12)*sum)-pi)
end

function pie6(n)
# Gauss-Legendre algorithm
z=BigFloat(2)
a=BigFloat(1)
b=1/sqrt(z)
t=1/(z*z)
p=BigFloat(1)
for k in 1:n
o=a
a=(o+b)/z
b=sqrt(o*b)
t=t-p*(a-o)*(a-o)
p=z*p
end
println((a+b)*(a+b)/(4*t)-pi)
end

fact(f,160)
@time pie1(10000000)
@time pie2(54)
@time pie3(22)
@time pie4(9)
@time pie5(5)
@time pie6(5)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-9.999999999999975000000000000312499999999990468750000000541015625007904590936847e-08
5.580711 seconds (60.04 M allocations: 3.131 GiB, 4.69% gc time)
3.454467422037777850154540745120159828446400145774512554009481388067436721264971e-77
0.055656 seconds (76.18 k allocations: 3.601 MiB, 7.32% gc time)
-6.908934844075555700309081490240319656892800291549025108018962776134873442529942e-77
0.071662 seconds (123.15 k allocations: 5.740 MiB)
-6.908934844075555700309081490240319656892800291549025108018962776134873442529942e-77
0.057192 seconds (93.11 k allocations: 4.357 MiB)
-3.454467422037777850154540745120159828446400145774512554009481388067436721264971e-77
0.068310 seconds (106.84 k allocations: 5.084 MiB)
0.0
0.036871 seconds (53.77 k allocations: 2.544 MiB)

（ネイピア数の計算に続く）

2019年1月29日火曜日

Juliaでパズル（５）

（Juliaでパズル（４）から続く）

結論：失敗でした。７時間以上かかるプログラムを高速化するために，重複計算部分を配列に落とし込もうとしたが，プログラムの見通しが非常に悪くなって，しかも時間も短縮できないという残念な結果に終わってしまった。

30年くらい前，スーパーコンピュータを使って原子核の殻模型計算を高速化するための試みについて，大阪大学の大型計算機センターニュースに投稿したことがある（スーパーコンピュータと原子核殻模型計算 1987）。導入されて間もなかったNECのベクトル型計算機の最初のモデルのSX-1を使い，$sd$殻の11粒子系である${}^{29}{\rm Si}$の基底状態（$J=\frac{1}{2}, T=\frac{1}{2}$，2360次元）を求める計算を高速化しようというものである（西村道明さんとの共同研究のための予備計算，大坪久夫先生にRCNPの計算費を確保していただいた）。

結果は，ベクトル化よりも重複計算を間接アドレッシングやメモリに蓄えて呼び出すことで避ける効果が高いというなさけないものであったが，その時から重複計算を避けようというのがスローガンとしてたたき込まれていたのであった（しかし，ベクトルマシンそのものが陳腐化して消えてしまうわけで，世の中の栄枯盛衰，諸行無常なり）。

今回の結果は，そもそも重複計算がそんなにあったのかというところからして問題であり，コードが複雑になって保守性にも欠けるという結末で終わってしまった。もともと重複だと思っている関数呼び出しも大したコストがかかっておらず，多次元の配列呼び出しとあまりかわらなかったのかもしれない。とりあえず，悪戦苦闘した（残念な）プログラムは以下の通り

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function sqr(x)
return sqrt(abs(x))
end

function cbr(x)
return cbrt(abs(x))
end

function rog(x)
return log(abs(x))
end

function cc(a,b)
if (a==1 && b==1) || (a==11 && b==1)
return 10*a+b
else
return 0
end
end

function ne(op, op1)
expr = Expr(:call, op, op1)
return expr
end

function me(op, op1, op2)
expr = Expr(:call, op, op1, op2)
return expr
end

op1=[:+,:sqr,:cbr,:exp,:rog]
op2=[:+,:-,:*,:/,:cc]
ne1=Array{Any}(undef, 5)
ne2=Array{Any}(undef, 5, 5, 5)
ne11=Array{Any}(undef, 5, 5, 5, 5)
nex1=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5)
ne1x=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5)
mexy11=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
mey1x1=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
mey11x=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
me1yx1=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
me1y1x=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)

function stor(ne1,ne11,ne1x,nex1,mexy11,mey1x1,mey11x,me1yx1,me1y1x,op1,op2)
for i in 1:5
ne1[i]=Expr(:call,op1[i],1)
end
for i in 1:5, j in 1:5, k in 1:5
ne2[i,j,k]=Expr(:call,op2[i],ne1[j],ne1[k])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5
ne11[i,j,k,l]=Expr(:call,op1[i],ne2[j,k,l])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5
nex1[i,j,k,l,m,n,o]=Expr(:call,op1[i],Expr(:call,op2[j],ne11[k,l,m,n],ne1[o]))
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5
ne1x[i,j,k,l,m,n,o]=Expr(:call,op1[i],Expr(:call,op2[j],ne1[k],ne11[l,m,n,o]))
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
mexy11[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne11[j,k,l,m],ne11[n,o,p,q])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
mey1x1[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],nex1[j,k,l,m,n,o,p],ne1[q])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
mey11x[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne1x[j,k,l,m,n,o,p],ne1[q])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
me1yx1[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne1[j],nex1[k,l,m,n,o,p,q])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
me1y1x[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne1[j],ne1x[k,l,m,n,o,p,q])
end
end

@time stor(ne1,ne11,ne1x,nex1,mexy11,mey1x1,mey11x,me1yx1,me1y1x,op1,op2)

function doe1(a,b,c,d)
# ((1,1)_a,(1,1)_c)_b
# x=ne11[d[5],a,d[1],d[2]]
# y=ne11[d[6],c,d[3],d[4]]
z=ne(op1[d[7]],mexy11[b,d[5],a,d[1],d[2],d[6],c,d[3],d[4]])
if b==5
z=:0
end
return (z,eval(z))
end

function doe2(a,b,c,d)
# (((1,1)_a,1)_b,1)_c
# x=ne11[d[3],a,d[1],d[2]]
# y=nex1[d[5],b,d[3],a,d[1],d[2],d[4]]
z=ne(op1[d[7]],mey1x1[c,d[5],b,d[3],a,d[1],d[2],d[4],d[6]])
if (a!=5 && b==5) || ((a!=5 || b!=5) && c==5)
z=:0
end
return (z,eval(z))
end

function doe3(a,b,c,d)
# ((a,(b,c)),d)
# x=me(b,ne(d[1],1),ne(d[2],1))
# y=ne1x[d[5],a,d[3],d[4],b,d[1],d[2]]
z=ne(op1[d[7]],mey11x[c,d[5],a,d[3],d[4],b,d[1],d[2],d[6]])
if a==5 || c==5
z=:0
end
return (z,eval(z))
end

function doe4(a,b,c,d)
# (a,((b,c),d))
# x=me(b,ne(d[1],1),ne(d[2],1))
# y=nex1[d[5],c,d[3],b,d[1],d[2],d[4]]
z=ne(op1[d[7]],me1yx1[a,d[6],d[5],c,d[3],b,d[1],d[2],d[4]])
if a==5 || (b!=5 && c==5)
z=:0
end
return (z,eval(z))
end

function doe5(a,b,c,d)
# (a,(b,(c,d)))
# x=me(c,ne(d[1],1),ne(d[2],1))
# y=ne1x[d[5],b,d[3],d[4],c,d[1],d[2]]
z=ne(op1[d[7]],me1y1x[a,d[6],d[5],b,d[3],d[4],c,d[1],d[2]])
if a==5 || b==5
z=:0
end
return (z,eval(z))
end

function main(pz)
d=[1,1,1,1,1,1,1]

for a in 1:5, b in 1:5, c in 1:5
for d[1] in 1:5, d[2] in 1:5, d[3] in 1:5, d[4] in 1:5, d[5] in 1:1, d[6] in 1:1, d[7] in 1:5
e,f = doe1(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe2(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe3(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe4(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe5(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
end
end
end

function fine(pz)
println(length(pz))
qz=[]
sort!(pz, by=x->x[2])
pl=length(pz)
for i in 1:pl
push!(qz,Int(ceil(pz[i][2])))
end

for i in 1:121
if in(i,qz)==true
j=findfirst(isequal(i),qz)
rz=repr(pz[j][1])
rz=replace(rz,":"=>"")
rz=replace(rz,"cc(+1, +1)"=>"11")
rz=replace(rz,"cc(+(11), +1)"=>"111")
rz=replace(rz,"+1"=>"1")
println(i," -> ",rz," = ",pz[j][2])
end
end
end

pz=[]
@time main(pz)
@time fine(pz)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
18.116161 seconds (20.53 M allocations: 1.368 GiB, 85.01% gc time)
420.117500 seconds (115.60 M allocations: 6.676 GiB, 4.81% gc time)
740539
1 -> (exp(1 + +((rog(1) - rog(exp(1) + exp(1)))))) = 0.5000000000000001
2 -> (cbr(1 + +(1 / exp(11)))) = 1.000005567202603
3 -> (cbr(1 - +(exp(1) * sqr(11)))) = 2.0012925727345556
4 -> (+((1 + +((exp(1) - cbr(1 / exp(1))))))) = 3.001750517885256
5 -> (exp(1 + +((exp(1) - sqr(exp(1) + exp(1)))))) = 4.0013741789600505
6 -> (+(1 * +(exp(1) / rog(1 - exp(1))))) = 5.021535230269061
7 -> (exp(1 + +(1 / cbr(1 + 1)))) = 6.011657650957131
8 -> (exp(+((1 + cbr(exp(1) * exp(1)))) - 1)) = 7.012778894114402
9 -> (+((+(exp(1) * sqr(11)) - 1))) = 8.015520899439874
10 -> (exp(1 + +((rog(1) + cbr(1 - exp(1)))))) = 9.004695764350682
11 -> (rog(1 + +(exp(11) / exp(1)))) = 10.000045398899218
12 -> (exp(+((1 + rog(11))) - 1)) = 11.000000000000002
13 -> (rog(1 + +(exp(1) * exp(11)))) = 12.000006144193478
14 -> (exp(+(rog(1 + exp(1)) * exp(1)) - 1)) = 13.063412466658708
15 -> (exp(+((sqr(exp(1) + rog(1)) + rog(1))) + 1)) = 14.135951028495938
16 -> (sqr(1 + +((exp(1) - exp(exp(1) + exp(1)))))) = 15.031080541832813
17 -> (+((1 + +(exp(1) / exp(1 - exp(1)))))) = 16.15426224147926
18 -> (+((+((1 + exp(exp(1) + rog(1)))) + 1))) = 17.154262241479262
19 -> (exp(+((exp(1) + exp(1 - exp(1)))) * 1)) = 18.131593378402822
20 -> (exp(1 + +((rog(1) + cbr(exp(1) * exp(1)))))) = 19.062709434872296
21 -> (+(+((exp(1) + rog(1))) * +(exp(1) * exp(1)))) = 20.085536923187664
22 -> (exp(+(+(11) / exp(1)) - 1)) = 21.045228353448717
23 -> (exp(1 + +((exp(1) - sqr(1 / exp(1)))))) = 22.46034183113636
24 -> (exp(1 + +(exp(1) / cbr(1 + 1)))) = 23.511783406204675
25 -> (exp(1 / +((rog(1 + exp(1)) - 1)))) = 24.34239023161335
26 -> (sqr(1 + +(exp(exp(1) + exp(1)) * exp(1)))) = 25.005158374895853
27 -> (exp(1 / +((1 - rog(1 + 1))))) = 26.020673411168517
28 -> (exp(1 + +((1 + sqr(1 - exp(1)))))) = 27.40793292860369
29 -> (cbr(1 - +(exp(11) / exp(1)))) = 28.031200676839138
30 -> (+(+((exp(1) + exp(1))) * +((exp(1) + exp(1))))) = 29.556224395722598
31 -> (exp(+((exp(1) + rog(exp(1) + exp(1)))) - 1)) = 30.308524482958514
32 -> (exp(+((exp(1) + cbr(1 / exp(1)))) * 1)) = 31.025614547492058
33 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) - 1)) = 32.350947201581
35 -> (exp(1 / +((1 - cbr(1 / exp(1)))))) = 34.046473983702036
36 -> (exp(+(sqr(1 - exp(1)) * exp(1)) * 1)) = 35.27632820034533
37 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) - 1)) = 36.68805458104296
38 -> (exp(+(11) - +(exp(1) * exp(1)))) = 37.00096158422464
39 -> (exp(+((sqr(1 + exp(1)) + exp(1))) - 1)) = 38.34279034771416
40 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) - 1)) = 39.09583226076508
41 -> (+(+((exp(1) + exp(1))) * +(exp(1) * exp(1)))) = 40.17107384637533
42 -> (exp(+((+((1 + 1)) + exp(1))) - 1)) = 41.1935556747161
43 -> (+((+((exp(1 + exp(1)) + rog(1))) + 1))) = 42.193555674716116
44 -> (+((+((1 + exp(1 + exp(1)))) + 1))) = 43.193555674716116
45 -> (exp(1 + +((exp(1) + exp(rog(1) - exp(1)))))) = 44.00353027860415
47 -> (exp(+(sqr(1 + 1) * exp(1)) * 1)) = 46.72274206040535
48 -> (exp(1 + +((exp(exp(1) - 1) - exp(1))))) = 47.30706718593521
50 -> (exp(1 + +((exp(1) + exp(1 - exp(1)))))) = 49.286780801520734
51 -> (exp(+((exp(1) + cbr(1 - exp(1)))) * 1)) = 50.20065233451701
52 -> (exp(+((exp(1) + cbr(11))) - 1)) = 51.5350376483723
54 -> (exp(+((cbr(1 + 1) + exp(1))) * 1)) = 53.42094397564407
55 -> (cbr(1 - +(exp(1) * exp(11)))) = 54.598038212039256
56 -> (exp(+(exp(1) / rog(1 - exp(1))) - 1)) = 55.78668552594455
57 -> (exp(+((exp(1) + sqr(1 - exp(1)))) * 1)) = 56.21110430560197
58 -> (exp(+(1 / exp(1)) * +(11))) = 57.20686180895067
60 -> (exp(1 + +((exp(1) + exp(rog(1) - 1))))) = 59.511005963978846
62 -> (exp(+((cbr(exp(1) + rog(1)) + exp(1))) * 1)) = 61.18452227596258
63 -> (exp(+((sqr(1 + 1) + exp(1))) * 1)) = 62.33327490494041
65 -> (exp(+((exp(1) + exp(1 / exp(1)))) * 1)) = 64.2607927021827
67 -> (sqr(1 - +(exp(exp(1) * exp(1)) * exp(1)))) = 66.31488392984271
68 -> (exp(+(cbr(1 + exp(1)) * exp(1)) * 1)) = 67.43919204479707
69 -> (exp(1 + +((1 + cbr(11))))) = 68.30480329758417
70 -> (exp(+(sqr(1 + exp(1)) * exp(1)) - 1)) = 69.52046855397886
71 -> (exp(1 + +(cbr(1 - exp(1)) * exp(1)))) = 70.51403099154743
72 -> (exp(+((cbr(1 + exp(1)) + exp(1))) * 1)) = 71.34344142876819
74 -> (exp(+(cbr(exp(1) * exp(1)) * exp(1)) - 1)) = 73.2948282500825
75 -> (exp(1 + +((rog(1) + sqr(11))))) = 74.93527870314952
76 -> (exp(1 + +((exp(1) + sqr(1 / exp(1)))))) = 75.55134476063684
77 -> (exp(+((1 + sqr(exp(1) + exp(1)))) + 1)) = 76.06924026135655
79 -> (exp(+((sqr(exp(1) + rog(1)) + exp(1))) * 1)) = 78.80710038074831
82 -> (+(+(exp(1) * exp(1)) * +(11))) = 81.27961708823715
83 -> (exp(+((exp(1) + rog(1 + 1))) + 1)) = 82.3871113494322
84 -> (+((+(exp(exp(1) + exp(1)) / exp(1)) - 1))) = 83.48412584713863
85 -> (exp(1 + +((exp(1) + cbr(1 / exp(1)))))) = 84.33636424122227
86 -> (+((1 + +(exp(exp(1) + exp(1)) / exp(1))))) = 85.48412584713863
88 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) * 1)) = 87.93899191149563
89 -> (exp(+(sqr(exp(1) + rog(1)) * exp(1)) * 1)) = 88.38383317988601
96 -> (exp(1 + +(sqr(1 - exp(1)) * exp(1)))) = 95.89100192175613
97 -> (exp(1 + +(rog(1 + exp(1)) * exp(1)))) = 96.52628755961122
98 -> (exp(+(1 / exp(1 - exp(1))) - 1)) = 97.02236556502673
100 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) * 1)) = 99.72847208916271
105 -> (exp(+((sqr(1 + exp(1)) + exp(1))) * 1)) = 104.22651025460627
107 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) * 1)) = 106.27349040292063
110 -> (+((+(111) - 1))) = 110
111 -> (+((+(exp(1 + exp(1)) * exp(1)) - 1))) = 110.9756938401968
112 -> (exp(+((1 + 1)) + +((exp(1) + rog(1))))) = 111.97569384019675
113 -> (+((1 + +(exp(1 + exp(1)) * exp(1))))) = 112.9756938401968
120 -> (exp(+(cbr(exp(1) + exp(1)) * exp(1)) * 1)) = 119.07124802191112
121 -> (exp(1 + +(cbr(exp(1) + rog(1)) * exp(1)))) = 120.74343163832631
10.825410 seconds (63.96 M allocations: 993.083 MiB, 36.17% gc time)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

2019年1月27日日曜日

Juliaでパズル（４）

（Juliaでパズル（３）からの続き）

目的はJuliaプログラミングの習得ということなのでこの話題がだらだたと続いている。高速が話題のJuliaで初めて実行時間が１時間以上のプログラムを作ったが，アルゴリズムがポンコツだといくらでも遅いプログラムができるのであった・・・orz

昨日の２項演算だけの組み合わせで整数を作るものから，単項演算を加えた組み合わせに拡張する。$5 \times 5^3 = 625$回から$5 \times 5^3 \times 5^7 ≒ 5000$万回に増えて７時間で終わらなかったので，次は時間短縮の方法を考えなければならない。

単項演算は１つのエントリーポイントで１回だけ作用するものとしている。
以下は，７重ループを５重ループに落として７分程度で実行したもの。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function sqr(x)
return sqrt(abs(x))
end

function cbr(x)
return cbrt(abs(x))
end

function rog(x)
return log(abs(x))
end

function cc(a,b)
if (a==1 && b==1) || (a==11 && b==1)
return 10*a+b
else
return 0
end
end

function ne(op, op1)
expr = Expr(:call,Symbol(op), op1)
return expr
end

function me(op, op1, op2)
expr = Expr(:call,Symbol(op), op1, op2)
return expr
end

function doe1(a,b,c,d)
# ((1,1)_a,(1,1)_c)_b
x=me(a,ne(d[1],1),ne(d[2],1))
y=me(c,ne(d[3],1),ne(d[4],1))
z=ne(d[7],me(b,ne(d[5],x),ne(d[6],y)))
if b=="cc"
z=:0
end
return (z,eval(z))
end

function doe2(a,b,c,d)
# (((1,1)_a,1)_b,1)_c
x=me(a,ne(d[1],1),ne(d[2],1))
y=me(b,ne(d[3],x),ne(d[4],1))
z=ne(d[7],me(c,ne(d[5],y),ne(d[6],1)))
if (a!="cc" && b=="cc") || ((a!="cc" || b!="cc") && c=="cc")
z=:0
end
return (z,eval(z))
end

function doe3(a,b,c,d)
# ((1,(1,1)_b)_a,1)_c
x=me(b,ne(d[1],1),ne(d[2],1))
y=me(a,ne(d[3],1),ne(d[4],x))
z=ne(d[7],me(c,ne(d[5],y),ne(d[6],1)))
if a=="cc" || c=="cc"
z=:0
end
return (z,eval(z))
end

function doe4(a,b,c,d)
# (1,((1,1)_b,1)_c)_a
x=me(b,ne(d[1],1),ne(d[2],1))
y=me(c,ne(d[3],x),ne(d[4],1))
z=ne(d[7],me(a,ne(d[6],1),ne(d[5],y)))
if a=="cc" || (b!="cc" && c=="cc")
z=:0
end
return (z,eval(z))
end

function doe5(a,b,c,d)
# (1,(1,(1,1)_c)_b)_a
x=me(c,ne(d[1],1),ne(d[2],1))
y=me(b,ne(d[3],1),ne(d[4],x))
z=ne(d[7],me(a,ne(d[6],1),ne(d[5],y)))
if a=="cc" || b=="cc"
z=:0
end
return (z,eval(z))
end

function main(pz)
uno=["+","-","*","/","cc"]
dos=["+","sqr","cbr","exp","rog"]
d=["+","+","+","+","+","+","+"]
for a in uno, b in uno, c in uno
#a="*"; b="-"; c="/";
for d[1] in dos, d[2] in dos, d[3] in dos, d[4] in dos, d[7] in dos
#, d[5] in dos, d[6] in dos, d[7] in dos
e,f = doe1(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe2(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe3(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe4(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe5(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
end
end
end

function fine(pz)
println(length(pz))
qz=[]
sort!(pz, by=x->x[2])
pl=length(pz)
for i in 1:pl
push!(qz,Int(ceil(pz[i][2])))
end

for i in 1:121
if in(i,qz)==true
j=findfirst(isequal(i),qz)
rz=repr(pz[j][1])
rz=replace(rz,":"=>"")
rz=replace(rz,"cc(+1, +1)"=>"11")
rz=replace(rz,"cc(+(11), +1)"=>"111")
rz=replace(rz,"+1"=>"1")
println(i," -> ",rz," = ",pz[j][2])
end
end
end

pz=[]
@time main(pz)
@time fine(pz)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

408.936226 seconds (144.22 M allocations: 8.591 GiB, 5.96% gc time)

751191

1 -> (exp(+((rog(1) - rog(exp(1) + exp(1)))) + 1)) = 0.5000000000000001

2 -> (cbr(+(1 / exp(11)) + 1)) = 1.000005567202603

3 -> (cbr(+(exp(1) * sqr(11)) - 1)) = 2.0012925727345556

4 -> (+((+((exp(1) - cbr(1 / exp(1)))) + 1))) = 3.001750517885256

5 -> (exp(+((exp(1) - sqr(exp(1) + exp(1)))) + 1)) = 4.0013741789600505

6 -> (+(+(exp(1) / rog(1 - exp(1))) * 1)) = 5.021535230269061

7 -> (exp(+(1 / cbr(1 + 1)) + 1)) = 6.011657650957131

8 -> (exp(+((1 + cbr(exp(1) * exp(1)))) - 1)) = 7.012778894114402

9 -> (+((+(exp(1) * sqr(11)) - 1))) = 8.015520899439874

10 -> (exp(+((rog(1) + cbr(1 - exp(1)))) + 1)) = 9.004695764350682

11 -> (rog(+(exp(11) / exp(1)) + 1)) = 10.000045398899218

12 -> (exp(+((1 + rog(11))) - 1)) = 11.000000000000002

13 -> (rog(+(exp(1) * exp(11)) + 1)) = 12.000006144193478

14 -> (exp(+(rog(1 + exp(1)) * exp(1)) - 1)) = 13.063412466658708

15 -> (exp(+((sqr(exp(1) + rog(1)) + rog(1))) + 1)) = 14.135951028495938

16 -> (sqr(+((exp(1) - exp(exp(1) + exp(1)))) + 1)) = 15.031080541832813

17 -> (+((+(exp(1) / exp(1 - exp(1))) + 1))) = 16.15426224147926

18 -> (+((+((1 + exp(exp(1) + rog(1)))) + 1))) = 17.154262241479262

19 -> (exp(+((exp(1) + exp(1 - exp(1)))) * 1)) = 18.131593378402822

20 -> (exp(+((rog(1) + cbr(exp(1) * exp(1)))) + 1)) = 19.062709434872296

21 -> (+(+((exp(1) + rog(1))) * +(exp(1) * exp(1)))) = 20.085536923187664

22 -> (exp(+(+(11) / exp(1)) - 1)) = 21.045228353448717

23 -> (exp(+((exp(1) - sqr(1 / exp(1)))) + 1)) = 22.46034183113636

24 -> (exp(+(exp(1) / cbr(1 + 1)) + 1)) = 23.511783406204675

25 -> (sqr(+(exp(exp(1) * exp(1)) / exp(1)) - 1)) = 24.37815447036041

26 -> (sqr(+(exp(exp(1) + exp(1)) * exp(1)) + 1)) = 25.005158374895853

27 -> (exp(+((exp(1) + rog(1 - exp(1)))) * 1)) = 26.039293433236857

28 -> (exp(+((1 + sqr(1 - exp(1)))) + 1)) = 27.40793292860369

29 -> (cbr(+(exp(11) / exp(1)) - 1)) = 28.031200676839138

30 -> (+(+((exp(1) + exp(1))) * +((exp(1) + exp(1))))) = 29.556224395722598

31 -> (exp(+((exp(1) + rog(exp(1) + exp(1)))) - 1)) = 30.308524482958514

32 -> (exp(+((exp(1) + cbr(1 / exp(1)))) * 1)) = 31.025614547492058

33 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) - 1)) = 32.350947201581

35 -> (exp(+((exp(1) - exp(1 - exp(1)))) + 1)) = 34.42929324111139

36 -> (exp(+(sqr(1 - exp(1)) * exp(1)) * 1)) = 35.27632820034533

37 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) - 1)) = 36.68805458104296

38 -> (exp(+(11) - +(exp(1) * exp(1)))) = 37.00096158422464

39 -> (exp(+((sqr(1 + exp(1)) + exp(1))) - 1)) = 38.34279034771416

40 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) - 1)) = 39.09583226076508

41 -> (+(+((exp(1) + exp(1))) * +(exp(1) * exp(1)))) = 40.17107384637533

42 -> (exp(+((+((1 + 1)) + exp(1))) - 1)) = 41.1935556747161

43 -> (+((+((exp(1 + exp(1)) + rog(1))) + 1))) = 42.193555674716116

44 -> (+((+((1 + exp(1 + exp(1)))) + 1))) = 43.193555674716116

45 -> (exp(+((exp(1) + exp(rog(1) - exp(1)))) + 1)) = 44.00353027860415

47 -> (exp(+(sqr(1 + 1) * exp(1)) * 1)) = 46.72274206040535

48 -> (exp(+((exp(exp(1) - 1) - exp(1))) + 1)) = 47.30706718593521

50 -> (exp(+((exp(1) + exp(1 - exp(1)))) + 1)) = 49.286780801520734

51 -> (exp(+((exp(1) + cbr(1 - exp(1)))) * 1)) = 50.20065233451701

52 -> (exp(+((exp(1) + cbr(11))) - 1)) = 51.5350376483723

54 -> (exp(+((cbr(1 + 1) + exp(1))) * 1)) = 53.42094397564407

55 -> (cbr(+(exp(1) * exp(11)) - 1)) = 54.598038212039256

56 -> (exp(+(exp(1) / rog(1 - exp(1))) - 1)) = 55.78668552594455

57 -> (exp(+((exp(1) + sqr(1 - exp(1)))) * 1)) = 56.21110430560197

58 -> (exp(+(1 / exp(1)) * +(11))) = 57.20686180895067

60 -> (exp(+((exp(1) + exp(rog(1) - 1))) + 1)) = 59.511005963978846

62 -> (exp(+((cbr(exp(1) + rog(1)) + exp(1))) * 1)) = 61.18452227596258

63 -> (exp(+((sqr(1 + 1) + exp(1))) * 1)) = 62.33327490494041

65 -> (exp(+((exp(1) + exp(1 / exp(1)))) * 1)) = 64.2607927021827

67 -> (sqr(+(exp(exp(1) * exp(1)) * exp(1)) - 1)) = 66.31488392984271

68 -> (exp(+(cbr(1 + exp(1)) * exp(1)) * 1)) = 67.43919204479707

69 -> (exp(+((1 + cbr(11))) + 1)) = 68.30480329758417

70 -> (exp(+(sqr(1 + exp(1)) * exp(1)) - 1)) = 69.52046855397886

71 -> (exp(+(cbr(1 - exp(1)) * exp(1)) + 1)) = 70.51403099154743

72 -> (exp(+((cbr(1 + exp(1)) + exp(1))) * 1)) = 71.34344142876819

74 -> (exp(+(cbr(exp(1) * exp(1)) * exp(1)) - 1)) = 73.2948282500825

75 -> (exp(+((rog(1) + sqr(11))) + 1)) = 74.93527870314952

76 -> (exp(+((exp(1) + sqr(1 / exp(1)))) + 1)) = 75.55134476063684

77 -> (exp(+((1 + sqr(exp(1) + exp(1)))) + 1)) = 76.06924026135655

79 -> (exp(+((sqr(exp(1) + rog(1)) + exp(1))) * 1)) = 78.80710038074831

82 -> (+(+(exp(1) * exp(1)) * +(11))) = 81.27961708823715

83 -> (exp(+((exp(1) + rog(1 + 1))) + 1)) = 82.3871113494322

84 -> (+((+(exp(exp(1) + exp(1)) / exp(1)) - 1))) = 83.48412584713863

85 -> (exp(+((exp(1) + cbr(1 / exp(1)))) + 1)) = 84.33636424122227

86 -> (+((+(exp(exp(1) + exp(1)) / exp(1)) + 1))) = 85.48412584713863

88 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) * 1)) = 87.93899191149563

89 -> (exp(+(sqr(exp(1) + rog(1)) * exp(1)) * 1)) = 88.38383317988601

96 -> (exp(+(sqr(1 - exp(1)) * exp(1)) + 1)) = 95.89100192175613

97 -> (exp(+(rog(1 + exp(1)) * exp(1)) + 1)) = 96.52628755961122

98 -> (exp(+(1 / exp(1 - exp(1))) - 1)) = 97.02236556502673

100 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) * 1)) = 99.72847208916271

105 -> (exp(+((sqr(1 + exp(1)) + exp(1))) * 1)) = 104.22651025460627

107 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) * 1)) = 106.27349040292063

110 -> (+((+(cc(+(11), 1)) - 1))) = 110

111 -> (+((+(exp(1 + exp(1)) * exp(1)) - 1))) = 110.9756938401968

112 -> (exp(+((1 + 1)) + +((exp(1) + rog(1))))) = 111.97569384019675

113 -> (+((+(exp(1 + exp(1)) * exp(1)) + 1))) = 112.9756938401968

120 -> (exp(+(cbr(exp(1) + exp(1)) * exp(1)) * 1)) = 119.07124802191112

121 -> (exp(+(cbr(exp(1) + rog(1)) * exp(1)) + 1)) = 120.74343163832631

6.644247 seconds (64.72 M allocations: 1004.683 MiB, 1.93% gc time)

2019年1月26日土曜日

Juliaでパズル（３）

（Juliaでパズル（２）からの続き）

世の中で流布している数学パズルだと，４つの数字と加減乗除を組み合わせて10をつくる Make 10 が有名らしい。もちろん，この手の問題では計算結果はすべて整数となるのが普通であるが，今回は，実数の範囲で計算して最後に整数化できるものを探している。

とりあえず，２項演算だけを組み合わせた結果をすべて探すプログラムを作ってみた。

(1) findfirst()のように，ある目的の数に一致する数をコレクション（配列またはタプル）の中から探すとともに，その位置の添字を返すことができれば有り難い。いまのところ「タプルの配列に対して１コマンドである数字を見つける方法」がわからない。

(2) ２次元配列ならば融通がききそうなのであるが，「タプルの配列を１コマンドで多次元配列に展開する方法」がわからない。もちろん地道にやればできることなのだが。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function cc(a,b)
return 10*a+b
end

function me(op, op1, op2)
expr = Expr(:call,Symbol(op), op1, op2)
return expr
end

function dop1(a,b,c)
# ((a,b),(cd))
x=me(a,1,1)
y=me(c,1,1)
z=me(b,x,y)
if b=="cc"
z=:0
end
return (z,eval(z))
end

function dop2(a,b,c)
# (((a,b),c),d)
x=me(a,1,1)
y=me(b,x,1)
z=me(c,y,1)
if (a!="cc" && b=="cc") || ((a!="cc" || b!="cc") && c=="cc")
z=:0
end
return (z,eval(z))
end

function dop3(a,b,c)
# ((a,(b,c)),d)
x=me(b,1,1)
y=me(a,1,x)
z=me(c,y,1)
if a=="cc" || c=="cc"
z=:0
end
return (z,eval(z))
end

function dop4(a,b,c)
# (a,((b,c),d))
x=me(b,1,1)
y=me(c,x,1)
z=me(a,1,y)
if a=="cc" || (b!="cc" && c=="cc")
z=:0
end
return (z,eval(z))
end

function dop5(a,b,c)
# (a,(b,(c,d)))
x=me(c,1,1)
y=me(b,1,x)
z=me(a,1,y)
if a=="cc" || b=="cc"
z=:0
end
return (z,eval(z))
end

pz=[]
uno=["+","-","*","/","cc"]
for a in uno, b in uno, c in uno
d,e = dop1(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
d,e = dop2(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
d,e = dop3(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
d,e = dop4(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
d,e = dop5(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
end

qz=[]
sort!(pz, by=x->x[2])
pl=length(pz)
for i in 1:pl
push!(qz,Int(ceil(pz[i][2])))
end

for i in 1:121
if in(i,qz)==true
j=findfirst(isequal(i),qz)
rz=repr(pz[j][1])
rz=replace(rz,":"=>"")
rz=replace(rz,"cc(1, 1)"=>"11")
rz=replace(rz,"cc(11, 1)"=>"111")
println(i," -> ",rz," = ",pz[j][2])
end
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

1 -> (1 - 1 / 11) = 0.9090909090909091

2 -> (1 + 1 / 11) = 1.0909090909090908

3 -> ((1 + 1) + 1 * 1) = 3

4 -> ((1 + 1) + (1 + 1)) = 4

6 -> (11 / (1 + 1)) = 5.5

9 -> (11 - (1 + 1)) = 9

10 -> (1 * 11 - 1) = 10

11 -> ((1 + 11) - 1) = 11

12 -> (1 + 1 * 11) = 12

13 -> ((1 + 1) + 11) = 13

22 -> ((1 + 1) * 11) = 22

110 -> (111 - 1) = 110

111 -> (1 * 111) = 111

112 -> (1 + 111) = 112

121 -> (11 * 11) = 121

「冬麗の不思議をにぎる赤ン坊」（野澤節子　1920-1995）

（Juliaでパズル（４）に続く）

2019年1月25日金曜日

Juliaでパズル（２）

（Juliaでパズル（１）からの続き）

４つの数字から整数をつくるパズルをコンピュータで処理させるためには，式に対応する文字列を，プログラム内で式として評価する必要がある。これは，JuliaドキュメントのMetaprogrammingのところに詳しいが，ここでのポイントは，Juliaでは，Juliaプログラムの中で，Juliaコードを生成したり操作したりすることができるということである。具体例としては次のような例文が示されている。

julia> function math_expr(op, op1, op2)
           expr = Expr(:call, op, op1, op2)
           return expr
       end
math_expr (generic function with 1 method)

julia>  ex = math_expr(:+, 1, Expr(:call, :*, 4, 5))
:(1 + 4 * 5)

julia> eval(ex)
21

そこで（1, 1, 1, 1）パズルで上記を用い，２項演算子を３つ与えた場合の５つのパターンの構文とその値を出力する関数を書いてみる。cc(a, b) = 10*a + bは２つの数字を並べて２桁の数（例えば a=1, b=1 なら cc(1,1) = 11である）を作る関数である。funcion me() においてSymbol(op)は文字列として与えた２項演算子を Expr型に変換する関数である（注：Ver 0.6以前の資料では symbol と小文字で始まっているので，Ver. 1.0 以降では大文字にすること）。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function cc(a,b)
return 10*a+b
end

function me(op, op1, op2)
expr = Expr(:call,Symbol(op), op1, op2)
return expr
end

function dop1(a,b,c)
# ((a,b),(cd))
x=me(a,1,1)
y=me(c,1,1)
z=me(b,x,y)
println((z,eval(z)))
end

function dop2(a,b,c)
# (((a,b),c),d)
x=me(a,1,1)
y=me(b,x,1)
z=me(c,y,1)
println((z,eval(z)))
end

function dop3(a,b,c)
# ((a,(b,c)),d)
x=me(b,1,1)
y=me(a,1,x)
z=me(c,y,1)
println((z,eval(z)))
end

function dop4(a,b,c)
# (a,((b,c),d))
x=me(b,1,1)
y=me(c,x,1)
z=me(a,1,y)
println((z,eval(z)))
end

function dop5(a,b,c)
# (a,(b,(c,d)))
x=me(c,1,1)
y=me(b,1,x)
z=me(a,1,y)
println((z,eval(z)))
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

dop1("cc","/","+")
dop2("cc","+","*")
dop3("-","*","/")
dop4("+","\\","^")
dop5("+","-","*")

(:(cc(1, 1) / (1 + 1)), 5.5)
(:((cc(1, 1) + 1) * 1), 12)
(:((1 - 1 * 1) / 1), 0.0)
(:(1 + (1 \ 1) ^ 1), 2.0)
(:(1 + (1 - 1 * 1)), 1)

（Juliaでパズル（３）に続く）

2019年1月24日木曜日

Juliaでパズル（１）

４つの同じ数字に演算をほどこして，１から順に整数を作るというパズルがある。
例えば，6÷6+6-6 = 1,　 6÷6+6÷6=2,　 (6+6+6)/6 = 3,　などなど。
これを，Juliaの基本演算子や基本関数だけを使って実行してみると，

1　1+1-1*1
2　1+1-1+1
3　1+1*1+1
4　1+1+1+1
5　Int(floor(11/(1+1)))
6　Int(ceil(11/(1+1)))
7　Int(floor(exp(1*1+1*1)))
8　Int(ceil(exp(1*1+1*1)))
9　11-1-1
10　11-1*1
11　11+1-1
12　11+1*1
13　11+1+1
14　Int(ceil(sqrt(sqrt(exp(11)))-1-1))
15　Int(floor(exp(1+sqrt(1+1+1))))
16　Int(ceil(exp(1+sqrt(1+1+1))))
17　Int(floor(exp((1+1)*sqrt(1+1))))
18　Int(ceil(exp((1+1)*sqrt(1+1))))
19　Int(floor(exp(1+1+1)-1))
20　Int(floor(exp(1+1*1+1)))
21　Int(floor(exp(1+1+1)+1))
22　Int(ceil(exp(1+1+1)+1))
23　Int(floor(exp(sqrt(1+1)*cbrt(11))))
24　Int(ceil(exp(sqrt(1+1)*cbrt(11))))
25　Int(floor(exp(sqrt(11))-1-1))
26　Int(floor(exp(sqrt(11))-1*1))
27　Int(floor(exp(sqrt(11))+1-1))
28　Int(ceil(exp(sqrt(11))+1-1))
29　Int(ceil(exp(sqrt(11))+1*1))
30　Int(ceil(exp(sqrt(11))+1+1))
31　Int(floor(sqrt(sqrt(exp(11)))*(1+1)))
32　Int(floor(cbrt(exp(sqrt(111)))-1))
33　Int(floor(cbrt(exp(sqrt(111)))*1))
34　Int(floor(cbrt(exp(sqrt(111)))+1))
35　Int(ceil(cbrt(exp(sqrt(111)))+1))
36　Int(floor(exp(sqrt(11+1+1))))
37　Int(floor(cbrt(exp(11))-1-1))
38　Int(ceil(cbrt(exp(11))-1-1))
39　Int(ceil(cbrt(exp(11))-1*1))
40　Int(ceil(cbrt(exp(11))+1-1))
41　Int(ceil(cbrt(exp(11))+1*1))
42　Int(ceil(cbrt(exp(11))+1+1))
43　Int(floor(cbrt(sqrt(1+1)*exp(11))))
44　Int(ceil(cbrt(sqrt(1+1)*exp(11))))
45　Int(ceil(exp(cbrt(111)-1)))
46　Int(ceil(exp(sqrt(sqrt(11))+1+1)))
47　Int(ceil(sqrt(exp(11-sqrt(11)))))
などなど

（１）もし，ユーザ定義関数が使えるならば，任意の n に対して，すべての整数を作り出すことが簡単にできるのであった・・・

f(n1,n2,n3,n4)=Int(n1/n2+n3-n4)
f(n)=Int(ceil(n+1/n))

f(1,1,1,1)
1
f(f(1,1,1,1))
2
f(f(f(1,1,1,1)))
3

（２）juliaドキュメントの算術的演算と基本関数から主な演算と関数を選んで使うことにする。例えば，２項演算ならば，+, -, *, / の４種とする（ここでは，整数ゼロ除算をさけるために÷と%は含めない）。数字を２〜４桁並べる操作が必要なのでこれに対応する２項演算を１つ付け加えると５種になる。この２項演算は cc(a,b)=10*a+b で実現できる（例えば３桁並べるには，cc(cc(a, b), c)とすればよい），ただし適応できる位置に制限があることに注意する。次に単項演算を選ぶ。ここでは，+, sqrt(abs()), qbrt(abs()), exp(), log(abs()) の５種とする（+は恒等演算子）。

（３）４つの数字（a, b, c, d, a = b = c = d）に対応する２項演算のパターンは５つある。
２項演算を（a, b）などで表し，ふつうのカッコの優先順位で計算するものとして，
（（a, b）,（c, d）），（（（a, b）, c）, d），（a,（b,（c, d））），
（（a,（b, c））, d），（a,（（b, c）, d））の５パターンである。１パターンあたりの２項演算の数は３であるから，パターンあたりの総数は$5^3$通りである。

（４）上記の各数字や２項演算の結果に対して５種の単項演算を任意の階数だけ作用することができる。単項演算のエントリーポイントは１パターンにつき$n$個あるとする。単項演算子の種類を$m$個とし，エントリーポイントに$k$個まで入るならば，パターンあたりの配位の数は，$\sum_{j=0}^k {}_{n+j-1}C_j m^k$となる。

したがって，この他に条件をつけなければ，総数で $5 \times 5^3 \times (5^k \dfrac{n^k}{k!}+\cdots) $個の式が得られる。これをしらみつぶしに調べるのだろうか・・・

（Juliaでパズル（２）に続く）

2019年1月23日水曜日

円周率が22/7より小さいことの証明

というタイトルのWikipedia記事がある（円周率が22/7より小さいことの証明）。

証明は次の通りである。
\begin{equation}
\begin{aligned}
0 &< \int_0^1 \dfrac{x^4(1-x)^4}{1+x^2} dx \\
& = \int_0^1\dfrac{x^4-4x^5+6x^6-4x^7+x^8}{1+x^2} dx\\
&= \int_0^1x^6 -4x^5 + 5x^4 -4x^2 +4 -\dfrac{4}{1+x^2} dx\\
&=\frac{x^7}{7}-\frac{2 x^6}{3}+x^5-\frac{4 x^3}{3}+4 x-4 \tan^{-1} x\ \Big|_0^1\\
&= \frac{1}{7}-\frac{2}{3}+1-\frac{4}{3}+4-\pi\\
&=\frac{22}{7}-\pi
\end{aligned}
\end{equation}

同様にして円周率が355/113より小さいことも証明される。
\begin{equation}
\begin{aligned}
0 &< \int_0^1 \dfrac{x^8(1-x)^8(25+816x^2)}{3164(1+x^2)} dx \\
&= \cdots\\
&= \frac{12 x^{17}}{791}
-\frac{102 x^{16}}{791}
+\frac{3151 x^{15}}{6780}
-\frac{703 x^{14}}{791}
+\frac{393 x^{13}}{452}
-\frac{23 x^{12}}{113}\\
&-\frac{145 x^{11}}{452}
-\frac{5 x^{10}}{791}
+\frac{1409 x^9}{3164}
-\frac{4 x^7}{7}
+\frac{4 x^5}{5}
-\frac{4 x^3}{3}
+4 x−4 \tan^{-1}x \Bigr|_0^1 \\
&= \cdots\\
&=\frac{355}{113}-\pi
\end{aligned}
\end{equation}

えーっ，3164と816と25はどこから連れてきたのだ・・・という方は上記記事から参考文献をたどってみよう。この世界にはまだまだ自分が知らないことが沢山埋もれている。

2018年12月27日木曜日

素数鉛筆

京都大学生協限定素数価格577円で販売している素数ものさしは皆様よくご存知のことと思う。それでは素数えんぴつ（A Prime Pencil）はどうだろう。このあたりで販売しているようだ。本日twitterでWolfram Japanさんが宣伝していたのでさっそく試してみた。いや，購入したのではなくて，「Juliaでプログラミング教育をはじめよう（仮題）」の練習問題を作ったのだ。本家のほうはバリバリのMathematica関数プログラミングなので，スマートすぎておじぃさんには手に負えない。

注意点：今日も型の整合性のところで難渋した。適当な演算では自動型変換は行われないので，明示的に型を指定する必要があった。10^mが Int128 (-2^{127} 〜 2^{127} - 1)を必要とする場面で，10^m → Int128(10)^mとすべきところに嵌まった。また，zeros()とかones() で配列を準備するところで，型指定が必要なのに気付かなかった。p = zeros(Int128,N)とか q = ones(Int,N) などとしなければならない。

結果は以下の通りで出力は省略するが，23桁には３つの素数（ 96686312646216567629137， 57686312646216567629137，95918918997653319693967）があり，最大の24桁では１つの素数（357686312646216567629137）が得られた。これ以上はありません。これを鉛筆に焼いて左から削っていけばよい。そうそう，素数鉛筆に印刷されている切り捨て可能素数の定義を説明するのを忘れていた。与えられた素数の上位桁（または下位桁）から１桁づつとりのぞいた数が再び素数であり，最後に１桁の素数が残るものである。

P. S. @ksasaoさんが右から削る場合にも言及していたので，こちらも加えた。いずれも終了判定はしていないので人力でがんばってほしい。

using Primes

#左から削っていく場合の切り捨て可能素数
function pencil1(N)
p=zeros(Int128,N)
q=ones(Int,N)
r=Int128(10)
p[1]=3
p[2]=7
nb = 0
nt = 2
while max(nb,nt) < N-1
nb=nb+1
for j in 1:9
k= (r^q[nb])*j+p[nb]
if isprime(k)
nt=nt+1
p[nt]=k
q[nt]=q[nb]+1
end
end
end
println(q)
println(p)
end

#右から削っていく場合の切り捨て可能素数
function pencil2(N)
p=zeros(Int128,N)
q=ones(Int,N)
r=Int128(10)
p[1]=2
p[2]=3
p[3]=5
p[4]=7
nb = 0
nt = 4
while max(nb,nt) < N-1
nb=nb+1
for j in 1:9
k= 10*p[nb]+j
if isprime(k)
nt=nt+1
p[nt]=k
q[nt]=q[nb]+1
end
end
end
println((q,p))

end

登録: 投稿 (Atom)