《Mathematica科学计算与程序设计》-CH07

第七章Mathematica程序设计7.4模块7.5程序设计实例7.6本章小结Mathematica软件平台上用于科学计算的指令和语句统称为Wolfram语言，Wolfram语言属于高级计算机语言，也是全球顶级的科学计算语言，是物理学等自然科学研究和数学与密码学等信息科学研究的利器。本章将详细介绍使用Wolfram语言进行程序设计的方法，主要内容包括函数定义与应用方法、常用程序设计控制语句、模块定义与调用方法和一些程序设计实例。7.1函数Wolfram语言中的函数包括内置函数、包函数和自定义函数三种，其中，内置函数名以大写字母开头，可以直接调用；包函数是指存储在Wolfram资源库中的线上函数，使用时需动态地装入函数所在的包文件；自定义函数为用户编写的函数，可用于扩展内置函数的功能。Wolfram语言中，函数的调用方法非常灵活，设函数名为f，其参数为x，则基本的调用方法为f[x]，即使用中括号将参数包括其中；另一种更简洁的方法为f@x或x//f，均等价于f[x]。如果多个函数f、g和h嵌套调用，基本用法为f[g[h[x]]]，也可使用f@g@h@x或x//h//g//f。函数调用的另一种常见表示为“/@”，例如，f/@x，表示函数f作用于列表x的每一个元素，因此，f/@{1,2,3}等价于{f[1],f[2],f[3]}。7.1函数Wolfram语言内置了常用的科学计算相关的函数。这里以伪随机数相关的函数为例，介绍内置函数的用法。伪随机数相关的常用函数及其基本用法如表7.1所示。7.1.1内置函数7.1.1内置函数7.1.1内置函数7.1.1内置函数例7.1伪随机数函数应用举例（1）生成长度为10的0或1伪随机序列，并统计其中1的个数，见表7.2序号1。（2）生成长度为12的10至99间的伪随机整数序列，见表7.2序号2。（3）生成3×4的伪随机实数矩阵，每个元素的取值在1和2之间，见表7.2序号3。（4）生成位于直角坐标系中顶点为(0,0)、(1,0)、(0,3)和(1,3)的矩形内的伪随机复数序列，序列长度为20，见表7.2序号4。（5）随机生成一个小于100的素数，见表7.2序号5。（6）从1至100中随机选择10个数（有重复采样），见表7.2序号6。（7）从1至100中随机选择10个数（不重复采样），见表7.2序号7。（8）生成服从正态分布的10个伪随机数序列，见表7.2序号8。（9）给定伪随机数种子为20200705，生成长度为10的伪随机整数序列，每个元素在10至99间取值，然后，从小至大进行排序。见表7.2序号9。表7.2伪随机数函数典型应用实例表7.2中序号1至8的执行结果与时间有关，读者每次执行时得到的结果是不同的，所以，这里不给出这些语句的执行结果。但是表7.2中序号9的执行结果是不随时间变化的，执行后，变量a的结果一定为列表{81,67,94,51,22,95,65,17,92,41}，而排序后的结果为{17,22,41,51,65,67,81,92,94,95}。7.1函数7.1.1内置函数7.1函数Mathematica不但是全球最先进的科学计算软件，而且也是全球最先进的科学数据库，实时地记录着全球乃至可观测宇宙部分的动态信息，供科研人员参考和使用。Mathematica线上资源特别丰富，其中包含了大量的软件包，每个软件包以文件的形式存储在服务器上。这些软件包中的函数称为包函数，这里以有限域算术包为例介绍包函数的用法。在线上资源库中，有限域算术包的名称为FiniteField，在计算机联网的情况下，可在Notebook中调入有限域算术包，即输入

