Corollaryは必然に。

このブログは「コロちゃんぬ」の提供でお送りします

恵方を表す関数を求めてみた

恵方巻といえば、節分の日に決まった方角(恵方)を向いて無言で食べると良いとされる巻き寿司のことです。

恵方は毎年変わり、以下のようにして決まるそうです:

西暦年の1の位 恵方
24方位 方位角*1 16方位 東基準反時計周り
4・9 075° 東北東やや東 015°
0・5 255° 西南西やや西 195°
1・6
3・8
165° 南南東やや南 285°
2・7 345° 北北西やや北 105°
歳徳神 - Wikipediaより一部引用)


例えば2021年の恵方「南南東やや南」ですね。
f:id:corollary2525:20210129043415p:plain

16方位だと「南寄りの南東のやや南(7.5°)」と聞こえてややこしいので、個人的には「南むいて15°左に回転」が分かりやすいと思います(他の恵方も同様)。



それはさておき、今年の恵方を計算で求められたら便利ですよね。毎年毎年「恵方 方角」で検索せずに済みますし。

ということで今回は、西暦nに対してその年の恵方(角度)を返す関数を恵方関数」と呼ぶことにし、その恵方関数\theta(n)を求めていこうと思います。

素朴な恵方関数

xy平面で考えることにして、東は0^{\circ},北は90^{\circ},西は180^{\circ},南は270^{\circ}で表すことにします(東基準反時計回り)。弧度法がよければ90^{\circ}=\frac{\pi}{2}としてください。

恵方はそれぞれ東,西,南,北から反時計回りに15^{\circ}ずれているので、まずは恵方をざっくりと「東,西,南,北」で考えます。

表をみると、西暦nを5で割った余りが0,1,2,3,4のとき,恵方は順にざっくりと西,南,北,南,東と変化します。よって

\begin{equation}
a_n=
\begin{cases}
2&(n\equiv0\pmod5)\\
3&(n\equiv1\pmod5)\\
1&(n\equiv2\pmod5)\\
3&(n\equiv3\pmod5)\\
0&(n\equiv4\pmod5)
\end{cases}\label{a}
\end{equation}

とおけば,a_n\times90^{\circ}で「ざっくりと恵方関数」ができます。これに15^{\circ}(弧度法だと\frac{\pi}{12})を足した角度の列
\begin{equation*}
\theta(n) = a_n\times90^{\circ} + 15^{\circ}
\end{equation*}はn=0から整数を順に代入したときに西南西やや西,南南東やや南,北北西やや北,南南東やや南,東北東やや東を繰り返します。ということで、求めていた恵方関数が得られました。





