matlab

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

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分类: 

R 的 vectorization, again

思考的方式,大概会受一个人正在说的语言影响,写 codes 似乎也是一样。如果是写 C++,思考的习惯大概会和 Matlab 不一样——不是说 syntax 或者是 concept,而是解决问题的抽象方式。

一个简单的题目,模拟一个 Markov Pure Jump Process (MPJP), transition matrix 是一个 8x8 的矩阵,模拟一百万次,先在 R 中写,习惯的先写一个1:8的loop 模拟一条路径,然后重复一百万次:

   1:  exp_sim=function(n){          
   2:    sim=array(0,c(n,64)) 
   3:    for (k in 1:n){ 
   4:      for(i in 1:8){ 
   5:        j=i 
   6:        t=rexp(1,rate=lam[j])           
   7:        while (t<5){ 
   8:          j=sample(1:8,1,prob=P[j,])    
   9:          t=t+rexp(1,rate=lam[j])       
  10:          }                             
  11:        sim[k,((i-1)*8+j)]=1 
  12:        } 
  13:      } 
  14:    result=array(0,c(8,8)) 
  15:    for (i in 1:8){ 
  16:      for (j in 1:8){ 
  17:      result[i,j]=mean(sim[,((i-1)*8+j)]) 
  18:      } 
  19:      } 
  20:      result 
  21:  }

结果消耗的时间可想而知,耗费了大约 340 秒才完成,实在太慢了,以为是 R 的问题,又在 matlab 里重新写过,运行大约耗费了 23 秒,正感叹 Matlab 速度快,优化作的好。。。可马上也意识到在 matlab 里写 codes 的思路竟然完全不同——根本没有去作那个 1:10^6 的 loop,而是按照习惯,一开始就生成了一个一百万行的 vector,避开了这个耗时的 loop。。。所以又参照 matlab 的 codes,重新写 R 的 codes:

   1:  exp_simv2=function(n){ 
   2:    mysample=function(vec){ 
   3:      m=length(vec) 
   4:      temp=numeric(m) 
   5:      for (i in 1:8){ 
   6:     temp=temp+(vec==i)*sample(1:8,m,prob=P[i,], replace=TRUE) 
   7:      } 
   8:      temp 
   9:      } 
  10:      
  11:      result=array(0,c(8,8)) 
  12:      for (i in 1:8) { 
  13:          t=array(0,n) 
  14:          state=i*array(1,n) 
  15:          t=t+rexp(n,rate=lam[state]) 
  16:          while (min(t)<5){ 
  17:              newstate=mysample(state) 
  18:              state=newstate*(t<=5)+state*(t>5) 
  19:              t=t+rexp(n,rate=lam[state]) 
  20:              } 
  21:              for (j in 1:8){ 
  22:                  result[i,j]=sum(state==j) 
  23:              } 
  24:          } 
  25:          result=result/n 
  26:          result 
  27:  }

像 matlab 一样思考,避开了那个一百万次的 loop,虽然写的时候有些别扭,但是运行以后发现,速度果然快了不少,效率大约提高了10倍,最后只耗时30秒左右,基本上和 matlab 的速度相仿。

习惯的力量真可怕 –_-#

P.S. 不过也发现了一个新的让人欣慰的习惯,不管是写 R 还是 matlab,已经不自觉习惯于 encapsulation,这肯定是受 C++ 的影响吧。看来习惯也是可以改的:)

Blog分类: