neorobin

google:
Date calculation algorithm
有些有价值的参考

neorobin

set /a "m+=9, m%%=12, y-=m/10, index=365*y + y/4 - y/100 + y/400 + (m*306 + 5)/10 + d - 1
复制代码

计算年: 从历法年 x.3.1 日到 (x+1).2.(28/29) 日, 定义为一个计算年,
一个计算年的长度不固定:
365: 当 x+1 是平年份
366: 当 x+1 是闰年份

一个计算年中, 3.1 日称为年始日, 2.(28/29) 日称为年终日
以 0 年 3 月 1 日为基准日期, 序号设为 0, 其后所有日期的序号依次加 1.
基准日期就是第 1 个计算年的年始日.
前 400 个计算年的始日, 终日, 长度, 始日序号, 终日序号如下:

p	first	last	len	i(first)	i(last)
1	0.3.1	1.2.28	365	0	364
2	1.3.1	2.2.28	365	365	729
3	2.3.1	3.2.28	365	730	1094
4	3.3.1	4.2.29	366	1095	1460
...
8	7.3.1	8.2.29	366	2556	2921
...
99	98.3.1	99.2.28	365	35794	36158
100	99.3.1	100.2.28	365	36159	36523
...
400	399.3.1	400.2.29	366	145731	146096

从基准日期开始, p 个计算年的总天数公式为(全作平年算, 加上 4 的倍数的闰年的个数, 减去整百年份的个数, 再加上 400 的倍数闰年份的个数):

365*p + p/4 - p/100 + p/400

其中包含的最后的年终日和年始日的序号为:

最后一个年终日的序号 = 365*p + p/4 - p/100 + p/400 - 1

最后一个年始日的序号 = 365*q + q/4 - q/100 + q/400      (其中, q=p-1)

对于历法日期 y.m.d, 首先假定它不是一个年终日,
从基准日期到 小于这个日期的最后一个年终日, 就得到 t 个计算年:

t = y:  当 m >= 3 时; 或者 t = y-1:  当 m < 3 时.

计算 t 的代码:

set /a n=m, n+=9, n%=12, t=y, t-=n/10
复制代码

第 t 个年终日的序号 = 365*t + t/4 - t/100 + t/400 - 1
第 t+1 个年始日的序号 = 365*t + t/4 - t/100 + t/400

上述计算中的

set /a n=m, n+=9, n%=12
复制代码

将 m 线性映射到一个新的变量 n, 映射关系为:
m -> n : {3, 4, 5, ... 12, 1, 2} -> {0, 1, 2, ... 9, 10, 11}

n 不再是历法月份的意义, 而是一个计算年中每一个月份对起始月份 3 月的偏移数.
这个 n 将有助于我们构造一个结构简单的离散函数来计算每个月的起始日对 3.1 日的偏移量.

于是:

y.m.d 的序号 = 第 t+1 个年始日(3.1 日)的序号 + n.1 日对 3.1 日的偏移 + d 日对 1 日的偏移

m.1 日对 3.1 日的偏移 用一个关于 n 的近似线性的离散函数来实现:

int(f(n)) = int((n*306 + 5)/10)

这里 int() 是截尾式取整, 即把参数的小数部分直接去掉, 而不是四舍五入的方式.

Month	n	f(n)*10	int(f(n))	int(f(n+1))-int(f(n))
Mar	0	5	0	31
Apr	1	311	31	30
May	2	617	61	31
Jun	3	923	92	30
Jul.	4	1229	122	31
Aug	5	1535	153	31
Sep	6	1841	184	30
Oct	7	2147	214	31
Nov	8	2453	245	30
Dec	9	2759	275	31
Jan	10	3065	306	31
Feb	11	3371	337

第 n 个计算月的天数: Δint(f(n)) = int(f(n+1))-int(f(n))   (n < 11)
Δint(f(n)) 的平均值 = 30.6 ≈ (31+30+31+30+31+31+30+31+30+31+31)/11

上面表中 f(n)*10, int(f(n)), int(f(n+1))-int(f(n)) 这 3 列可用如下代码在命令行运行得到

for /l %a in (0 1 11) do >nul set /a "fn=%a*306+5" & set fn
for /l %a in (0 1 11) do >nul set /a "fn=%a*306+5, ifn=fn/10" & set ifn
for /l %a in (0 1 10) do >nul set /a "delta=((%a+1)*306+5)/10-(%a*306+5)/10" & set delta
复制代码

