Perl 与数学

一直都有人问讨论有关于 Perl 与数学的问题。有时候一些 perl 玩家问如何使用 perl 做一些高级数学的工作。另一方面，一些数学家又反过来问如何利用 perl 来帮助完成他们本身的工作。所以，现在我提供一些方便的参考文献，比较和说明一些常用的 perl 数学模块，以及对 perl 有用的软件。它并不是完整的Perl数学编程手册，而只是对于一些常用的数学模块和软件的简要综述。我忽略了 bioperl 以及生物信息学的内容，因为他们涵盖范围太广，无法简单地称之为“数学”。

一般来说，如果你在 CPAN 上搜索与数学相关的模块，那么你应该从以下关键字入手： Math::×， Statistics::×， 以及 AI::× ，Algorithm::×， Cript::×， Date::， Graph::， GraphViz::， Inline::， 等等。 GNU项目也是
寻找 perl 扩展模块的好地方。

Perl数学模块及其相关软件
以下列举了一些有关于数学的 perl 模块及软件：

 

 

模块/软件	类型	可用性	评论
Inline::C	通用模块，perl 与 C 语言的接口	任何操作系统	可能不是你用 perl 做数学的第一选择
Inline::Octave & Octave	矩阵代数以及数值分析	Inline::Octave 无法从ActiveState获得。 Octave 可以在任何操作系统运行	Matlab的开源版本，运行快，语法简单。Octave 可以交互式的使用，或用做脚本语言。其他类似的商业软件有：Gauss，APL等。
Math::Cephes	通用工程数学库	任何操作系统	特征系统（Eigensystem）仅适用与实对称矩阵。
Math::LP	线性编程	ActiveState perl 没有，LP solve则有 DOS 的版本。	免费，但 LP Solve 可能不再持续维护了。类似的商用软件有：CPLEX，Lindo，Minos，AMPL等。
Math::Pari	数论	任何操作系统	很适合密码学分析
Math::MatrixReal	矩阵代数（仅适用用于实数）	任何操作系统	纯 Perl 代码实现，特征系统仅用于实对称矩阵
PDL	Perl 数据语言（Perl Data Language）	最新版本为2.4.0，windows预编译版2.3.1。	为 perl 提供更多的数学语法
R	Statistical software	任何操作系统	可交互式使用或作为脚本语言使用。在作图方面非常出色（可以参考这里）最新版本1.7.1（译者注：本文翻译的时候的最新版本为2.1.0）是Splus的开源版本。可作为 windows 系统的双向“nix & COM”服务器界面（从R v.1.7.0开始）引自 Omegahat。类似的商业软件：SPSS，SAS等。

　　语法比较

　　我们来比较一下这些模块和软件的语法。我们以运算一个2×2矩阵与一个向量的乘积为例，以下是它们各自的语法。

 

模块/软件	代码
Math::Cephes	use Math::Cephes::Matrix qw(mat); $a = mat [[1, 2], [3, 4]]; $b = [1, 1]; print "@{$a->mul($b)}"; # print "3 7"
Math::RealMatrix	use Math::MatrixReal; $a = Math::MatrixReal->new_from_rows( [[1, 2], [3, 4]] ); $b = Math::MatrixReal->new_from_cols( [ [1, 1] ] ); print $a*$b; # print "[ 3.000000000000E+000 ] # [ 7.000000000000E+000 ]"
PDL	use PDL; $a = pdl [[1, 2],[3, 4]]; $b = pdl [[1], [1]]; print $a x $b; # print "[ # [3] # [7] # ]"
Octave	>> a = [1, 2; 3, 4]; >> b = [1; 1]; >> a*b ans = 3 7
R	> a = matrix(c(1,2,3,4), ncol=2, byrow=T) > b = matrix(c(1,1), ncol=1) > a%*%b [,1] [1,] 3 [2,] 7

　　仅就数学上来说，以上这些模块和软件的语法看起来都相当简洁。不过 Octave 和 R 比 perl 还是要好懂的多了。许多数学语言的一个突出特点是它们的“向量操作”以及“下标操作”。

具个例子来讲，例如现在我们要从一个矩阵中，提取一个子矩阵。

实例：提取一个子矩阵

 

