妖精現実

[ 新着記事 | 数学・プログラミング | 天文・暦
| シリア語・Unicode・詩 | ジョーク | 漫画・アニメ
| 字幕 | 哲学・ファンタジー | 全記事 ]

チラ裏コーナー（過去のチラ裏 | 野蛮な解法）

2020-07-02　野蛮な解法（その28の追記）

（♡）の方法は、右辺が +2 の場合にも同様に成立、（＊）が x = A² − 1 になる。右辺 ±2 の形の（♡）を適宜使うことで、d が4k+3型素数のときの直接検索は必要なくなった。d = 331, 379 などはもはや問題でなく、4k+3型素数については、恐らく700台まで進むことができる。d が合成数でも、それを法として ±2 が平方剰余なら、同じ手を試みることが可能。その結果、速度的に遅かった d = 214, 334 もほぼ一瞬になった。新しい壁は、4k+1型素数の d = 421。

d = 421 以外の d については、525 まで全部対応できるが、合成数のケースで少し時間がかかるものが4個ある。d = 301（依然として10秒単位の時間がかかる）、d = 436（分単位の時間がかかる）など。d = 421 についても、ブルートフォースで解決できるめどが立った。

2020-07-01　地球の自転が一時的に速くなっている　1日の長さ24時間を切る

精密に測ると地球の自転は日々変動しているが、長期的にはだんだん遅くなっている。「自転による1日は、定義上の24時間＝86400秒より1ミリ秒（1000分の1秒）くらい長い」というのが、これまでの大ざっぱな感覚だった。だから1000日くらい（数年）のオーダーで地球は時計より1秒遅れ、時計合わせのため「うるう秒」が挿入される…と。

ところが実測で、2020年6月5日ごろから、86400秒かからないで自転が終わる状態が続いている。「自転と時計の差」（UT1−UTC）が普通ならだんだん減るはずなのに（自転より時計の方が進みが速いので、マイナスが大きくなる）、逆に少しずつ増えている。この傾向は9月頃まで続き、その後も時々自転が速めになると予想されている（[1]）。このこと自体は、決して異常ではない。2004年ごろにも同様のことが起こり、その影響で7年間くらい一度もうるう秒がないことがあった（[2]）。24時間を切ったと言っても違いは、せいぜい1ミリ秒。地球の自転と無関係に1日は通常86400秒に固定されているので、生活に影響が出るわけではない。大昔の地球はもっと速く回っていて、1日が23時間だったときもあるというのだから、1ミリ秒ずれたくらいで、大騒ぎすることもないだろう。

しかし年オーダーの短期的な話として、「うるう秒が入らない最長期間」の新記録になる可能性がある。「うるう秒」制度が始まって以降、これまでの最高記録は「7年間うるう秒なし」。協定世界時（UTC）1998年12月31日の末尾（日本時間1999年1月1日8時59分59秒の後ろ）にうるう秒が挿入されてから、UTC2005年の末尾に再びうるう秒が入るまで7年あいた。前回のうるう秒はUTC2016年末。2020年末のうるう秒の有無が公式発表されるのは数日後だが、UT1−UTCマイナス0.2秒台で「うるう秒」という先例は、21世紀に入ってから1度もない。マイナス0.5秒になりそうな辺りで挿入してプラス0.5秒側に振るのが、常識的に無難な線だし、最近の運用は実際そんな感じ。

前回の「7年間なし」記録のときは、うるう秒を入れてから4年で、またマイナス0.3秒くらいまで地球が遅れた。一方、2020年末にうるう秒がない場合、結果は「ずれマイナス0.2秒くらい」＝前回の同じ条件より、ずれのペースが小さめ。この調子だと「8年間うるう秒なし」、場合によっては「10年間なし」といった新記録もあり得る。

地球の自転は、大ざっぱに10～20年の周期で速くなったり遅くなったりしているらしい。グラフを見ての漠然とした印象だが、現在入りつつある「谷」（1日の長さが短くなる）は前回よりちょっと深いかも…。地球の自転は気象の影響も受けているので、気象の乱れも関係しているのかもしれない。地球の自転が速めになるのは、UTCがずれにくいという点では良いこと。これまでの遅れの「貯金」がたっぷりあるので、「負のうる秒」などという事態にはならない。長期的には、どっちにしても誤差の範囲内の揺らぎだろう。

[1] https://datacenter.iers.org/data/latestVersion/6_BULLETIN_A_V2013_016.txt

[2] https://en.wikipedia.org/wiki/File:Deviation_of_day_length_from_SI_day.svg

