SymPy で二重根号を外す

以前の記事で，「SymPy では二重根号を外せない」と書いたが，そんなことはなかった。
sqrtdenest() があった。Maxima にもあるようだが，SymPy の sqrtdenest() に言及した日本語の記事は 1 つしか見つからなかった。
 

ちなみに，そこには「二重根号外しを二重にやらなければいけないケースにどのくらい対応できているか試してみた」とあり，うまく行かないという以下の例が示されていた。
 
Python の場合

import sympy as S
S.init_printing()
a = S.sqrt(S.sqrt(49 + 20 * S.sqrt(6)))
a

S.sqrtdenest(a)

うまくやるためには sqrtdenest を 2 回使うことである（Python の sqrtdenest のソースドキュメントにヒントが書いてあった）。

S.sqrtdenest(S.sqrt(S.sqrtdenest(S.sqrt(49+20*S.sqrt(6)))))

Julia の場合

Julia では，ルート中の整数は Sym() でシンボリックナンバーにすることもできる（そのようにすれば sqrt は自動的に対応する）。しかし，sqrtdenest() はインポートされていないので sympy.sqrtdenest() として使用する。

using SymPy
a = sqrt(sqrt(49 + 20sqrt(Sym(6))))

a  |>  N

3.1462643699419723423291350657155704455124771291873287012324867174426654953709

sympy.sqrtdenest(sqrt(sqrt(49 + 20sqrt(Sym(6)))))

`sqrtdenest` を 2 回使う。

c = sympy.sqrtdenest(sqrt(sympy.sqrtdenest(sqrt(49 + 20sqrt(Sym(6))))))

c |> N

3.1462643699419723423291350657155704455124771291873287012324867174426654953709

---------------------------------

なお，以下のような場合には 1 回 `sqrtdenest` を使うだけでちゃんと動く。

Python の場合

b = S.sqrt(1 + S.sqrt(40 + 16 * S.sqrt(6)))
b


N(b)

3.1462643699419723423291350657155704455124771291873287012324867174426654953709

S.sqrtdenest(b)

Julia の場合

b = sqrt(1 + sqrt(40 + 16 * sqrt(Sym(6))))

sympy.sqrtdenest(b)

大学入試問題を SymPy で解く

以下にある問題を SymPy で解いてみる。
大学入試数学の問題
http://mikiotaniguchi.com/main/center/main.htm

22092501剰余の定理22京都産業大

22082401定積分01芝浦工業大

14121806円と直線05大阪市立大

バーゼル問題

1/n^2 を n = 1, 2, ..., ∞ として，和を取ると π^2/6 に収束する。

using SymPy

summation(1/n^2, (n, 1, oo)) |> string

"pi^2/6"

ゲルフォントの定数

julia> ℯ^π

23.140692632779267

julia> big(ℯ)^π

23.14069263277926900572908636794854738026610624260021199344504640952434235069032

有効桁をほぼ 300 桁になるように設定する。

julia> setprecision(BigFloat, 1000)

1000

julia> big(ℯ)^π

23.1406926327792690057290863679485473802661062426002119934450464095243423506904527835169719970675492196759527048010877731444280444146938358447174458796098493653279658636692422302689910137417646844014103951838684772430680595881624498444914309667784136716319634147840382165112876377314703473538331628212899

 

年齢の計算

年齢の計算

生年月日 year, month, day として，ある時点 year2, month2, day2 での年齢を求めるのに以下のようなプログラムを書いていた。

関数定義
age(year, month, day, year2, month2, day2) =
    year2 - year - (month > month2 || (month == month2 && day > day2))

実行
age(2020, 9, 19, 2022,  8, 29),
age(2020, 9, 19, 2022,  9, 18),
age(2020, 9, 19, 2022, 10, 18),
age(2020, 9, 19, 2022,  9, 19)

実行結果
(1, 1, 2, 2)

日付を yyyymmdd のような 8 桁の整数で受け渡しすると，以下のように簡単なプログラムになる。

ちゃんと上のプログラムと同じように月と日の大小順を考慮することになっている。

関数定義
age(ymd, ymd2) = (ymd2 - ymd) ÷ 10000

実行
age(20200919, 20220829),
age(20200919, 20220918),
age(20200919, 20221018),
age(20200919, 20220919)

実行結果
(1, 1, 2, 2)

 

Julia で閏年判定

Dates パッケージの isleapyear() を使う

julia> using Dates
julia> isleapyear.([2020, 2022, 2100, 2400])
4-element BitVector:
 1
 0
 0
 1

パッケージを使わないなら以下のような関数を書く。

isleapyear2(y) = y % 400 == 0 || y % 4 == 0 && y % 100 != 0

Rump の例題--その2 Julia の SymPy で解く

Julia の SymPy は，Python の sympy を利用している。Python の整数型は必要に応じて長精度整数になる。

求める式には実数が含まれるが，この式を 4 倍することで整式になる。

Rump の例題

Rump の例題
https://www.tuhh.de/ti3/paper/rump/Ru88a.pdf

a = 77617, b = 33096 として，以下の計算を行う

