コグノスケ

経過日を年月日に変換 - その3

前回（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

実装する関数の仕様は下記の通りです。

year

グレゴリオ歴の西暦を渡します。値域は1以上です。

month

グレゴリオ歴の月を渡します。値域は1〜12です。ただし1年の場合は3〜12です。

days

グレゴリオ歴の日を渡します。値域は1〜31です。グレゴリオ暦にない日（閏年でない年の2月29日など）を渡した場合はエラーになります。

mjd

修正ユリウス日を返します。値域は -678516（1年3月1日）以上 です。

返り値

成功ならば0を返します。エラーが起きた場合は -1を返します。

この処理を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月と考えます。これで経過日、年月日ともにフリーゲルの公式と揃えることができました。

実装する関数の仕様は下記の通りです。

mjd

修正ユリウス日を渡します。値域は -94493（1600年3月1日）以上 です。

year

グレゴリオ歴の西暦を返します。値域は1600以上です。

month

グレゴリオ歴の月を返します。値域は1〜12です。

days

グレゴリオ歴の日を返します。値域は1〜31です。

返り値

成功ならば0を返します。エラーが起きた場合は -1を返します。

この処理を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回のループがそんなに遅いだろうか…？計算一発で出す方式を考えて、ループ方式と速度比較しようと思います。