Q

素数·派

  • Posted on: 17 March 2015
  • By: kzeng

π 节(3月14日)已经过去了,所以现在写 π 就当是向魔幻圆周率致敬吧。其实数学的发展最需要牛人,就像祖冲之,冷不丁的从历史里给你扔出一个 355/113 (=3.1415929204) 来,等你回过神来,这个 355/113 是如何得到的已经湮没在历史的尘埃中了,所以牛人之后还需要不太牛的传承者。也是因为中国的数学研究一直缺乏传承(后来很大一部分因为科举),在元朝时数学大跃进之后,到了明朝,就大踏步的倒退到了简单算术(虽然珠算得到了很大的发展)。

相比较而言,欧洲的数学一直到十六世纪初都不过尔尔(古希腊的数学传统到那时已经断掉),但是在突然冒出了笛卡尔,莱布尼茨,牛顿,伯努利等人以后,数学的发展风驰电掣,而到了欧拉,简直是神一样的人物了。

莱布尼茨给出了一个用无穷级数求π的方法:

pi= 4*sum (-1)^n/(2n+1)

在 Q 里快速的实现一下:

 {sum raze 4%(1 -1)*/:(0N,2)#1+2*til x} 1000000
 3.1415916536

上例一共用了一百万项,精度略低于祖冲之的密率。

然后是欧拉用 p-级数的方法:

pi = sqrt( 6 * sum (1/n^2))

在Q里快速实现一下:

` {sqrt 6 * sum 1 % {x*x} 1_til x} 1000000 3.1415916987

同样一百万项,收敛速度和莱布尼茨的无穷级数差不多,但是作为大神,欧拉进一步证明了上述p-级数与素数乘积的关系:

sum 1/n^s = prd 1/(1- p^-s)   p 是所有素数

当 s = 2 时, 我们就可以得到

prd 1/(1-p^-2) = sum 1/n^2 = pi^2 /6 

这等同于

如果我们从 1 到 N 的N个自然数中随意抽取两个数,那么这两个数互质(最大公约数为一)的概率在当 N 趋近于正无穷的时候趋近于 6 / pi ^ 2

以上推论的证明很简单,如果两个数互质,那么它们没有共同的素数因子。任意一个自然数被一个素数 p整除的概率是 1/p,譬如,能被5 整除的数字每间隔5个出现一次。任意两个自然数能被 p 整除的概率是 1/p^2. 至少一个不能被整除的概率是 1-1/p^2。把这个结论扩展到任意的 p,那么这个概率就是 prd (1- 1/p^2) for all p。 简单的变形,我们知道这个概率是 6/pi^2

简单的蒙特卡洛一下:

{sqrt 6% avg 1={$[y=0;:x;:.z.s[y;x mod y]];}./:(0N;2)#(1+ x?x)} 1000000
3.1431107401

从一百万个自然数中随机选取50万对,然后算互质的概率。上述模拟中求最大公约数时用的是欧几里得算法:{$[y=0;:x;:.z.s[y;x mod y]];}

转了一圈,终于回到了古希腊:)

Blog分类: 

脍炙《通鉴》

  • Posted on: 23 September 2012
  • By: kzeng

(这是一篇关于很枯燥的技术,很枯燥的历史文本,和不太枯燥的统计的 blog)

看过一篇关于《全宋词》词频统计文章,挺有趣的,想用类似的方法处理一下《资治通鉴》,所以就趁周末花了几个小时作了一下。

词是长短句,统计两个字组成的词频比较合适,《通鉴》是古文,文字结构不同,所以我统计了单字频,两字词词频,三字词词频,四字词词频,和五字词词频。同时也记录各个统计单位(字或词)出现的卷数。《通鉴》294卷,从三家分晋到五代结束共共1362年,所以卷数可以作为时间的度量。

《全宋词》的词频是用 R 作的。R 虽然是不错的统计软件,也是我的最爱之一,但是 R 并不适合作文本分析,更不适合来作数据库操作。所以就用了 C# 和 Kdb +3.0。 C# 用来分析文本,.Net 是懒人的福音,并且多线程运算非常简单,能够大大提升文本处理速度,Kdb+用来储存数据,它差不多是性能最好的 in-memory 数据库了,从它的网站上能下载到免费版本。这个分析里数据库是重头戏,因为需要查询数百万行的数据 row,如果用 MySQL,估计会龟速到死。另外 Kdb + 本身只有 300多K,不用安装,很方便。还有就是 Kdb+ 的 Q 语言也能满足编程需要。

Kdb+ 的网站提供了各种语言 API 的源码,C# 的 API 不支持多线程,所以需要在适当的地方加锁。Kdb 唯一的问题是不支持 UTF-8。它用的是 UTF-7,所以在注入中文文字数据的时候可能会出现乱码,为了省事,从 C# 里 publish 数据的时候,直接 publish 为三字节的 int[] 了。query kdb 时用了一个免费的 GUI QPad。QPad 似乎是用 Java 写的,它的编码默认是 UTF-8,所以在 query  Kdb 的时候直接把三字节的 int vector  cast 成 char,在 QPad 里显示的就是中文了,所以也很方便。

产生数据的 C# 代码非常简单,发布数据的时候自动生成 Kdb 的 schema。使用的《资治通鉴》的文本是网上广为流传的国学网简体版,在生成数据前,先用 C# 作了预处理,主要是用正则表达式替换掉了现代语言的“污染”(譬如:“后一页”,公元xxx年 等)

下表是各个字、词频的数据量:

类别 数据行数
单字            2,586,329
双字            2,102,023
三字            1,633,875
四字            1,221,713
五字                851,403

从上表看,《资治通鉴》应该有近两百六十万字。

单字的字频统计如下:

排名 次数 百分比 累积百分比 
1 66087 2.56% 2.56%
2 39874 1.54% 4.10%
3 35677 1.38% 5.48%
4 34376 1.33% 6.81%
5 21578 0.83% 7.64%
6 21279 0.82% 8.46%
7 20182 0.78% 9.24%
8 20100 0.78% 10.02%
9 19035 0.74% 10.76%
10 18209 0.70% 11.46%
11 18083 0.70% 12.16%
12 使 17160 0.66% 12.82%
13 16116 0.62% 13.45%
14 16031 0.62% 14.07%
15 15600 0.60% 14.67%
16 15558 0.60% 15.27%
17 15252 0.59% 15.86%
18 14746 0.57% 16.43%
19 12826 0.50% 16.93%
20 12536 0.48% 17.41%

“之”字当之无愧的排在了第一位。第一个非虚词是“王”,它包含了姓和爵位,第一个动词是“曰”。“人”的频率也很高,“将”,“军” 在双字词频中也会遇到。“帝”字排名32,“后”字排名33。但是因为是简体字,“后”并不专指皇\王后。

下面是价值观念的排名:

排名 次数 百分比
152 3507 0.14%
181 3004 0.12%
223 2475 0.10%
240 2287 0.09%
253 2190 0.08%
294 1935 0.07%
767 694 0.03%

义、忠、孝排名在前,智排名最后,倒正印证了司马温公那句话:“凡取人之术,苟不得圣人、君子而与之,与其得小人,不若得愚人。”当然这个数据里噪音很多。

另外还有很多有趣的东西,就不一一叙述了。下面看一下两个字的词频:

排名 最早卷数 次数 百分比
            1 将军 1             6,176 0.29%
            2 刺史 21             4,790 0.23%
            3 州刺 21             4,110 0.20%
            4 节度 29             3,698 0.18%
            5 以为 1             3,479 0.17%
            6 度使 203             3,202 0.15%
            7 天下 1             2,972 0.14%
            8 尚书 20             2,742 0.13%
            9 太子 1             2,584 0.12%
          10 陛下 6             2,492 0.12%
          11 不能 1             2,375 0.11%
          12 不可 1             2,351 0.11%
          13 太后 3             2,165 0.10%
          14 皇帝 6             2,050 0.10%
          15 太守 5             2,010 0.10%
          16 大将 6             1,813 0.09%
          17 遣使 4             1,501 0.07%
          18 司马 1             1,480 0.07%
          19 二月 4             1,477 0.07%
          20 馀人 2             1,463 0.07%

“将军”出现的频率最高,在第一卷里就出现了,“度史”显然是“节度使”里出现的,虽然在203卷才出现,但是它居然出现了3202次,唉,唐朝啊!“节度”一次出现的要比“节度使”早。“皇帝”一次最早在第6卷出现,其实那时还是昭襄王元年,但是因为文本中出现了“秦始皇帝上”。

两字地名出现的最多的是“长安”,排名43,最早出现在第5卷,不过那里的“长安”并不是长安城,而是赵国的长安君。“洛阳”其次,排名81,最早出现在第2卷,三家分晋不久,洛阳附近就成了三晋与秦国的战场。

三个字的词频:

排名 最早卷数 次数 百分比
1 州刺史 21 4102 0.25%
2 节度使 210 3195 0.20%
3 大将军 6 1547 0.09%
4 平章事 203 933 0.06%
5 同平章 203 901 0.06%
6 十二月 4 704 0.04%
7 之子也 13 700 0.04%
8 十一月 7 686 0.04%
9 部尚书 70 655 0.04%
10 指挥使 254 578 0.04%

比较有趣的是“之子也”,老子英雄儿好汉。

四字字频:

排名 最早卷数 次数 百分比
1 同平章事 203 900 0.07%
2 仪同三司 49 403 0.03%
3 都指挥使 254 374 0.03%
4 日有食之 1 368 0.03%
5 中书侍郎 73 325 0.03%
6 节度使李 217 312 0.03%
7 散骑常侍 69 308 0.03%
8 开府仪同 79 300 0.02%
9 府仪同三 79 283 0.02%
10 御史大夫 9 281 0.02%

鉴于“州刺史”在三字字频中频繁出现,所以频率出现比较高的各个州刺史的频率单列出来:

排名 最早卷数 次数 百分比
15 豫州刺史 58 220 0.02%
19 荆州刺史 49 195 0.02%
26 兗州刺史 37 172 0.01%
32 徐州刺史 30 167 0.01%
38 扬州刺史 24 152 0.01%
42 雍州刺史 65 147 0.01%
52 江州刺史 86 122 0.01%
64 益州刺史 39 113 0.01%
68 二州刺史 49 111 0.01%
74 秦州刺史 79 109 0.01%
101 青州刺史 21 91 0.01%
105 梁州刺史 84 90 0.01%
107 冀州刺史 27 88 0.01%
114 并州刺史 52 84 0.01%
135 凉州刺史 31 71 0.01%
165 幽州刺史 50 62 0.01%
196 广州刺史 80 55 0.00%

豫州刺史第58卷时才登场,但是雄踞第一,而豫、荆、兖、徐、扬也勾勒出了中国政治地理的热点。顺便提一句,最早登场的豫州刺史是王允,而最早登场的荆州刺史是杨震。

五字词已经没有太大的意义:

排名 最早卷数 次数 百分比
               1 府仪同三司 79 283 0.03%
               2 开府仪同三 79 283 0.03%
               3 尚书左仆射 77 167 0.02%
               4 皇帝上之下 10 140 0.02%
               5 为中书侍郎 84 118 0.01%
               6 尚书右仆射 81 115 0.01%
               7 军都指挥使 256 111 0.01%
               8 骠骑大将军 39 107 0.01%
               9 河东节度使 214 101 0.01%
            10 督中外诸军 74 97 0.01%

最后看看慕容家的英杰们谁的全名被提到的次数最多:

次数
慕容彦超 27
慕容垂 26
慕容廆 24
慕容绍宗 15
慕容恪 18
慕容评 18
慕容皝 17
慕容农 15
慕容翰 15
慕容仁 12
慕容白曜 10

似乎是慕容彦超险胜慕容垂……慢着!慕容垂最初的名字是慕容霸,而慕容霸被提及了10次,所以慕容垂以 36 次远远胜出 (慕容缺这个全名并没有出现在《通鉴》中:))。

Blog分类: 

(-1)^(1/3):从 C++ 到傅立叶变换

  • Posted on: 17 January 2012
  • By: kzeng

image 一直好奇一个关于 power function 算法的问题,直到自己动手写 power function。

在不同的编程语言里,遇到过同一个问题:(-1)^(1/3) 是多少。 很显然,在实数集里,这个表达式是有意义的,正如右图 Google Calculator 给的结果,在实数集里,它等于 -1。但是在很多编程语言的实数集运算中,这个表达式是无意义的。

譬如在 C++ 里,用 cmath library, 当你计算 (-1)^(1/3) 时,你得到输出结果是:

 -1.#IND 

C#里,用 System.Math, 作同样的计算,得到的输出结果:

NaN

R 里,不使用复数集,得到的结果也是一样的:

[1] NaN

同样的在 VBA 里,如果不 call Excel的 power function ( i.e. Application.WorksheetFunction.Power) ,而是直接使用 ^, 得到的结果仍然是 run-time error. 其它的编程语言也类似。

复数集中,如果用 MATLAB, 得到的结果:

 0.5000 + 0.8660i

Maple,结果:

.5000000001+.8660254037*I

R 的复数集(i.e. as.complex(-1)^(1/3)),得到的结果:

[1] 0.5+0.8660254i

很显然,所有的语言的 power function 用的是同一种算法。这种算法无法得到实数解,而复数解得到是同一个数值。因为 power function 太基本,虽然有疑问,但是也没有过多的想这个问题,直到后来用 Q

Q 的语法和 C++/C#/JAVA,或者 MATLAB/R/MAPLE 都不太一样,"^" 符号的定义和 C++/C# 相同,不是 power function。刚开始用的时候,不知道 Q  的 power function 是 xexp,觉得 power function 又不难,自己写一个吧,但是真正开始写,却又卡壳了:(int,int) 的函数好写,那 (double, double) 的呢?

拿出算法圣经《 Numerical Recipe 》(第三版),但是却发现它没有给出 power function 的算法,大概是太基础了吧,所以自己又想了一下,幸亏 Q  里的 log 和 exp 还是 logexp,后来就想到用

exp(y*log(x)) = exp(log(x^y)) = x^y

来定义 power function,解决非整数的问题。这样的以来,一般的问题都解决了,但是因为用到了 log(x) , x 的值必须非负(0 的问题可以很简单的处理),所以一旦 x<0,这个算法就不适用了 —— 这时才突然的想到莫非上面的那些问题的症结正在此?实数集的问题是由于 log(x), x<0 在实数集里无定义,那复数集呢?在 MATLAB 里试了

>> exp(1/3*log(-1))
ans =
   0.5000 + 0.8660i

果然是这样的。所有的算法都依赖于 log 函数来获得 power 函数的值,这导致了上述问题在实数集无定义,而在复数集因为 log(-1)  =  3.1416i  这个默认值,导致了 0.5000 + 0.8660i  这个结果。

但是问题还没有结束。看到 log(-1),自然想到了 2*log(i), 然后自然而然的想到傅立叶变换里常用的 trick 可以解出 log(-1)的一般表达式(为了省事儿,下面用 LaTex写了):

image

有了 log(-1)的通解,我们可以让 power function 获得任意 x^y, x<0 的所有解。譬如 (-1)^(1/3)简单测试一下,在 MATLAB 里,

x=2*i*(pi/2+2*(-3:3)*pi)
exp(1/3*x)

得到:

x =
        0 -34.5575i        0 -21.9911i        0 - 9.4248i        0 + 3.1416i        0 +15.7080i        0 +28.2743i        0 +40.8407i
ans =
   0.5000 + 0.8660i   0.5000 - 0.8660i  -1.0000 - 0.0000i   0.5000 + 0.8660i   0.5000 - 0.8660i  -1.0000 + 0.0000i   0.5000 + 0.8660i

结果里面包含了它的所有三个解。因为当 x<0, x^y = (-1)^y*(|x|)^y,所以只要有 (-1)^y ,就可以得到任何负数的 power function ( exp(a+bi)也可以用上面的方法转化成三角函数来解)。

上面的长篇累牍都起源于一开始的时候不知道 Q 的 power function 是 xexp,但是如果不是自己去写 power function,恐怕也没有机会搞明白 (-1)^(1/3) 这个简单的问题,俗谚云:“看人挑担不吃力,事非经过不知难。”诚哉:)

Blog分类: