今，テレビでやってた。クロワッサンを作るときに，最適なバターの割合をx%からはじめてx+1%にしてとやって，20年かかってやっと最適な割合を見つけたんですって。
ご苦労様でした。二分法を知っていれば，もっと速く最適な割合が見つかったんでしょうにねえ。
数学は偉大だよ。
-2 の平方根を求めてみる
Python は (-2)**0.5 で求めることができる。しかし，誤差が表示される。虚数単位は j で表される。
>>> (-2)**0.5
(8.659560562354934e-17+1.4142135623730951j)
numpy.sqrt(-2) ではエラーになる。
>>> import numpy as np
>>> np.sqrt(-2)
<stdin>:1: RuntimeWarning: invalid value encountered in sqrt
nan
虚数として与えるとちゃんと計算する。
>>> np.sqrt(-2+0j)
1.4142135623730951j
R ではどうか。
> options(digits=16)
> (-2)^0.5
[1] NaN
> (-2.0+0i)^0.5
[1] 0+1.414213562373095i
Octave はどうか。
>> (-2)^0.5
ans = 8.659560562354932e-17 + 1.414213562373095e+00i
sqrt() は誤差を含まない。
>> sqrt(-2)
ans = 0 + 1.414213562373095i
Julia では?
julia> (-2)^0.5
ERROR: DomainError with -2.0:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
単にエラーを出すだけではなく，解決方法を提案してくれる。
julia> (-2+0im)^0.5
0.0 + 1.4142135623730951im
sqrt() も，解決方法を提案してくれる。
julia> sqrt(-2)
ERROR: DomainError with -2.0:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
julia> sqrt(Complex(-2))
0.0 + 1.4142135623730951im
以下のプログラムは Julia なのだけど，他の言語でも同じ。
function calc(s1, s2)
# 消費支出 s1 円と飲食費 s2 円からエンゲル係数を計算する。
return s2*100/s1
end
function engel()
s1 = Base.prompt("消費支出? ")
s2 = Base.prompt("飲食費? ")
s1 = parse(Float64, s1)
s2 = parse(Float64, s2)
x1 = calc(s1, s2)
if x1 >= 80
println("エンゲル係数は ", round(Int, x1),
" です。飲食費を使いすぎです。")
else
println("エンゲル係数は ", round(Int, x1), " です。")
end
end
まあそれでもいいのですが，こなれたプログラムならこう書くべき。
つまり，同じプログラム片を何回も(2回でも)書かない。
似たプログラム片，変数名などが出てくると，前に出てきたのと同じなのか，少し違うのかが気にかかる。
同じ部分は繰り返さないようにするし，違うところはそれがはっきりわかるように書く。
ついでに，金額は整数と思うのでそのようにしておく。
function engel()
s1 = Base.prompt("消費支出? ")
s2 = Base.prompt("飲食費? ")
s1 = parse(Int, s1)
s2 = parse(Int, s2)
x1 = calc(s1, s2)
print("エンゲル係数は $(round(Int, x1)) です。")
println(x1 >= 80 ? "飲食費を使いすぎです。" : "")
end
なお，一般的に言えば，オンラインの実行時にコンソールからデータを入力するプログラムは書くべきではない。
更に蛇足
本来のエンゲル係数は，貧困度の測定にある。収入が少ないために，必須である食費支出の割合が多くなることを問題とする。飽食の時代では，有り余る収入に贅沢な食生活の割合を測定するものになっているのかなあ?
そんなこともあるので，ソースは明示しない。
まずは、Julia, JupyterLab desktop のバージョンを上げる
次に、バージョンアップ後の Julia を立ち上げ、
pkg > add IJulia
pkg > build IJulia
そのあと、Juila を抜ける（その状態で作業したければどうぞ）
JupyterLab desktop を立ち上げ、Select Kernel で新しいバージョンを選択できるようになっているはず。
IEEE754 - IEEE Standard for Floating-Point Arithmetic では、
atan2( 0, -0) = π
atan2(-0, -0) = -π
と定義されている。
R はこれに従っている。
> options(digits=16)
> atan2(0, 0)
[1] 0
> atan2(0, -0)
[1] 3.141592653589793
> atan2(-0, -0)
[1] -3.141592653589793
Octave でも R と同じである。
>> format long
>> atan2(0, 0)
ans = 0
>> atan2(0, -0)
ans = 3.141592653589793
>> atan2(-0, -0)
ans = -3.141592653589793
AWK/GAWK も R と同じである。
$ awk "BEGIN{print atan2(0, 0)}"
0
$ awk "BEGIN{print atan2(0, -0)}"
3.14159
$ awk "BEGIN{print atan2(-0, -0)}"
-3.14159
Julia は atan2(x, y) は atan(x, y) である。結果はいずれも 0.0 とされる。
julia> atan(0, 0)
0.0
julia> atan(0, -0)
0.0
julia> atan(-0, -0)
0.0
Python でも Julia と同じく，結果はいずれも 0.0 とされる。
>>> import math
>>> math.atan2(0, 0)
0.0
>>> math.atan2(0, -0)
0.0
>>> math.atan2(-0, -0)
0.0
C++ (Apple clang version 14.0.0 (clang-1400.0.29.202)) でも Julia と同じく，結果はいずれも 0.0 とされる。
#include <stdio.h>
#include <math.h>
int main()
{
printf("%f\n", atan2(0, 0));
printf("%f\n", atan2(0, -0));
printf("%f\n", atan2(-0, -0));
return 0;
}
$ ./a.out
0.000000
0.000000
0.000000
WolframAlpha では未定義とされる。
atan2(0, 0), atan2(0, -0), atan2(-0, -0)
いずれも，結果は "(未定義)" と表示される
2.675 を小数点以下3桁目を四捨五入せよ。
2.68 となることを予想しているのだが，
★ Python
>>> round(2.675, 2)
2.67
で，decimal パッケージの Decimal を使いましょう...ってさ
★ R
> round(2.675, 2)
[1] 2.67
★ Julia
julia> round(2.675, digits=2)
2.68
★ AWK には round 関数はないのだが int(2.675*100+0.5)/100 というような関数を定義すればよい。
$ awk 'BEGIN{print int(2.675*100+0.5)/100}'
2.68
まあ，完全かどうかは虱潰しにでも調べるしかないか?
macOS Ventura バージョン 13.0 (22A380)
になったよ！！！(10月25日）
- iPhoneをWebカメラ代わりに使える「連係カメラ」
- 複数のアプリやウィンドウをまとめてグループ化する「ステージマネージャ」
- メールアプリでは、一度送信操作を行なった直後の取り消しや送信予約
- 「FaceTime」がHandoff機能に対応し、MacからiPhoneやiPadへシームレスに通話を引き継げるようになった
- 「Safari」に共有タブグループ機能を追加。ほかのユーザーと閲覧中のタブを共有できる
- WebサイトへのログインにFace IDやTouch IDが利用できる「パスキー」機能を追加(Webサイト側の対応が必要)
- 「メッセージ」アプリでSharePlay連携が可能になったほか、メッセージの編集や送信取り消し、削除したメッセージの復元、未読状態への復元機能などを追加 ]
- 「iCloud」に共有写真ライブラリ機能を追加。5人まで招待でき、個々のライブラリを統合できる
- 「Spotlight」がクイックルックに対応。検索結果の上でスペースキーを押すことでファイルをプレビューできる様になった。また画像検索がWebとローカル、メッセージ、メモ、Finderを横断するようになった
- 一時停止した動画の画面内にあるテキストを認識する機能が日本語と韓国語に対応
- 画像認識が人間に加え新たに動物、鳥、昆虫、彫像、ランドマークを認識するようになった ・ライブキャプションなどアクセシビリティ関連ツールを強化
- iOSに搭載している「天気」「時計」アプリがMacにも追加
- 標的型攻撃を受けた際にシステムをロックダウンする「ロックダウンモード」を追加
Python 3.11.0 がダウンロード可能になった
Python.org からダウンロードしよう
> Python を使って開発や実験を行うときは、用途に応じて専用の実行環境を作成し、切り替えて使用するのが一般的です。
などというフェークニュースには惑わされないで...
まあ，一般ユーザが複数の実行環境なんで必要ないですわ。
0^0 はいくつ？
インターネット上にも多くのページがある。
多くのプログラミング言語では 0^0 = 1 とされている。Julia も同じである。
0^0
1
図を描いてみよう。
using Plots
x = 0:0.001:2
plot(x, x.^x, label="", size=(300, 200))
正の方向から 0 に近づけていくと確かに 1 に近づいていくようだ。
x = 1e-15
x^x
0.9999999999999655
x = big"1" / big"100000000000000000000000000000000000000000000000000000000000000"
x^x
0.999999999999999999999999999999999999999999999999999999999998572397242343691675
SymPy の limit() で見てみよう。
using SymPy
@syms x
limit(x^x, x, 0, dir="+") # dir="+" はデフォルト
1
確かに 1 になる。
ちなみに，x^x が最小値を取るときの x はなんだろう。
f = diff(x^x) # 微分して
x^x(log(x) + 1)
x0 = solve(f)[1] # 傾きが 0 になる x0 を求める。
e^(-1)
x0.evalf()
0.367879441171442
x0^x0.evalf()
0.692200627555346
x が exp(-1) ≒ 0.36788 のとき最小値 0.69220 となる。
以前の記事で，「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)
年齢の計算
生年月日 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)
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
Julia の SymPy は，Python の sympy を利用している。Python の整数型は必要に応じて長精度整数になる。
求める式には実数が含まれるが，この式を 4 倍することで整式になる。
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
PVアクセスランキング にほんブログ村