<<FiniteFields`或Needs["FiniteFields`"]然后，可以使用有限域算术包中的全部函数。这里重点介绍有限域算术中创建有限域对象的函数，如表7.3所示。7.1.2包函数7.1函数函数是Mathematica程序设计的基本要素。在Mathematica中，程序由函数组成，而函数一般以模块的形式供用户调用。模块将在第7.3节中介绍，这里重点介绍自定义函数及其用法。由于Mathematica系统内置函数和包函数均以大写字母开头，因此，自定义函数尽可能以小写字母开头。自定义函数的格式为：函数名[参数1_,参数2_,…,参数n_]=表达式。需要注意的是，在参数表中，每个参数后有一个下划线“_”，而表达式中的参数不再具有下划线。这里的参数以所谓的符号或变量的形式出现，其对应的实参可以为数值、列表或函数名。7.1.3自定义函数7.2控制语句程序语句的执行只有三种方式，即顺序、分支和循环。在Mathematica的Notebook中输入的计算语句按照顺序执行的方式执行。这里重点介绍分支和循环控制语句。7.2.1分支控制（1）If[条件，语句组1，语句组2，语句组3]如果条件为逻辑真，则执行语句组1；如果条件为逻辑假，则执行语句组2；如果条件在逻辑上非真非假，则执行语句组3。每个语句组可以包含多个语句，语句间用分号“;”分隔。If语句的其他形式简化形式有：(i)If[条件，语句组1]当条件为逻辑真时，执行语句组1；否则，无操作。(ii)If[条件，，语句组2]当条件为逻辑假时，执行语句组2；否则，无操作。注意，这里中间的逗号不可缺少。(iii)If[条件，语句组1，语句组2]当条件为逻辑真时，执行语句组1；当条件为逻辑假时，执行语句组2。7.2控制语句程序语句的执行只有三种方式，即顺序、分支和循环。在Mathematica的Notebook中输入的计算语句按照顺序执行的方式执行。这里重点介绍分支和循环控制语句。7.2.1分支控制（2）Switch[表达式，情况1，语句组1，情况2，语句组2，……]如果表达式为情况1，则执行语句组1；如果表达式为情况2，则执行语句组2。每个语句组可以有多个语句，使用分号“;”分隔。（3）Which[逻辑表达式1，值1，逻辑表达式2，值2，……]Which语句依次计算各个逻辑表达式的值，返回第一个为真的逻辑表达式i对应的值i。（4）条件控制符（“/;”）后面接一个逻辑表达式test当逻辑表达式test为真时，条件控制起作用。现在，使用分支控制语句实现下述分段函数，7.2控制语句Mathematica提供了三种循环控制语句，即Do语句、While语句和For语句。这三种语句的语法如表7.10所示。7.2.2循环控制7.2控制语句Wolfram语言是一种高级程序设计语言，具有高级循环控制语句，如Table、Array、Nest和NestList等。在大多数情况下，While、Do和For循环常被这些高级循环控制语句替代。这些高级循环控制语句的语法如表7.12所示。7.2.3高级循环控制7.2.3高级循环控制7.2.3高级循环控制7.2.3高级循环控制7.3纯函数纯函数又称纯匿名函数。纯函数的格式有两种，其一，“Function[函数体]”；其二，“函数体&”。这里“函数体”中，#表示形式参数；当有多个形式参数时，#1表示第一个形式参数，#2表示第二个形式参数，依此类推。在第一种情况下，还可以指定形式参数，如Function[形式参数列表，函数体]。下面列举几个典型语句，如表7.14所示。7.3纯函数纯函数一个重要的作用体现在表7.14序号7的应用中。Wolfram语言的基本数据结构为列表，而借助于纯函数可以把函数作用施加到列表的各个元素。注意，有些函数具有Listable属性（例如，Sin函数，使用“??Sin”可查看Sin函数的属性），这类函数作用于列表时，将作用于列表的每个元素。例如，Sin[{0,Pi/2,Pi,3Pi/2,2Pi}]将返回{0,1,0,-1,0}。7.4模块在Wolfram语言中，模块是程序的基本单位，模块对应着程序设计实现的子功能，类似于C语言中的函数。常用的模块有四种，即With、Block、Module和Compile。这里借助于Hénon映射生成伪随机序列为例，介绍模块的应用方法。Hénon映射是一种二维离散混沌，其吸引子的方程为

其中，a=1.4，b=0.3。给定初始值x0和y0，生成指定长度的伪随机序列，序列的每个元素为0至255间取值的整数。With语句的基本语法为With[初始化变量列表,语句组]7.4.1With模块7.4.1With模块其中，“语句组”可包含多个语句，各个语句间用分号“；”分隔。“语句组”中使用的变量可以为全局变量，也可以为“初始化变量列表”中的变量。如果“语句组”中使用了“初始化变量列表”中的变量，这些变量的初始值直接代入语句组中。特别注意：“初始化变量列表”中的变量均为With语句内部可见的局部变量。With语句的典型实例如表7.15所示。With语句的基本语法为With[初始化变量列表,语句组]7.4.1With模块With语句的处理速度比第7.4.3节将要介绍的Module语句快。With语句的典型用法如表7.15的序号3所示。注意：在With模块中出现的非“初始化变量列表”中的变量均为全局变量。借助于With语句使用Hénon映射生成伪随机序列的程序如下：henon:={#[[2]],1.0-1.4#[[2]]^2+0.3#[[1]]}&;f[x0_,y0_,m_]:=With[{x1=x0,y1=y0,m1=m},dat1=NestList[henon,{x1,y1},m1];dat2=Flatten[dat1][[3;;-1;;2]];dat3=Mod[Floor[(2+dat2)*10^6],256]]f[0.3,0.23,30]上述代码中，定义了两个函数henon和f，其中，henon使用纯函数定义；函数f具有三个参数x0、y0和m，x0和y0作为Hénon映射的迭代初值，m为生成的伪随机序列的长度。调用f[0.3,0.23,30]的执行结果如下：{240,4,204,52,14,106,206,79,232,189,113,50,255,186,203,163,254,170,182,123,147,207,33,81,174,249,115,33,211,113}需要特别注意的是，在函数f的定义中，dat1、dat2和dat3在函数f被调用后，将成为Notebook中的全局变量。但是大多数情况下，只希望函数返回值，并不产生任何全局变量。Block和Module模块可实现这种变量局部化的要求。7.4模块类似于With模块，Block模块可以直接使用全局变量，并可以生成全局变量。而全局变量的存在破坏了算法的模块化，因此，Block模块还实现了全部变量局部化方法。Block语句的语法为：

Block[局部变量列表,语句组]或Block[局部变量初始化列表,语句组]这里的“语句组”可包含多条语句，各条语句间用分号“；”分隔，语句组的最后一条语句的执行结果为返回值。Block语句的典型用法实例如表7.16所示。7.4.2Block模块借助于Block语句使用Hénon映射生成伪随机序列的程序如下：henon:={#[[2]],1.0-1.4#[[2]]^2+0.3#[[1]]}&;f[x0_,y0_,m_]:=Block[{x1=x0,y1=y0,m1=m,dat1,dat2,dat3},dat1=NestList[henon,{x1,y1},m1];dat2=Flatten[dat1][[3;;-1;;2]];dat3=Mod[Floor[(2+dat2)*10^6],256]]f[0.3,0.23,30]上述代码中，先定义了henon函数和f函数，在f函数中使用了Block语句，将Block语句中使用了变量均定义为局部变量，运行结果由“语句组”的最后一条语句返回。执行函数“f[0.3,0.23,30]”返回：{240,4,204,52,14,106,206,79,232,189,113,50,255,186,203,163,254,170,182,123,147,207,33,81,174,249,115,33,211,113}在使用Block模块时，如果“语句组”中的表达式包含了局部变量，将使用局部变量定义的初始值，这种方式称为动态计算“语句组”中的表达式（这是和Module模块唯一的区别，下文的Module模块是“静态”处理的）。例如：y=x^2+2x+1Block[{x=a},x+y]返回1+3a+a^2。而Module[{x=a},x+y]返回1+a+2x+x^2。7.4.2Block模块7.4模块在绝大多数情况下，使用Module语句实现模块的功能，最主要的原因，可能是因为“Module”英文有“模块”的含义。由7.4.2节可知，当全部使用局部变量时，Block模块与Module模块完全通用，事实上，Wolfram语言中Block模块的处理速度比Module模块更快。Module语句的语法如下：

Module[局部变量列表,语句组]或Module[局部变量初始化列表,语句组]这里的“语句组”可包含多条语句，各条语句间用分号“；”分隔，语句组的最后一条语句的执行结果为返回值。Module语句的典型用法实例如表7.17所示。7.4.3Module模块由Module语句使用Hénon映射生成伪随机序列的程序如下Clear["`*"]henon:={#[[2]],1.0-1.4#[[2]]^2+0.3#[[1]]}&;f[x0_,y0_,m_]:=Module[{x1=x0,y1=y0,m1=m,dat1,dat2,dat3},dat1=NestList[henon,{x1,y1},m1];dat2=Flatten[dat1][[3;;-1;;2]];dat3=Mod[Floor[(2+dat2)*10^6],256];{dat1[[2;;-1]],dat3}]{ps1,ps2}=f[0.7,0.13,3000]上述程序代码首先定义了henon函数和f函数，f函数的参数x0和y0为Hénon映射的迭代初值，m设定产生的混沌伪随机序列的长度。函数f的返回值为列表{dat1[[2;;-1]],dat3}，列表中有两个元素dat1[[2;;-1]]和dat3，依次为Hénon映射的状态值序列和伪随机序列。执行“{ps1,ps2}=f[0.7,0.13,3000]”获得长度为3000的状态序列和伪随机序列。下面使用ps1借助于ListPlot[Reverse/@ps1,AxesLabel->{"x","y"},AspectRatio->Automatic,ImageSize->Large]绘制Hénon映射的相图如下：建议在Module模块中全部使用局部变量，并将Module模块定义为函数的实现部分。由Module语句使用Hénon映射生成伪随机序列的程序如下7.4模块Mathematica软件使用C语言实现的，用户编写的Mathematica程序是顺序调用Mathematica系统函数实现的。一般地，用户程序的执行效率远远低于C语言可执行程序。为了提高Mathematica用户程序的执行效率，可以使用Compile模块，设计经过编译的用户函数。这种编译后的函数的执行效率非常接近于C语言可执行程序的效率，但是这类函数只能使用常规的数据类型，即整型（_Integer）、浮点型（_Real）、复数类型（_Complex）和逻辑变量（True或False），使用的列表必须为数值型数组。可见，虽然Compile模块可以编译为机器代码，极大地提高了执行速度，但是远没有Module模块灵活。7.4.4Compile模块Compile模块的典型语法为Compile[{{变量名1,变量类型1},{变量名2,变量类型2},…},语句组]或Compile[{{变量名1,变量类型1},{变量名2,变量类型2},…},语句组,属性]或Compile[{{变量名1,变量类型1,变量1维数},{变量名2,变量类型2,变量2维数},…},语句组]上述语法中，“语句组”可以包含多条语句，各条语句间用分号“；”隔开。Compile模块的典型用法实例如表7.18所示。由Complile语句使用Hénon映射生成伪随机序列的程序如下：Clear["`*"]henon:={#[[2]],1.0-1.4#[[2]]^2+0.3#[[1]]}&;f=Compile[{{x0,_Real},{y0,_Real},{m,_Integer}},Module[{x1=x0,y1=y0,m1=m,dat1,dat2,dat3},dat1=NestList[henon,{x1,y1},m1];dat2=Flatten[dat1][[3;;-1;;2]];dat3=Mod[Floor[(2+dat2)*10^6],256]]]ps=f[0.3,0.23,30]返回值为：{240,4,204,52,14,106,206,79,232,189,113,50,255,186,203,163,254,170,182,123,147,207,33,81,174,249,115,33,211,113}。上述程序中在Compile模块内部嵌入了Module模块。Compile模块的最大用处在于可以借助于Compile模块比较不同算法的运算速度，Compile模块本质上是机器语言程序，可以准确地反映算法的运行速度。本节以两个常用的对称密码算法（RC4和SM4）为例，介绍Mathematica程序设计方法。这里使用了Module模块作为程序设计基本单元。对于加密算法而言，输入为明文和密钥，输出为密文；对于解密算法而言，输入为密文和密钥，输出为明文。RC4是一种常用的流密码，而SM4是我国的一项文本加密标准。7.5程序设计实例7.5.1RC4加密原理与实现RC4密码，全称为“RivestCipher4”，是一种典型的分组密钥，习惯上称之为流密码，因为RC4可用于互联网中的实时数据传输。RC4的密钥长度可为1至256个字节，建议实际保密通信应用中使用128字节以上的密钥。这里，设p表示明文，k表示密钥，c表示密文，均为基于字节的向量。结合表7.6可知，对于RC4加密过程，输入为密钥k和长度为n个字节的明文p，输出为长度为n个字节的密文c。具体的加密步骤如下：（1）密码流初始化第1步：将密钥k扩展为长度为256字节的key。设密钥k的长度为m个字节，则key[i++]=k[(i++)modm],i=0,1,2,…,255(7.4)第2步：初始化长度为256字节的数组sbox，即sbox=[0,1,2,…,255]。第3步：循环变量i从0至255，循环执行以下两条语句：(i)j=(j+sbox[i]+key[i])mod256RC4加密过程如图7-6所示(ii)互换sbox[i]与sbox[j]的值。经过上述3步得到的sbox称为初始密码流。（2）加密算法已知明文p的长度为n。初始化变量i=0、j=0。变量u从0至n-1，循环执行以下语句：(i)i=(i++)mod256；(ii)j=(j+sbox[i])mod256；(iii)互换sbox[i]与sbox[j]的值；(iv)t=(sbox[i]+sbox[j])mod256；(v)c[u]=sbox[t]异或p[u]。最后得到的c即为密文。RC4加密过程如图7-6所示需要注意的是RC4密码不是一次一密算法，使用RC4密码的通信双方在“密码流初始化”之后，将随着图7-6中k的增加持续加密过程。RC4可能的不安全性在于密码流的重复（或循环再现）。因此，RC4密码不宜长期使用，在使用一段时间（加密了足够长的数据）后，应借助于公钥技术替换RC4密码的密钥k。此外，RC4不宜加密大量的重复性内容，这种情况下即使密码流是变化的，仍然有信息泄露的危险。RC4密码的解密过程与加密过程相似，除了有两点不同：（1）输入为密钥k和密文c，输出为还原后的明文p；（2）图7-6中有灰色填充的方框中的内容由原来的“c[u]=sbox[t]异或p[u]”变为“p[u]=sbox[t]异或c[u]”。现在，借助于Mathematica软件实现RC4算法。首先，设密钥k为71个字节长的字符串：“Donotforonerepulsegiveupthepurposethatyouresolvedtoeffect.”（仅作为示例，选自莎士比亚先生的名言，由于这个密钥中存在大量重复的字符（包括空格），这并不是一个优秀的密钥）。然后，设明文p取为莎士比亚先生的另一段名言：“Ignoranceisthecurseofnature,knowledgethewingwherewithweflytoheaven.”（长80个字节）。7.5程序设计实例7.5.2SM4加密原理与实现商用密码SM4是我国的一项文本数据加密标准，输入密钥长度为128比特，输入的明文长度为128比特，输出的密文长度也为128比特。如果用SM4加密长文本，需要使用CBC（密文分组链接）模式。SM4和DES（数据