裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

職人の技も,数学を応用すべし。

2023年03月16日 | ブログラミング

今,テレビでやってた。クロワッサンを作るときに,最適なバターの割合をx%からはじめてx+1%にしてとやって,20年かかってやっと最適な割合を見つけたんですって。

ご苦労様でした。二分法を知っていれば,もっと速く最適な割合が見つかったんでしょうにねえ。

数学は偉大だよ。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

負の数の平方根(虚数の冪乗)

2023年03月15日 | ブログラミング

-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

 

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

プログラムのお作法

2023年01月03日 | ブログラミング

以下のプログラムは 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

なお,一般的に言えば,オンラインの実行時にコンソールからデータを入力するプログラムは書くべきではない。

 

更に蛇足

本来のエンゲル係数は,貧困度の測定にある。収入が少ないために,必須である食費支出の割合が多くなることを問題とする。飽食の時代では,有り余る収入に贅沢な食生活の割合を測定するものになっているのかなあ?

 

そんなこともあるので,ソースは明示しない。

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia 1.8.4

2022年12月25日 | ブログラミング

Julia 1.8.4 が出ました

Current stable release: v1.8.4 (December 23, 2022)

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Windows の JupyterLab desktop で使用する Julia のバージョンを上げる

2022年12月17日 | ブログラミング

まずは、Julia, JupyterLab desktop のバージョンを上げる

次に、バージョンアップ後の Julia を立ち上げ、

pkg > add IJulia

pkg > build IJulia

そのあと、Juila を抜ける(その状態で作業したければどうぞ)

JupyterLab desktop を立ち上げ、Select Kernel で新しいバージョンを選択できるようになっているはず。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

atan2(0, 0) はどのように評価されるか

2022年12月14日 | ブログラミング

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)

いずれも,結果は "(未定義)" と表示される

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

2.675 を四捨五入して,小数点以下2桁で表示せよ

2022年11月14日 | ブログラミング

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

まあ,完全かどうかは虱潰しにでも調べるしかないか?

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

新 macOS

2022年10月27日 | ブログラミング

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にも追加
  • 標的型攻撃を受けた際にシステムをロックダウンする「ロックダウンモード」を追加

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Python 3.11.0 が出たよ!

2022年10月27日 | ブログラミング

Python 3.11.0 がダウンロード可能になった

Python.org からダウンロードしよう

> Python を使って開発や実験を行うときは、用途に応じて専用の実行環境を作成し、切り替えて使用するのが一般的です。

などというフェークニュースには惑わされないで...

まあ,一般ユーザが複数の実行環境なんで必要ないですわ。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

 0^0 はいくつ?

2022年10月03日 | ブログラミング

 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 となる。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

SymPy で二重根号を外す

2022年09月30日 | ブログラミング

以前の記事で,「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)


コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

年齢の計算

2022年09月19日 | ブログラミング

年齢の計算

生年月日 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)

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で閏年判定

2022年09月19日 | ブログラミング

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

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

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

2022年09月13日 | ブログラミング

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

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

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rump の例題

2022年09月12日 | ブログラミング

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

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村