ウーン…(´・ω・`)

しかし、a_nは半ば都合のいいように値を指定しており、ぶっちゃけ「東、北、西、南」を「ぜろ、いち、に、さん」に置き換えた数列に過ぎないとも思えるので、やや面白みに欠けます。

場合分けしないで済むもっといい表示はないだろうか…?


ラグランジュ補間で

解決策の一つとして,\eqref{a}の代わりに
\begin{equation}
\begin{cases}
f(0)=2\\
f(1)=3\\
f(2)=1\\
f(3)=3\\
f(4)=0
\end{cases}\label{f}
\end{equation}をみたす関数を見つける方法があります。そのような関数が見つかれば、西暦nを5で割った余りを\operatorname{Mod}(n,5)とおいて

\begin{equation}
\theta(n) = f(\operatorname{Mod}(n,5))\times90^{\circ} + 15^{\circ}
\end{equation}

とおけば恵方関数が得られます(ちなみに\operatorname{Mod}(x,m)はDesmos,GeoGebraなどの計算アプリで用意されている関数です)。

これなら場合分けして定義する必要がなくなりますね!

\eqref{f}をみたす関数はいくらでもありますが、今回は4次多項式で見つけようと思います。多項式「多項」は「多幸」とも読めるので縁起がいいですし(?)
ということで
\begin{equation*}
f(x) = c_0 + c_1x + c_2x^2 + c_3x^3 + c_4x^4
\end{equation*}とおくと\eqref{f}より連立方程式

\begin{alignat*}{5}
c_0 &&&&&&&&&=2\\
c_0 &+& c_1 &+& c_2 &+& c_3 &+& c_4 &=3\\
c_0 &+& 2c_1 &+& 4c_2 &+& 8c_3 &+& 16c_4 &=1\\
c_0 &+& 3c_1 &+& 9c_2 &+& 27c_3 &+& 81c_4 &=3\\
c_0 &+& 4c_1 &+& 16c_2 &+& 64c_3 &+& 256c_4 &=0
\end{alignat*}

を解けばいいのですが、今日はこれを解かずにラグランジュ補間」で求めてみようと思います。ラグランジュ「ジュ」は「寿」とも読めるので縁起がいいですし(?)

ラグランジュ補間とは

ラグランジュ補間とはx_1, \ldots, x_kたちは互いに異なるものとしたときに、k+1個の点(x_0,y_0),\ldots,(x_k, y_k)をすべて通るようなk次以下の多項式を求める方法のことです。その多項式L(x)は以下で与えられます:

\begin{equation}
L(x)=\sum_{i=0}^k y_i\prod_{\substack{0\le m\le k \\ m\neq i}}\frac{x-x_m}{x_i - x_m}
\end{equation}

一応言っておくと、x_iたちは互いに異なるものとしているので分母にゼロは現れません。

見た目はごついですが、意味が分かると上手いことできてて面白い関数です。例えばこの関数にx_0を代入してみます:

\begin{align*}
L(\color{red}{x_0}) &=\sum_{i=0}^k y_i\prod_{\substack{0\le m\le k \\ m\neq i}}\frac{\color{red}{x_0}-x_m}{x_i - x_m}\\
&=y_{\color{blue}{0}}\prod_{\substack{0\le m\le k \\ m\neq \color{blue}{0}}}\frac{\color{red}{x_0}-x_m}{x_{\color{blue}{0}} - x_m} + y_{\color{blue}{1}}\prod_{\substack{0\le m\le k \\ m\neq \color{blue}{1}}}\frac{\color{red}{x_0}-x_{m}}{x_{\color{blue}{1}} - x_{m}} + \cdots + y_{\color{blue}{k}}\prod_{\substack{0\le m\le k \\ m\neq \color{blue}{k}}}\frac{\color{red}{x_0}-x_{m}}{x_{\color{blue}{k}} - x_{m}}\\
&=y_0 + 0 + \cdots +0\\
&=y_0
\end{align*}

何が起きたか分かりました?

まず、2行目の最初の項は
\begin{align*}
&y_{\color{blue}{0}}\prod_{\substack{0\le m\le k \\ m\neq \color{blue}{0}}}\frac{\color{red}{x_0}-x_m}{x_{\color{blue}{0}} - x_m}\\
=& y_{\color{blue}{0}}\frac{\color{red}{x_0}-x_1}{x_{\color{blue}{0}} - x_1}\cdots\frac{\color{red}{x_0}-x_k}{x_{\color{blue}{0}} - x_k}
\end{align*}なので積をとる部分はすべて約分されることが分かると思います。次の項は
\begin{align*}
&y_{\color{blue}{1}}\prod_{\substack{0\le m\le k \\ m\neq \color{blue}{1}}}\frac{\color{red}{x_0}-x_m}{x_{\color{blue}{1}} - x_m}\\
=&y_{\color{blue}{1}}\frac{\color{red}{x_0}-x_0}{x_{\color{blue}{1}} - x_0}\frac{\color{red}{x_0}-x_2}{x_{\color{blue}{1}} - x_2}\cdots\frac{\color{red}{x_0}-x_k}{x_{\color{blue}{1}} - x_k}
\end{align*}です。ここで、m0\le m\le kのうちm= \color{blue}{1}除外した範囲を動くことに注意してください。ということは当然m=\color{red}{0}も範囲内なので、積の中に\color{red}{x_0}-x_{\color{red}{0}}=0が現れていますね。よって積をとる部分はゼロになります。

それ以降の項も同様にゼロになってしまうため、結果y_0だけが生き残る、という仕組みになっています。x_1,…,x_kを代入した場合でもちゃんとL(x_1)=y_1,…,L(x_k)=y_kとなっています。いや~よくできてますね。

ラグランジュ補間を使ってみる

それでは本題に戻ります。\eqref{f}は「関数f(x)が5つの点(0,2)(1,3)(2,1)(3,3)(4,0)を通る」と言い換えることができるのでラグランジュ補間が使えます。その5つの点を通るような4次多項式f(x)

\begin{align*}
f_0(x)&=\frac{(x-1)(x-2)(x-3)(x-4)}{(0-1)(0-2)(0-3)(0-4)}\\
&=\frac{1}{24}(x-1)(x-2)(x-3)(x-4),\\
f_1(x)&=\frac{x(x-2)(x-3)(x-4)}{1(1-2)(1-3)(1-4)}\\
&=-\frac{1}{6}x(x-2)(x-3)(x-4),\\
f_2(x)&=\frac{x(x-1)(x-3)(x-4)}{2(2-1)(2-3)(2-4)}\\
&=\frac{1}{4}x(x-1)(x-3)(x-4),\\
f_3(x)&=\frac{x(x-1)(x-2)(x-4)}{3(3-1)(3-2)(3-4)}\\
&=-\frac{1}{6}x(x-1)(x-2)(x-4),\\
f_4(x)&=\frac{x(x-1)(x-2)(x-3)}{4(4-1)(4-2)(4-3)}\\
&=\frac{1}{24}x(x-1)(x-2)(x-3)
\end{align*}

とおいたとき

\begin{equation*}
f(x)=2f_0(x)+3f_1(x)+f_2(x)+3f_3(x)+0f_4(x)
\end{equation*}

で表せます(0かけたからf_4(x) は用なしじゃったのぉ)。あとは気合で(or 機械で)計算すれば

\begin{equation}
f(x)=-\frac{2 x^{4}}{3}+\frac{31 x^{3}}{6}-\frac{37 x^{2}}{3}+\frac{53 x}{6}+2\label{fx}
\end{equation}

となります。これは確かに\eqref{f}をみたします!



以上より恵方関数は

\begin{align*}
\theta(n) &= f(\operatorname{Mod}(n,5))\times90^{\circ} + 15^{\circ}\\
&=\left(-\frac{2 \operatorname{Mod}(n,5)^{4}}{3}+\frac{31 \operatorname{Mod}(n,5)^{3}}{6}-\frac{37 \operatorname{Mod}(n,5)^{2}}{3}+\frac{53 \operatorname{Mod}(n,5)}{6}+2\right)\times90^{\circ} + 15^{\circ}
\end{align*}

と表せました。


やったぁ~~~





ウーン…(´・ω・`)


