コグノスケ

子供の頃の夢

社会人になってすっかり忘れていたのですが、子供の頃、ラジコンが欲しかったのです。特に車タイプと空飛ぶタイプが欲しかったのです。

車タイプは大学生の時に買いました。部屋の中でも遊べるくらいの小さいやつを買ったのですが、モーター音があまりにうるさくて近所迷惑だったのですぐ止めました…。

今日、ヘリコプターのラジコンを、しかも「室内用」と書いてあるラジコンを見つけたので思わず買ってしまいました。単三電池がたくさん必要と書いてあったので、充電式エボルタ8本と一緒にまとめ買いです。

いざ寮の廊下で飛ばしてみると、結構楽しいのですが、やはりモーター音がうるさくて、かなり近所迷惑な音がします。

寮の部屋でもやってみましたが、やはりうるさいのと、狭すぎてすぐに事故って全然楽しくありません。いったんぶつかると吹っ飛んで、ピンボールのようにあちこちにぶつかっります。今にも壊れそうです。

「空を自由に飛びたいな〜。」「えっ？何！？聞こえないよ！！」

飛行機もヘリコプターも、ものすごい音たてて飛んでいます。エンジン近くでは会話が不可能なほどの爆音です。

空を自由に飛ぶと言えば、ドラえもんに登場するタケコプターですが、作中では人間一人を軽々と持ち上げ、それでもかなりの加速力を発揮しています。

飛行機やヘリの騒音から想像するに、タケコプターの騒音も相当なもので、飛行中の会話なんて不可能ではないでしょうか？

でもドラえもんやのび太はタケコプターで飛行中でも、地上に居るときと同じように会話しています。未来の道具はとても静かなのでしょう、羨ましいですね…。

誕生日？焼き肉だろ！

同期の誕生祝いに焼き肉屋に行きました。

同期で誕生祝いをするときはケーキを買って行って、食事の最後に出すのですが、今回は焼き肉とケーキというヘビーすぎる組み合わせになりました。

店のチョイスからしておかしくね？というツッコミはその通りな気がしますが、こまけぇこたぁいいんだよ。みんな食ってたからいいんです。

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

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

前回（2009年11月16日の日記参照）は経過日から月日への計算を行いました。残るは年の計算です。

グレゴリオ歴の周期性

前回も紹介したとおり、グレゴリオ暦の平年、閏年には下記の法則があります。

1. 平年は365日、閏年は366日
2. 西暦が4で割り切れる年は閏年(1年が366日)
3. 西暦が100で割り切れる年は平年(1年が365日)
4. 西暦が400で割り切れる年は閏年(1年が366日)

4年周期は平年、平年、平年、閏年のパターンです。

100年周期は4年周期を25回繰り返すだけですが、最後の100年目だけは閏年ではありません。

400年周期は100年周期を4回繰り返すだけですが、最後の400年目だけは閏年になります。400年以上を扱うルールはありませんので、以降400年周期で同じパターンが続きます。

図示すると下記のようになります。オレンジの四角が閏年を表しています。

閏年のパターン

各周期の日数の計算式は下記の通りです。

1. 1年は365日
2. 4で割り切れるときは閏年、つまり最後の年だけ +1日
   365+365+365+366
   =365*4+1
   =1461
3. 100で割り切れるときは平年、つまり最後の年だけ -1日
   1461+1461+…++1461+1460
   =1461*25-1
   =36524
4. 400で割り切れるときは閏年、つまり最後の年だけ +1日
    36524+36524+36524+36525
   =36524*4+1
    =146097

パターンがわかってしまえば、与えられた経過日にこのパターンがいくつ含まれているか？を計算するのみです。

経過日から年への変換

では実際に経過日から年を計算してみます。しつこいですがパターンは下記4つです。経過日にこれらのパターンがいくつ含まれているか調べます。

• パターンA: 400年 = 146097日
• パターンB: 100年 = 36524日
• パターンC: 4年 = 1461日
• パターンD: 1年 = 365日

なぜ上記パターンの数を求めれば年数が出るのか？がわかる人は次の章を飛ばして読んでください。

桁の話

ここでは桁の概念についてと、年への換算にどう使うか？を補足させていただきます。

私たちにおなじみの10進数は1234のように書きますが、これってどういう意味でしょうか？

10進数の1桁目は10^0がいくつ含まれているか、2桁目は10^1がいくつ含まれているかを示します。n桁目は10^(n-1) がいくつ含まれているか？を示します。

ですから10進数で1234は（1 * 10^3 + 2 * 10^2 + 3 * 10^1 + 4 * 10^0 = 1234）という数（10進数表記）を表します。

当たり前？でもこれは m進数でも同様で、n桁目はm^(n-1) がいくつ含まれているか？を示します。

