[新手上路]批处理新手入门导读[视频教程]批处理基础视频教程[视频教程]VBS基础视频教程[批处理精品]批处理版照片整理器
[批处理精品]纯批处理备份&还原驱动[批处理精品]CMD命令50条不能说的秘密[在线下载]第三方命令行工具[在线帮助]VBScript / JScript 在线参考
返回列表 发帖
回复 15# neorobin


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

TOP

本帖最后由 plp626 于 2012-4-7 19:54 编辑

回复 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)";
  1. static double cdm[] = {0,31,59,90,120,151,181,212,243,273,304,334};
  2.         if (mon < 1) {
  3.             mon = 1;
  4.         }
  5.         if (mon > 12) {
  6.             y += (mon-1)/12;
  7.             mon = ((mon-1) % 12) + 1;
  8.         }
  9.         *t = 365.*y + ceil(y/4.) - ceil(y/100.) + ceil(y/400.) + cdm[mon-1] + *d;
  10.         if (mon > 2) {
  11.             iy = (int) y;
  12.             if ((iy%4 == 0) && (iy%100 != 0) || (iy%400 == 0)) {
  13.                 *t += 1.;
  14.             }
  15.         }
复制代码
对于月份的转换,没有像儒略日那样,把1,2月份当做上年的末月;而是传统思路,把其当做平年处理,再判断,当所在年份为闰年,且月份大于2时,再+1;

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

比较了下,matlab里的datenummx.c的源代码兼容公元前,月份为任意整数,天数为任意整数的情形,但思想在cmd中实现挺繁琐;
附件: 您需要登录才可以下载或查看附件。没有帐号?注册

TOP

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

TOP

本帖最后由 plp626 于 2012-4-7 19:03 编辑

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

TOP

本帖最后由 terse 于 2012-4-7 23:55 编辑

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日是空白的

TOP

本帖最后由 neorobin 于 2012-4-9 01:55 编辑

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

TOP

本帖最后由 plp626 于 2012-4-9 12:13 编辑

回复 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字节,致精致简;

TOP

回复 22# plp626


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

TOP

发出后 看到P兄已发
  1. (153*M+2)/5
复制代码

TOP

本帖最后由 neorobin 于 2012-4-9 15:44 编辑

上接 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,
再重新计算序号误差:
  1. set /a "y=(i-1)*250/91311+1, dd=i-(365*y + y/4 - y/100 + y/400)"
  2. 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 日的偏移
monthlenindsef(s)*1000f(e)*1000trunc(f(s))trunc(f(e))
Mar3100301399300
Apr30131601026197311
May31261912006298622
Jun303921213019396733
Jul3141221524000498044
Aug3151531835013599355
Sep3061842136026697366
Oct3172142447006798677
Nov3082452748019896788
Dec3192753059000998099
Jan311030633610013109931010
Feb281133736411026119081111
Feb(leap)291133736511026119411111


上面表格可用代码得到:
  1. @echo off & setlocal enabledelayedexpansion
  2. echo mon len ind    s     e     f(s)    f(e)  int(f(s))   int(f(e))
  3. for %%i in ("Feb 28" "Feb 29") do (
  4.   set /a "i=100, s=1000, e=1000-1"
  5.   for %%a in ("Mar 31" "Apr 30" "May 31" "Jun 30" "Jul 31" "Aug 31"
  6.   "Sep 30" "Oct 31" "Nov 30" "Dec 31" "Jan 31" %%i) do (
  7.     set "l=%%~a"
  8.     set /a "e+=!l:~-2!, fs=(5*(s-1000)+2)*1000/153+100000,fe=(5*(e-1000)+2)*1000/153+100000"
  9.     set /a "ifs=(5*(s-1000)+2)/153+100,ife=(5*(e-1000)+2)/153+100"
  10.     set "lin=!l!   !i:~-2!   !s:~-3!   !e:~-3!   !fs:~-5!   !fe:~-5!     !ifs:~-2!          !ife:~-2!"
  11.     set "lin=!lin: 0=  !"
  12.     if %%i=="Feb 28" (echo !lin!) else if %%a=="Feb 29" echo !lin!
  13.     set /a "s+=!l:~-2!, i+=1"
  14. ))