模块/软件	代码
Math::Cephes	use Math::Cephes::Matrix qw(mat); $mat = mat [[1, 2, 3], [4, 5, 6], [7, 8, 9]]; $mat = $mat->coef; for my $i (1..2) { # 0 first index print "@{$mat->[$i]}[1..2]\n"; } # print "5 6 # 8 9"
Math::RealMatrix	use Math::MatrixReal; $mat = Math::MatrixReal->new_from_rows( [[1, 2, 3] [4, 5, 6] [7, 8, 9] ) for my $i (2..3 { firs ndex fo ..3) {
PDL	use PDL; $mat = pdl [[1, 2, 3], [4, 5, 6], [7, 8, 9]]; print $mat->slice("1:2,1:2"); # 0 first index # print "[ # [5 6] # [8 9] # ]"
Octave	mat = [1,2,3; 4,5,6; 7,8,9]; mat(2:3, 2:3) ans = 5 6 8 9
R	> mat = matrix(1:9, ncol=3, byrow=T) > mat[2:3, 2:3] [,1] [,2] [1,] 5 6 [2,] 8 9

　　向量的串行运算是一个非常不错的功能. 考虑一下下面的有关于 R 的源码:

　　> vec = 1:10
　　> vec
　　[1] 1 2 3 4 5 6 7 8 9 10
　　> vec %% 2
　　[1] 1 0 1 0 1 0 1 0 1 0
　　> vec[vec %% 2 == 1]
　　[1] 1 3 5 7 9
　　> vec[vec %% 2 == 1] + 1
　　[1] 2 4 6 8 10

　　我们刚才在前面也提到了 Pari，但它只是一个自成体系的模块。我们来看一个 Pari 的例子。

　　#! /usr/local/bin/perl -w
　　use strict;
　　# -----------------------------------------------------------
　　# RSA algorithm -- assymetrical\public-key cryptography
　　# -----------------------------------------------------------
　　use Math::Pari qw(gcd PARI) ;
　　# -----------------------------------------------------------

　　# m -- message
　　my $m = 'Perl' ;
　　print "original: $m\n" ;
　　my $tmpl = 'C*' ;
　　my @m = unpack($tmpl, $m) ; # string -> unsigned char values
　　print "coded: @m\n" ;

　　# n = pq -- in RSA, p & q = prime, each 1024 bits/308 digits long
　　my $p = PARI("prime(".int(rand 50).")") ;
　　my $q = PARI("prime(".int(rand 50).")") ;
　　my $n = $p*$q ; # $n = Pari's obj

　　# choose a random number r, s.t.
　　# 1 < r < (p-1)(q-1) = b
　　# gcd(r, b) = 1 -- relative prime
　　my $b = ($p-1)*($q-1) ;
　　my $r ;
　　do {$r = int rand $b ; } until (gcd($r,$b) == 1) ;
　　$r = PARI $r ;

　　# rk = 1 mod (p-1)(q-1) -- k = private key; (n, r) public
　　my $k = (1/$r)%$b ; # Pari's math operators, since vars = Pari

　　# encrypt -- c = (m ^ r) mod n
　　my @c ;
　　map { $c[$_] = ($m[$_]**$r)%$n } 0..$#m ; # Perl: ** for power
　　print "ciphered: @c\n" ;

　　# decrypt -- m = (c ^ k) mod n
　　my @d ;
　　map { $d[$_] = PARI("($c[$_]^$k)%$n") } 0..$#c ; # Pari: ^ for power
　　print "deciphered: @d\n" ;
　　print "decoded: " . pack($tmpl, @d) . "\n" ;

　　__END__
　　original: Perl
　　coded: 80 101 114 108
　　ciphered: 18431 6512 5843 7236
　　deciphered: 80 101 114 108
　　decoded: Perl
　　</font></tt>

　　有时侯 perl 和数学软件之间并没有直接的相互接口，但是那些软件可以通过命令行来运行，于是我们就可以利用命令行界面来操作他们。以下是 perl 与 R 之间的一个例子：

　　<tt class="code"><font size="-1">#! /usr/local/bin/perl -w
　　use strict ;

　　R("getwd()");

　　sub R {
　　　　my $Rpath = "C:\R\rw\bin\" ;
　　　　my $Rcmd = $Rpath . "rterm --vanilla --quiet --slave" ;
　　　　my $Rscript = shift ;
　　　　$Rscript =~ s/(\r|;\r)/ ;/gm ;
　　　　$Rscript =~ s/<-/=/gm ; # \r or <- will break "echo"
　　　　return `echo $Rscript | $Rcmd` ;
　　}

　　如果你只有 DOS LP Solce 的可执行程序，那么你也可以这样操作：

　　my $dir = 'D:\tmp\prog\math\lp32';
　　my $lp_solve = "d:\mydir\lp_solve.exe";
　　my $lpfile = "d:\mydir\model.lp";
　　(my $lp = << " EOF") =~ s/^\s+//gm;
　　min: 8 x1 + 15 x2 ;
　　c1: 10 x1 + 21 x2 > 156 ;
　　c2: 2 x1 + x2 > 22 ;
　　EOF
　　open OUTFILE, "+>$lpfile";
　　print OUTFILE $lp;
　　close OUTFILE;
　　$output = qx($lp_solve < $lpfile);
　　$output =~ s/\r//gm;
　　print $lp;
　　print $output;
　　system("del $lpfile");

　　这也是为什么 Perl 被成为“胶水语言”的原因。

　　评测

　　人们做数学的时候往往很关心速度，所以我们来做些评测。

　　use strict;
　　use warnings;

　　my ($a1,$a2,$a3,$ref);

　　for my $x (0..20) {
　　　　for my $y (0..20) {
　　　　　　$ref->[$x][$y] = rand 100;
　　　　}
　　}

　　use Math::Cephes::Matrix qw(mat);
　　$a1 = mat $ref;

　　use Math::MatrixReal;
　　$a2 = Math::MatrixReal->new_from_rows( $ref );

　　use PDL;
　　use PDL::Slatec;
　　$a3 = pdl $ref;

　　use Benchmark qw(cmpthese);
　　cmpthese(100,
　　　　{
　　　　　　cephes=>sub{$a1->inv()},
　　　　　　matrixreal=>sub{$a2->inverse()},
　　　　　　pdl=>sub{matinv($a3)}
　　　　}
　　);

　　__END__
　　Rate matrixreal cephes pdl
　　matrixreal 5.28/s -- -94% -99%
　　cephes 83.3/s 1479% -- -87%
　　pdl 625/s 11744% 650% --

　　注意，对于任何编译的数学模块，编译的方式将对计算速度有着非常大的影响。关于此，你可以在这里找到Octave 和 R 以及其他一些数学工具的评测。

　　假如你不太清楚 Perl 和其他数学软件与模块如何写作来完成一个项目，以下是一些提示（但不是规定）。

流程	Perl	数学模块/软件	评论
数据存储	N/A	N/A	很多人用电子表格存储数据。不过可能用数据库会更好些，比如Oracle，MySQL，MS SQL等等。这部分工作和 Perl 和数学工具都关系不太大。
数据存取	是	可能	这部分工作是“插入”，“更新”和“删除”数据。你可能会用DBI 来做这件事。如果你的数学软件可以直接完成这些数据操作的工作，那部分就不再需要 Perl 了。
数据操作	YES	可能	数据操作的意思是以正确的格式获取原始数据以用于数据分析。比如从 log 文件中提取一系列的时间，把一些格式混乱的文本文件转换为 CSV 文件，等等。当这项工作变得复杂的时候，Perl 可能比直接操作电子表格或者数学软件方便得多，同时你可能会用到很多正则表达式。
数据分析	可能	是	在做复杂数学运算的时候， Perl 没有 C 那么快的速度，也没有 Octave 或其他数学模块那样简单的语法，所以数据分析的部分常常是由专业的数学软件完成的。当然，如果你大部分工作都是用的 Perl，而你需要的数学又相对简单，那么再去求助专业软件就可能是杀鸡用牛刀了。
生成报表	是	可能	一些数据库（如Oracle）或数学/统计学软件（比如SAS）提供“报表生成器”。如果这些能够满足你的要求，那么你就用它好了。但是如果你需要更复杂的报告，Perl 可能是你的选择。 你经常会用 Perl 来创建 HTML 格式的报告（比其他格式如 PDF 要简单得多），然后将 HTML 再转换成其他格式，比如 WORD 或 PDF。有些人会更关心 XML。XML实际上是展示数据的一种方法，但不是用来做视觉上的数据展示。用设计与开发的术语来说，XML 可能更接近于数据存储，而离数据报表更远。数据存储的工作是你应该在整个项目早期考虑的事情，而不是在数据汇报的时候。

　　越来越多的人，通过互联网来获取数据，例如 XML 等格式，然后对这些数据进行分析（股票数据，DNA 数据，等等）。

	Perl	数学模块/软件	可能用到的工具
XML 获取	是	可能	LWP, XML::libXML
数据处理	是	可能	正则表达式
数据分析	可能	是	前面提到的软件或工具
汇总报告	是	可能	Text::Template, HTML::Template, XML::libXSLT