ですから8進数の1234は（1 * 8^3 + 2 * 8^2 + 3 * 8^1 + 4 * 8^0 = 668）という数（10進数表記）を表します。

経過日から年を求める場合も、これと似た考え方ができます。各パターンの数を桁と見なして（400年の桁、100年の桁、4年の桁、1年の桁）、10進数へと換算してやればよいのです。

各パターンの数（400年、100年、4年、1年）が1, 2, 3, 3であれば（1 * 400 + 2 * 100 + 3 * 4 + 3 * 1 = 615）という年（10進数表記）を表します。

このように各パターンの数を求めてあげることで、年数が計算可能なのです。

経過日から年への計算式

パターンAの数は（経過日 / 146097）です。余りの経過日は400年未満のどこか（0年3月1日〜399年2月29日）を表します。これはパターンBの計算に回します。

パターンBの数は（経過日 / 36524）です。余りの経過日は100年未満のどこか（0年3月1日〜99年2月28日）を表します。余りは同様にパターンCの計算に回します。

しかし146096日（399年2月29日）の場合に3となるべきところが4になり間違えてしまうため、特別に3とします。その際の余り経過日は99年2月29日ですから36525 - 1 = 36524日（開始が0日のため1引く）となります。

パターンCの数は（経過日 / 1461）です。余りの経過日は4年未満のどこか（0年3月1日〜3年2月29日）を表します。余りは同様にパターンDの計算に回します。

パターンDの数は（経過日 / 365）です。余りの経過日は1年未満のどこか（0年3月1日〜0年2月28日）を表します。ここの余りは経過日から月日を求める計算に回します。

しかし1460日（3年2月29日）の場合に3となるべきところが4になり間違えてしまうため、特別に3とします。その際の余り経過日は0年2月29日ですから366 - 1 = 365日（開始が0日のため1引く）となります。

実装例

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

date

400で割り切れる西暦（1600年、2000年など）の3月1日からの経過日を渡します。値域は正の整数です。

year

年を返します。値域は正の整数です。

mod

1年に満たない経過日を返します。値域は0〜365です。

返り値

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

この処理をC言語で書くと下記のようになります。

経過日から年への変換、実装例

int date_to_year(int date, int *year, int *mod)
{
        int a400, a100, a4, a;
        int m400, m100, m4, m;

        if (date < 0) {
                return -1;
        }

        a400 = date / 146097;
        m400 = date % 146097;
        if (m400 == 146096) {
                a100 = 3;
                m100 = 36524;
        } else {
                a100 = m400 / 36524;
                m100 = m400 % 36524;
        }
        a4 = m100 / 1461;
        m4 = m100 % 1461;
        if (m4 == 1460) {
                a = 3;
                m = 365;
        } else {
                a = m4 / 365;
                m = m4 % 365;
        }

        if (year) {
                *year = a400 * 400 + a100 * 100 + a4 * 4 + a;
        }
        if (mod) {
                *mod = m;
        }

        return 0;
}

これで経過日から年への変換は完成ですが、まだテストが終わっていません。

月日の変換のときは、入力に対する答えがたかだか366通りと少なかったため、目でチェックしました。しかし今回は400年つまり146097日もの数をチェックしなければなりません。こんな数を目で一々チェックしていたら頭がおかしくなります。

テストの方法についてはまた次回に。

経過日を年月日に変換

ある時点からの経過日（0基点）を年月日にする、なんていかにも学校の授業で出てきそうなアルゴリズムですが、世間では需要がないのか、名前が付いているようなアルゴリズムがありません。

練習問題にちょうど良かったので、いっちょ考えてみました。鈴木アルゴリズム完成！なーんて言えれば良かったんですけど、大した物はできませんでした。

その代わりといってはなんですが、考え方を細かくメモっていきたいと思います。今日は月日の変換、後日に年の変換を扱います。

文章をうだうだ読むのが面倒くせえ！って人のために、実装例も載せる予定です。ライセンスを特に明示していなければ、修正BSDライセンスを適用します。

グレゴリオ歴の性質

年月日のルールを決めるのは暦です。現在採用されている暦の主流は1582年頃に採用が始まった「グレゴリオ歴」です。下記の性質があります。

• 1年は365日（平年）か366日（閏年）
• 1年は12ヶ月
• 1月は28日（平年2月）29日（閏年2月）30日（4,6,9,11月）31日（1,3,5,7,8,10,12月）のどれか
• 1日は24時間
• 1時間は60分
• 1分は60秒

閏年は暦と地球の季節がずれないようにする仕組み（※1）です。もっと簡単に言うと、平年より1日長い年のことです。平年/閏年は下記のルールで決まります。

1. 西暦が4で割り切れる年は閏年(1年が366日)
2. 西暦が100で割り切れる年は平年(1年が365日)
3. 西暦が400で割り切れる年は閏年(1年が366日)