以上, 对 年始日序号 和 月始日序号偏移 的计算, 为了方便而引入了变量 t 和 n, 现在我们去掉这两个变量,
给出计算日期 y.m.d 的序号 i 的伪代码:

i(y, m, d) {
  set /a "m+=9, m%=12, y-=int(m/10)"
  return 365*y + y/4 - y/100 + y/400 + int((m*306 + 5)/10) + d - 1
}

仅就此函数而言, 返回值也可不要最后的 "- 1", 只是把任意日期的序号增长了 1, 或者认为将日期序号轴的原点左移了一天.

cmd 的代码在本文开头已给出.

下接 25# 序号求日期的算法分析

参考原文:
http://alcor.concordia.ca/~gpkatch/gdate-algorithm.html
http://alcor.concordia.ca/~gpkatch/gdate-method.html

neorobin

回复 17# plp626

在 cmd 下不能用 set/a "t=a/b+!!(a%b)" 实现 ceil(a/b)

在 cmd 下, 当 a, b 正负同号时, a/b 采取的是截尾式取整 (floor() 的效果); 但当 a, b 正负异号时, 其实和 ceil(a/b) 是等同的.
或者换一种理解方式: 将 商的绝对值按 floor 的方式取整, 再附上符号, 同取正, 异取负.

set /a 16/7
2
set /a -16/-7
2
set /a -16/7
-2
set /a 16/-7
-2

对于 floor((m*306 + 5)/10) 中的 5, 并非唯一, 只要几个常数配合上能让最后的结果得到 0,31,61,92,122,153,184,214,245,275,306,337 这样一个序列就是理想的组合, 或者说对这个序列给出了一个 完美拟合 并且 形式上足够简洁 而且易于实现 的 解析式.

m*306 + 5 得到的序列:
5,311,617,923,1229,1535,1841,2147,2453,2759,3065,3371 中 个位数最小为 1, 最大达到 9, 所以 5  换成 4, 仍然可以完美拟合, 但能和 306, 10 理想配合的除了 5 和 4 外再没有别的数字了.

neorobin

上接 15# 日期序号算法分析

序号求日期的算法分析

设算法函数为 date(i), i 为待求日期的序号, 返回值为 y.m.d, 其中包含 p 个计算年, d1st 为第 p + 1 个计算年的年始日序号.
首先对 p 作估值, 以

(d1st-1) * 250 / 91311

计算初始估值, 在 32位 cmd 下估值溢出下限为 23519, 故而 p ∈ [1, 20000] 估值不会发生溢出错误.

在第 p + 1 个计算年中:
估值误差范围: [(d1st-1) * 250 / 91311 - p, ((d1st-1) + 365) * 250 / 91311 - p]
在 p ∈ [1, 20000] 内计算得出此估值误差的最小值, 最大值为:

[-0.088522, 0.998171] (-0.088522 : 当 p = 19903;  0.998171 : 当 p = 96)

以 trunc 取整再加 1 后, 估值最小值, 最大值为:
  [trunc(p-0.088522)+1, trunc(p+0.998171)+1] = [p, p+1]

估值后, 计算序号误差, 若超过原序号, 表明估值计算得 p+1, 将估值减 1,
再重新计算序号误差:

set /a "y=(i-1)*250/91311+1, dd=i-(365*y + y/4 - y/100 + y/400)"
set /a "y+=dd>>31, dd=i-(365*y + y/4 - y/100 + y/400)"
复制代码

上面的计算中, 最后得到的 dd 即为第 p+1 个计算年中 m.d 日 对 3.1 日 的偏移值, 从这个偏移值计算 m, 采用如下近似线性的拟合函数:

ind(m) = trunc(f(dd)) = trunc((5*dd + 2)/153)

下表呈现了该函数是完美拟合的, s, e 为每一月的始日和终日对 3.1 日的偏移

