儒略日-题解

前言

时隔两年，这个极为经典的题目终于被我 AC 了。经过诸多优化改良，最终得到了这个个人认为比较优美的做法，写篇题解纪念一下，也供参考。

首先，建议读者先对照无注释代码自行理解一下大致过程。

无注释代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mon[13]={0,31,28,31,30,31,30,31,31,30,31,30,31};
const int _1=365,_100=_1*100+24,_400=_100*4+1;
struct date{int y,m,d,ru;} t[4000010];
bool isrun(int year) {
	if(year<0) return (year+1)%4==0;
	return year%4==0&&(year<=1582||year%100!=0||year%400==0);
}
date nxtday(date a) {
	a.d++,a.ru++;
	if(a.y==1582&&a.m==10&&a.d==5) return a.d=15,a;
	if(a.m==2&&a.d==29&&isrun(a.y));
	else if(a.d>mon[a.m]) a.m++,a.d=1;
	if(a.m>12) a.y++,a.y+=!a.y,a.m=1;
	return a;
}
void print(date x) {
	cout<<x.d<<" "<<x.m<<" "<<abs(x.y)<<(x.y<0?" BC":"")<<endl;
}
signed main() {
	t[0]={-4713,1,1,0};
	int _20000101;
	for(int i=1;t[i-1].y<=2000;i++) {
		t[i]=nxtday(t[i-1]);
		if(t[i].y==2000&&t[i].m==1&&t[i].d==1)
			_20000101=t[i].ru;
	}
	int _0=_20000101-_400*5,T,x;
	cin>>T;
	while(T--) {
		cin>>x;
		if(x<=_20000101) print(t[x]);
		else {
			date ans{(x-_0)/_400*400,1,1,(x-_0)/_400*_400+_0};
			while(ans.ru+_100+isrun(ans.y)<=x)
				ans.ru+=_100+isrun(ans.y),ans.y+=100;
			while(ans.ru+_1+isrun(ans.y)<=x)
				ans.ru+=_1+isrun(ans.y),ans.y++;
			while(ans.ru<x) ans=nxtday(ans);
			print(ans);
		}
	}
	return 0;
}

解析

大体思路是这样的：\(1582\) 年之前的暴力 nxtday 跑出来，而 \(1582\) 年之后的想办法逐渐逼近计算。

这里我们暴力跑到了 \(2000\) 年（以下讲解只考虑 \(2000\) 年以后），并保存了 _20000101 这个常数，表示 \(2000\) 年 \(1\) 月 \(1\) 日的儒略日。原因之一是避免考虑边界情况，而原因之二，就是本做法的灵魂所在：

重点在于 _20000101 这个常数。将这个常数减去五个 \(400\) 年（代码中的 _400），就可以反推出“理想状态下”公元 \(0\) 年 \(1\) 月 \(1\) 日的儒略日 _0。这样，我们就可以 \(O(1)\) 计算出所求儒略日位于哪一个“\(400\) 年周期”中，并得到该周期的第一天的确切 date：

• year=(x-_0)/_400*400
• ru=(x-_0)/_400*_400+_0

这两行为本 trick 的精华，建议仔细理解（见程序第 36 行）

求出了本周期第一天的儒略日，接下来就是简单的逼近答案了：先尝试每次加上 \(100\) 年（不超过 \(4\) 次），再每次加上 \(1\) 年（不超过 \(100\) 次），最后暴力 nxtday 跑到答案（不超过 \(365\) 次）。时间复杂度在可接受范围之内。

详细代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mon[13]={0,31,28,31,30,31,30,31,31,30,31,30,31};
const int _1=365,_100=_1*100+24,_400=_100*4+1;
//方便计算，预处理出格里高利历1年、100年和400年的天数
struct date{int y,m,d,ru;} t[4000010];
//公元前的y为负数
bool isrun(int year) { //闰年判断
	if(year<0) return (year+1)%4==0;
	return year%4==0&&(year<=1582||year%100!=0||year%400==0);
}
date nxtday(date a) { //暴力计算下一天
	a.d++,a.ru++; //天数和儒略日首先增加
	if(a.y==1582&&a.m==10&&a.d==5) return a.d=15,a;
	//特判中间消失的10天
	if(a.m==2&&a.d==29&&isrun(a.y)); //闰年的2月29无需进位
	else if(a.d>mon[a.m]) a.m++,a.d=1; //由日向月进位
	if(a.m>12) a.y++,a.y+=!a.y,a.m=1; //由月向年进位（特判了公元0年不存在）
	return a;
}
void print(date x) { //输出
	cout<<x.d<<" "<<x.m<<" "<<abs(x.y)<<(x.y<0?" BC":"")<<endl;
}
signed main() {
	t[0]={-4713,1,1,0};
	int _20000101;
	for(int i=1;t[i-1].y<=2000;i++) { //暴力打表到2000年
		t[i]=nxtday(t[i-1]);
		if(t[i].y==2000&&t[i].m==1&&t[i].d==1)
			_20000101=t[i].ru; //保存常数2000.1.1儒略日
	}
	int _0=_20000101-_400*5,T,x; //计算“理想状态”下0.1.1儒略日
	cin>>T;
	while(T--) {
		cin>>x;
		if(x<=_20000101) print(t[x]); //直接查表
		else {
			date ans{(x-_0)/_400*400,1,1,(x-_0)/_400*_400+_0};
			//计算400周期的第一天（格式为y,m,d,ru）
			while(ans.ru+_100+isrun(ans.y)<=x)
				ans.ru+=_100+isrun(ans.y),ans.y+=100;
			//100年为步长逼近
			while(ans.ru+_1+isrun(ans.y)<=x)
				ans.ru+=_1+isrun(ans.y),ans.y++;
			//1年为步长逼近
			while(ans.ru<x) ans=nxtday(ans);
			//暴力nxtday得到答案
			print(ans);
		}
	}
	return 0;
}