2020-11-09

Redlog 使ってみた

数学

Quantifier Elimination (QE，限量子消去) を解いてくれるREDUCEのパッケージRedlogを使って大学の入試問題を解いてみます．

東工大 1969前期
実数 $a,b,c,x,y,z,p$ が次の次の4条件を満たしている： $a^2-b^2-c^2>0, ax+by+cz=p, ap<0, x>0,$ このとき $x^2-y^2-z^2$ の符号を調べよ．

名古屋大 1963前期
正方形 $ABCD$ とその内部の $1$ 点 $P$ がある．線分 $AP,BP,CP,$ の長さがそれぞれ $7,5,1$ であるとき，この正方形の面積を求めよ．

Quantifier Elimination (QE，限量子消去) とは

実数の理論において\begin{align}(a\neq0) \wedge\bigl( \forall x (ax^2+bx+c>0)\bigr)\end{align} は次の論理式と同値です．\begin{align} (a>0) \wedge (b^2-4ac<0)\end{align} 上の論理式は ∀(任意記号) とそれに縛られた束縛変数 $x$ がありますが，下の論理式は自由変数 $a,b,c$ だけの式になっています．∃ (存在記号) と ∀(任意記号) をあわせて限量子 (quantifier) と呼びます．このように論理式を，限量子が含まない同値な論理式に変形することをQuantifier Elimination (QE，限量子消去) と呼びます．実数の理論では，代数的な一階述語論理式を限量子が含まない同値な論理式に変形するアルゴリズムが存在することが知られています（タルスキの定理）．特に閉論理式（自由変数を含まない論理式）は，真か偽かをアルゴリズムで判定できます．

Redlog の使い方

REDUCE を Obtaining REDUCE から入手してインストールしましょう．
REDUCE を起動して，Redlog（reduce logic system）パッケージが読み込まれていないならロードしておきます．

load_package "redlog";

次にQEのためのおまじないをします（実数体 $\mathbb{R}$ の理論ですという宣言）．

rlset r$

Redlog には QE 用の関数 rlqe (Virtual Substitution)，rlcad (Cylindrical Algebraic Decomposition) などが用意されています．論理式は次の表くらいが分かっていればとりあえずは十分でしょう．

記号論理	Redlog
$\neq, \le, \ge$	<>,<=, >=
[式1] $\wedge$ [式2]	[式1] and [式2]
[式1] $\vee$ [式2]	[式1] or [式2]
[式1] $\rightarrow$ [式2]	[式1] impl [式2]
$\lnot$ [式]	not [式]
$\forall x_1,\forall x_2,\dots, \forall x_n$ [述語]	all({ $x_1,x_2,\dots,x_n$ }, [述語])
$\exists x_1,\exists x_2,\dots, \exists x_n$ [述語]	ex({ $x_1,x_2,\dots,x_n$ }, [述語])

問題をQEに翻訳

東工大 1969前期
実数 $a,b,c,x,y,z,p$ が次の次の4条件を満たしている： $a^2-b^2-c^2>0, ax+by+cz=p, ap<0, x>0,$ このとき $x^2-y^2-z^2$ の符号を調べよ．

次の論理式を考えてみましょう．\begin{align}\exists a\exists b\exists c\exists x\exists y\exists z\exists p\bigl(&(a^2-b^2-c^2>0)\wedge (ax+by+cz=p)\\&\wedge (ap<0)\wedge (x>0) \wedge (r=x^2-y^2-z^2)\bigr) \end{align}これと同値な自由変数 $r$ だけに関する論理式が，条件 \begin{align} a^2-b^2-c^2>0, ax+by+cz=p, ap<0, x>0\end{align}のもとでの $r:=x^2-y^2-z^2$ の取りうる値の条件となります．

名古屋大 1963前期
正方形 $ABCD$ とその内部の $1$ 点 $P$ がある．線分 $AP,BP,CP,$ の長さがそれぞれ $7,5,1$ であるとき，この正方形の面積を求めよ．