month	len	ind	s	e	f(s)*1000	f(e)*1000	trunc(f(s))	trunc(f(e))
Mar	31	0	0	30	13	993	0	0
Apr	30	1	31	60	1026	1973	1	1
May	31	2	61	91	2006	2986	2	2
Jun	30	3	92	121	3019	3967	3	3
Jul	31	4	122	152	4000	4980	4	4
Aug	31	5	153	183	5013	5993	5	5
Sep	30	6	184	213	6026	6973	6	6
Oct	31	7	214	244	7006	7986	7	7
Nov	30	8	245	274	8019	8967	8	8
Dec	31	9	275	305	9000	9980	9	9
Jan	31	10	306	336	10013	10993	10	10
Feb	28	11	337	364	11026	11908	11	11
Feb(leap)	29	11	337	365	11026	11941	11	11

上面表格可用代码得到:

@echo off & setlocal enabledelayedexpansion
echo mon len ind    s     e     f(s)    f(e)  int(f(s))   int(f(e))
for %%i in ("Feb 28" "Feb 29") do (
  set /a "i=100, s=1000, e=1000-1"
  for %%a in ("Mar 31" "Apr 30" "May 31" "Jun 30" "Jul 31" "Aug 31"
  "Sep 30" "Oct 31" "Nov 30" "Dec 31" "Jan 31" %%i) do (
    set "l=%%~a"
    set /a "e+=!l:~-2!, fs=(5*(s-1000)+2)*1000/153+100000,fe=(5*(e-1000)+2)*1000/153+100000"
    set /a "ifs=(5*(s-1000)+2)/153+100,ife=(5*(e-1000)+2)/153+100"
    set "lin=!l!   !i:~-2!   !s:~-3!   !e:~-3!   !fs:~-5!   !fe:~-5!     !ifs:~-2!          !ife:~-2!"
    set "lin=!lin: 0=  !"
    if %%i=="Feb 28" (echo !lin!) else if %%a=="Feb 29" echo !lin!
    set /a "s+=!l:~-2!, i+=1"
))
复制代码

综上, 由序号求日期的完全代码:

set /a "y=(i-1)*250/91311+1,dd=i-(365*y+y/4-y/100+y/400)"
set /a "y+=dd>>31,dd=i-(365*y+y/4-y/100+y/400)"
set /a "mi=(5*dd+2)/153"
set /a "m=(mi+2)%%12+1" & rem {0,1,2,..9,10,11} -> {3,4,5,..12,1,2}
set /a "y-=m-3>>31" & rem 1,2月到了下一年
set /a "d=1+dd-(153*mi+2)/5"
复制代码

PS:
计算序号的代码在开始可以将 y 加上一个 400 的倍数, 在序号转为日期的最后将 y 再减去同一个 400 的倍数, 这样就可以将年份 y 的可计算范围向负数调整. 上面 [1, 20000] 的范围分作正负数各一半仍足够大.

和 plp626, terse不约而同的都想到了 153, 2, 5 这个组合.

下面是一个算法互逆性测试代码, 但不能对序号生成的日期正确性作验证测试:

@echo off & setlocal enabledelayedexpansion
REM 7305155=20000.12.31
for /l %%i in (0 1 7305155) do (
  set /a "i=%%i"
  set /a "y=(i-1)*250/91311+1,dd=i-(365*y+y/4-y/100+y/400)"
  set /a "y+=dd>>31,dd=i-(365*y+y/4-y/100+y/400)"
  set /a "mi=(5*dd+2)/153"
  set /a "m=(mi+2)%%12+1" & rem {0,1,2,..9,10,11} -> {3,4,5,..12,1,2}
  set /a "y-=m-3>>31" & rem 1,2月到了下一年
  set /a "d=1+dd-(153*mi+2)/5"
  set "test=%%i:   !y!.!m!.!d!"

  set /a "m+=9,m%%=12,y-=m/10,i=365*y+y/4-y/100+y/400+(m*306+5)/10+d-1"
  set "test=!test!     !i!"
  if "!i!" neq "%%i" set "test=!test!     error"
  echo !test!
  if "!test:error=!" neq "!test!" pause
)
复制代码

neorobin

回复 26# plp626


    http://alcor.concordia.ca/~gpkatch/gdate-method.html

y = (10000*g + 14780)/3652425

10000 / 3652425  (相当于 除以 400 年的年平均天数 365.2425)
除以 40
为了 20000] 的界限不溢出
23519 * 91311 = 2147543409 上溢了

neorobin

回复 28# plp626


    常数再小的话, 估值误差的范围会增大, 后面的代码不好简洁实现了