R

f = function(a, b) 333.75*b^6 + a^2 *(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b)
library(gmp)
a = as.bigq(77617)
b = as.bigq(33096)
ans = f(a, b)
as.numeric(ans) # -0.8273960599468214

Julia

f(a, b) = 333.75b^6 + a^2*(11a^2*b^2 - b^6 - 121b^4 - 2) + 5.5b^8 + a/2b
a = big(77617)
b = big(33096)
f(a, b) # -0.8273960599468213681411650954798162919990331157843848199178148416727096930142628

Julia 1.8.1 が公開されました

Julia 1.8.1 が公開されました

 Version 1.8.1 (2022-09-06)

Official https://julialang.org/ release

AWK と R のコラボレーション 定式化された入力ファイルからデータファイルへ

一般解は

https://stackoverflow.com/questions/69385888/parse-multi-line-result-text-file-with-awk-bash-or-r

にあった。AWK スクリプトを平易に書き変えた。

★ 以下のような定式化された入力ファイル melting.txt 

  The MELTING results are :

 Enthalpy : -181,100 cal/mol ( -756,998 J /mol)

 Entropy : -467.3 cal/mol-K ( -1,953.31 J /mol-K)

 Melting temperature : 75.13 degrees C.

  The MELTING results are :

 Enthalpy : -170,800 cal/mol ( -713,944 J /mol)

 Entropy : -444 cal/mol-K ( -1,855.92 J /mol-K)

 Melting temperature : 70.6 degrees C.

からデータを取り出し，R で

  Enthalpy Entropy MeltingTemperature
1  -181100  -467.3              75.13
2  -170800  -444.0              70.60

のようなデータフレームを作る。

★ 用意する AWK スクリプト melting.awk

/Enthalpy/ || /Entropy/ { # "Enthalpy" または "Entropy" を含む行について
    n = split($0, a, " ") # 空白でフィールドに区切る
    gsub(",", "", a[3])   # a[3] が求めるデータである。数値区切りの「,」を消去する
    printf "%g, ", a[3]   # 書き出す
}
/Melting temperature/ {   # "Melting temperature" を含む行について
    n = split($0, a, " ") # 空白でフィールドに区切る
    print a[4]            # 求めるデータ a[4] を書き出す
}

★ R コンソールで以下を実行する。

CSV ファイルに結果を書き出す。
以下で |> を使っているが，tidyverse ではなく Base R のパイプである。
最初の文字列がコンソールで入力されるコマンドである。（脚注参照）

"awk -f melting.awk melting.txt" |>
    pipe() |>
    readLines() |>
    write("results.csv")

後は，ありふれた処理。CSV ファイルをデータフレームに読み込むだけである。

df = read.csv("results.csv", header=FALSE,
         col.names=c("Enthalpy", "Entropy", "MeltingTemperature"))

出来上がり。

  Enthalpy Entropy MeltingTemperature
1  -181100  -467.3              75.13
2  -170800  -444.0              70.60

なお，入力ファイルがちゃんと同じフォーマットであることを仮定している。
その仮定ができない場合の処理を AWK スクリプトに含めるかどうかは，程度問題。

＝＝＝＝＝＝＝＝＝＝
脚注

以下の一行をコンソールで入力すれば

$ > awk -f melting.awk melting.txt 

以下が表示される

-181100, -467.3, 75.13
-170800, -444, 70.6

結果をファイルに入れるには出力を results.csv へリダイレクトすればよい。

$ > awk -f melting.awk melting.txt > results.csv

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

なお，R だけでやろうとすると以下のようにもなるだろうか。

txt = readLines("melting.txt")
df = NULL
x = NULL
for (i in 1:length(txt)) {
	if (grepl("Enthalpy", txt[i]) || grepl("Entropy", txt[i])) {
		target = sub(",", "", unlist(strsplit(txt[i], " "))[4])
		x = append(x, target)
	} else if (grepl("Melting temperature", txt[i])) {
		x = append(x, unlist(strsplit(txt[i], " "))[5])
		df = rbind(df, x)
		x = NULL
	}
}
df = as.data.frame(df, col.names=c("Enthalpy", "Entropy", "MeltingTemperature"))
print(df)

Julia でやろうとすると以下のようになる（書式に関して別の仮定を置く）

using DataFrames
txt = readlines("melting.txt")
x = Float64[]
for s in txt
    pos = findfirst(r"[-+]?[0-9,]+\.?[0-9]*", s)
    isnothing(pos) == false && append!(x, parse(Float64, replace(s[pos], "," => "")))
