あちらの「はじめての Maxima プログラミング」に若干の補足説明を追加したもの。Maxima には構造化プログラミングの機能もあり,条件分岐や繰り返し処理などを利用することで,プログラミング言語としても利用できます。
エンゲル係数を求める例
エンゲル係数とは,家計の消費支出に占める飲食費(食料費,食費)の割合(パーセント単位)です。
消費支出 110,000円,飲食費 37,000円の人のエンゲル係数を求める例。Maxima では以下のような電卓的な使い方ですぐ答えが出ます。
本稿のように Jupyter Notebook 環境では,以下の1行を書いたあとに,Shift キーを押しながら Enter キー(または Return キー)を押して実行します。
Maxima の文末は,;
か $
で終わります。$
で終わる場合は結果が表示されません。
37000/110000 * 100;
Maxima では割り算の結果は厳密な有理数として出力されます。(近似的な)実数(小数点数)で表示させたい場合は,とりあえずは数値のどれか1つにでも小数点をつけます。(あとで float()
関数で実数に変換することも説明します。)
37000/110000 * 100.0;
ここからは,単に電卓的に使うだけではなく,応用・発展も見据えて以下のようなプログラム文を作成して実行することにします。
以下の例では,
/*
と*/
で囲まれる部分はコメント。コメントはそれを読む人間のための注釈であり,プログラムの実行時には無視されます。- 3行目で変数
spending
に消費支出である110000
(円)という数値を代入し, - 4行目で変数
food
に飲食費である37000
(円)という数値を代入し, - 5行目でエンゲル係数を計算して,変数
eng
に代入し,表示させています。
Maxima では,変数への代入は :
です。
/* エンゲル係数を求める
消費支出 spending, 飲食費 food, エンゲル係数 eng */
spending: 110000$
food: 37000$
eng: food / spending * 100.0;
エラーとデバグ
意に反して,プログラムがうまく動作しない場合があります。たとえば,上記の 3行目を以下のように書いた場合。
掛け算を示す演算子は *
であるところを誤って,x
と書くと,これらの記号は Maxima では演算子として定義されていないのでエラーとなります。このようなエラーは文法エラー(incorrect syntax) と呼ばれ,Jupyter Notebook では,以下のように誤りやその箇所を指摘してくれます。
spending: 110000$
food: 37000$
eng: food x 100 / spending;
別の種類の誤りもあります。たとえば,パーセントにするためにかける「100.0」を誤って「1000.0」と入力した場合…
spending: 110000$
food: 37000$
eng: food / spending * 1000.0;
上記の例では incorrect syntax
(文法エラー)の警告が出ませんが,答えは明らかに変です。(エンゲル係数の最大値は100だから。)
このような場合はコンピュータにとって誤りではなく,指示通りに計算し,その結果を表示します。
このようなプログラム作成者の不注意や勘違いによって生じるエラーには注意が必要です。
プログラムが正しく動作しない場合は,誤りの箇所を探し,修正して再度動作を確認します。正しく動作するようになるまで,プログラムを修正し,実行(Shift + Enter)します。この作業をデバグといいます。
データの型と変数
上記の例の 1行目では,spending という名前の変数に 110000 という 数値を代入しています。変数名はアルファベットで始まる英数字列にします。大文字と小文字は区別されます。
データ型には数値(整数・有理数・浮動小数点数)と文字列などがあります。
変数 moji に文字列を代入する場合は,以下のように ” で囲みます。
moji: "エンゲル係数"$
print(moji)$
式と演算
主な演算子
足し算は +
,
引き算は -
,掛け算は *
,割り算は /
です。
優先的に計算したい箇所は丸括弧 ( )
で グループ化します。
2 + 50 - 5 * 6 /2;
2 + (50 - 5 * 6)/2;
文字列の連結
文字列の連結演算は concat()
を使います。以下の例を参照。
mo: "エンゲル"$
ji: "係数"$
concat(mo, ji);
文字列と数値を直接連結するときも,concat()
を使います。
eng: 37000 / 110000 * 100.0$
output: concat("エンゲル係数は ", eng, " です。")$
print(output)$
文字列と連結しなくても,以下のようにも書くことができます。
print()
文は,カンマ (,
) で区切ることによって複数の項目(変数)を表示することができます。
print("エンゲル係数は ", eng, " です。")$
参考:整数を3桁0詰めで表示する例
いわゆる %03d
のように,整数を3桁,0詰めで表示する需要があるかもしれません。後述する書式指定でできそうにないので,以下のような奇策を考えてみました。
concat(i)
で整数i
を文字列にし,substring( , 2)
で2文字目からの文字列にしてprint()
for i:1000 thru 1003 do(
print(substring(concat(i), 2))
)$
割り算の商と余り
Maxima の割り算は以下のように厳密な値(有理数のまま)を返します。
17/3;
34/6; /* 約分してくれる */
割り算の「商」と「余り」を求める場合は,quotient()
関数や mod()
関数を使います。
/* 商 */
quotient(17, 3);
/* 余り */
mod(17, 3);
/* 検算 */
quotient(17, 3) * 3 + mod(17, 3);
/* 別途述べる,切り捨てする関数 floor() を使う例 */
/* 商 */
floor(17/3);
/* 余り */
17 - floor(17/3) * 3;
浮動小数点数表示する場合は,すでに述べたように数値に小数点をつけて整数でないことをアピールするか,浮動小数点数に変換する関数 float()
を使います。
17.0/3;
float(17/3);
べき乗は **
(または ^
)です。以下の例では $2^4 = 2 \times 2 \times 2 \times 2$ と $2^{1/2} = \sqrt{2}$ の値を表示します。
2**4;
2**(1/2);
float(%);
練習 1
1 天文単位(1 au)の距離を光速 c で進むのに要する時間は何分何秒か。ただし,
- 1 au は 12 桁の定義定数(誤差無し)で
149597870700
m - c は 9 桁の定義定数(誤差無し)で
299792458
m/s である。
参考:
対話型プログラム
上記の例では,プログラムの中で消費支出や飲食費を決めていましたが,今度はプログラムを実行する人が入力できるようにしてみます。
read()
関数は入力された値を返します。
注意: Maxima では read()
に対して値をキーボードから入力する場合も,最後は ;
をつけてから Enter キーを押します。
spending: read("消費支出? (数値の後に ; をつけて Enter キー)")$
food: read("飲食費? (数値の後に ; をつけて Enter キー)")$
eng: float(food / spending * 100)$ /* float() で変換 */
print("エンゲル係数は ", eng, " です。")$
関数の定義と呼び出し
(これまでの例では1回しか計算していませんが)よく使う計算を「関数」として定義する例です。
以下の例では,1行目でエンゲル係数を計算する関数 calc()
を定義し,5行目でその関数を呼び出しています。
忘れないで!: Maxima では read()
に対して値をキーボードから入力する場合も,最後は ;
をつけてから Enter キーを押します。
calc(spending, food):= float(food / spending * 100)$
spending: read("消費支出? ")$
food: read("飲食費? ")$
/* 上で定義した関数 calc() を使う */
eng: calc(spending, food)$
print("エンゲル係数は ", eng, " です。")$
いったん定義された関数は,以下のように何回でも利用できます。
calc(120000, 53000);
練習 2
消費支出を入力すると,エンゲル係数が 20 ~ 40 となる飲食費を表示するプログラムをつくりなさい。(以下の例はエンゲル係数が 20 以上となる飲食費を表示します。編集して完成させなさい。)
spending: read("消費支出? ")$
food20: spending * 0.2$
print("エンゲル係数が 20 ~ 40 となる飲食費は...",
food20, "円〜")$
定義定数や数学関数を使う
Maxima には,$\sin x, \cos x, \tan x$などの数学関数が組み込まれています。また,数学定数としては円周率 %pi
,虚数単位 %i
,ネイピア数 %e
,無限大 inf
などが使えます。
/* 円周率 π */
%pi;
sin(%pi/2);
cos(%pi/3);
tan(%pi/4);
/* 虚数単位の2乗は... */
%i**2;
/* log は自然対数 */
log(%e);
以下の例では,x
以下の最大の整数を与える関数 floor(x)
を使って,小数点以下を切り捨てたエンゲル係数の値を表示させています。
また,エンゲル係数を計算する部分はすでに関数 calc()
として定義していますが,プログラム全体を関数 engel()
として定義してしまいます。
以下の例のように,block()
は複数の関数や式を ,
(カンマ)で区切って記述します。
engel():= block(
spending: read("消費支出? "),
food: read("飲食費? "),
/* 上で定義された calc() を使用 */
eng: calc(spending, food),
print("エンゲル係数は ", floor(eng),
" です。(小数点以下切り捨て)")
)$
engel()$
Maxima で使える数学関数
Maxima で使える数学関数や数学定数については,以下のマニュアルのページを参照してください。
練習 3
小数点以下を切り上げたエンゲル係数の値を表示させるように上記のプログラムを変更しなさい。
ヒント: Maxima 5.42.2 Manual: 10. Mathematical Functions の ceiling()
を参照。
参考:数値の丸めと四捨五入
Maxima の組み込み関数に round()
があります。数値を丸める関数ですが,いわゆる四捨五入とはちょっと違うようです。
以下の例からわかるように,小数点以下が 0.5
未満の場合と 0.5
を超える場合は四捨五入と同様の丸め方ですが,ちょうど 0.5
のときは,偶数に丸められていることがわかります。
これを偶数への丸めといい,端数が 0.5
より小さいなら切り捨て,端数が0.5
より大きいなら切り上げ,端数がちょうど 0.5
なら切り捨てと切り上げのうち結果が偶数となる方へ丸めます。
print(round(12.49), " ", round(12.50), " ", round(12.51))$
print(round(13.49), " ", round(13.50), " ", round(13.51))$
条件分岐
条件分岐とは,例えばエンゲル係数の値によって実行文を変えることです。例を示します。
関数 engel()
を以下のように変更します。
/* 忘れてるかもしれないのであらためて定義しておく */
calc(spending, food):= float(food / spending * 100)$
engel():= block(
spending : read("消費支出? "),
food : read("飲食費? "),
/* 上で定義した関数 calc() を使う */
eng : calc(spending, food),
if eng >= 40 then
print("エンゲル係数は ", round(eng),
" です。高い値です。")
else
print("エンゲル係数は ", round(eng), " です。")
)$
engel()$
engel()$
if
文の書式は以下の通りです。
if 条件式1 then
条件式1が満たされる場合に実行する文
elseif 条件式2 then
条件式2が満たされる場合に実行する文
else
上記以外の場合に実行する文 $
関係演算子
上の例では,eng >= 40
つまり eng
が 40
以上のとき,という関係演算子を使いました。他の例は以下のとおりです。
a < b
a
が b
より小さい
a <= b
a
が b
以下
a > b
a
が b
より大きい
a >= b
a
が b
以上
is(equal(a, b))
は a
と b
が等しいとき true, そうでないとき false を返します。equal(a, b)
の否定は notequal(a, b)
です。
練習 4
上記の条件分岐を,エンゲル係数が
- 40 以上の場合には「高い値です。」
- 20 未満の場合には「低い値です。」
と表示させるように変更しなさい。
繰り返し処理
while
以下の例では,1行目で代入された消費支出に対してエンゲル係数を計算します。
1回だけ計算して表示するのではなく,2行目の飲食費の初期値(ここでは 30000
円)から,5行目の条件式(ここでは 50000
円以下)が成り立つ場合に,8行目にあるように 5000
円ずつ増やしながらエンゲル係数を表示します。
spending : 110000 $
food : 30000 $
print("消費支出", spending, "円に対するエンゲル係数の値")$
while food <= 50000 do(
eng : calc(spending, food),
print(" 飲食費", food, "円の場合: ", eng),
food : food + 5000
)$
while
文の書式は以下の通りです。
while 条件式 do(
条件式が成り立つ場合に実行する式1,
条件式が成り立つ場合に実行する式2, ... ,
条件式が成り立つ場合に実行する式
)
for
同様の繰り返し処理を for
文を使って行うこともできます。
spending: 110000 $
food: 30000 $
print("消費支出", spending, "円に対するエンゲル係数の値")$
for i: 30000 thru 50000 step 5000 do(
eng: calc(spending, food),
print(" 飲食費", food, "円の場合: ", eng),
food: food + 5000
)$
for
文の書式は以下の通りです。
for ivar: istart thru iend step istep
do(
ivar が istart から iend の間,実行される文
)$
ivar
は整数変数で,その初期値は istart
。文が実行されるごとに,ivar
の値は istep
だけ加算され,iend
の値になるとループが終了します。
step
の部分を省略すると,1
ずつ加算されます。
例:数列の和
$\displaystyle \sum_{n = 1}^{10} n$ を繰り返し処理の練習として求めてみます。
/* Maxima の関数 sum() を使えばこの通りだが... */
sum(n, n, 1, 10);
/* 一般の n までの和を求める例 */
nusum(i, i, 1, n);
/* 繰り返し処理 while の練習として */
answer: 0$
i: 1$
while i <= 10 do(
answer: answer + i,
i: i + 1
)$
answer;
/* 繰り返し処理 for の練習として */
answer: 0$
for i: 1 thru 10 do(
answer: answer + i
)$
answer;
参考:再帰的定義関数と繰り返し処理
参考までに,再帰的定義関数を使った方法について。
$i=1$ から $i = n$ までの $i$ の和 $\displaystyle \sum_{i=1}^n i$ を求める関数 mysum(n)
を以下のように定義する。
\begin{eqnarray}
\mbox{mysum}(1) &=& 1 \\
\mbox{mysum}(2) &=& 1 + 2 = \mbox{mysum}(1) + 2 \\
\mbox{mysum}(3) &=& 1 + 2 + 3 = \mbox{mysum}(2) + 3 \\
&\vdots& \\
\mbox{mysum}(n) &=& \mbox{mysum}(n-1) + n
\end{eqnarray}
/* 自分自身の一つ前の値を使って自分自身を定義する */
/* このように定義された関数を再帰的定義関数という */
mysum(n):= if n = 1 then 1 else mysum(n-1) + n$
mysum(10);
練習 5
以下の値を求める関数を mysum3()
として作成し,$n = 10$ および $n = 100$ の場合を出力せよ。
$$\sum_{i = 1}^{n} i^3$$
/* Maxima の関数を使うと... */
sum(i**3, i, 1, 10);
sum(i**3, i, 1, 100);
/* だけど,自分で作ってみるんですよ */
リストとファイルの読み書き
リスト
要素が5個のリストの各成分に数値を入れたり表示したりする例です。
a: [5, 4, 3, 2, 1]$
print("全成分を一挙に表示")$
print(a)$
少し表示を工夫してみます。sprint()
関数は(改行せずに)横に続けて出力します。
for c in a do()$
はリスト a
の要素を c
として順番にとります。
print("全成分を横に表示")$
for i in a do(sprint(i))$
print("全成分を横に表示")$
for i in a do(print(i))$
5個の成分を全てゼロにしてリストを初期化する例です。
/* リスト b の5つの要素を全て 0 に*/
b : makelist(0, 5)$
print(b)$
/* 3番目の要素の値を設定 */
b[3] : 3$
print(b)$
ファイルからリストへの読み込み
以下のようなテキストファイル mydata.txt
があるとします。ない場合は「新規」の「テキストファイル」で作成します。
$ cat mydata.txt
1
2
3
4
5
6
ファイル mydata.txt
の数値をリスト a
に読み込みます。
a: read_list("mydata.txt");
/* 要素の総和,平均,最大値,最小値を出力します。*/
print("要素の総和は", lsum(i, i, a))$
print("要素の平均は", float(lsum(i, i, a)/length(a)))$
print("最大値は ", lmax(a))$
print("最小値は ", lmin(a))$
次は,以下のような「6行2列」のデータが入ったテキストファイル mydata2.txt
を読み込む例です。
$ cat mydata2.txt
1 1
2 4
3 9
4 16
5 25
6 36
b: read_nested_list("mydata2.txt")$
print(b)$
/* 入れ子になっているリストの成分を取り出す方法 */
b[2];
b[2][2];
リストではなく,行列として読み込む場合は以下のようにします。
M: read_matrix("mydata2.txt");
/* 行列として読み込んだ場合の成分の取り出し方 */
M[2, 2];
ファイルへの書き込み
ファイルへの書き込み例です。この例では,$𝑥$ と $\sin 𝑥$ の値を $𝑥=0.1$ から $𝑥=1.0$ まで書き込んでいます。
x: makelist([0.1*i, sin(0.1*i)], i, 1, 10)$
write_data(x, "myoutput.txt")$
JupyterHub のホームに「myoutput.txt」が作成されていることを確認します。ファイル名をクリックすると,内容を確認できます。
0.1 0.09983341664682815
0.2 0.1986693307950612
0.3 0.2955202066613396
0.4 0.3894183423086505
0.5 0.479425538604203
0.6000000000000001 0.5646424733950355
0.7000000000000001 0.6442176872376911
0.8 0.7173560908995228
0.9 0.7833269096274834
1.0 0.8414709848078965
書式を指定してファイルへ書き込み
上記の write_data()
では数値の桁数に凸凹があります。後で説明する printf()
を使って書式設定してファイルへ書き込む例です。
/* 書き込み用に myoutput2.txt をオープンする */
s: openw("myoutput2.txt")$
for i:1 thru 10 do(
/* 小数点以下3桁 小数点以下15桁 */
printf(s, "~,3f ~,15f~%", 0.1*i, sin(0.1*i))
)$
/* 書き込みが終わったらクローズする */
close(s)$
JupyterHub のホームに「myoutput2.txt」が作成されていることを確認します。ファイル名をクリックすると,内容を確認できます。
0.100 0.099833416646828
0.200 0.198669330795061
0.300 0.295520206661340
0.400 0.389418342308651
0.500 0.479425538604203
0.600 0.564642473395035
0.700 0.644217687237691
0.800 0.717356090899523
0.900 0.783326909627483
1.000 0.841470984807897
readline()
でファイルから1行ずつ読み込み
以下のようにして,ファイルから1行ずつ読み込んで表示させることができます。ファイルの最後まで readline()
で1行ずつ読み込み,読み込む文字列がなくなったら終了します。
/* 読み込み用に myoutput2.txt をオープンする */
s: openr("myoutput2.txt")$
/* ファイルの最後まで1行ずつ読み込み表示する */
while stringp(tmp:readline(s)) do print(tmp)$
/* 読み込みが終わったらクローズする */
close(s)$
リスト・行列の相互変換
行列からリストへの変換は args()
関数で行います。
M;
m: args(M);
また,リストから行列への変換は以下のようにします。
b;
M1: apply('matrix, b);
書式指定
浮動小数点数の書式指定 "~f"
実数を表示する際,全体の文字数や小数点以下の桁数を指定して表示させる例です。
例えば "~8,5f"
は全体で 8 桁の表示範囲をとり,小数点以下は 5 桁表示します。"~%"
は改行を表します。
(小数点が ,
になっているのはなんとなく… と思って調べると,Wikipedia 日本語版では .
「ピリオド」(イギリス式)で ,
「コンマ」(フランス式)とある。フランスだけでなく非英語圏では「コンマ」がよく使われるらしい。)
sqrtwo : float(sqrt(2))$
printf(true, "12345678901234567890 桁数~%")$
printf(true, "~f~%", sqrtwo)$
printf(true, "~8,5f~%", sqrtwo)$
printf(true, "~18,15f~%", sqrtwo)$
printf(true, "2の平方根は ~18,15f です~%", sqrtwo)$
全体の桁数のみを指定することもできます。例えば,"~7f"
のように全体の桁数 7 桁のみを指定します。また,小数点以下の桁数のみを"~,2f"
のように指定することもできます。
printf(true, "12345678901234567890 桁数~%")$
printf(true, "~7f~%", sqrtwo)$
printf(true, "~,2f~%", sqrtwo)$
fpprintprec
による表示桁数設定
また,以上のように個別に指定しなくても,
fpprintprec: 5$
のようにすると,その後の小数表示の桁数が(この場合であれば)5桁になります。表示桁数のみの変更であり,計算の精度等には影響しません。既定に戻すには
fpprintprec: 0$
fpprintprec: 5$
sqrtwo;
fpprintprec: 0$
sqrtwo;
練習 6
以下の出力を,エンゲル係数の値を(四捨五入して)小数点以下1桁まで表示するように変更しなさい。
spending: 110000$
food: 37000$
eng: float(food * 100 / spending)$
print("エンゲル係数は ", eng, " です。")$
指数表式の指定 "~e"
非常に大きい数や小さい数を表示するときには, $1.9891 \times 10^{30}$
や $6.67430 \times 10^{-11}$
などのような指数表示をします。
例えば "~9,2e"
では,全体で 9 桁の表示分を確保し,小数点以下は 2 桁表示します。
c: 299792458$
printf(true, "12345678901234567890 桁数~%")$
printf(true, "~9,2e~%", c)$
printf(true, "~15,8e~%", c)$
printf(true, "光速は ~15,8e m/s です。~%", c)$
全体の桁数のみを指定したり,小数点以下の桁数のみを指定することもできます。
printf(true, "12345678901234567890 桁数~%")$
printf(true, "~,2e~%", c)$
printf(true, "~15e~%", c)$
"~f"
または "~e"
となる "~g"
"~g"
は数値の大きさによって,"~f"
または "~e"
に(だいたい)なります。
printf(true, "~,2g~%", sqrtwo)$ /* 総桁数が 2 になる? */
printf(true, "~,2f~%", sqrtwo)$
printf(true, "~,2g~%", c)$
printf(true, "~,2e~%", c)$
文字型の書式指定と substring()
文字列の書式指定は "~a"
です。"~a"
では文字数の設定ができないようなので,以下のように substring()
を使ってみます。
substring()
で文字列の一部を取り出します。
substring(str, i, j)
は文字列 str
の i
番目から始まり,j
番目以降を切り取った(つまり,i
番目から j-1
番目までの)文字列を与えます。
ten: "1234567890"$
printf(true, "~a ~a~%", substring(ten, 1, 6),
substring(ten, 1, 6))$
printf(true, "~a ~a~%", substring(ten, 1, 7),
substring(ten, 1, 5))$
整数の書式指定
整数の書式指定は "~d"
です。"~5d"
のように桁数を指定することもできます。
v: [1, 23, 456, 7890]$
/* v の各要素を整数5桁右詰めで表示 */
for c in v do printf(true, "~5d~%", c)$
参考:整数を 0 詰めで表示する策
いわゆる %05d
のように 0 詰めで表示する機能はないようなので,苦肉の策:
/* v の各要素を整数5桁 0 詰めで表示する */
for c in v do(
printf(true, "~a~%", substring(concat(c+100000), 2))
)$
m: [[1, 0, 0, 0],
[0, 23, 0, 0],
[0, 0, 456, 0],
[0, 0, 0, 7890]]$
/* sprint() は改行せずに出力。newline() は改行のみ。*/
for i in m do(
for j in i do sprint(j),
newline()
)$
/* 書式指定で表示桁数を6桁とり,右寄せで各成分を表示 */
for i in m do(
for j in i do printf(true, "~6d", j),
newline()
)$
参考:printf(true, )
の true
以外の使い方
書式指定で出力する際,printf(true, )
としているが,true
以外にも例えば false
とかもあるのか,について。
すでに,以下のように true
のところに書き込み用のファイルの設定を書く例は示している。
s: openw("myoutput2.txt")$
for i:1 thru 10 do(
printf(s, "~,3f ~,15f~%", 0.1*i, sin(0.1*i))
)$
close(s)$
以下は,printf(false, )
の例。false
と書くと printf(false, )
は出力を含む文字列を返すので,文字列を別の変数(たとえば tmp
)に代入したいときに使ったらどうでしょうか。
参考:0詰め連番ファイル名の作成
以下は 0 詰め連番ファイル名の作成例。
for i: 0 thru 3 do(
tmp: printf(false, "~4d", 1000+i),
filename: concat("fig", substring(tmp, 2), ".png"),
print(filename)
)$