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 和其他数学软件与模块如何写作来完成一个项目，以下是一些提示（但不是规定）。