前回(2009年11月18日の日記参照)は経過日から年の計算を行いました。今回はそのテストを行います。
そもそもこの問題を解き始めた動機は、年月日から経過日への変換関数はあるけど、逆はないよね?でした。要するに世の中には年月日から経過日を計算する、信頼できる関数があるってことです。
テストの方法は下記の通り、
経過日 --(今回作成の関数)-> 年月日 --(実績ある関数)-> 経過日
として、同じ経過日が出てきたらOK、違っていたらNGです。
年月日から経過日(修正ユリウス日)への変換関数にはフリーゲルの公式を選択しました。
(Wikipediaユリウス通日の項から引用)
グレゴリオ暦y年m月d日午前0時の修正ユリウス日は、x以下で最大の整数をfloor(x)で表すと、
floor(365.25 * y) + floor(y / 400) - floor(y / 100) + floor(30.59 * (m - 2)) + d - 678912
実装する関数の仕様は下記の通りです。
この処理をC言語で書くと下記のようになります。floatへのキャストがかなりウザいですが、気にしないでください…。
int cal_to_mjd(int year, int month, int day, int *mjd)
{
int chk[] = {
31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31
};
int u, d;
if (month < 1 || 12 < month ||
day < 1 || 31 < day) {
return -1;
}
if (year % 400 == 0) {
u = 1;
} else if (year % 100 == 0) {
u = 0;
} else if (year % 4 == 0) {
u = 1;
} else {
u = 0;
}
if (chk[month - 1] + u < day) {
return -1;
}
if (month == 1 || month == 2) {
year -= 1;
month += 12;
}
if (year < 1) {
return -1;
}
d = (int)(
floor(365.25 * (float)year) +
floor((float)year / 400.0) +
- floor((float)year / 100.0) +
floor(30.59 * ((float)month - 2.0)) +
(float)day - 678912.0 );
if (mjd) {
*mjd = d;
}
return 0;
}
テストを行う前に、私の作った関数とフリーゲルの公式が取る「引数の差」について考える必要があります。
引数 | 私の作った関数 | フリーゲルの公式 |
---|---|---|
経過日 | 400で割り切れる年(基準年とする) の3月1日を0日とする経過日 |
修正ユリウス日 (1858年11月27日を0日とする経過日) |
年 | 基準年からの経過年 | 西暦(グレゴリオ歴) |
月 | 月(3月〜14月) | 月(1月〜12月、グレゴリオ歴) |
日 | 日(1日〜31日、グレゴリオ歴) | 日(1日〜31日、グレゴリオ歴) |
以上から、私の作った関数とフリーゲルの公式の引数では、経過日、年、月の3つが異なっています。まずはこの差を埋める補助関数を作ります。
経過日については、修正ユリウス日から適当な値(でたらめという意味ではありません)を引いて、3月1日を0とするような日に直してあげます。
基準とする日は400年で割り切れる年ならどこでも良いです。しかしグレゴリオ暦の採用年が1582年以降であることを考えると、そこ付近のグレゴリオ暦にはあまり意味がありません。1600年あたりを開始日とするのが妥当でしょう。
経過日0日が表す1600年3月1日の修正ユリウス日は -94493日です(フリーゲルの公式より、計算省略)。修正ユリウス日に94493を足せば経過日へ変換できます。
年については、経過日をグレゴリオ暦の1600年を基準としたので、結果に1600を足せばグレゴリオ暦の年に変換できます。
月については、3月〜12月まではそのまま、13月、14月を翌年の1月、2月と考えます。これで経過日、年月日ともにフリーゲルの公式と揃えることができました。
実装する関数の仕様は下記の通りです。
この処理をC言語で書くと下記のようになります。
int mjd_to_cal(int mjd, int *year, int *month, int *day)
{
int base, date, y, m, d, mod;
int result;
base = 1600;
date = mjd + 94493; //date 0 is 1600/3/1
date += (1600 - base) / 400 * 146097; //date 0 is base/3/1
result = date_to_year(date, &y, &mod);
if (result == -1) {
return result;
}
result = date_to_month(mod, &m, &d);
if (result == -1) {
return result;
}
y += base;
if (m > 12) {
y += 1;
m -= 12;
}
if (year) {
*year = y;
}
if (month) {
*month = m;
}
if (day) {
*day = d;
}
return 0;
}
上記で説明したこと以外に、mjd_to_cal関数では基準年を変更できるようになっています。base = 1200にすれば1200年3月1日以降の修正ユリウス日を扱えます。ただし負の数、0年には対応していないため、最小値は400年3月1日(base = 400)です。
きちんと変換できるか、異常値に対してエラーを返すかどうか、-94495(1600年3月1日の2日前)から944929(4446年1月2日)までの経過日を与えてテストします。
int main()
{
int y, m, d, mjd;
int f, result, i;
int tmin, tmax;
f = 0;
tmin = -94495;
tmax = 944930;
for (i = tmin; i < tmax; i++) {
result = mjd_to_cal(i, &y, &m, &d);
if (result == -1) {
printf("mjd->date:%3d -> error\n", i);
continue;
}
result = cal_to_mjd(y, m, d, &mjd);
if (result == -1) {
printf("date->mjd:%4d/%2d/%2d -> error\n", y, m, d);
continue;
}
if (i != mjd) {
printf("mismatch!!\n");
printf("mjd->date:%5d -> %4d/%2d/%2d\n", i, y, m, d);
printf("date->mjd:%4d/%2d/%2d -> %5d\n", y, m, d, mjd);
f = 1;
}
}
if (f == 0) {
printf("Test Passed, mjd [%d, %d]\n", tmin, tmax);
}
return 0;
}
実行結果は下記の通りです。
mjd->date:-94495 -> error mjd->date:-94494 -> error Test Passed, mjd [-94495, 944930]
1600年3月1日より前の日付はエラーになり、それ以外はテストにパスしました。どうやら正しく変換できているようです。
ネットで調べていたら、経過日から月日へ変換する際にループを回すのは遅くてナンセンスという記述を見つけてしまいました。な、なんだってぇー!?
高々12回、平均6回のループがそんなに遅いだろうか…?計算一発で出す方式を考えて、ループ方式と速度比較しようと思います。
< | 2009 | > | ||||
<< | < | 11 | > | >> | ||
日 | 月 | 火 | 水 | 木 | 金 | 土 |
1 | 2 | 3 | 4 | 5 | 6 | 7 |
8 | 9 | 10 | 11 | 12 | 13 | 14 |
15 | 16 | 17 | 18 | 19 | 20 | 21 |
22 | 23 | 24 | 25 | 26 | 27 | 28 |
29 | 30 | - | - | - | - | - |
合計:
本日:
管理者: Katsuhiro Suzuki(katsuhiro( a t )katsuster.net)
This is Simple Diary 1.0
Copyright(C) Katsuhiro Suzuki 2006-2023.
Powered by PHP 8.2.20.
using GD 2.3.3(png support.)