読者です 読者をやめる 読者になる 読者になる

kanetaiの二次記憶装置

プログラミングに関するやってみた、調べた系のものをQitaに移して、それ以外をはてブでやる運用にしようと思います。http://qiita.com/kanetai

日付関連

リポジトリ

閏年(leap year)

ユリウス暦(Julian calendar)は、紀元前46年、古代ローマで採用された。4年に1回、西暦年が4で割り切れる年を閏年としていた。

現行のグレゴリオ暦(Gregorian calendar)は1852年10月15日を金曜日として施行された。
グレゴリオ暦では、

  • 西暦年が4で割り切れる年は閏年
  • ただし、西暦年が100で割り切れる年は平年
  • ただし、西暦年が400で割り切れる年は閏年
コード 閏年判定
public static boolean isLeapYear(int y){
	//if( y < 0 ) y = -y; //グレゴリオ暦のまま紀元前までさかのぼって考えるとき 1 BD-> y=0, 2 BD-> y=-1, ... ,n BD-> y=1-n
	return y % 4 == 0 && y % 100 != 0 || y % 400 == 0;
}

Fairfieldの公式(Fairfield's congruence)

1年1月1日年から y m d日までの日数 Dは、
 D =  365y + \lfloor \frac{y}{4} \rfloor - \lfloor \frac{y}{100} \rfloor + \lfloor \frac{y}{400} \rfloor + \lfloor \frac{306(m+1)}{10} \rfloor +d - 428
特に、 h = D  mod 7とすると、曜日が次のように求まる。
 h = 0 → Sun, 1 → Mon,  2 → Tue,  3 → Wed,  4 → Thu,  5 → Fri,  6 → Sat
ただし、求めたい日の月が1月、2月の場合はそれぞれ前年の13月、14月とする
( m = 1, 2 の場合は、 y \leftarrow  y - 1, m = \leftarrow m+12とし、1年を、3月1日 〜 14月28日(閏年は29日)と仮定する)。

コード Fairfieldの公式
public static final int Fairfield(int y, int m, int d){
	if (m < 3) { --y; m += 12; }
	return 365*y+y/4-y/100+y/400+306*(m+1)/10+d-428;
}
証明 Fairfieldの公式

0年13月1日(1年1月1日)〜0月14月28日(1年2月28日)の日数は
31+28
1年3月1日 〜  y - 1 年14月末日 の日数は ※閏年を考慮しなければ
 365(y - 1)
0年13月1日〜 y - 1年14月末日の閏年の回数は
 \begin{array}{rl} &\lfloor \frac{1+(y-1)}{4} \rfloor - \lfloor \frac{1+(y-1)}{100} \rfloor + \lfloor \frac{1+(y-1)}{400} \rfloor \\ =& \lfloor \frac{y}{4} \rfloor - \lfloor \frac{y}{100} \rfloor + \lfloor \frac{y}{400} \rfloor \end{array}
 \sum d_{m-1} y年3月1日〜 y m - 1月末日までの日数とすると
 \begin{array}{|c|c|c|c|}\hline m & m-1 & \sum d_{m-1} & \lfloor \frac{306(m+1)}{10} \rfloor - 122 \\\hline \\\hline 3	&	&0	&0\\\hline 4	&3	&31	&31\\\hline 5	&4	&61	&61\\\hline 6 &5	&92	&92\\\hline 7	&6	&122	&122\\\hline 8	&7	&153	&153\\\hline 9	&8	&184	&184\\\hline 10	&9	&214	&214\\\hline 11	&10	&245	&245\\\hline 12	&11	&275	&275\\\hline 13	&12	&306	&306\\\hline 14	&13	&337	&337\\\hline \end{array}
 y年3月1日〜 y m - 1月末日の日数は、
 \lfloor \frac{306(m+1)}{10} \rfloor - 122 \cdots (1)
 y年3月1日〜 y m - 1月末日の日数の計算式は、式(1)以外にも考えられる。
 y m月1日〜 y m d日までの日数は
 d
以上より、
 \begin{array}{rl}D =& 31+28 + 365(y - 1) + \lfloor \frac{y}{4} \rfloor - \lfloor \frac{y}{100} \rfloor + \lfloor \frac{y}{400} \rfloor + \lfloor \frac{306(m+1)}{10} \rfloor + d - 122\\ =& 365y + \lfloor \frac{y}{4} \rfloor - \lfloor \frac{y}{100} \rfloor + \lfloor \frac{y}{400} \rfloor + \lfloor \frac{306(m+1)}{10} \rfloor + d - 428 \end{array}
1582年10月15日を金曜日としてグレゴリオ暦が施行されている。
 h= D|_{y=1582, m=10, d=15} mod 7 = 5となるので、
 h = D  mod 7 = 0 → Sun, 1 → Mon,  2 → Tue,  3 → Wed,  4 → Thu,  5 → Fri,  6 → Sat
(証明終わり)

Zellerの公式(Zeller's congruence)

 h =  ( y +  \lfloor \frac{y}{4}  \rfloor -  \lfloor \frac{y}{100}  \rfloor +  \lfloor \frac{y}{400}  \rfloor +  \lfloor \frac{13m+8}{5}  \rfloor +d )  mod 7  \cdots (2)
 h = 0 → Sun, 1 → Mon,  2 → Tue,  3 → Wed,  4 → Thu,  5 → Fri,  6 → Sat
ただし、求めたい日の月が1月、2月の場合はそれぞれ前年の13月、14月とする
( m = 1, 2 の場合は、 y \leftarrow  y - 1,  m = \leftarrow m+12とし、1年を、3月1日 〜 14月28日(閏年は29日)と仮定する)。

また、 y=100J+Kとしたとき、
 h = ( 5J + K + \lfloor \frac{K}{4} \rfloor + \lfloor \frac{J}{4} \rfloor + \lfloor \frac{26(m+1)}{10} \rfloor + d + 6 )  mod 7  \cdots (3)
とも表せる。
ただし、 (0\leq K \leq 99), K,Jは整数。

コード Zellerの公式
enum WeekDay {
	Sun("Sunday"), Mon("Monday"), Tue("Tuesday"), Wed("Wednesday"), Thu("Thursday"), Fri("Friday"), Sat("Saturday");
	private static final WeekDay[] Table = WeekDay.values();
	public static final WeekDay getWeekDay(int i){ return Table[i]; }
	public final String description;
	private WeekDay(String description){ this.description = description; }
	@Override public String toString(){ return description; }
}
/**
 * Calculates a day of the week from (y,m,d) via Zeller's congruence. <br>
 * @param y	year(>=-4800?) 1 BD-> y=0, 2 BD-> y=-1, ...,  n BD-> y=1-n
 * @param m	month
 * @param d	day
 * @return 	a day of the week
 */
public static final WeekDay Zeller(int y, int m, int d){
	y += 4800; //y<=0で使いたい場合
	if (m < 3) { --y; m += 12; }
	return WeekDay.getWeekDay((y + y/4 - y/100 + y/400 + (13*m + 8)/5 + d) % 7);
}
/**
 * Calculates a day of the week from (100J+K, m, d) via Zelle's congruence. ※y > 0<br>
 * @param J	year/100
 * @param K	year mod 100
 * @param m	month
 * @param d	day
 * @return 	a day of the week
 */
public static final WeekDay Zeller(int J, int K, int m, int d){
	if (m < 3) { 
		int y=100*J+K-1;
		J=y/100; K=y%100; m += 12;
	}
	return WeekDay.getWeekDay((5*J+K+K/4+J/4+26*(m+1)/10 + d +6) % 7);
}
証明 Zellerの公式

Fairfieldの公式より
 h
 =\left( 365y +  \lfloor \frac{y}{4}  \rfloor -  \lfloor \frac{y}{100}  \rfloor +  \lfloor \frac{y}{400}  \rfloor +  \lfloor \frac{306(m+1)}{10}  \rfloor +d - 428 \right) mod 7
 = \left( 7\cdot 52y + y  +  \lfloor \frac{y}{4}  \rfloor -  \lfloor \frac{y}{100}  \rfloor +  \lfloor \frac{y}{400}  \rfloor +  \lfloor \frac{306(m+1)}{10}  \rfloor +d - (61 \cdot 7+1) \right)mod 7
 = \left( y  +  \lfloor \frac{y}{4}  \rfloor -  \lfloor \frac{y}{100}  \rfloor +  \lfloor \frac{y}{400}  \rfloor +  \lfloor \frac{140(m+1)}{5} + \frac{13(m+1)}{5}  \rfloor +d - 1 \right)mod 7
 = \left( y  +  \lfloor \frac{y}{4}  \rfloor -  \lfloor \frac{y}{100}  \rfloor +  \lfloor \frac{y}{400}  \rfloor +  \lfloor \frac{7\cdot 20(m+1)}{5} + \frac{13m+13}{5}  \rfloor +d - 1 \right)mod 7
 = \left( y  +  \lfloor \frac{y}{4}  \rfloor -  \lfloor \frac{y}{100}  \rfloor +  \lfloor \frac{y}{400}  \rfloor +  \lfloor 7\cdot 4(m+1) + \frac{13m+13}{5}  \rfloor +d - 1 \right)mod 7
 = \left( y  +  \lfloor \frac{y}{4}  \rfloor -  \lfloor \frac{y}{100}  \rfloor +  \lfloor \frac{y}{400}  \rfloor +  \lfloor \frac{13m+13}{5}  \rfloor +d - 1 \right)mod 7
 kが整数なら \lfloor k + x \rfloor = k + \lfloor x \rfloorなので、
 h
 = \left( y  + \lfloor \frac{y}{4} \rfloor - \lfloor \frac{y}{100} \rfloor + \lfloor \frac{y}{400} \rfloor + \lfloor \frac{13m+13}{5} - \frac{5}{5} \rfloor +d \right) mod 7
 = \left( y  + \lfloor \frac{y}{4} \rfloor - \lfloor \frac{y}{100} \rfloor + \lfloor \frac{y}{400} \rfloor + \lfloor \frac{13m+8}{5} \rfloor +d \right) mod 7  \equiv 式(2)

ここで y=100J+Kを代入すると、
 h
 = \left( 100J+K  + \lfloor \frac{100J+K}{4} \rfloor - \lfloor \frac{100J+K}{100} \rfloor + \lfloor \frac{100J+K}{400} \rfloor + \lfloor \frac{13m+8}{5} \rfloor +d \right) mod 7
 = \left( 100J+K  + \lfloor 25J+\frac{K}{4} \rfloor - \lfloor J+\frac{K}{100} \rfloor + \lfloor \frac{J}{4} + \frac{K}{400} \rfloor + \lfloor \frac{26m+16+10}{10}-\frac{10}{10} \rfloor +d \right) mod 7
 = \left( 100J+K  + \lfloor 25J+\frac{K}{4} \rfloor - \lfloor J+\frac{K}{100} \rfloor + \lfloor \frac{J}{4} + \frac{K}{400} \rfloor + \lfloor \frac{26m+26}{10}-1 \rfloor +d \right) mod 7
 kが整数なら \lfloor k + x \rfloor = k + \lfloor x \rfloorなので、
 h
 = \left( 100J+25J+K  + \lfloor \frac{K}{4} \rfloor - J - \lfloor \frac{K}{100} \rfloor + \lfloor \frac{J}{4} + \frac{K}{400} \rfloor +  \lfloor \frac{26m+26}{10}\rfloor +d -1 \right) mod 7
 = \left( 124J+K  + \lfloor \frac{K}{4} \rfloor - \lfloor \frac{K}{100} \rfloor + \lfloor \frac{J}{4} + \frac{K}{400} \rfloor + \lfloor \frac{26(m+1)}{10}\rfloor +d -1 + 7\right) mod 7
 = \left( (7\cdot 17 + 5)J+K  + \lfloor \frac{K}{4} \rfloor - \lfloor \frac{K}{100} \rfloor + \lfloor \frac{J}{4} + \frac{K}{400} \rfloor + \lfloor \frac{26(m+1)}{10} \rfloor +d +6 \right) mod 7
 = \left( 5J+K  + \lfloor \frac{K}{4} \rfloor - \lfloor \frac{K}{100} \rfloor + \lfloor \frac{J}{4} + \frac{K}{400}\rfloor +\lfloor \frac{26(m+1)}{10}\rfloor +d +6 \right) mod 7
 0\leq K \leq 99, 0\leq \frac{K}{100}\leq 0.99<1, 0\leq \frac{K}{400}\leq 0.2475<1より、 \lfloor \frac{K}{100} \rfloor =0
また、 \frac{J}{4}の小数部は、0.025, 0.5, 0.75のいずれかをとるので、 \frac{J}{4}+\frac{K}{100}の小数部は高々0.75+0.2475 = 0.9975<1
従って、 \lfloor \frac{J}{4} + \frac{K}{100} \rfloor = \lfloor \frac{J}{4} \rfloor
これらより、
 h = \left( 5J+K  + \lfloor \frac{K}{4} \rfloor + \lfloor \frac{J}{4} \rfloor + \lfloor \frac{26(m+1)}{10} \rfloor +d + 6 \right) mod 7  \equiv 式(3)
(証明終わり)

ユリウス(通)日(Julian day)

ユリウス暦BC4713年1月1日の正午(世界時)を0日としたときの通算の日数。
またユリウス日から2400000.5を引いた数を修正(準)ユリウス日(modified Julian day)と呼ぶ。

ユリウス日⇔暦の求め方を探してみると、いろいろ見つかったけど、
負の数に対する除算、剰余演算結果によって計算結果が異なる恐れのなさそうな、
Julian day - Wikipedia, the free encyclopediaSpaghetti Source - 日付関連を参考にさせてもらった(というか殆ど写しただけです)。
この辺の話って、天文学の本とかに乗ってんのかなぁ...
※紀元前の時代は、グレゴリオ暦ではありませんが、グレゴリオ暦施行からさかのぼって、計数した値になってます。
※紀元前n年を1-n年と表現しています(紀元前(後)0年は存在しない)。

ユリウス日グレゴリオ暦

C++のBoost.GregorianとPHPのcal_from_jdでユリウス日 = 0 〜5373485日(10000/01/01)の範囲で確かめました。
※cal_from_jdでは紀元前n年を-nと表すが、下のプログラムは1-n年と表す。
Boost.Gregorianは1400/01/01〜10000/01/01の間で正常に動作するはずだが、
修正ユリウス日に関しては、値が負になる時、boost::gregorian::date::modjulian_day()(Boost. 1.47.0)は正しく求まらなかった。おそらくunsignedで計算しているため。ダメじゃん...。
※以降C++

/* transform hour, minute and sec into day */ 
inline double hms2day(int hour, int minute, int sec){
	return (hour + (minute + sec/60.0)/60.0)/24.0;
}
/*****************************************************
 * Gregorian Calendar -> (modified) Julian day
 * 1 BD-> y=0, 2 BD-> y=-1, ...,  n BD-> y=1-n
 * (Julian day + 1)%7 = 0->Sun, 1->Mon, .., 6->Sat
 * ((modified Julian day)%7 + 10 )% 7 = 0->Sun, 1->Mon, .., 6->Sat
 * ※負の数に対する除算や剰余の計算結果は、言語あるいは処理系によって異なるので注意
 *****************************************************/
// http://www.prefield.com/algorithm/misc/date.html
int greg2jd(int y, int m, int d, bool modflag = false) { //Ver. 1
    y += 4800;
	if (m < 3){ --y; m += 12; }
    return 365*y + y/4 - y/100 + y/400 + (153*m - 457)/5 + d-32045 - (modflag ?  2400001 : 0);
}
double greg2jd_double(int y, int m, int d, bool modflag = false, int hour = 0, int minute = 0, int sec = 0){
	return greg2jd(y, m, d, modflag) + hms2day(hour, minute, sec) - (modflag ? 0.0 : 0.5);
}
// http://en.wikipedia.org/wiki/Julian_day
int gregtojd(int y, int m, int d, bool modflag = false){ // Ver. 2
	int a  = (14-m)/12;
	y = y + 4800 - a;
	m = m + 12*a -3;
	return d + (153*m+2)/5 + 365 * y + y/4 - y/100 + y/400 - 32045 - (modflag ?  2400001 : 0);
}
double gregtojd_double(int y, int m, int d, bool modflag = false, int hour = 0, int minute = 0, int sec = 0){
	return gregtojd(y,m,d, modflag)  + hms2day(hour, minute, sec) - (modflag ? 0.0 : 0.5);
}

/*****************************************************
 * Julian day -> Gregorian Calendar
 * 1 BD-> y=0, 2 BD-> y=-1, ...,  n BD-> y=1-n
 *****************************************************/
// http://en.wikipedia.org/wiki/Julian_day
void jd2greg(int jd, int& y, int& m, int& d){
	int temp = jd + 32044; // +0.5
	y = -4800;
	
	jd = temp/146097;		temp %= 146097;			y += (jd * 400);
	jd = (temp/36524 + 1)*3/4;	temp -= (jd * 36524);	        y += (jd * 100);
	jd = temp/1461;			temp %= 1461;			y += (jd * 4);
	jd = (temp/365 + 1)*3/4;	        temp -= (jd * 365);		y += jd;

	m = (temp*5 + 308)/153 -2;
	d = temp - (m + 4)*153/5 + 122;
	y += ((m + 2)/12 );

	m = (m+2) % 12 + 1;
	d++;
}
ユリウス暦ユリウス日

こっちは確かめてません

/***************************************************
 * Julian Calendar-> (modified) Julian day 
 * 1 BD-> y=0, 2 BD-> y=-1, ...,  n BD-> y=1-n
 ***************************************************/
//http://www.prefield.com/algorithm/misc/date.html
long julian2jd(int y, int m, int d, bool modflag = false) {
    y += 4716;
    if (m < 3){ --y; m += 12; }
    return 365L*y + y/4 + (153*m - 457)/5 + d-1402 - (modflag ? 2400001 : 0);
}
double julian2jd_double(int y, int m, int d, bool modflag = false, int hour = 0, int minute = 0, int sec = 0){
	return julian2jd(y, m, d, modflag) + hms2day(hour, minute, sec);
}