mを法とする連立合同式

ラグランジュ補間で求めた式、無理矢理1行にするとModがしつこくて長いですね(やる前から知ってた)(Modは悪くない)(2行じゃダメなんですか)。あと、分数が出てくるのもちょっとかっこ悪いかも(だいたいそうなるもんだぞ)。整数の方がほら…なんか縁起がよさそうだし(せめて何かこじつけろ)

そこで、こういうアイデアはどうでしょう?
\begin{equation}
\begin{cases}
g(0)\equiv2 \pmod5\\
g(1)\equiv3 \pmod5\\
g(2)\equiv1 \pmod5\\
g(3)\equiv3 \pmod5\\
g(4)\equiv0 \pmod5
\end{cases}\label{g}
\end{equation}となるような整数係数多項式g(x)=c_0+c_1x+\cdots+c_k x^kを見つけるんです。そのような多項式g(x)が見つかれば、g(5)以降もmod5で「2, 3, 1, 3, 0」を繰り返すので
\begin{equation*}
\theta(n) = \operatorname{Mod}(g(n),5)\times90^{\circ} + 15^{\circ}\label{2}
\end{equation*}として恵方関数を定義できます。これならModは一回しか現れないのでだいぶスッキリします。

それでは、\eqref{g}をみたす整数係数多項式を求めていきましょう!