neorobin

回复 28# plp626


    [1,20000]

pOffs := ((d1st - 1) * 250) / 91311 - p;

The pOffsMin is: -0.088522;  the pOffsMax is: -0.001161
The pMin is: 19903;  the pMax is: 96
pOffsMin = -0.088522,  pOffsMax + 365*250 / 91311 = 0.998171

pOffs := ((d1st - 1) * 25) / 9131 - p;

The pOffsMin is: -0.004709;  the pOffsMax is: 0.134158
The pMin is: 303;  the pMax is: 20000
pOffsMin = -0.004709,  pOffsMax + 365*25 / 9131 = 1.133501

这个误差范围, 还可以调整到 [p-1, p] 了
下面更甚

pOffs := ((d1st - 1) * 5) / 1825 - p;

The pOffsMin is: -0.002740;  the pOffsMax is: 13.284932
The pMin is: 1;  the pMax is: 20000
pOffsMin = -0.002740,  pOffsMax + 365*5 / 1825 = 14.284932

neorobin

回复 28# plp626


    刚看到你回复的 30 楼, 突然不见了
测试结果
pOffs := (d1st * 10000.0) / 3652425 - p;

The pOffsMin is: -0.004045;  the pOffsMax is: 0.001971
The pMin is: 8303;  the pMax is: 8496
pOffsMin = -0.004045,  pOffsMax + 365*10000 / 3652425 = 1.001307

-0.004045 * 3652425 = -14774.059125
14780 比较恰好的调整误差范围到 [0.000000, 1.005354]:

pOffs := (d1st * 10000.0 + 14780) / 3652425 - p;

The pOffsMin is: 0.000000;  the pOffsMax is: 0.006018
The pMin is: 4;  the pMax is: 1296
pOffsMin = 0.000000,  pOffsMax + 365*10000 / 3652425 = 1.005354


pOffs := (d1st * 10000.0 -7300) / 3652425 - p;

The pOffsMin is: -0.006044;  the pOffsMax is: -0.000027
The pMin is: 8303;  the pMax is: 2096
pOffsMin = -0.006044,  pOffsMax + 365*10000 / 3652425 = 0.999309

Pascal 测试代码

var
  d1st, p, pMin, pMax: longint;
  pOffs, pOffsMin, pOffsMax: real;
begin
  pOffsMin := 0;
  pOffsMax := -1;
  for p := 1 to 20000 do
  begin
    d1st := 365 * p + p div 4 - p div 100 + p div 400;
    pOffs := (d1st * 10000.0) / 3652425 - p;
    if pOffs < pOffsMin then begin pOffsMin := pOffs; pMin := p; end
    else if pOffs > pOffsMax then begin pOffsMax := pOffs; pMax := p; end;
    WriteLn(p, '  ', d1st, '  ', pOffs: 0: 6);
  end;
  WriteLn('The pOffsMin is: ', pOffsMin: 0: 6, ';  the pOffsMax is: ', pOffsMax: 0: 6);
  WriteLn('The pMin is: ', pMin: 0, ';  the pMax is: ', pMax: 0);
  WriteLn('pOffsMin = ', pOffsMin: 0: 6, ',  pOffsMax + 365*10000 / 3652425 = ', pOffsMax + 365 * 10000 / 3652425: 0: 6);
  Readln;
end.
复制代码

neorobin

回复 31# neorobin
   
http://alcor.concordia.ca/~gpkatch/gdate-method.html

. In the intervening years of this 400 year cycle, the error maxima are 1.4775 days in year 303 of the cycle, and -0.72 days in year 96 of the cycle. Because of this error, finding the year given a number of days (the inverse function) is not exact, but we can find a very close approximation. Given d days, the year number can be approximated as
y = d / 365.2425     (3)

查看 400 年周期内, 平均数估算天数的误差峰值, 在命令行运行:

>nul (set /a "mi=0,ma=0,ii=1,ia=1"&(for /l %y in (1 1 400) do set /a "y=%y,d=y*365+y/4-y/100+y/400,d*=10000,d=y*3652425-d,ii^=(d-mi>>31)&(y^ii),ia^=(ma-d>>31)&(y^ia),mi^=(d-mi>>31)&(d^mi),ma^=(ma-d>>31)&(d^ma)"))&set m&set i
复制代码