end
df = DataFrame(Matrix(reshape(x, 3, :)'), ["Enthalpy", "Entropy", "MeltingTemperature"])

AWK で製表

以下のような入力ファイル CodeFormat

Code
Format
Example
0
dd-mmm-yyyy HH:MM:SS
07-Sep-2000 15:38:09
1
dd-mmm-yyyy
07-Sep-2000
2
mm/dd/yy
09/07/00
3
mmm
Sep
4
m
S

からデータを読み込み

Code | Format | Example
:--|:--|:--
0 | dd-mmm-yyyy HH:MM:SS | 07-Sep-2000 15:38:09
1 | dd-mmm-yyyy | 07-Sep-2000
2 | mm/dd/yy | 09/07/00
3 | mmm | Sep
4 | m | S

のように整形して出力ファイル table に書き込む。

つまり，最終形の表は3列からなり，各列は | で区切られている。元のデータは3行で表の1行になる。

awk で書くとこんなスクリプト

{
	getline a
	getline b
	printf "%s | %s | %s\n", $0, a, b
}
FNR == 3 {
	print ":--|:--|:--"
}

こんなふうに実行する。

awk -f temp.awk < CodeFormat >table

RStudio で awk を使う

前の記事で ```{awk} ... ``` が動かないと書いた。

https://gedevan-aleksizde.github.io/rmarkdown-cookbook/custom-engine.html

の「15.1 カスタム言語エンジンを登録する」を読んで，以下のように解決した（のかな?）。

最初に一回だけ，以下の R スクリプトを実行する。

knitr::knit_engines$set(awk = function(options) {
  code <- paste(options$code, collapse = '\n')
  out  <- system2('awk', shQuote(code), stdout = TRUE)
  knitr::engine_output(options, code, out)
})

この後，awk によるチャンクを書く。

```{awk}
BEGIN {
    for(i = 1; i <= 10; i++) {
        print i, sqrt(i)
    }
}
```

実行すると，以下の結果を得る。

## 1 1
## 2 1.41421
## 3 1.73205
## 4 2
## 5 2.23607
## 6 2.44949
## 7 2.64575
## 8 2.82843
## 9 3
## 10 3.16228

```{awk}
BEGIN {print "Hello, world!"}
```

## Hello, world!

RStudio で Julia を使う

RStudio で Julia を使う方法を述べる。

とは言っても，RStudio のノートブック作成機能を使うだけである。通常は R 言語を選択するために {r} となっているところを {julia} にするだけだ。

1. まず，RStudio のファイルメニューから，File --> New File --> R Notebook でノートブックウインドウを出す。

2. 以下の例のように，コードブロック ```{julia} ... ``` で挟んで Julia プログラムを書く。

3. 右上にある Run ボタンをクリックし，Run Current Chunk をクリックする。
実行すべきチャンクを選択して shift + command + return でも良い。

4. 実行結果はチャンクの下に書かれる。

{ } 内に書けるものは，R, Python, Julia, Perl, Octave ... 

https://bookdown.org/yihui/rmarkdown/language-engines.html

によれば，以下のものどもが使えるそうだ。

awk, gawk はリストにあるけど動かなかったなあ...

> names(knitr::knit_engines$get())
 [1] "awk"       "bash"      "coffee"    "gawk"      "groovy"    "haskell"  
 [7] "lein"      "mysql"     "node"      "octave"    "perl"      "psql"     
[13] "Rscript"   "ruby"      "sas"       "scala"     "sed"       "sh"       
[19] "stata"     "zsh"       "asis"      "asy"       "block"     "block2"   
[25] "bslib"     "c"         "cat"       "cc"        "comment"   "css"      
[31] "ditaa"     "dot"       "embed"     "exec"      "fortran"   "fortran95"
[37] "go"        "highlight" "js"        "julia"     "python"    "R"        
[43] "Rcpp"      "sass"      "scss"      "sql"       "stan"      "targets"  
[49] "tikz"      "verbatim" 

もっとも，それらが別途インストールされている必要がある。

Octave でトーラスを描いてみた

Octave でトーラスを描いてみた。

Julia に翻訳--219　トーラス，3次元プロット 参照

Julia の range(), repeat() に相応する関数がないので，自前の関数を書く。

function [rng] = range_(begin_, end_, n)
    step = (end_ - begin_) / (n - 1);
    rng = begin_:step:end_;
end

function [res] = repeat(x, n)
    res = x;
    for i = 2:n
        res = [res, x];
    end
end

function torus(n=100, a=2, b=9)
    u = reshape(repeat(range_(0, 2*pi, length=n), n), n, n);
    v = u';
    x = (b + a * cos(u)) * cos(v);
    y = (b + a * cos(u)) * sin(v);
    z = a * sin(u);
    colormap(pink)
    surf(x, y, z)
    pbaspect([4, 4, 1])
end

torus(50, 5, 10)

matgeom パッケージにトーラスを描く関数があった。

標準では緑色なのだが，赤にしてもやはり浮き輪にしか見えない。赤道の位置（?）に明るい線があるからかな。

pkg load matgeom

figure;
drawTorus([50 50 50 30 10 30 45], "facecolor", "red");
axis equal; view([95 10]); light;

「三次元リング」という説明もあったので，金色にしてみたが，安っぽい色だ。真鍮の指輪みたいだな。

pkg load matgeom

figure;
drawTorus([0 0 0 30 10 0 0], "facecolor", [0xff 0xd3 0x26]/255, "nPhi", 240, "nTheta", 120);
axis equal; view([60 30]); light;