mを法とした連立合同式の解

ここで合同式における基本的な命題を振り返ります。

命題
mabを整数とする.maが互いに素のとき,1次合同式
\begin{equation}
ax\equiv b\pmod m
\end{equation}は整数解xをもつ.また,解はmを法として一意に定まる.

この命題をmを法とする連立合同式

\begin{equation}
\left\{\begin{array}{cc}
a_{11} x_{1}+a_{12} x_{2}+\cdots+a_{1 n} x_{n} \equiv b_{1} & \pmod m \\
a_{21} x_{1}+a_{22} x_{2}+\cdots+a_{2 n} x_{n} \equiv b_{2} & \pmod m \\
\vdots & \vdots \\
a_{n 1} x_{1}+a_{n 2} x_{2}+\cdots+a_{n n} x_{n} \equiv b_{n} & \pmod m
\end{array}\right.\label{eqs}
\end{equation}

拡張することを考えます。ここで、線形代数で習ったように\eqref{eqs}をn次正方行列A=(a_{ij})と,n次元の列ベクトル\boldsymbol{b}=(b_i)\boldsymbol{x}=(x_i)を用いて
\begin{equation*}
A\boldsymbol{x}\equiv\boldsymbol{b}\pmod m
\end{equation*}と表記することにします。すると次の定理が成り立ちます(!)
定理1(連立合同式行列式
mを整数,A整数成分をもつn次正方行列,\boldsymbol{b}整数成分n次列ベクトルとする.m\det Aが互いに素のとき,連立合同式
\begin{equation}
A\boldsymbol{x}\equiv\boldsymbol{b}\pmod m\label{mat}
\end{equation}は整数成分の解\boldsymbol{x}をもつ.また,解はmを法として一意に定まる.
証明
\tilde{A}Aの余因子行列とすると\tilde{A}は整数を成分にもち,
\begin{equation*}
A\tilde{A}=\tilde{A}A=(\det A) E
\end{equation*}が成り立つ(E単位行列).ここでm\det Aは互いに素なのでc\det A\equiv 1\pmod mとなる整数cがとれる.そこで\boldsymbol{x}=c\tilde{A}\boldsymbol{b}とおくと\boldsymbol{x}は整数を成分にもち,
\begin{align*}
A\boldsymbol{x}&=A(c\tilde{A}\boldsymbol{b})=cA\tilde{A}\boldsymbol{b}=c(\det A)\boldsymbol{b}\\
&\equiv \boldsymbol{b}\pmod m
\end{align*}をみたすので\eqref{mat}の解である.

次に一意性を示す.\boldsymbol{x}_1\boldsymbol{x}_2を\eqref{mat}の解とするとA(\boldsymbol{x}_1 - \boldsymbol{x}_2)\equiv\boldsymbol{0}\pmod mとなるので左から\tilde{A}をかけると
\begin{equation*}
\det A(\boldsymbol{x}_1 - \boldsymbol{x}_2)\equiv\boldsymbol{0}\pmod m .
\end{equation*}左から整数cをかければ\boldsymbol{x}_1 \equiv \boldsymbol{x}_2\pmod m が得られる.■


連立合同式を解く

準備ができたのでg(x)を求めていきましょう。ラグランジュ補間のときと同様に4次以下の多項式

\begin{equation*}
g(x)=c_0+c_1x+c_2x^2+c_3x^3+c_4x^4
\end{equation*}

を見つけていこうと思います。\eqref{g}より

\begin{alignat*}{5}
c_0 && && && && &\equiv2 \pmod5\\
c_0 &+& c_1 &+& c_2 &+& c_3 &+& c_4 &\equiv3 \pmod5\\
c_0 &+& 2c_1 &+& 4c_2 &+& 3c_3 &+& c_4 &\equiv1 \pmod5\\
c_0 &+& 3c_1 &+& 4c_2 &+& 2c_3 &+& c_4 &\equiv3 \pmod5\\
c_0 &+& 4c_1 &+& c_2 &+& 4c_3 &+& c_4 &\equiv0 \pmod5
\end{alignat*}

となりますが、これを行列で表すと

\begin{equation*}
\begin{pmatrix}
1 & 0 & 0 & 0 & 0\\
1 & 1 & 1 & 1 & 1\\
1 & 2 & 4 & 3 & 1\\
1 & 3 & 4 & 2 & 1\\
1 & 4 & 1 & 4 & 1
\end{pmatrix}
\begin{pmatrix}
c_0\\
c_1\\
c_2\\
c_3\\
c_4
\end{pmatrix}\equiv
\begin{pmatrix}
2\\
3\\
1\\
3\\
0
\end{pmatrix} \pmod5
\end{equation*}

となります。
\begin{equation*}
\begin{vmatrix}
1 & 0 & 0 & 0 & 0\\
1 & 1 & 1 & 1 & 1\\
1 & 2 & 4 & 3 & 1\\
1 & 3 & 4 & 2 & 1\\
1 & 4 & 1 & 4 & 1
\end{vmatrix}=9
\end{equation*}は5と互いに素なので、定理1により解をもつことが保証されました!定理1の証明中の解\boldsymbol{x}=c\tilde{A}\boldsymbol{b}は余因子行列を求めるのがだるそうなので、拡大係数行列の基本変形mod5バージョンでこねくり回して解きます。

\begin{align*}
& \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 2\\
1 & 1 & 1 & 1 & 1 & 3\\
1 & 2 & 4 & 3 & 1 & 1\\
1 & 3 & 4 & 2 & 1 & 3\\
1 & 4 & 1 & 4 & 1 & 0
\end{pmatrix}\\
\overset{\text{2行以降}\atop-\text{1行}}{\longrightarrow} & \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 2\\
0 & 1 & 1 & 1 & 1 & 1\\
0 & 2 & 4 & 3 & 1 & -1\\
0 & 3 & 4 & 2 & 1 & 1\\
0 & 4 & 1 & 4 & 1 & -2
\end{pmatrix}\\
\overset{\text{3行以降}\atop -\text{2行}\times n}{\longrightarrow} & \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 2\\
0 & 1 & 1 & 1 & 1 & 1\\
0 & 0 & 2 & 1 & -1 & -3\\
0 & 0 & 1 & -1 & -2 & -2\\
0 & 0 & -3 & 0 & -3 & -6
\end{pmatrix}\\
\overset{\text{2行} -\text{4行},\atop\text{3行} -\text{4行}}{\underset{\color{blue}{\text{5行}\div(-3)}}{\longrightarrow}} & \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 2\\
0 & 1 & 0 & 2 & 3 & 3\\
0 & 0 & 1 & 2 & 1 & -1\\
0 & 0 & 1 & -1 & -2 & -2\\
0 & 0 & 1 & 0 & 1 & 2
\end{pmatrix}\\
\overset{\text{4行} -\text{3行},\atop\text{5行} -\text{3行}}{\longrightarrow} & \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 2\\
0 & 1 & 0 & 2 & 3 & 3\\
0 & 0 & 1 & 2 & 1 & -1\\
0 & 0 & 0 & -3 & -3 & -1\\
0 & 0 & 0 & -2 & 0 & 3
\end{pmatrix}\\
\overset{\text{2行} +\text{5行},\atop\text{3行} +\text{5行}}{\underset{\text{4行}\times(-1)+\text{5行}}{\longrightarrow}} & \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 2\\
0 & 1 & 0 & 0 & 3 & 6\\
0 & 0 & 1 & 0 & 1 & 2\\
0 & 0 & 0 & 1 & 3 & 4\\
0 & 0 & 0 & -2 & 0 & 3
\end{pmatrix}\\
\overset{\text{5行} +\text{4行}\times2}{\longrightarrow} & \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 2\\
0 & 1 & 0 & 0 & 3 & 6\\
0 & 0 & 1 & 0 & 1 & 2\\
0 & 0 & 0 & 1 & 3 & 4\\
0 & 0 & 0 & 0 & 6 & 11
\end{pmatrix}\\
\overset{\bmod5}{\equiv} & \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 2\\
0 & 1 & 0 & 0 & 3 & 6\\
0 & 0 & 1 & 0 & 1 & 2\\
0 & 0 & 0 & 1 & 3 & 4\\
0 & 0 & 0 & 0 & \color{blue}{1} & \color{blue}{1}
\end{pmatrix}\\
\overset{\text{5行で攻める}}{\longrightarrow} & \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 2\\
0 & 1 & 0 & 0 & 0 & 3\\
0 & 0 & 1 & 0 & 0 & 1\\
0 & 0 & 0 & 1 & 0 & 1\\
0 & 0 & 0 & 0 & 1 & 1
\end{pmatrix}\\
\end{align*}

ふぅ、ということで
\begin{equation*}
\begin{pmatrix}
c_0\\
c_1\\
c_2\\
c_3\\
c_4
\end{pmatrix}\equiv
\begin{pmatrix}
2\\
3\\
1\\
1\\
1
\end{pmatrix} \pmod5
\end{equation*}が求まったので
\begin{equation}
g(x)=x^4+x^3+x^2+3x+2\label{gx}
\end{equation}が得られました。以上より

\begin{align*}
\theta(n) &= \operatorname{Mod}(g(n),5)\times90^{\circ} + 15^{\circ}\\
&= \operatorname{Mod}(n^4+n^3+n^2+3n+2, 5)\times90^{\circ} + 15^{\circ}
\end{align*}

恵方関数です。これは縁起がいい!


ベクトル値にして遊ぶ

最後に、恵方関数\theta(n)があれば
\begin{equation*}
(\cos(\theta(n)), \sin(\theta(n)))
\end{equation*}として恵方をベクトルで表すことができるので、計算アプリでチェックします。あるいは複素数「複」は「福」なので、複素平面で考えて
\begin{equation*}
\cos(\theta(n))+i\sin(\theta(n))
\end{equation*}で計算した方が縁起がいいかもしれません。

ということで、作ってみました。
恵方関数 - Desmos

やっぱり思い通りに動いてくれると嬉しいですね~



また、せっかく多項式で補間したので、nを連続的に動かして変化を観察しましょう。
まずはfの方から。

f:id:corollary2525:20210202162849g:plain

役目を終えた後のヤケクソ感。

まあ4次関数ですし、f(x)のグラフをみれば原因は明らかなんですけどね。
ヤケクソ感を解消したかったらf(x)5次多項式にしてf(5)=2という条件も課してラグランジュ補間をしたらいいでしょう。



次にgで作った方を見ましょう。

f:id:corollary2525:20210202162746g:plain

この曲を聴きたくなりました。

まぁね、先にg(n)を計算してからmod5してますからね。変化の激しさはg微分して2020付近を代入したら分かるわな。



最後になりましたが、個人的には定理1という、合同式行列式が関係する定理が面白かったです。できれば合同式ヴァンデルモンドの行列式との関係まで深堀りしたかったのですが、長くなったのでまた今度にします。


*1:北基準時計回り