[3] https://en.wikipedia.org/wiki/File:Leapsecond.ut1-utc.svg

全然関係ないが、変な夢を見たのでメモ。なぜか首藤氏が登場して「花の子ルンルンは、日本に二つだけの珍しいタイプの魔法少女物ですね」みたいな訳の分からないことを言っている。それで自分は「ようこそようこも珍しい魔法少女ですよね」と言おうとしていた…。起きた直後は「そもそも、ルンルンも、ようこも魔法少女じゃないし。訳の分からん夢だな」と思った。が…今よく考えると「魔法を使わない魔法少女だなんて、珍しいですね…」って、首藤にぴったりの、深い意味のある言葉だったのかもしれない。あの世で首藤に会ったら、ようこ→ホッキョンのノリで「脚本家って、すぐ交通事故を起こしたがる人のことだったんですかぁ？」って言うぞ…とか時々くだらないことを考えているが、夢の中ではそんなことは思わなかった。

2020-07-01　野蛮な解法（その28）　d=331の壁を突破

ペル方程式 x² − dy² = 1 について、普通に解かず、ゲーム感覚で全数検索をしてきた。d が4k+1型の素数なら便利なショートカットが存在するのだが、4k+3型の素数のときには同じショートカットを使えず、4k+3型の d = 331 が壁となっていた。何時間も計算すれば解は得られるだろうが、4k+3型素数に使えるうまいショートカットは無いものか…？

4k+1型のショートカットは、与式の右辺を −1 に置き換えた「負のペル」からの間接攻略だった。同様に、与式の右辺を −2 に置き換えた次の方程式を解くことで、d = 331 に対応できる:
A² − dB² = −2　（♡）

仮に（♡）を満たす正の整数 A, B が得られた場合、（♡）の両辺を2乗して
(A² − dB²)² = 4
(A²)² − 2(A²)(dB²) + (dB²)² = 4
(A²)² + 2(A²)(dB²) + (dB²)² − 4(A²)(dB²) = 4
(A² + dB²)² − d[4(AB)²] = 4　（♡♡）
になるが、（♡♡）の右辺と、左辺第2項・角かっこ内は、どちらも4の倍数なので、当然、左辺第1項も4の倍数。だから各項を4で割ることができて、
C − d(AB)² = 1
の形になる。ここで C は、（♡♡）の左辺第1項を4で割った結果の整数だが、もしそれが平方数 C = x² になってくれれば、後は y = AB と置くだけで、ペル方程式が解けたことになる。このアイデアは実を結ぶ。実際、（♡）より
dB² = A² + 2
であり、これを使うと、（♡♡）の左辺第1項は
(A² + A² + 2)² = (2A² + 2)² = [2(A² + 1)]² = 4(A² + 1)²
すなわち、それを4で割った結果 C とは (A² + 1)²。結局、（♡）が解けるなら
x = A² + 1, y = AB　（＊）
が、対応するペル方程式の解。

d = 331 の場合、（♡）は A = 52778687, B = 2900979 という解を持つ（これは7～8桁なので、容易に全数検索可能）。これらを（＊）に入れると:
x = 2785589801443970,
y = 153109862634573
が、ペル方程式 x² − 331y² = 1 の解。検算すると、確かに合っている！

たった16桁、1000兆のオーダーだが、直接検索は難しい。d が4k+1型素数のとき、すぐに解が京や垓のオーダーになるが、あれは一種の「バブル」。内部的な検索は10桁やそこらなのだが、最終結果ではそれが大ざっぱに平方されるため、見掛け上、20桁（1000京）オーバーも珍しくない。現状、ブルートフォースで直接検索できているのは、秒単位では12～13桁くらいまで。分単位では14桁～ギリギリ15桁。あとたった1桁とはいえ、ブルートフォースでは、それは計算量的に10倍。「ギリギリできていることの10倍」は、きつい。並列処理とか、強力なCPUとか、高速なコンパイル言語とか、そういう手もあるけど…

d = 331 については
x = 8278, y = 455 ⇒ x² − 331y² = 9
という関係も見つかったが、これは各項を9で割れないので、役立ちそうにない。それと、一つ宿題が残っている。d = 331 のペル方程式については、上記の通り、間接的な方法で15～16桁の解を得たが、pell4 で直接検索したのは12桁まで。上記が基本解（最小の正の解）であることを言うためには「12～16桁のまだ探してない範囲に、これより小さい他の解がない」ことを示す必要がある。4k+1型関連の「巨大バブル解」にも、同様の問題が…。「12桁未満の解がなく16桁の解があるなら、それが最小解」というのは、ペル方程式の基本性質から分かることだが、このメモでは、まだそういう考察をしていない。

