CrLf

回复 15# neorobin


    叹为观止，m*306 真是绝妙的一招

plp626

回复 15# neorobin

算法思想体现出了数学之美，深刻，很欣赏；之前用（365.2422±x）常数*year 试着矫正year2day， 对这类整数域上的曲线拟合缺少资料。。。该问题搁置。。。

看了下matlab里面的datenummx.c源代码，year2day用了一个函数，ceil(x)，向右取整；
ceil函数在cmd中实现相对繁琐；t=ceil(a/b)  同 set/a "t=a/b+!!(a%b)"；

static double cdm[] = {0,31,59,90,120,151,181,212,243,273,304,334};
        if (mon < 1) {
            mon = 1;
        }
        if (mon > 12) {
            y += (mon-1)/12;
            mon = ((mon-1) % 12) + 1;
        }
        *t = 365.*y + ceil(y/4.) - ceil(y/100.) + ceil(y/400.) + cdm[mon-1] + *d;
        if (mon > 2) {
            iy = (int) y;
            if ((iy%4 == 0) && (iy%100 != 0) || (iy%400 == 0)) {
                *t += 1.;
            }
        }
复制代码

对于月份的转换，没有像儒略日那样，把1,2月份当做上年的末月；而是传统思路，把其当做平年处理，再判断，当所在年份为闰年，且月份大于2时，再+1；

neorobin给出的资料把1，2月份当做上年月份处理，可以将ceil函数转换为fix函数
（fix(a/b)=set/a "a/b"，向零取整；15楼的int()取整函数就是向零取整）
也使得月份转换为天数之间建立了向整数近似的线性关系（如果把1，2月份当做当年月份处理，映射关系不是单调的）
从而大大简化了代码；

比较了下，matlab里的datenummx.c的源代码兼容公元前，月份为任意整数，天数为任意整数的情形，但思想在cmd中实现挺繁琐；

plp626

大家 讨论下 儒略日 的思想，为何 “参考日期（第零天）是 1858 年 11 月 17 日”

plp626

http://alcor.concordia.ca/~gpkatch/gdate-algorithm.html
月份转当年天数：(306*m+5)/10是如何构造的？常数5很像是凑出来的，换成4也可以；

terse

1858 年 11 月 17 日 应该是简化后的吧
儒略日应该推前-4712年1月1日 是儒略日的第0日
这里的 306 是否可看作是30.5833333的近似值  也就是 30.5833333 * M  
这里 (306*m+5)/10 可看作 30.6*M+0.5

这里的 0.5 应该是一个修正 就是+0.5天 因为儒略日简化前 是以中午12点起始计算的

现在计算儒略日公式  譬如（SET /A "M=(M+9)%%12+3,Y=Y-M/13,JDB=365*(4712+Y)+(4712+Y)/4+153*(M+1)/5-Y/100+Y/400+D-61"）

只是计算 1582年10月15后的 之前的要+10天  也就是 1582年10月4日至1582年10月14日是空白的

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

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字节，致精致简；

fatcat

回复 22# plp626


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

terse

发出后 看到P兄已发

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

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

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

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

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

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