$xy$ -平面上で $A=(0,a), B=(0,0), C=(a,0), P=(b,c)$ とおくと条件は \begin{align} b^2+(a-c)^2=7^2, b^2+c^2=5^2, (a-b)^2+c^2=1^2, a > b > 0, a > c > 0 \end{align} と書けます．この条件のもとでの $a^2$ の値を求めたいので \begin{align}\exists a\exists b\exists c\bigl(&(b^2+(a-c)^2=7^2)\wedge (b^2+c^2=5^2)\wedge ((a-b)^2+c^2=1^2) \\ &\wedge (a > b)\wedge(b > 0) \wedge (a > c) \wedge(c > 0)\wedge (s=a^2)\bigr) \end{align} から限量子を消去して， $s:=a^2$ の取りうる値を求めればよいことになります．

Redlogで解いてみる

Redlog の QE 用の関数 rlqe，rlcad を使って実際に解いてみましょう．showtime で直前の計算でどれくらい時間がかかったのかを表示できます．

f:id:egory_cat:20201109023040p:plain
というわけで東工大の問題の答えは「負」，名古屋大の問題の答えは「 $32$ 」となります．

2020-08-15

gcd(nᵖ+a,(n+1)ᵖ+a) の突然の裏切り予想

数学

前回の記事で $\mathrm{gcd}(n^5+5,(n+1)^5+5)$ や $\mathrm{gcd}(n^{17}+9,(n+1)^{17}+9)$ はずっと $1$ のままなのに，大きな $n$ で突然 $1$ でなくなるという「突然の裏切り」が起こることを紹介しました．
n⁵+5 と (n+1)⁵+5 の最大公約数 - 現実と数学の区別が付かない
他にも探してみると似た性質を持つ多項式がいくつか見つかります．それらを観察して，次の予想を立てました．

予想
$p$ を素数とし， $x,t$ を変数とする． $x^p +t$ と $(x+1)^p+t$ の $x$ に関する終結式を $r_p(t)$ とする．このとき次が成り立つ．

$r_p(t)$ は $\mathbb{Z}[t]$ の元として既約．
$r_p(a)$ が素数となる自然数 $a$ が無限個存在する.
$q:=r_p(a)$ が素数のとき， $x^p+a$ は $\mathbb{F}_q[x]$ において相異なる $p$ 個の１次式の積に因数分解される．
さらに $\beta^p +a=(\beta+1)^p+a=0$ となる $\beta\in \mathbb{F}_q$ がただ1つ存在する．

この予想が正しければ $r_p(a)$ が素数になる $a$ はたくさん存在し，そのとき $\mathrm{gcd}(n^p+a,(n+1)^p+a)$ で「突然の裏切り」が起こることになります;
$r_p(a)$ が素数のとき $\mathbb{F}_q$ での像が $\beta$ となる $b\in\mathbb{Z},~0\le b\le q-1$ を取ると，\begin{align*}\mathrm{gcd}(n^{p}+a,(n+1)^{p}+a)= \begin{cases}
q &(n\equiv b \mod q)\\
1 &(\mbox{otherwise})
\end{cases}\end{align*} となります．

この予想を解けた方には「裏切り」に関する何かしらの商品の進呈します．

状況証拠

$r_p(t)$ を小さい $p$ に対して計算してみると \begin{align*}
r_2(t)=&~4t+1\\
r_{3}(t)=&~27t^2+1\\
r_{5}(t)=&~3125t^4+625t^2+1\\
r_{7}(t)=&~823543t^6+6000099t^4+12005t^2+1\\
r_{11}(t)=&~285311670611t^{10}+86346568369165434t^8\\&~+15483521175490489t^6+1782585424103t^4+3879865t^2+1\\
\vdots
\end{align*}となります． $p\le 67$ までは既約であることを計算機で確認しました（でも小さいところで成り立っているからといってずっと成り立つとは限らないという話を最近聞いたような……）

追記： $r_p(t)$ は定数項が $1$ なので原始的で， $\mathbb{Q}[t]$ で既約であることを示せば十分です． $t^{p-1} r_p(\frac{1}{tp})$ が既約であることを Eisensteinの判定法で証明できそうです．

各 $p$ に対し， $r_p(a)$ が素数となる $50$ 以下の $a$ を列挙すると次のようになります．

