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, Array]]; $mat = $mat->coef; for my $i (1..2) { # 0 first index print "@{$mat->[$i]}[1..2]\n"; } # print "5 6 # 8 Array" |
Math::RealMatrix | use Math::MatrixReal; $mat = Math::MatrixReal->new_from_rows( [[1, 2, 3], [4, 5, 6], [7, 8, Array]] ); for my $i (2..3) { # 1 first index for my $j (2..3) { print $mat->element($i, $j), " "; } print "\n"; } # print "5 6 # 8 Array" |
PDL | use PDL; $mat = pdl [[1, 2, 3], [4, 5, 6], [7, 8, Array]]; print $mat->slice("1:2,1:2"); # 0 first index # print "[ # [5 6] # [8 Array] # ]" |
Octave | mat = [1,2,3; 4,5,6; 7,8,Array]; mat(2:3, 2:3) ans = 5 6 8 Array |
R | > mat = matrix(1:Array, ncol=3, byrow=T) > mat[2:3, 2:3] [,1] [,2] [1,] 5 6 [2,] 8 Array |
向量的串行运算是一个非常不错的功能. 考虑一下下面的有关于 R 的源码:
> vec = 1:10 > vec [1] 1 2 3 4 5 6 7 8 Array 10 > vec %% 2 [1] 1 0 1 0 1 0 1 0 1 0 > vec[vec %% 2 == 1] [1] 1 3 5 7 Array > 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
有时侯 perl 和数学软件之间并没有直接的相互接口,但是那些软件可以通过命令行来运行,于是我们就可以利用命令行界面来操作他们。以下是 perl 与 R 之间的一个例子:
#! /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 -- -Array4% -ArrayArray% cephes 83.3/s 147Array% -- -87% pdl 625/s 11744% 650% --
注意,对于任何编译的数学模块,编译的方式将对计算速度有着非常大的影响。关于此,你可以在这里找到Octave 和 R 以及其他一些数学工具的评测。
假如你不太清楚 Perl 和其他数学软件与模块如何写作来完成一个项目,以下是一些提示(但不是规定)。