复制代码
综上, 由序号求日期的完全代码:
  1. set /a "y=(i-1)*250/91311+1,dd=i-(365*y+y/4-y/100+y/400)"
  2. set /a "y+=dd>>31,dd=i-(365*y+y/4-y/100+y/400)"
  3. set /a "mi=(5*dd+2)/153"
  4. set /a "m=(mi+2)%%12+1" & rem {0,1,2,..9,10,11} -> {3,4,5,..12,1,2}
  5. set /a "y-=m-3>>31" & rem 1,2月到了下一年
  6. set /a "d=1+dd-(153*mi+2)/5"
复制代码
PS:
计算序号的代码在开始可以将 y 加上一个 400 的倍数, 在序号转为日期的最后将 y 再减去同一个 400 的倍数, 这样就可以将年份 y 的可计算范围向负数调整. 上面 [1, 20000] 的范围分作正负数各一半仍足够大.

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

下面是一个算法互逆性测试代码, 但不能对序号生成的日期正确性作验证测试:
  1. @echo off & setlocal enabledelayedexpansion
  2. REM 7305155=20000.12.31
  3. for /l %%i in (0 1 7305155) do (
  4.   set /a "i=%%i"
  5.   set /a "y=(i-1)*250/91311+1,dd=i-(365*y+y/4-y/100+y/400)"
  6.   set /a "y+=dd>>31,dd=i-(365*y+y/4-y/100+y/400)"
  7.   set /a "mi=(5*dd+2)/153"
  8.   set /a "m=(mi+2)%%12+1" & rem {0,1,2,..9,10,11} -> {3,4,5,..12,1,2}
  9.   set /a "y-=m-3>>31" & rem 1,2月到了下一年
  10.   set /a "d=1+dd-(153*mi+2)/5"
  11.   set "test=%%i:   !y!.!m!.!d!"
  12.   set /a "m+=9,m%%=12,y-=m/10,i=365*y+y/4-y/100+y/400+(m*306+5)/10+d-1"
  13.   set "test=!test!     !i!"
  14.   if "!i!" neq "%%i" set "test=!test!     error"
  15.   echo !test!
  16.   if "!test:error=!" neq "!test!" pause
  17. )
复制代码
1

评分人数

    • plp626: 很喜欢这样的讨论,收益匪浅PB + 10 技术 + 1

TOP

本帖最后由 plp626 于 2012-4-9 17:25 编辑

回复 25# neorobin
  1. struct sdate dtf(long d) { /* convert day number to y,m,d format */
  2.   struct sdate pd;
  3.   long y, ddd, mm, dd, mi;
  4.   y = (10000*d + 14780)/3652425;
  5.   ddd = d - (y*365 + y/4 - y/100 + y/400);
  6.   if (ddd < 0) {
  7.     y--;
  8.     ddd = d - (y*365 + y/4 - y/100 + y/400);
  9.     }
  10.   mi = (52 + 100*ddd)/3060;
  11.   pd.y = y + (mi + 2)/12;
  12.   pd.m = (mi + 2)%12 + 1;
  13.   pd.d = ddd - (mi*306 + 5)/10 + 1;
  14.   return pd;
  15.   }
复制代码
原作者的C代码;我精简后,除了开头部分,后面的和你基本一致;

你代码中的开始一行中的常数250,91311,这些如何得来的?

TOP

本帖最后由 neorobin 于 2012-4-9 17:45 编辑

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

TOP

本帖最后由 plp626 于 2012-4-9 17:54 编辑

回复 27# neorobin


    这样啊,那你的代码可以再化简了
  1. ->m 1/365.2425
  2. 0.00273790701
  3. ->m 250/91311
  4. 0.00273789576
  5.  
  6. ->m 250/91310
  7. 0.00273792575
复制代码
可以看出,用91310代替91311更接近真实值;这样一来250,91310又可以化简了;同除以10;得25,9131;
----------
上面的接近分析有误,实际是:
|0.00273790701-0.00273789576| < |0.00273790701-0.00273792575|

TOP

回复 28# plp626


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

TOP

本帖最后由 neorobin 于 2012-4-9 18:03 编辑

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

TOP

返回列表