変なルールに見えるかもしれませんが、人類が苦労して「時間とはなんぞや？暦とはなんぞや？」を定義してきた証なのでしょう。

（※1）地球の公転周期（1年）は自転周期（1日）の365倍と1/4日くらいです。単純に暦を1年 = 365日としてしまうと4年で1日分、季節と暦がずれてしまいます。その差を補正するために閏年が作られました。

邪魔な2月

初めに経過日→月日の変換を考えます。その際に邪魔なのが年によって長さが変わる2月です。

例として「1/1から60日後は何月何日か？」を考えてみます。平年（2月が28日間）ならば答えは3/2、閏年（2月が29日間）ならば答えは3/1です。2月の存在によって、経過日と月の対応関係がずれてしまうのです。

年によって計算方法を変えるのは面倒です。以下のように経過日の基準点を変えて一通りの計算方法で、経過日→月日を計算できるようにします。

3/1を経過日0日とし、2月を最後に持って行く

つまり3月を一年の始まり（正確には3/1を経過日0日とする）とし、来年の2月を当年の最終月14月として考えます。こうすると経過日と各月の対応関係が固定されるため、2月の長さが変化しても計算に影響が出ません（後ほどまた説明します）

経過日から月日への変換

では実際に経過日から月日を計算してみます。使うのは下記のテーブル（※2）です。

テーブルの2列目は「その月の1日を表す経過日」を表しています。調べたい経過日と、2列目で大小比較すれば、調べたい経過日が何月なのかがわかります。何月なのかわかれば何日か？も簡単に分かります。

例として100を考えましょう。
まず92（6/1）< 100 < 122（7/1）から6月であることがわかります。
さらに100 - 92（6/1）= 8ですから6/1から8日後、つまり6/9だとわかります。

この計算方法の利点は、各月の1日に対応する経過日だけ考えれば計算できることです。さきほど2月の長さは気にしなくて良いと言った理由はここにあります。

2月は最終月のため、長さが変化しても各月の1日に対応する経過日は変化しません。というより各月の1日に対応する経過日を変化させないようにするために、わざわざ2月を最後に回したのです。

実装例

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

date

3/1からの経過日を渡します。値域は0〜365です。

month

変換後の月を返します。値域は3〜14（翌年2月）です。

days

変換後の日を返します。値域は1〜31です。

返り値

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

この処理をC言語で書くと下記のようになります。

経過日から月日への変換、実装例

int date_to_month(int date, int *month, int *days)
{
        static int m_date[] = {
                0, 31, 61, 92, 122, 153, 184, 214, 245, 275, 306, 337,
        };
        int m, d;
        int i;

        if (date < 0 || 365 < date) {
                return -1;
        }

        for (i = 0; i < 12; i++) {
                if (date < m_date[i]) {
                        break;
                }
        }
        m = i + 2;
        d = date - m_date[i - 1] + 1;

        if (month) {
                *month = m;
        }
        if (days) {
                *days = d;
        }

        return 0;
}

負の値や、1年（最長366日）以上の日数に対しては正常に動作しません。エラー処理として366以上の値や負の値を渡したときにエラーを返すこととします。

きちんと変換できるか、異常値に対してエラーを返すかどうか、-3から369までの経過日を与えてテストします。

経過日から月日への変換、テスト関数

int main()
{
        int m, d;
        int result, i;

        for (i = -3; i < 370; i++) {
                result = date_to_month(i, &m, &d);
                if (result == -1) {
                        printf("date:%3d -> error\n", i);
                        continue;
                }
                printf("date:%3d -> %2d/%2d\n", i, m, d);
        }

        return 0;
}

実行結果は下記の通りです。適当に端折ってあります。

経過日から月日への変換、実行結果

date: -3 -> error
date: -2 -> error
date: -1 -> error
date:  0 ->  3/ 1
date:  1 ->  3/ 2
...（略）...
date:100 ->  6/ 9
date:101 ->  6/10
date:102 ->  6/11
...（略）...
date:364 -> 14/28
date:365 -> 14/29
date:366 -> error
date:367 -> error
date:368 -> error
date:369 -> error

この関数は月のパラメータに13月やら14月を返します。しかし後ほど翌年の1月、2月へ変換してつじつまを合わせますので、ここでは何も変換しません。

ある月の1日を表す経過日を求めるには、下記の漸化式を用います。
（ある月の1日を表す経過日） = （前月の1日を表す経過日） + （前月の長さ）
要はテーブルの2列目（各月1日の経過日）と3列目（月の長さ）を足すと、次の月の1日を表す経過日が計算できるということです。

また今度

経過日から月日へ変換できたところで、また今度。次は経過日から年への変換を書く予定です。