儒略日-题解

前言

时隔两年，这个极为经典的题目终于被我 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;
}