互いに素な数で割った余りがそれぞれ、であるようなnはを法として一意に存在します。(中国の剰余定理)解nは以下の式で求めることができます。ここで、、、はを法とするの逆元です。すなわち、を満たします。 これをMaximaで素直に書き下してみましょう。 p :…

Newton法を使って逆数を求めるには、以下の漸化式を使います。これは多項式にも適用することができます。具体的な例として、Bernouli数を計算してみます。Bernouli数は以下の式で定義される数列です。微分してt=0を代入すると不定になるので代数的に求めるこ…

例えば、8ビット*8ビットの乗算器を持つマイコンを想像してください。このマイコンで2つのmバイトの整数の掛け算を行うには、普通、回、乗算器を使います。 この掛け算を、のオーダーの乗算回数で行うアルゴリズムを考えます。 一般化して整数x, yを基数wで…

Maximaには素数テストの組み込み関数primepがあります。これはMiller-Rabinの確率的アルゴリズムを用いたものです。 決定的アルゴリズムであるLucasの素数テストを実装して、の素数性を証明しましょう。 /* 4.3.15m */ p : (10^23 - 1)/9$ factors : ifactor…

6番目のFermat数の因数を楕円曲線法で探してみましょう。 /* 4.3.13m */ /* 因数発見を示すグローバルスイッチ */ notFoundYet : true$ /* Maximaの組み込み関数gcdexは、最大公約数の符号が不定なので、正となるように規格化する関数を用意。 */ extendedGC…

Maximaは素因数分解を実行する関数が用意されています。使ってみましょう。 /* 4.3.1m */ factor(54444439); 使用した機能 factor マニュアルによると、factorはデフォルトでは、Berlekampのアルゴリズムというものを使うそうです。 組み込み関数を使う代わ…

複数の極を持つがそれ以外では正則な函数の単純閉曲線上の積分は留数定理を使って計算することができます。 留数定理を活用した複素積分の計算をMaximaでやってみましょう。 例は、0" align="middle" />を以下の経路で積分した場合です。 外周の積分は極限で…

で表される曲面の各点の平均曲率Cは以下のように定義することができます。 球面の平均曲率を計算してみましょう。 /* 4.2.13m */ assume(r > 0)$ f(x, y) := sqrt(r^2 - x^2 - y^2)$ fx : diff(f(x, y), x); fy : diff(f(x, y), y)$ den : sqrt(1 + fx^2 + f…

を導出します。 以下のように周期2の函数を定義します。は実定数です。 のフーリエ変換を求めてみましょう。 /* 4.2.6m, 4.2.7m, 4.2.8m, 4.2.9m */ load(fourie) f(x) := a*x^2 + b*x^4$ f(x) = totalfourier(f(x), x, 1); 使用した機能 totalfourier を代…

Jacobiの三番目のテータ定数は以下のように定義されます。(はqの函数ですが、テータ函数の原点の値であることから定数と呼ばれます。) この函数に関して、以下の恒等式が成り立つことが知られています。これを有限項に限定して確認してみましょう。まず、右…

円周等分多項式は、の複素数解のうち、n乗すると初めて1になる解を根に持つ既約多項式です。正確な定義は以下の式になります。円周等分多項式をMaximaでいくつか書き出してみましょう。 /* 4.1.24m */ cyclotomic(k, z) := trigrat(product(if gcd(j, k) = 1…

整数全体に渡る以下の和を計算します。函数に以下のPoisson和公式を適用するところから始めます。 最後の積分をMaximaに解かせてみましょう。 /* 4.1.10m */ term : integrate(exp(-y)*cos(m*y), y, 0, inf); 使用した機能 integrate inf 従って、以下が成り…

Legendre多項式は、以下の母函数を展開することで得られます。この性質から、Legendre多項式は双極子ポテンシャルを距離の逆数で展開する際に現れます。早速展開してみましょう。 /* 4.1.2m */ legendreGen : taylor((1 - 2*z*t + t^2)^(-1/2), t, 0, 4); 使…

三つ葉結び目を描いてみましょう。 まずは2次元平面への投影図です。 /* 3.6.4m */ a : 1.0$ b : 0.4$ r(t) := a + b*cos(3*t/2)$ plot2d([parametric, r(t)*cos(t), r(t)*sin(t), [t, 0, 4*%pi + 0.1], [nticks, 200]], [x, -2.0, 2.0], [y, -1.5, 1.5])$ …

パラメトリックプロットの例として、等角螺旋を描いてみましょう。 /* 3.5.1m */ r(t) := exp(-0.1*t)$ plot2d([parametric, r(t)*cos(t), r(t)*sin(t), [t, 0, 30], [nticks, 500]], [x, -1.5, 1.5], [y, -1, 1]); 螺旋状の点と原点を結んだ線とその点の接…

正五角錐を回転させるアニメーションを製作しましょう。 /* 3.4.1m.mac */ load(draw)$ xyz2polar(v) := block([r : sqrt(v[1]^2 + v[2]^2 + v[3]^2)], [r, acos(v[3]/r), mod(atan2(v[2], v[1]), 2*%pi)])$ t : makelist(ev([cos(a), sin(a), 0], a = 2*%pi…

の主値について3Dプロットしてみましょう。arctanはtanの逆函数です。tanは周期函数ですからarctanは多価函数になります。arctanの主値はからまでと定義されます。 /* 3.3.2m */ plot3d(atan(y/x), [x, -1.5, 1.5], [y, -1.5, 1.5]); 使用した機能 atan は座…

函数をプロットしてその様子を見る時には、プロットする解像度が重要です。 シンク函数を例に取ります。これは、矩形函数(パルス函数)をFourier変換すると得られる函数です。一次元のFraunhofer回折像の波動関数でもあります。 まずシンク函数を定義しましょ…

2次元、3次元グラフィックスの利用例として、ガンマ函数を可視化してみましょう。 ガンマ函数は階乗の複素数への拡張です。以下のように定義される解析函数です。が成り立ちます。 ガンマ函数は実数上では実数値を取ります。グラフを表示してみましょう。 /*…

東北地方太平洋沖地震に被災された皆様へお見舞い申し上げます。 やるべきことが見えるまでは日常をしっかり続けることと自分に言い聞かせて。 Mathematica―理工系ツールとしての (アジソン ウェスレイ・トッパン情報科学シリーズ)作者: R.E.クランドール,伊…