terse

没完全测试 翻版一下网上资料及万年历

@echo off&setlocal enabledelayedexpansion
FOR /l %%i in (0 1 7305155) DO (
    set JD=%%i
    IF !JD! GEQ 2299161 set /a "JD+=1+(JD*100-186721625)/3652425-(JD*100-186721625)/3652425/4"
    set /a "B=JD+1524,Y=(B*100-12210)/36525,D=36525*Y/100"
    set /a "M=(B-D)*10000/306001,D=B-D-306001*M/10000,M=(M-2)%%12+1,Y-=4715+^!(2/M)"
    set "test=%%i:   !y!.!m!.!d!"
    set /a "M=(M+9)%%12+3,Y-=M/13,JD=36525*(Y+4716)/100+306001*(M+1)/10000+D-1524"
    set /a "JD+=(2-Y/100+Y/400)*^!(2299161/JD)"
    set "test=!test!     !JD!"
    if "!JD!" neq "%%i" set "test=!test!     error"
    echo !test!
    if "!test:error=!" neq "!test!" pause
)
PASUE
复制代码

Batcher

回复 31# neorobin


    现在还有学校在教Pascal啊，好神奇。

plp626

回复 30# neorobin
理论上的分析没经得起实践测试；漏掉了一个变量，分析有误，便删了；
---------------------------
date2index的反函数index2date；
应该可以精简优化到两行150个字节以内
很期待这个函数。。。

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

回复 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

回复 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


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

plp626

回复 27# neorobin


    这样啊，那你的代码可以再化简了

->m 1/365.2425
0.00273790701

->m 250/91311
0.00273789576
　
->m 250/91310
0.00273792575
复制代码

可以看出，用９１３１０代替９１３１１更接近真实值；这样一来２５０，９１３１０又可以化简了；同除以１０；得２５，９１３１；
－－－－－－－－－－
上面的接近分析有误，实际是：
｜0.00273790701－0.00273789576｜　＜　｜0.00273790701－0.00273792575｜

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 上溢了

plp626

回复 25# neorobin

struct sdate dtf(long d) { /* convert day number to y,m,d format */
  struct sdate pd;
  long y, ddd, mm, dd, mi;

  y = (10000*d + 14780)/3652425;
  ddd = d - (y*365 + y/4 - y/100 + y/400);
  if (ddd < 0) {
    y--;
    ddd = d - (y*365 + y/4 - y/100 + y/400);
    }
  mi = (52 + 100*ddd)/3060;
  pd.y = y + (mi + 2)/12;
  pd.m = (mi + 2)%12 + 1;
  pd.d = ddd - (mi*306 + 5)/10 + 1;
  return pd;
  }
复制代码

原作者的Ｃ代码；我精简后，除了开头部分，后面的和你基本一致；

你代码中的开始一行中的常数２５０，９１３１１，这些如何得来的？

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
)
复制代码

terse

发出后 看到P兄已发

(153*M+2)/5
复制代码

fatcat

回复 22# plp626


    几种取整方式: floor(), ceil(), trunc()(或表示为 fix()), round() cmd 中最易于实现的是 trunc(): 直接截去小数部分取整.

plp626

回复 21# neorobin

是我不够严密；应该是 当a/b>0时； set/a a/b+!!(a%b) 等同于 ceil(a/b)
【set/a a/b+!!(a%b) 可进一步化简为 set/a a/b-!(a%b)；】

概括来讲，set/a "a/b" 等价于fix(a/b)；而ab同号时，等价于floor(a/b)；ab异号时等价于ceil(a/b)；
我记忆fix,floor,ceil的方式是 ---- floor  --- (fix) --- ceil --->
即floor向左取整，fix向零(原点)取整，ceil向右取整；

对于x为正的情况，用fix如何表示ceil，需要根据具体的x来变换，画函数图像可看出 ceil(x)=N+fix(x-N) (-N<x<N，N为正整数) ；

比如，当a/b所表示的实数的绝对值小于10000时； ceil(a/b)=set/a (a-10000*b)/b+10000
------------------------------------
对于那个常数；
昨天思考了其一般问题即：如何判断两组数据在整数范围上是否线性相关，若存在如何求解关系式；已有初步结果；
（更一般的问题是 寻求在整数范围上有理函数建立的关系；这个数学问题很有意义）

这些整数常数（y=(kx+b)/c ，k,b,c为整数常数） 可以通过给出的数据直接求出k/c, 和修正数b的范围；
在给定c或者k的情况下，另外的两个常数便可确定范围；需要一些代数知识（像是整数范围上的最小二乘法）；
昨天贴了几张图都是和数学有关，有点离题，便又删了；

其实 (306*m+5)/10 能换成换成(306*m+4)/10，就可以化简了， (153*m+2)/5，这样代码少了1字节，致精致简；

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 外再没有别的数字了.