2020-06-30　野蛮な解法（その27）　普通に解いてもつまらない！　数学の縛りプレイ

x² − 214y² = 1 のような形の2変数方程式（「ペル方程式」と呼ばれる）の、正の整数解を求めたい。パズル的に面白いだけでなく、数学的にもいろいろと重要な意味を持つのだが、とりあえず意味なんか気にせず、遊んでみる。

x² − dy² = 1 の d（平方数ではない自然数）が小さい場合、大抵は簡単に解を見つけることができる。例えば
x² − 11y² = 1
の最小解は、容易に暗算可能。左辺に y = 1, 2, 3, … を入れると y = 3 のとき第2項は −99 になるが、そのとき x² = 100 となり、それを満たす正の整数 x は、もちろん10。次の例題も、同様にして、暗算で簡単に解けるので、やってみよう。
x² − 12y² = 1

ペル方程式の解法としては、平方根の近似分数が利用されることが多い。ここでは普通の道を歩まず、上記の暗算のような「全数検索」的アルゴリズム pell4(d) をテスト実装した。pell4 は、まず小さい y について、上の例そのままの素朴な検索を試み、それでは解が見つからない場合、d が 4k+1型の素数なら pell4a(d)、d が 4k+3型の素数なら pell4b(d)、それ以外の場合、pell4c(d) をそれぞれ呼び出す。

pell4a は「その26」の pell3a の改良版。解きたい方程式の右辺を −1 に置き換えたもの（負のペル方程式）の解を見つけることで、間接的に与式を解く。まず 64d を法として解 x になり得る数を決定。それらの数について、1ステップ 64d の「ホイールシーブ」のような検索を1000万程度まで行う。これが不成功なら、今度は 64⋅81⋅25⋅49d を法として1200億まで検索。…「ホイールシーブ」の準備では、使いたい法の因子（素数べき）について、それぞれを法として「負のペル」の解 x になり得る数を決定し、得られた複数の条件に中国剰余定理を適用して、大きな法における1種類の条件にまとめておく。法 d に関しては、q² ≡ −1 (mod d) を満たす q を求めることに相当。これも全数検索可能だが、数論的には、原始根の (d − 1)/4 乗を q として {q, q³} が求めるもの。

pell4b も pell4a とほぼ同じ実装。こちらは「負のペル」の解を探さず、直接与式を解く。まず法 256d で1億まで。それで駄目なら、法 256⋅81⋅25⋅49d で2800億まで。q² ≡ 1 (mod d) を満たす q が必要だが、それは ±1 (mod d) に他ならない。

最後の pell4c は、合成数を扱うため、ややこしい。現在のテスト実装では、まず 128⋅81q を法として、5億まで検索。ここで q は、d を2と3でそれぞれ割れるだけ割った結果の数（言い換えれば、2とも3とも互いに素であるような q の約数うち、最大のもの）。それで駄目なら、今度は対応する「負のペル」が解を持ち得るか調べ、持ち得るならそれを探す。それも駄目な場合、法 256⋅81⋅25⋅49q で検索。この q は、d を2、3、5、7で割れるだけ割った結果の数。

以上によって、330以下の全ての d について、ブルートフォースでペル方程式を解けるようになった。所要時間は、易しいケースなら一瞬、難しいケースでもおおむね1秒以内だが、d = 214 は数秒かかり、d = 301 はさらに遅い。この実装（特に pell4c）はまだ試作段階。うまくやれば、もっと速くなるかもしれない。楽に解きたいだけなら、もっと速いやり方が複数あるのだが、あえてブルートフォースで。

pell4 がうまく対応できない最初のケース（時間がかかり過ぎる）は d = 331, 379（「その24」参照）。壁になっている d = 331 では、探したい数が16桁。今の速度（想定ターゲット最大13桁くらい）では、きつい。なにしろ全数検索なので、1桁増えるごとに探す場所が10倍、3桁増えれば1000倍…それだけ時間がかかる。「どこかで探し切れなくなる」という結末を承知の上で、あの手この手の高速化を行い、どこまで追いすがれるか…

pell4c の d = 301, 334 も時間がかかるが（13～14桁）、数十秒～5分くらい待てば解が得られるので、一応、実用になる（速度的には改善の余地あり）。pell4a が実用的でなくなる最初の数は d = 421 で、探したい数は17桁（「その27参照」）。今はまだ330だし、この全数検索はただの遊びなので、多分そのうち飽きるだろう。本当の目標は、とりあえず一般化ペル方程式の基本解の上限について、すっきりさせること。