$p=2$ のとき $a={\small 1,3,4,7,9,10,13,15,18,22,24,25,27,28,34,37,39,43,45,48,49}$
$p=3$ のとき $a= 2,4,12,26,28,40,42$
$p=5$ のとき $a=2,5,7,8,13,17,18,22,27,31,40,46$
$p=7$ のとき $a=4,18,38,40,42$
$p=11$ のとき $a=4,5,8,14,17,18,23,26,49,50$
$p=13$ のとき $a=36$
$p=17$ のとき $a=9,15,34,38$
$p=19$ のとき $a=6,28,33,49$
$p=23$ のとき $a=45$
$p=29$ のとき $a=14,15$
$p=31$ のとき $a=50$
$p=37$ のとき $a=14$

そしてこの表のすべての $p,a$ で予想3，4 が成り立っています．

$p=2$ のとき

$r_2(a)=4a+1$ が素数となる自然数 $a$ は無限個存在します．
ピタゴラス素数 - Wikipedia
また， $q=4a+1$ が素数のとき $\mathbb{F}_q[x]$ において\begin{align} (x+2a)(x+2a+1)=x^2+(4a+1)x+4a^2+2a=x^2-a+2a=x^2+a\end{align}となり， $(2a+1)-2a=1$ なので，予想3，4も成り立っています．

$p=3$ のとき

$r_3(a)=27a^2+1$ が素数となる自然数 $a$ 無限個存在するかどうかは私には分かりません．しかし， $\mathrm{gcd}(r_3(1), r_3(2))=\mathrm{gcd}(28,109)=1$ なので，ブニャコフスキー予想が正しいと仮定すれば， $r_3(a)=27a^2+1$ が素数となる自然数 $a$ が無限個存在することになります．

ブニャコフスキー予想 - Wikipedia

また， $q=27a^2+1$ が素数のとき $\mathbb{F}_q[x]$ において \begin{align}&\left(x+\cfrac{3a+1}{2}\right)\left(2x+\cfrac{3a-1}{2}\right)(x-3a)=x^3-\cfrac{27a^2+1}{4} x+\cfrac{-27a^3+3a}{4}\\
&=x^3+\cfrac{a+3a}{4}=x^3+a
\end{align} が成り立ちます． $q$ は $28$ 以上の素数なので $2$ が可逆元であることに注意してください． $\cfrac{3a+1}{2}-\cfrac{3a-1}{2}=1$ なので，予想3，4が成り立っています．

$p=5$ のとき

$\mathrm{gcd}(r_5(1),r_5(2))=\mathrm{gcd}(3751,52501)=1$ なので，ブニャコフスキー予想が正しいと仮定すれば， $r_5(a)=3125a^4+625a^2+1$ が素数となる自然数 $a$ が無限個存在することになります．

$q=r_5(a)$ が素数のとき， $\mathbb{F}_q[x]$ において \begin{align}
x^5+a=
&~\left(x+\cfrac{625a^3+90a+11}{22}\right)\\
&~\left(x+\cfrac{625a^3+90a-11}{22}\right)\\
&~\left(x+\cfrac{625a^3+145a}{11}\right)\\
&~\left(x-\cfrac{1250a^3+125a^2+235a+7}{22}\right)\\
&~\left(x-\cfrac{1250a^3-125a^2+235a-7}{22}\right)
\end{align} が成り立ちます． $\cfrac{625a^3+90a+11}{22}-\cfrac{625a^3+90a-11}{22}=1$ であることから，予想3，４が成り立つことが分かります．

$p=37, a=14$ のときを鑑賞

最後に $p=37, a=14$ のときの「突然の裏切り」を鑑賞してみましょう．\begin{align*} &b=\\
&114639163986838022433471493033582231487059843696971501729723\\
&838399007961015906196300583570493634776987376087643269329057\\
&161223105482392984049807969081020156915883248950085087033826\\
&60163326302581676327242502043057172204427\\
&q= \\
&104310066642138695684562957075767539994224919231022785321022\\
&998638201111250690582693240070811357134844988196278442827230\\
&538607105018619585940730548217113231821013606082099727905440\\
&865311886180325006446369415049503294784509\\
\end{align*}とおくと，\begin{align*}\mathrm{gcd}(n^{37}+14,(n+1)^{37}+14)= \begin{cases}
q &(n\equiv b \mod q)\\
1 &(\mbox{otherwise})
\end{cases}\end{align*} となります．ここまでくると面白いかどうかも分からないですが何か凄いです．

2020-08-13

n⁵+5 と (n+1)⁵+5 の最大公約数

数学

Twitterで見かけた答えが意外過ぎる問題．

多項式の公約数と言えば、昔どこかに投稿したんだけど、nが自然数の時の n^5+5 と (n+1)^5+5 の正の公約数としてあり得る整数が、おそらく見た目からは予想できない結果で、面白い。
— nishimura (@icqk3) 2020年8月10日

自然数 $a,b$ の最大公約数 (greatest common divisor) を $\mathrm{gcd}(a,b)$ で表します．

$n$ を自然数とする． $\mathrm{gcd}(n^5+5,(n+1)^5+5)$ を求めよ．

この手の問題は，小さい $n$ で試してみるのが常套手段です． $n\le 100$ くらいまで試してみると，すべて $\mathrm{gcd}(n^5+5,(n+1)^5+5)=1$ となります．その後， $n$ を 1万，10万と増やしていっても，ずっと gcd は $1$ のままです．こうなると，はいはいパターン見えてきたよと $\mathrm{gcd}(n^5+5,(n+1)^5+5)=1$ であると予想を立て，数学的帰納法で証明しようという気になります．しかしこれはうまくいきません．実は $n<533360$ のときは $\mathrm{gcd}(n^5+5,(n+1)^5+5)=1$ なのに， $n=533360$ で急に\begin{align}\mathrm{gcd}(533360^5+5,533361^5+5)= 1968751\end{align}となります．こんなことが起こることも面白いですし，これに気が付いたことも凄いと思います．この $n$ が約53万でフリーザの戦闘力とほぼ同じなことも芸術点が高いです．

Why Japanese people?

なぜこのようなことが起こるのか考えてみましょう．

$a_1,\dots,a_k \in \mathbb{Z}$ で生成される $\mathbb{Z}$ のイデアルを $\langle a_1,\dots,a_k \rangle$ で表すことにします．すると次が成り立ちます．

\begin{align}\langle a, b\rangle=\langle \mathrm{gcd}(a,b)\rangle\end{align}

最大公約数をこの等式で理解することは結構大切です．例えばユークリッドの互除法はイデアルの生成系の取り換えをしているだけと見ることがでいます．また， $\langle a, b\rangle$ の要素は $\mathrm{gcd}(a,b)$ の倍数になっているので，具体的な整数 $m\in \langle a, b\rangle$ を $1$ つ見つけることができると， $\mathrm{gcd}(a,b)$ の候補が有限個に絞られます．

$\mathrm{gcd}(n^5+5,(n+1)^5+5)$ がどうなるか考えていきましょう．
\begin{align}
f_0(n):=&~(n+1)^5+5\\
f_1(n):=&~n^5+5\\
f_2(n):=&~f_0(n)-f_1(n)=5n^4+10n^3+10n^2+5n+1\\
f_3(n):=&~5f_1(n)-(n-2)f_2(n)=10n^3+15n^2+9n+27\\
f_4(n):=&~4f_2(n)-(2n+1)f_3(n)=7n^2-43n-23\\
f_5(n):=&~49f_3(n)-(70n+535)f_4(n)=25056n+13628\\
f_6(n):=&~156950784 f_4(n)-(43848n-293201)f_5(n)=385875196
\end{align} なので， $385875196\in \langle n^5+5,(n+1)^5+5\rangle$ となります．ちょっと注意をしておくと，この計算はユークリッドの互除法そのものではありませんので，最後に出てきた $385875196$ は最大公約数とは限りません．しかし，定義から $k\ge 2$ のとき $f_k(n) \in \langle f_{k-1}(n), f_{k-2}(n) \rangle$ であるので，\begin{align} 385875196&=f_6(n)\in \langle f_{5}(n), f_{4}(n) \rangle \subset\langle f_{4}(n), f_{3}(n) \rangle \subset\langle f_{3}(n), f_{2}(n) \rangle \\
& \subset\langle f_{2}(n), f_{1}(n) \rangle \subset\langle f_{1}(n), f_{0}(n) \rangle = \langle n^5+5,(n+1)^5+5\rangle \end{align} が言えるのです．

