maple

一直好奇一个关于 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 还是 log 和 exp，后来就想到用

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写了）：

有了 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) 这个简单的问题，俗谚云：“看人挑担不吃力，事非经过不知难。”诚哉：）

天天和各种程序的编辑器打交道，有一个恶习，就是苛求编辑器的颜色搭配和字体。譬如用 Lyx 编辑 Tex 文档之前，会不断的去修改它的 preferences 里的各项颜色设定（路径在 "Documents and Settings\Username\Application Data\Lyx"），以便使用习惯的配色。

也是因为配色的缘故，我舍弃了 WinEdt，舍弃了 Tinn-R，舍弃了 Visual Studio 去使用 Sublime Text 编辑 Tex，R，或者是 VB，VC++；也是因为不支持 TT 字体（并且界面太丑），舍弃了 Texnicenter —— 当然，从功能上讲，Sublime Text 并不算强大，譬如对 LaTex 的支持远不如 WinEdt 或者 Texnicenter，以至于编译复杂的 LaTex 文档时 （譬如使用 pstricks package 画图，必须使用 LaTex –> DVIPS->PS2PDF 的时候），就不得不先在 Sublime Text 里面编辑好，然后用 WinEdt 或者 Texnicenter 来编译，因为 Sublime Text 只支持简单的 PDFLatex （顺便感慨一下世界如此之小， Sublime Text 的 LaTex 插件的开发者，竟然是我的微观经济学老师 。）

其实与其说喜欢 Sublime Text，倒不如说喜欢 Monokai 这款颜色搭配（如上图），不仅纯粹从视觉效果上来说，这样的灰底彩色搭配更容易减缓视觉疲劳（现在绝对不用白底黑字，看一个小时以后，眼睛肯定要花），并且这个色彩的名字跟我很有缘分（迷信！），所以在无可替代的时候，譬如 Lyx，譬如 Maple，就只好改这些程序的 style 设定了，好在现在一般的程序都支持客制化。为了方便，特提取了 Monokai 所包含 11 种颜色的 Hex 码如下，方便使用。

要搁在战国，我肯定是那个买椟还珠的郑人：“楚人有卖其珠于郑者，为木兰之柜，薰以桂、椒，缀以珠玉，饰以玫瑰，辑以翠羽。郑人买其椟而还其珠。”：）

收到 Maple 的电子邮件，Maple 14 上市了，不知道新的版本会有多大的更新，希望很快能就能拿到来测试一下。Maple 是我使用频率最高的数学软件之一，几乎所有的 paper 都是先在 Maple 里演算并且简要写完，然后再在 LaTex 里完工。事实上，写 paper 最 enjoy 的就是在 Maple 中尝试各种新的想法，满足自己的好奇心，真的到了需要把 paper 彻底完成、在 LaTex 里编译出来时，却觉得无比的枯燥，因为最新奇的部分已经结束，而一点点的写 paper 的过程就像吃自己嚼过的馒头 ，毫无新意…… 当然从客观而不是情绪的角度讲，写的部分却是最重要的。

不过 Maple 这样的工具也让人变得越来越笨了，譬如用惯了 dsolve()，结果现在手算解微分方程的水平大大的降低，同样的因为 diff(), int(), limit()，彻底的把常用的微分、积分和极限公式忘到九霄云外； Statistics package 大大简化了涉及随机变量的符号运算，特别是复杂的 conditional probability 和 expectation 的计算，但结果是在某次的 talk 里，我想举一个简单的贝叶斯更新的例子，但是卡在那里半天算不出 conditional on  signal 的 expectation。

呵呵，依赖工具对于人脑来说，不知道是进步还是退步；不过念 Ph.D. 确实会造成 permanent head damage，这点是肯定的。

P.S. Office 2010 RTM 也出来了，从订阅的 MSDN 那里得到了注册码和下载链接。用了几天，最直观的感觉就是比 beta 版更漂亮了，其它的新功能，唯一值得一提的是 OneNote 总算支持数学符号输入了，并且还能作简单计算。

终于搞到了这个最新的版本，正在试用中，外观上和 Maple 12 没有太大的差别，据说效率提高了很多，正在试验几个小程序，看看结果到底如何。

收到 Maplesoft 的邮件， Maple 13 发布了。有一些新的特色，不过最期待的，是能把它的 solver 做的好一些，特别是 Optimization 的，可能是因为没有用单独另外发售的 Opitmization Toolbox 而只是用了自带的 with(Optimization), 一个简单的 piecewise function 的极值它都搞错了，这个版本里据说对于 differential equations 的支持会增强，希望能惠及优化，虽然别的软件优化要好很多，但是符号运算确实不如 Maple 直观方便，不知道什么时候能拿到最新的 Maple 13。期待！

呵呵，不经意两者就扯上关系了。正在写一篇 paper，中间牵扯一系列复杂的符号微分，用 Maple 算好以后懒得手动的输入到 LaTex 环境中，于是就用 Maple 10 的 Export As... 直接导出为 LaTex 文件，不过这个文件却不通过 MikTex 的编译，于是就 google 一下看看是什么问题，用了" maple export latex "作为关键词，第一个结果就挺相关的，点击去一看，网站的界面蛮熟悉，心里想该不会是 Drupal 做的吧，虽然网站用了 url alias，但是点了一下注册，确认就是 Drupal 做的，呵呵， Drupal 的SEO 很自然的就把这个网站放在搜索结果的第一位了。反复浏览了一下，觉得一些设计挺好的，这个叫做 Maple Primes 的网站也使用了 User Points 模块，并且还根据 Points 给用户分了等级，并用不同的枫叶来表示，很有趣，所以就开始研究它是如何做的，因为正在给 Drupal China 做改造，也使用了 User Points 模块，所以就像学习一下。大致研究了一番，Maple Primes 应该是自己写了一个叫做 maplerate 的模块，这个模块可以根据 User Points模块生成的 Points 将用户分为不同的级别，显示不同的图标，这个想法挺好的，于是就开始琢磨如何写这个模块，又做了一下初步的测试，然后突然想起来，其实我到这个网站是有正经事情的。。。嗯，去找 Maple 输入 LaTex 的问题的解答。。。ft，跑了大半天的神，现在折腾到凌晨4点多了，sigh。。。