2020-06-27　現金でちょうど50円払う方法は37通り

1円、5円、10円、50円という通貨があるとする。ちょうど50円を払う方法は何通りあるか。1円を50個、1円を45個・5円を1個…といった全ての可能性を考える。

50円を1個という単純な方法（1通り）の他に、下記のように
11 + 9 + 7 + 5 + 3 + 1 = (11+1) + (9+3) + (7+5) = 12*3 = 36
通りの方法がある。一見たわいもない算数だが、通貨の種類と支払う額を一般化した場合（例: 63円、84円、94円、100円の切手を組み合わせてちょうど1万円払う方法は何通りあるか）、解けることは解けるものの、計算量的に困難。計算法（アルゴリズムの実装）もいろいろで、興味深い。この問題の有名なバージョンは「1セント、5セント、10セント、25セント、50セント、1ドルを使って1ドル払う方法」。50円バージョンは37通り、1ドルバージョンは293通りと、どちらも答えが素数になるところが、誘惑的。

  1 : [ 10*0 + 5*0 + 50 = 50 ] 10円0枚なら残り50円: 5円は0～10枚＝11通り
  2 : [ 10*0 + 5*1 + 45 = 50 ] （不足分は1円で払う。
  3 : [ 10*0 + 5*2 + 40 = 50 ] 1円の枚数は自動的に決まる）
  4 : [ 10*0 + 5*3 + 35 = 50 ]
  5 : [ 10*0 + 5*4 + 30 = 50 ]
  6 : [ 10*0 + 5*5 + 25 = 50 ]
  7 : [ 10*0 + 5*6 + 20 = 50 ]
  8 : [ 10*0 + 5*7 + 15 = 50 ]
  9 : [ 10*0 + 5*8 + 10 = 50 ]
 10 : [ 10*0 + 5*9 + 5 = 50 ]
 11 : [ 10*0 + 5*10 + 0 = 50 ]
 12 : [ 10*1 + 5*0 + 40 = 50 ] 10円1枚なら残り40円: 5円は0～8枚＝9通り
 13 : [ 10*1 + 5*1 + 35 = 50 ]
 14 : [ 10*1 + 5*2 + 30 = 50 ]
 15 : [ 10*1 + 5*3 + 25 = 50 ]
 16 : [ 10*1 + 5*4 + 20 = 50 ]
 17 : [ 10*1 + 5*5 + 15 = 50 ]
 18 : [ 10*1 + 5*6 + 10 = 50 ]
 19 : [ 10*1 + 5*7 + 5 = 50 ]
 20 : [ 10*1 + 5*8 + 0 = 50 ]
 21 : [ 10*2 + 5*0 + 30 = 50 ] 10円2枚なら残り30円: 5円は0～6枚＝7通り
 22 : [ 10*2 + 5*1 + 25 = 50 ]
 23 : [ 10*2 + 5*2 + 20 = 50 ]
 24 : [ 10*2 + 5*3 + 15 = 50 ]
 25 : [ 10*2 + 5*4 + 10 = 50 ]
 26 : [ 10*2 + 5*5 + 5 = 50 ]
 27 : [ 10*2 + 5*6 + 0 = 50 ]
 28 : [ 10*3 + 5*0 + 20 = 50 ] 10円3枚なら残り20円: 5円は0～4枚＝5通り
 29 : [ 10*3 + 5*1 + 15 = 50 ]
 30 : [ 10*3 + 5*2 + 10 = 50 ]
 31 : [ 10*3 + 5*3 + 5 = 50 ]
 32 : [ 10*3 + 5*4 + 0 = 50 ]
 33 : [ 10*4 + 5*0 + 10 = 50 ] 10円4枚なら残り10円: 5円は0～2枚＝3通り
 34 : [ 10*4 + 5*1 + 5 = 50 ]
 35 : [ 10*4 + 5*2 + 0 = 50 ]
 36 : [ 10*5 + 5*0 + 0 = 50 ] 10円5枚なら残り0円: 5円は0枚＝1通り

このタイプの問題は「危険」。入り口はスイートでエレガントで簡単そうだが、スッキリ解決させようとすると予想外に深く、最初に考えていたことは氷山の一角で、一般の問題が面白くてハマってしまう…
Lectures on Integer Partitions
http://www.math.upenn.edu/%7Ewilf/PIMS/PIMSLectures.pdf