つまり $\mathrm{gcd}(n^5+5,(n+1)^5+5)$ は $385875196$ の約数です． $385875196$ を素因数分解してみると\begin{align}
385875196=2^2\times 7^2\times 1968751
\end{align}となります．

素因数 $p=2,7,1968751$ ごとに， $p$ が $\mathrm{gcd}(n^5+5,(n+1)^5+5)$ を割り切るかどうか考えてみます．

$p$ が $\mathrm{gcd}(n^5+5,(n+1)^5+5)$ を割り切ることは $x^5+5\in \mathbb{F}_p[x]$ が $n, n+1$ の両方を根に持つことと同値です．

$p=2$ のとき．

$x^5+5=x^5+1 \in \mathbb{F}_2[x]$ は $x=1$ しか根を持ちません．
よって $\mathrm{gcd}(n^5+5,(n+1)^5+5)$ は $2$ で割り切れません．

$p=7$ のとき．

$x^5+5 \in \mathbb{F}_7[x]$ は $x=4$ しか根を持ちません．
よって $\mathrm{gcd}(n^5+5,(n+1)^5+5)$ は $7$ で割り切れません．

$p=1968751$ のとき．

\begin{align} x^5+5=(x-1066696)(x-1707583)(x-96502)(x-533360)(x-533361) \end{align} という因数分解が $\mathbb{F}_{1968751}[x]$ で成り立ちます．つまり $n^5+5$ は \begin{align} n\equiv 1066696, 1707583, 96502, 533360, 533361 \mod 1968751\end{align} のときに $1968751$ の倍数となります． $n^5+5, (n+1)^5+5$ が共に $1968751$ の倍数となるのは $n\equiv 533360 \mod 1968751$ のときのみです．

よって，冒頭の問題の答えは \begin{align} \mathrm{gcd}(n^5+5, (n+1)^5+5)= \begin{cases}
1968751 & (n\equiv 533360 \mod 1968751)\\
1 & (\mbox{otherwise})
\end{cases}\end{align} となります．

以上．

追記

世の中には整数係数グレブナー基底という超便利な道具があり，Sage などに実装されているそうです．
doc.sagemath.org
$\langle x^5+5,(x+1)^5+5\rangle \subset \mathbb{Z}[x]$ の整数係数グレブナー基底を計算すると $\{x + 1435391, 1968751\}$ となります（下記リンクで計算結果を見ることができます）．
sagecell.sagemath.org
これは $\mathbb{Z}[x]$ のイデアルとして \begin{align}\langle x^5+5,(x+1)^5+5\rangle =\langle x + 1435391, 1968751 \rangle\end{align} であることを示しています．このことから，自然数 $n$ に対して \begin{align} \mathrm{gcd}(n^5+5, (n+1)^5+5)=\mathrm{gcd}(n + 1435391, 1968751) \end{align} が成り立つことが分かります． $1968751$ が素数であることも計算機で確かめることができるので， $\mathrm{gcd}(n^5+5, (n+1)^5+5)$ は $1$ まはた $1968751$ であり， $\mathrm{gcd}(n^5+5, (n+1)^5+5)=1968751$ となるのは \begin{align} n\equiv - 1435391 \equiv 533360 \mod 1968751\end{align} のときということが分かります．

$x^{17}+9, (x+1)^{17}+9$ でも類似の現象が起こるそうです．
sagecell.sagemath.org
このグレブナー基底の計算によると \begin{align} &\langle x^{17}+9, (x+1)^{17}+9\rangle \\
&= \left\langle\begin{array}{r} x + 512149312322827330662764931050044963334032796143126, \\
8936582237915716659950962253358945635793453256935559~ \end{array} \right\rangle\end{align}が成り立ちます．この2番目の生成元である自然数が素数になっています．

さらに追記

くさだんごさんが記事を書いてくれました．
mochi-mochi61.hatenablog.com
イデアル $I=\langle f(x), g(x) \rangle \subset \mathbb{Z}[x]$ に対して， $I\cap \mathbb{Z}$ に含まれる元を得るために終結式を使っています．
終結式 - Wikipedia
私の記事ではなんちゃってユークリッドの互除法を使いましたが，余計な素因数 $2,7$ が出ていました．終結式を使う方法では余分な因子が出ません．