GitHub coq/coq/issues/21452: Termination check passed in Rocq 9.0.0 but not in Rocq 9.1.0
以前から気になっていたというか腑に落ちていなかったこととして、 グレゴリオ歴の年月日からユリウス通日への変換がある。
Wikipedia:ユリウス通日によれば、 1月と2月は前年の13月と14月として、 y年m月d日午前0時の修正ユリウス日は以下で表される、とある
mjd = floor(365.25*y) + floor(y/400) - floor(y/100) + floor(30.59 * (m-2)) + d - 678912
(ユリウス日は正午が整数となってわかりにくいので、深夜0時が整数となる修正ユリウス日を扱うことにする。)
ここで、この式がなにをやっているのかだいたいはわかるのだが、30.59 の 0.59 ってどこから出てきたのだ?
365.25 の 0.25 は 4年に一回の閏年に関連しているのだろう。 y/400 と y/100 も、閏年の規則の 100年に一回と 400年に一回と対応しているのだろう。 678912 はなんらかのオフセットだろうから気にならない。
でも、0.59 はどうして 0.59 なのか?
これど関連して、(結論からいえば非常に強く関連しているのだが) 月の日数は一定じゃないのにこんな式でちゃんと求まるのか?
という疑問が長年あったのだが、思い立って Rocq で確認してみた。
結論からいうと、0.59 というのがうまく月の日数が一定でないのを扱っている。
Rocq でこれを扱うとして、まずちょっと問題になるのが、 大きな整数を扱う必要があることである。 Rocq の nat はペアノの自然数の素朴な実装なので、自然数 n を表現するのに n に比例したメモリを使う。 コンストラクタ O にコンストラクタ S を n 回適用した項を使うので。 まぁ、現実的な話としては数千くらいになるのは避けたい。
ということで、まず 678912 は無視することにする。 これは起点をしめしているだけで、私の疑問とは関係ない。 (あと、-678912 とすると負数を扱わなければならないので nat では表現できなくて困るというのもある)
本当は、年は自然数じゃなくて整数にしたほうがいいのだが、 MathComp で大きい整数をコンパクトに (2進で) 扱う方法はみあたらなかったので nat でやることにする。
ただ、少し大きな値は 2進で表現した項を使う。 Rocq には BinNums というライブラリがあって、そういう表現も可能であり、 ssrnat には、それで表現した数を記述するための %num という記法が用意されている。
が、それが動かない。 いろいろ調べたところ、以下のように最初に書いて動かした。 (そしてバグレポートしした。)
(* workaround for %num. https://github.com/math-comp/math-comp/issues/1512 *) From Stdlib Require Import BinNatDef. Delimit Scope N_scope with num. From mathcomp Require Import all_ssreflect.
まず、閏年かどうかを判定する関数を定義する。 ここで d %| m は、m が d で割りきれることを示す (bool を返す) 述語である。 グレゴリオ歴における閏年の定義そのものであるが、 今回は 1月2月のかわりに 13月14月を扱い、 13月14月は翌年の 1月2月なので、最初に年を補正する必要がある。 MathComp には、bool -> nat の coercion が定義されていて、false が 0, true が 1 になるので、この補正は y + (12 < m) と書ける。
Definition is_leap_year y m := let y' := y + (12 < m) in (4 %| y') && ~~(100%num %| y') || (400%num %| y').
月の長さを求める関数を定義する。 まぁ、グレゴリオ歴のままのふつうである。
Definition monlen (y m : nat) : nat := match m with | 13 => 31 | 14 => if is_leap_year y m then 29 else 28 | 3 => 31 | 4 => 30 | 5 => 31 | 6 => 30 | 7 => 31 | 8 => 31 | 9 => 30 | 10 => 31 | 11 => 30 | 12 => 31 | _ => 0 end.
正しい年月日かどうかを判定する述語を定義する。 (nat つまり 0以上の整数を表現する型を使っているので、正しくない値も表現できてしまうため)
Definition valid_ymd (ymd : nat * nat * nat) : bool := let: (y, m, d) := ymd in (3 <= m <= 14) && (1 <= d <= monlen y m).
Wikipedia日本語版にあった mjd を求める式を (最後の 678912 を除いて) 関数として定義する。
Definition jd_of_ymd (ymd : nat * nat * nat) : nat := let: (y, m, d) := ymd in 365%num * y + (y %/ 4 + y %/ 400%num - y %/ 100%num) + 30 * (m - 2) + (59 * (m - 2) %/ 100%num) + d.
年月日から、翌日の年月日を求める関数を定義する。 日をインクリメントしてオーバーフローしたら月をインクリメントしてそれもオーバーフローしたら年をインクリメントする。 これは素朴に正しいと思う。
Definition succ_ymd_y (y m d : nat) : nat := if d < monlen y m then y else if m < 14 then y else y.+1. Definition succ_ymd_m (y m d : nat) : nat := if d < monlen y m then m else if m < 14 then m.+1 else 3. Definition succ_ymd_d (y m d : nat) : nat := if d < monlen y m then d.+1 else 1. Definition succ_ymd (ymd : nat * nat * nat) := let: (y, m, d) := ymd in (succ_ymd_y y m d, succ_ymd_m y m d, succ_ymd_d y m d).
とりあえず、月の長さは 28日以上、というのを証明する。 月 m を 14回場合分けして証明している。 こうやって場合分けすると、0月から13月までは自明に証明できて、 15月以降も自明に証明できる。 14月 (つまり 2月) は閏年かどうかの場合分けが必要で、場合分けすれば証明できる。
Lemma monlen28 (y m : nat) : 3 <= m <= 14 -> 28 <= monlen y m. Proof. do 14 (case: m => [|m]; first by []). case: m => [|m]; last by []. rewrite /monlen. by case: ifP. Qed.
valid な年月日の翌日の年月日は valid であることを証明する。
Lemma succ_ymd_mrange y m d :
valid_ymd (y, m, d) ->
3 <= succ_ymd_m y m d <= 14.
Proof.
rewrite /valid_ymd => /andP [Hm Hd].
rewrite /succ_ymd_m.
case: ifPn => [|Hd2]; first by [].
case: ifPn => [Hm2|]; last by [].
apply/andP.
split; last by [].
case/andP: Hm.
by move/ltnW.
Qed.
Lemma succ_ymd_drange y m d :
valid_ymd (y, m, d) ->
1 <= succ_ymd_d y m d <= monlen (succ_ymd_y y m d) (succ_ymd_m y m d).
Proof.
rewrite /valid_ymd => /andP [Hm Hd].
rewrite /succ_ymd_y /succ_ymd_m /succ_ymd_d.
case: ifPn => [|Hd2]; first by [].
case: ifPn => Hm2.
apply/andP; split; first by [].
apply (@leq_trans 28); first by [].
apply monlen28.
apply/andP.
split; last by [].
apply ltnW.
by case/andP: Hm.
apply/andP; split; first by [].
apply (@leq_trans 28); first by [].
by apply monlen28.
Qed.
Lemma valid_succ_ymd ymd : valid_ymd ymd -> valid_ymd (succ_ymd ymd).
Proof.
case: ymd => [[y m] d] H.
rewrite /succ_ymd /valid_ymd /=.
apply/andP; split.
by apply succ_ymd_mrange.
by apply succ_ymd_drange.
Qed.
グレゴリオ歴は400年周期であることを証明する。 400年の中には97回閏年があって、400*365+97=146097日が周期である。
Lemma cycle400 y m d n :
jd_of_ymd (y + 400%num * n, m, d)
= jd_of_ymd (y, m, d) + (400%num * 365%num + 97) * n.
Proof.
rewrite /jd_of_ymd.
rewrite [365%num * (y + 400%num * n)]mulnDr.
rewrite [LHS]addn.[ACl 3*4*5*6*2*1].
rewrite [RHS]addn.[ACl 2*3*4*5*6*1].
congr (_ + 365%num * y).
rewrite (_ : (400%num * 365%num + 97) * n = 97 * n + 365%num * (400%num * n)); last first.
by rewrite mulnDl addnC [400%num * 365%num]mulnC mulnA.
rewrite addnA.
congr (_ + 365%num * (400%num * n)).
rewrite [RHS]addn.[ACl 1*5*2*3*4].
congr (_ + 30 * (m - 2) + (59 * (m - 2)) %/ 100%num + d).
rewrite [in (y + 400%num * n) %/ 4](_ : 400%num * n = n * 100%num * 4); last first.
by rewrite -[RHS]mulnA mulnC.
rewrite divnDMl; last by [].
rewrite [in (y + 400%num * n) %/ 400%num]mulnC.
rewrite divnDMl; last by [].
rewrite [in (y + 400%num * n) %/ 100%num](_ : 400%num * n = n * 4 * 100%num); last first.
by rewrite -mulnA mulnC.
rewrite divnDMl; last by [].
rewrite [y %/ 4 + n * 100%num + (y %/ 400%num + n)]addnACA.
rewrite -mulnSr.
rewrite subnDAC.
rewrite -addnBA; last first.
by rewrite leq_mul2l orbT.
rewrite (_ : n * 100%num.+1 - n * 4 = 97 * n); last first.
by rewrite -mulnBr mulnC.
rewrite addnBAC; first by [].
apply (@leq_trans (y %/ 4)); last by apply leq_addr.
rewrite -[y %/ 4](@divnMr 25); last by [].
apply leq_div2r.
by apply leq_pmulr.
Qed.
ここが個人的に腑に落ちなかった部分の確認である。0.59 * (m - 2) に月の日数を加えると 0.59 * (m.+1 - 2) になるのである。 ここでも 14回場合分けして証明している。 この場合分け (と内部的な計算) による証明は、正しいという確信は得られるが、 なぜ正しいのかという理由は正直わからない。 でも、(わかりたければ) ここに集中して考えればいいということはわかった。 (たぶんグラフを描くともうちょっとわかる気がするのだが...)
Lemma add_monlen y m :
3 <= m < 14 ->
(59 * (m - 2)) %/ 100%num + monlen y m
= (59 * (m.+1 - 2)) %/ 100%num + 30.
Proof.
do 14 (case: m => [|m]; first by []).
by [].
Qed.
補題を定義する。 この補題が MathComp にあってもおかしくない気がするのだが、ないということはなんかうまい証明のやりかたがあるのだろうか。
Lemma dvdn_mull (d1 d2 m : nat) : d1 * d2 %| m -> d1 %| m. Proof. move/dvdnP => [] k. rewrite mulnA => H. apply/dvdnP. exists (k * d2). rewrite H. by rewrite mulnAC. Qed. Lemma dvdn_mulr (d1 d2 m : nat) : d1 * d2 %| m -> d2 %| m. Proof. rewrite mulnC. by apply dvdn_mull. Qed.
ここで、年月日から翌日の年月日を求めてユリウス日を求めるのと、もともとの年月日からユリウス日を求めて 1 を足すのは等しい、ということを証明する。 けっこう長い証明になった。
難しかったのは、閏年であるという前提 (4 %| y.+1) && ~~ (100%num %| y.+1) || (400%num %| y.+1) があるときに (4 %| y.+1) && ~~ (100%num %| y.+1) の場合と (400%num %| y.+1) の場合に場合分けするとうまくいかないというところがあった。 前者の場合に y.+1 が 400で割りきれるかどうかの情報がないので、証明が進まない。 これはまず (4 %| y.+1) && ~~ (100%num %| y.+1) && ~~ (400%num %| y.+1) || (400%num %| y.+1) に変形して、 (4 %| y.+1) && ~~ (100%num %| y.+1) && ~~ (400%num %| y.+1) と (400%num %| y.+1) に場合分けしないといけないのだった。
Lemma jd_of_ymd_succ y m d :
valid_ymd (y, m, d) ->
jd_of_ymd (succ_ymd (y, m, d)) = (jd_of_ymd (y, m, d)).+1.
Proof.
rewrite /valid_ymd => /andP [Hm Hd].
rewrite /succ_ymd /succ_ymd_y /succ_ymd_m /succ_ymd_d.
case: ifPn => Hd2.
rewrite [in LHS]/jd_of_ymd.
by rewrite addnS.
case: ifPn => Hm2.
have {Hd Hd2}Hdmax : d = monlen y m.
rewrite -leqNgt in Hd2.
move: Hd => /andP [Hd3 Hd4].
apply/eqP.
rewrite eqn_leq.
by rewrite Hd2 Hd4.
have {}Hm : 2 < m.
by case/andP: Hm.
rewrite /jd_of_ymd.
rewrite [LHS]addn1.
congr _.+1.
rewrite -[LHS]addnA.
rewrite [RHS]addn.[ACl 1*2*(3*(4*5))].
congr (365%num * y + (y %/ 4 + y %/ 400%num - y %/ 100%num) + _).
rewrite [in RHS]Hdmax.
rewrite (add_monlen y m); last by rewrite Hm2 Hm.
rewrite [(59 * (m.+1 - 2)) %/ 100%num + 30]addnC.
rewrite [RHS]addnA.
congr (_ + (59 * (m.+1 - 2)) %/ 100%num).
rewrite subSn; last first.
by apply ltnW.
by rewrite mulnSr.
have {Hd Hd2}Hdmax : d = monlen y m.
rewrite -leqNgt in Hd2.
move: Hd => /andP [Hd3 Hd4].
apply/eqP.
rewrite eqn_leq.
by rewrite Hd2 Hd4.
have {Hm Hm2}Hmmax : m = 14.
rewrite -leqNgt in Hm2.
move: Hm => /andP [Hm3 Hm4].
apply/eqP.
rewrite eqn_leq.
by rewrite Hm4 Hm2.
rewrite {}Hdmax {}Hmmax.
rewrite /jd_of_ymd.
rewrite addn1.
congr _.+1.
change (30 * (3 - 2)) with 30.
change ((59 * (3 - 2)) %/ 100%num) with 0.
change (30 * (14 - 2)) with (nat_of_bin 360%num).
change ((59 * (14 - 2)) %/ 100%num) with 7.
rewrite addn0.
rewrite -[_ + 360%num + 7]addnA.
change (360%num + 7) with (nat_of_bin 367%num).
rewrite mulnS.
rewrite [LHS]addn.[ACl 3*4*1*2].
rewrite [RHS]addn.[ACl 2*3*4*1].
congr (_ + 365%num * y).
rewrite -[_ + 30 + 365%num]addnA.
change (30 + 365%num) with (28 + 367%num).
rewrite addnA [RHS]addnAC.
congr (_ + 367%num).
rewrite /monlen.
case: ifPn => Hleap.
suff: (y.+1 %/ 4 + y.+1 %/ 400%num - y.+1 %/ 100%num) =
(y %/ 4 + y %/ 400%num - y %/ 100%num) + 1.
move=> ->.
by rewrite -addnA.
rewrite addnBAC; last first.
apply (@leq_trans (y %/ 4)); last by apply leq_addr.
by apply leq_div2l.
rewrite [y.+1 %/ 100%num]divnS; last by [].
rewrite subnDA.
congr (_ - y %/ 100%num).
rewrite divnS; last by [].
rewrite divnS; last by [].
move: Hleap.
rewrite /is_leap_year.
change (y + (12 < 14)) with (y + 1).
rewrite addn1.
rewrite (_ : (4 %| y.+1) && ~~ (100%num %| y.+1) || (400%num %| y.+1)
= (4 %| y.+1) && ~~ (100%num %| y.+1) && ~~ (400%num %| y.+1) || (400%num %| y.+1)); last first.
by case: (4 %| y.+1); case: (100%num %| y.+1); case: (400%num %| y.+1).
case/orP.
move=> /andP [] /andP [] H4 H100 H400.
rewrite H4 (negbTE H100) (negbTE H400) /=.
rewrite add0n subn0.
by rewrite addn.[ACl 2*3*1].
move=> H400.
have H4: 4 %| y.+1.
by apply (dvdn_mulr 100%num).
have H100: 100 %| y.+1.
by apply (dvdn_mulr 4).
rewrite H4 H100 H400 /=.
by rewrite 2!add1n addn1 addnS.
congr (_ + 28).
rewrite divnS; last by [].
rewrite divnS; last by [].
rewrite divnS; last by [].
move: Hleap.
rewrite /is_leap_year.
change (y + (12 < 14)) with (y + 1).
rewrite addn1.
rewrite negb_or negb_and negbK.
move/andP => [] H H400.
rewrite (negbTE H400) add0n {H400}.
move: H.
rewrite (_ : (~~ (4 %| y.+1) || (100%num %| y.+1)) =
((~~ (4 %| y.+1) && ~~(100%num %| y.+1)) || (100%num %| y.+1))); last first.
by case: (4 %| y.+1); case: (100%num %| y.+1).
case/orP.
case/andP=> H4 H100.
by rewrite (negbTE H4) (negbTE H100) 2!add0n.
move=> H100.
have H4: 4 %| y.+1.
by apply (dvdn_mulr 25).
rewrite H4 H100.
by rewrite add1n addSn add1n subSS.
Qed.
最後に、年月日の n日後を、翌日を求める関数を n回適用して求め、そこからユリウス日を求めるのと、もともとの年月日からユリウス日を求めて n を足すのは等しい、ということを証明する。 これは帰納法を使う。 上の 1日後の定理があるので難しくない。 (でも求まった年月日が valid であるという部分を追加して命題を強化しないといけないという話はあった)
Lemma jd_of_ymd_add y m d n :
valid_ymd (y, m, d) ->
jd_of_ymd (iter n succ_ymd (y, m, d)) = (jd_of_ymd (y, m, d)) + n.
Proof.
move=> Hymd.
suff: valid_ymd (iter n succ_ymd (y, m, d)) /\
jd_of_ymd (iter n succ_ymd (y, m, d)) = jd_of_ymd (y, m, d) + n.
by case.
elim: n y m d Hymd.
move=> y m d Hymd.
by rewrite addn0.
move=> n IH y m d Hymd.
rewrite addnS.
specialize (IH y m d Hymd).
move: IH.
rewrite iterS.
case: (iter n succ_ymd (y, m, d)) => -[] y2 m2 d2 [Hymd2 Hjd].
split.
by apply valid_succ_ymd.
rewrite -Hjd.
by apply jd_of_ymd_succ.
Qed.
ともあれ、グレゴリオ歴の年月日からユリウス日を求める式のが正しいという確信は得られた。
[latest]