ma=14775
mi=-7200
ia=303
ii=96

neorobin

回复 33# plp626

把 25 楼的代码改了一点点, 没有实质性的变化, 凑成 2 行 155 字节:

set /a "y=(i*99+145)/36159,y+=i-365*y-y/4+y/100-y/400>>31,dd=i-365*y-y/4+y/100-y/400"
set /a "mi=(5*dd+2)/153,m=(mi+2)%%12+1,y-=m-3>>31,d=1+dd-(153*mi+2)/5"
复制代码

neorobin

回复 35# terse


    不错, 我一直测试到 7305155 没有发现错误

neorobin

回复 34# Batcher

事实上我没有上过电脑课, Pascal: 端庄典雅秀丽的淑女; C: 自由奔放不羁的野马

neorobin

回复 39# plp626


    弄在一行了, 不必要的变量都省了, 137 字节

set /a "y=(i*99+145)/36159,y+=i-365*y-y/4+y/100-y/400>>31,d=i-365*y-y/4+y/100-y/400,m=(5*d+2)/153,d+=1-(153*m+2)/5,y+=m/10,m=(m+2)%%12+1"
复制代码

neorobin

回复 41# plp626

@echo off & setlocal enabledelayedexpansion
REM nextIndex := date2index(1.1.1);
for /l %%y in (1 1 20000) do (
  set "y=%%y"
  set /a "leap = (^!(%%y %% 4) & ^!^!(%%y %% 100)) | ^!(%%y %% 400)"
  if "!leap!"=="1" (set "eFeb=29") else set "eFeb=28"
  for %%t in ("1 31" "2 !eFeb!" "3 31" "4 30" "5 31" "6 30" "7 31" "8 31" "9 30" "10 31" "11 30" "12 31") do (
    for /f "tokens=1,2" %%m in ("%%~t") do (
      set "m=%%m"
      for /l %%d in (1 1 %%n) do (
        set "d=%%d"

        REM testing...
        REM if date2index(y.m.d) <> nextIndex (
          REM error
        REM )
        REM if index2date(date2index(y.m.d)) <> y.m.d (
          REM error
        REM )
        REM nextIndex += 1;
        echo !y!.!m!.!d!
      )
    )
    pause
  )
)
复制代码

和 13 楼说的差不多, 我没想到可以证明正确的方式:

设想一个并不聪明的一定范围内穷举式的测试方式:

公历的历法是明确的:
每年 12 个月: m∈[1..12]
4,6,9,11 月份固定 30 天, d ∈[1..30]
2 月根据年份平闰取 28 天或 29 天,  d∈[1..29]
其余月份固定 31 天, d∈[1..31]

闰年:
(被 4 整除 且 不被 100 整除) 或 被 400 整除 的年份
闰年外的年份都是平年

从这个公历历法的原型定义出发设计一个 尽可能直接符合 且 足够简单 的日期生成算法(假定它是正确的算法):
从一个初始日期开始, 依次生成每一个明天,
把这个算法生成的日期作为 date2index 的参数, 得到相应的 序号:

i_0: 初始日期的序号(值由算法自定);   (_0 表示 0 为下标)

i_p: 初始日期 或者 初始日期之后的 某个日期的序号;
i_q: i_p 对应的那个日期的下一日的序号;

那么, 如果对于任意 i_p 与 i_q, 都得到:

i_q = i_p + 1

就可以把 date2index 看作是测试范围内正确的算法.

在上面的测试中可以同时 作 index2date 测试:
如果始终都有:
index2date( date2index(y.m.d) ) = y.m.d

在 date2index 测试正确的情况下, 就可以认为 index2date 也是正确的.


date2index 和 index2date 在任意参数的情况下的输出都是唯一的, 且是确定的(对于一个确定的输入 I, 不会得到 O1, O2两种或更多可能的输出),

因为算法是基于 数字计算机 的, 输入参数 不具备不确定性(每次输入的值的集合都是确定的),
且算法中不涉及随机函数的调用, 所以 输出是 具有 唯一性和确定性的.

neorobin

回复 44# terse

惭愧, 对儒略日, 历法变化等问题完全没有了解或者了解甚少, 前面都完全是按一个固定的历法规则处理的, 如要适应相关的 规则约定 或者 实际情况, 肯定要作出进一步的修改.