第5章模块编程技术 Wolfram语言是函数式的语言,模块编程是自定义函数的主要方法。模块编程通过内置函数Module、Block、With或Compile实现,这些内置函数起到“容器”的作用,将其中的语句或语句组组织在一起,实现特定的自定义功能。除此之外,Module函数可定义局部变量,Module函数是定义自定义函数的常用“容器”,每次调用Module函数都将生成不同的局部变量; Block函数可以使全局变量的值局部化,即在Block函数中出现的全局变量的赋值操作不影响全局变量的值,例如,在笔记本中定义“x=5”,然后,执行Block函数“Block[{x, y}, y=x; x=10; {y, x}]”将得到“{10, 10}”,而在Block函数外部访问x的值,仍然为5。With函数在执行前用常量替换掉函数中相应的符号,例如,“With[{x=3}, x+1]”将返回4。一般地,由于With函数单纯地执行常量替换符号操作,所以With函数的执行速度最快; Block函数直接使用全局变量(或局部定义的变量)的名称,执行速度比Module函数快; Module函数每次执行都将创建新的局部变量,其执行速度最慢。但是,由于Module模块具有将变量局部化的特征,其应用最广泛。Compile函数将其中的语句组编译为可执行的机器代码,以提高函数的运行效率。 5.1Module模块 Module模块由Module函数实现。在Wolfram语言中,任一笔记本中定义的变量均为全局变量,在所有其他的笔记本中均可访问。为了使变量的作用域局部化,最常用的方法是借助于Module函数。在Module函数中,定义的变量均为局部变量,其作用域为整个Module函数。 5.1.1Module函数 Module函数的语法如下: (1) Module[{x, y, z, …}, 语句组],这里的“x, y, z, …”为定义的局部变量名,“语句组”为由分号“;”分隔的任意数量的语句,“语句组”的最后一条语句的执行结果为Module函数的返回值。 (2) Module[{x, y=y0, z, …}, 语句组],该语法与上述语法相似,该语法说明自定义的局部变量可以在定义时赋初始值。 一般地,使用Module函数实现自定义函数,例如计算一个数的绝对值,如图51所示。 图51使用Module函数自定义函数 使用Trace函数可以跟踪自定义函数的全部执行过程,如图52所示。 图52自定义函数abs的执行过程 在图52中,“In[4]”中的“Trace[abs[-3.14]]”将展示自定义函数的执行过程,如“Out[4]”所示。由“Out[4]”可知,“abs[-3.14]”的执行步骤如下: 第一步: “abs[-3.14]”,将函数调用发送到Mathematica内核。 第二步: “{NumericQ[-3.14],True}”,判断参数x(这里为-3.14)是否为数值,返回“True”。 第三步: “Module[{y=-3.14},If[y<0,-y,y]]”,解析函数体。 第四步: “{y$7371=-3.14,-3.14}”,在Module函数中重命名局部变量y为“y$7371”,这里的后缀“$7371”称为模块编号,该编号保存在一个系统全局变量“$ModuleNumber”中,随着模块的创建和调用其数值自动增加。Wolfram语言使用这种方法为Module函数分配专用的临时局部变量。这一步执行了自定义变量并初始化的语句“{y=x}”。 第五步: “{{{y$7371,-3.14}, -3.14<0, True}, If[True, -y$7371, y$7371], -y$7371, {y$7371, -3.14}, -(-3.14), 3.14}”,这一步首先解析语句“If[y<0,-y,y]”中的条件表达式,即将所有的局部变量和常数代入该条件表达式,判断结果为“True”; 然后,解析“If”函数,得到“-y$7371”; 最后,解析“-y$7371”得到结果“3.14”。 第六步: 返回结果“3.14”。 通过在上述第五步的解析过程可知,Wolfram语言使用回溯的方式,解析各个符号(变量),直到发现它们的数值。例如,解析“-y$7371”时,将再次解析“y$7371”。 Module函数中可以包括内置函数和自定义函数,如图53所示。 图53Module函数中可内嵌自定义函数 在图53中,“In[1]”使用Module函数定义了函数gcdlcm,该函数具有两个整型参数x和y,返回这两个参数的最大公约数和最小公倍数。在Module函数中,定义了局部变量z1、z2和符号lcm,其中,z1用于保存最大公约数,z2用于保存最小公倍数,符号lcm为自定义函数“lcm[a_,b_]:=a b/GCD[a,b]”,用于计算两个整数的最小公倍数,这里的“GCD”为计算最大公约数的内置函数。然后,语句“z1=GCD[x,y]”得到x和y的最大公约数; 语句“z2=lcm[x,y]”得到x和y的最小公倍数。最后,“{z1,z2}”作为Module函数的返回值,即为包含最大公约数z1和最小公倍数z2的列表。 在“In[2]”中调用“gcdlcm[36,64]”计算36和64的最大公约数和最小公倍数,其结果为“{4,576}”,如“Out[2]”所示。 Module函数中的局部变量可以赋初始值,如图54所示。 图54Module函数局部变量初始化 在图54中,“In[2]”使用Module函数定义了函数fact,具有一个整型参数x,计算x的阶乘。在Module函数中,定义了局部变量t和v,t用x初始化,v的初始值为1。由于参数x在自定义函数fact调用时将被赋予数值,所以,x在Module函数中以常数的形式存在,不能作为变量赋值,一般地,将参数x赋给一个局部变量,例如这里的“t”。在“While[t>1,v=v*t;t=t-1]”循环体中,计算t的阶乘,结果保存在局部变量v中。最后,“v”作为Module函数的输出。 在“In[3]”中,“{fact[5], fact[10], fact[20]}”依次计算了5、10和20的阶乘,如“Out[3]”所示。 5.1.2Module模块实例 这里借助Module函数解决一个实际问题,即房贷计算问题。假设某人为买一套住房向银行贷款100万元,年利率为4.5%,按月计算复利(月利率为0.045/12),计划15年还清全部贷款,且每月还款金额相同(按月等额还款方式),编程计算每月应还款多少元。这里使用循环搜索解的方式(不套用公式计算),开始时设置较大的步长,然后,设置小步长,精确度至少为0.01元。 这里设每月应还款x元,月利率为rate,设从银行发放贷款的下一个月开始还款,设剩余还款金额保存在变量pt中,于是,第0个月pt初始化为总贷款额度; 第1月后,pt = pt*(1+rate)-x; 第二个月后,pt = pt*(1+rate)-x(此处pt为第一个月后的pt),以此类推,直到pt为0,表示还款完成。按这种思路编写的程序如图55所示。 图55房贷计算问题 在图55中,定义了house函数,具有money、years和rateofyear三个参数,分别表示贷款总额、还款年限和年利率。在Module模块中,定义了局部变量“{pt, y=years, rate, x, n=108}”,这里pt表示剩余还款额度,y初始化为还款年限数,rate用于保存月利率,x为每月还款额,n设置为108,用作Do循环的最大循环次数。然后,执行语句“rate=rateofyear/12; pt=money; x=0;”设置月利率rate,设置剩余还款额pt为总贷款额,设置每月还款额x为0元。接着进入如下的Do循环: "Do[Table[pt=pt*(1+rate)-x,{i,1,y*12}]; If[pt>0,pt=money;x=x+100,Break[]], {j,1,n}];" 在上述Do循环中,首先使用给定的月还款额x,进行还款操作“Table[pt=pt*(1+rate)-x, {i,1,y*12}]”,如果执行后,“pt>0”,说明每月还款额x过少,于是“pt=money; x=x+100”,即令pt为总还款额,每月还款额x增加100元。在Do循环中重复这一过程,直到pt小于或等于0,说明每月还款额x已大于需要的真实每月还款额。由于步长为100元,所以,真实的每月还款额在区间[x-100,x]内。 在下一个Do循环前,先执行“x=x-100; pt=money;”,即将每月还款额x设置为刚好不够还款的额度(在步长为100元的情况下),并将待还款额pt设为总还款额money,再执行下面的Do循环: "Do[Table[pt=pt*(1+rate)-x,{i,1,y*12}]; If[pt>0,pt=money;x=x+10,Break[]],{j,1,n}];" 上述Do循环与前一个Do循环实现的功能类似,只是将每月还款额x的步长设为10,经过上述Do循环,真实的每月还款额将落在区间[x-10,x]内。 回到图55中,上述Do循环后面的Do循环操作与此类似,每月还款额x的步长取为1元,如下所示: "x=x-10;pt=money; Do[Table[pt=pt*(1+rate)-x,{i,1,y*12}]; If[pt>0,pt=money;x=x+1,Break[]],{j,1,n}];" 这里,首先将每月还款额x设置为刚好不够还款的额度(在步长为10元的情况下),并将待还款额pt设为总还款额money,接着执行Do循环直到pt小于或等于0,此时,真实的每月还款额将落在区间[x-1,x]内。 然后,进一步缩小步长至0.1元,执行如下代码: "x=x-1;pt=money; Do[Table[pt=pt*(1+rate)-x,{i,1,y*12}]; If[pt>0,pt=money;x=x+0.1,Break[]],{j,1,n}];" 上述代码的分析与前述Do循环类似,执行完该Do循环后,真实的每月还款额将落在区间[x-0.1,x]内。 之后,进一步缩小步长至0.01元,执行如下代码: "x=x-0.1;pt=money; Do[Table[pt=pt*(1+rate)-x,{i,1,y*12}]; If[pt>0,pt=money;x=x+0.01,Break[]],{j,1,n}];" 上述代码执行完成后,真实的每月还款额将落在区间[x-0.01,x]内。 接着,再次缩小步长至0.001元,执行如下代码: "x=x-0.01;pt=money; Do[Table[pt=pt*(1+rate)-x,{i,1,y*12}]; If[pt>0,pt=money;x=x+0.001,Break[]],{j,1,n}];" 上述代码执行完成后,真实的每月还款额将落在区间[x-0.001,x]内。 最后,将每月还款额的步长设为0.0001元,执行如下代码: "x=x-0.001;pt=money; Do[Table[pt=pt*(1+rate)-x,{i,1,y*12}]; If[pt>0,pt=money;x=x+0.0001,Break[]],{j,1,n}];" 执行完上述Do循环后,最后得到的每月真实还款额将位于区间[x-0.0001, x]内,而题目要求精确度至少为0.01元,因此,这时可取x为每月还款额。 在“In[2]”中调用“house[1000000,15,0.045]”,得到结果为“{7649.93, -0.0030425}”,即每月还款额应为7649.93元,如“Out[2]”所示。 为了方便介绍算法的实现过程,将图55中的自定义函数house代码写得比较冗长,图56中将house函数的代码作了简化整理。 图56简化后的house函数 对比图55,在图56自定义函数house的Module函数中,多定义了一个局部变量dx,用于保存每月还款额的步长,初始值为100元。在Module函数中,令“x=money/(12y)”,设置x的初始值为零利率的每月还款额。然后,使用两重Do循环,内层Do循环与图55中的每个Do循环类似,如下: "Do[pt=money; Table[pt=pt*(1+rate)-x,{i,1,y*12}]; If[pt>0,x=x+dx, Break[]], {j,1,n}];" 在这个Do循环中,首先将剩余还款额pt赋为总还款额money; 然后,使用当前的每月还款额x进行还款; 如果pt大于0,说明每月还款额x小于实际的每月应还款额,将x增加到x+dx; 循环执行这个操作,直到pt小于或等于0,说明实际的每月应还款额在区间[x-dx, x]上。 在外层的Do循环中,将x设为x-dx,即在步长为dx的情况下,刚好不够还款的每月还款额度; 将dx缩小10倍,即“dx=dx/10.0;”。然后,回到循环体的开头继续执行,直到剩余还款额pt的绝对值小于0.01元为止。最后,将x=x+dx作为实际的每月还款额。 在图56的“In[3]”中,调用“house[1000000,15,0.045]”,得到实际的每月还款额为7649.93元,如“Out[3]”所示。 图55和图56的程序具有通用性,输入任意的贷款额、贷款年限和年利率均可以计算出每月应还款额度。在图56的“In[4]”中,调用“house[2000000,20,0.050]”计算了贷款额200万元、贷款20年和年利率为5%时每月还款额,结果为13199.1元,如“Out[4]”所示。 图56的程序可以使用While循环实现,代码更加简洁,如图57所示。 图57使用While实现的house函数 对比图56可知,在图57中,使用While循环替代了Do循环,改进后无须关心Do循环的执行次数。 5.2Block模块 Block模块由Block函数实现。与Module函数类似,Block函数是变量局部化的重要方法。与Module函数不同之处在于,Block函数为每个变量(或符号)在Block函数内部分配临时的值,这些临时的值不影响每个变量(或符号)在Block函数外的取值; 而不是像Module函数那样创建新的局部变量。一定意义上,Block函数可以完全替代Module函数。Block函数的执行速度比Module函数更快,Block函数可用于所有使用Module函数的情况下; 此外,Block函数可以临时地调整系统全局变量的值,以达到这些全局变量作用局部化的效果。 5.2.1Block函数 Block函数的语法如下所示: (1) Block[{x, y, z, …}, 语句组],其中,“{x, y, z, …}”为局部变量列表,这些局部变量可以与笔记本中的全局变量同名; “语句组”为由分号“;”分隔的任意多条语句。 (2) Block[{x, y=y0, z, …}, 语句组],该语法说明局部变量列表中的各个局部变量可以赋初始值。 下面通过图58所示程序说明Block函数与Module函数的异同点。 图58Block函数和Module函数的异同点 图58中,在“In[2]”中定义了全局变量x和y,分别赋值为1和2。“In[3]”为Block函数,定义了局部变量x、u和v,其中,u初始化为10; 然后,调用“Print[x]”打印x,将得到“x”; 接着,对局部变量x赋值3,这不影响全局变量x; 之后,对全局变量y赋值5,令“v=u+x+y”,此时,v=10+3+5=18,最后“Print["v=",v]”打印v的值为“v=18”。在“In[4]”中显示“{x,y}”的值为“{1,5}”,即全局变量x不受同名的局部变量的影响,全局变量y在Block函数中重新赋值为5。 在“In[5]”中使用Module函数实现相同的功能,这里“Print[x]”得到了“x$7370”,即一个新的局部变量,命名规则为“局部变量名$模块编号”。在“In[6]”中显示“{x,y}”的值为“{1,5}”,即全局变量x不受Module函数中“同名”局部变量的影响,而全局变量y在Module函数中被赋值,这与Block函数实现的功能相同。 图59中的程序只能使用Block函数。 图59Block函数使系统变量作用域局部化 在图59中,“In[2]”使用Block函数计算精度为32的浮点数运算。这里只能使用Block函数,因为Module函数为每个局部变量创建新的变量形式。在Block函数中,设定了系统全局变量$MaxPrecision和$MinPrecision的新值,均为32,表示使用精度为32的浮点数运算,这并不影响系统全局变量$MaxPrecision和$MinPrecision的全局设定值。然后,计算了“Exp[0.84'32]”,这里的“'32”表示浮点数的精度为32,计算结果如“Out[2]”所示。在“In[3]”中函数“RealDigits”将“Out[2]”中的浮点数的各位提取出来,以列表的形式保存,如“Out[3]”所示,列表中的第1个元素“{2, 3, 1, 6, 3, 6, 6, 9, 7, 6, 7, 8, 1, 0, 9, 1, 7, 3, 3, 5, 0, 0, 2, 4, 7, 1, 9, 2, 8, 6, 5, 5}”为“Out[2]”各位上的数字,第2个元素“1”表示小数点的位置在第1个元素之后。在“In[4]”中“Length[First[%]]”显示了“Out[2]”中数位的个数为32,如“Out[4]”所示。 由图58和图59可知,Block函数可以取代Module函数,实现符号(或变量)的作用域局部化,只是Block函数为符号分配新的存储空间(而不改变符号的名称),而Module函数将为符号创建新的副本(用“符号名$模块编号”命名该副本),由于每次运行Module函数,其“模块编号”均增加,故每次调用Module函数为局部符号(或变量)创建的副本都不相同。但是Module函数却不能完全替代Block模块,特别是在如图59所示的情况下,只能使用Block函数,而不能使用Module函数。由于Block函数的运算速度快于Module函数,所以,应尽可能使用Block函数进行模块化编程。 5.2.2Block模块实例 假设在图510所示的边长为2的正方形(正方形的中心位于坐标原点)中随机撒入细沙粒,若沙粒在正方形中服从均匀分布,则位于单位圆中的沙粒数与落入正方形中的总沙粒数之比应为单位圆的面积与正方形的面积之比。通过这种方式可以近似计算圆周率,如图511所示。 图510蒙特卡实验 图511计算圆周率 在图511的“In[1]”中定义了函数pi,具有一个整型参数n,表示使用的沙粒的总数量。在Block函数中定义了两个局部变量sand和m,m初始化为0,用于记录落在单位圆内的沙粒个数。然后,调用RandomVariate内置函数生成在图510所示正方形内均匀分布的n个随机变量,赋给sand。接着,“Table[If[Norm[sand[[i]]]<1,m++],{i,n}]”统计落入单位圆中的沙粒数,赋给变量m。最后,计算“4m/n//N”得到圆周率的近似值。 在“In[2]”中,调用“{pi[2000], pi[20000], pi[200000]}”计算了当沙粒个数为2000、20000和200000时的圆周率的值,如“Out[2]”所示。可见,沙粒数越多,圆周率的计算结果越准确。 现在回到图57中,添加一条语句统计使用Module函数的house函数的执行时间,如图512所示。 图512Module函数实现的house函数执行时间 然后,将图57中的“Module”改为“Block”,即使用Block函数实现house函数,其余内容保持不变,如图513所示。 图513Block函数的house函数执行时间 在图513中,使用Block函数实现了自定义函数house,house函数的其余内容与图57中的内容相同。对于图512中“In[4]”和图513中“In[3]”的执行结果,可知,Block函数实现的house函数的执行时间为0.0149613秒,比Module函数实现的house函数的执行时间(0.0153025秒)短。当程序更加复杂时,Block函数将明显地快于Module函数。对比图57和图513可知,任意的Module函数均可直接将“Module”改为“Block”函数,程序仍然工作正常。 5.3With模块 With模块由With函数实现,With函数通过给局部变量赋初始值的方式使局部变量成为With函数中的“局部常数”。With函数体中首先执行这些局部变量的初始值替换,然后,对替换后的表达式进行计算,不需要创建新的局部变量存储空间,因此,With函数的执行效率比Block和Module函数更高。 5.3.1With函数 With函数的语法如下: With[{x=x0, y=y0, z=z0, …}, 语句组] 在With函数中,“{x=x0, y=y0, z=z0, …}”声明局部“变量”的同时必须赋初始值,而且,这里的局部“变量”被其赋予的值替代,而不再是通常意义上的“变量”,即不能再对它进行赋值操作。 With函数的基本用法如图514所示。 图514With函数基本用法 在图514中,“In[2]”用Trace函数跟踪了With函数的执行。由“Out[2]”可知,在解析了“With[{x=0},x+1]”之后,直接将x替换为常数0,最后返回1。因此,无法在With函数中对x进行赋值操作,因为With函数体中没有x了。可以将With函数理解为将“常量”局部化。 在“In[3]”中,使用With函数定义了函数sqr,这里的With函数直接将x替换为n,返回n2。在“In[4]”中调用“sqr[4]”,得到16,如“Out[4]”所示。 5.3.2With模块实例 With函数的语法为“With[{x=x0, y=y0, z=z0, …}, 语句组]”,其实现的功能为用局部“变量”列表中的初始值“常量”替换“语句组”中的局部“变量”,有些类似于“ReplaceAll”函数实现的功能。With函数的这种替换有两个主要的应用场合,一是生成复杂的纯函数,二是函数的替换。 图515展示了With函数的两个典型应用实例。 图515With函数的典型应用实例 在图515中,“In[2]”定义了函数f,具有a、b和c三个参数,在With函数中,令“y=a x2+b x+c”,代换纯函数中的y。在“In[3]”中,“{f[3,4,5], f[3,4,5][1], f[1,2,1][1]}”依次得到参数为3、4和5的纯函数“Function[x,5+4 x+3 x2]”和当这个纯函数作用在1上的函数值12以及参数为1、2和1的纯函数作用于1上的函数值4,如“Out[3]”所示。 在“In[4]”中,用With函数定义了函数g,具有两个参数fun和x,其中fun表示函数名,x表示“fun”的参数。在“In[5]”中调用了“{g[Sin,x], g[Sin,π/3], g[Cos,π/3]}”,结果如“Out[5]”所示,为“Sin[x],32,12”。 5.4Compile模块 Wolfram语言的内置函数均作了优化处理,具有与机器代码相当的执行效率。对于自定义函数,应尽可能调用内置函数实现所定义的功能,这样保证自定义函数具有较高的执行效率。对于执行数值计算的自定义函数,可以借助于Compile模块将其编译为机器代码执行。Wolfram语言提供了两种编译目标: ①编译为运行于Wolfram虚拟机上的机器代码; ②编译为C语言机器代码。Compile模块对自定义函数执行了编译处理,这种编译仅对整型、浮点型、复数类型和逻辑型(True和False)的数据有效,因此,仅能对一小部分内置函数进行编译,也就是说,只有使用整型、浮点型、复数类型和逻辑型的数据且包含了可以被编译的内置函数的自定义函数,才能被Compile模块编译。可以通过指令“Compile'CompilerFunctions[]”查看可以被编译的部分内置函数列表(一些内置函数可以被编译但是没有列于其中)。 5.4.1配置MinGW64编译器 这里使用MinGW64编译器实现自定义函数的代码编译。 登录网址http://mingww64.org/doku.php/download,如图516所示。 图516下载MinGWW64安装程序 在图516中,单击MingWW64builds,进入图517所示界面。 图517MinGWW64下载链接 在图517中单击“Sourceforge”,下载MinGWW64安装文件。下载后的文件名为mingww64install.exe,文件大小约为938KB。在Windows 10(64位)环境下运行文件mingww64install.exe,进入图518所示安装界面,其中,架构“Architecture”中选择“x86_64”,版本号为“8.1.0”。 在图518中单击“Next”按钮,进入图519所示窗口。 图519中默认安装目录为C:\Program Files\mingww64\x86_648.1.0posixsehrt_v6rev0。然后,单击“Next”进行安装(需联网)。安装完成后,MinGWW64所在的目录为“C:\Program Files\mingw64\mingw64”,目录和文件结构如图520所示。 图518安装MinGWW64 图519安装目录设置 图520MinGWW64目录结构 右击“我的电脑”,在弹出菜单中选择“属性”; 然后,进入“高级系统设置”; 在“高级”选项卡中,单击“环境变量(N)…”,在弹出的“系统变量(S)”界面“编辑”路径“Path”,在其列表的最后一行添加路径“C:\Program Files\mingw64\mingw64\bin”。 现在,在目录E:\ZYMaths\YongZhang\ZYCPrj中编写myhello.c文件,如图521所示。 图521myhello.c文件 打开“命令提示符”工作窗口,工作路径设为E:\ZYMaths\YongZhang\ZYCPrj,如图522所示。 图522编译myhello.c并执行 在图522中,调用了gcc和x86_64w64mingw32gcc编译了myhello.c文件,并分别生成了myhello1.exe和myhello2.exe文件。执行这两个可执行文件均显示“Hello World!”,说明MinGWW64安装成功。 现在,在计算机C盘根目录下建立子目录MinGWw64,然后,将图520中所示内容复制到目录C:\MinGWw64中。之后,编辑C:\ProgramData\Mathematica\Kernel目录下的init.m文件,设定其内容如下: Needs["CCompilerDriver'GenericCCompiler'"]; $CCompiler = {"Compiler"->GenericCCompiler, "CompilerInstallation"->"C:/MinGW-w64", "CompilerName"->"x86_64-w64-mingw32-gcc.exe"}; 文件init.m在Mathematica启动时自动被调用,将编译器配置为MinGWw64。 在Mathematica软件中,新建笔记本Y0513.nb,输入如图523所示的代码,在Block函数中指定编译目标为C语言可执行代码,即“$CompilationTarget="C"”,可实现对Compile函数代码的编译。 图523使用MinGWW64编译Block模块 在图523中,Compile函数有一个参数x,该参数将用作自定义函数ccp的参数,如“In[2]”所示,执行结果如“Out[2]”所示。此外,在Compile模块中使用选项“CompilationTarget->"C"”,可将Wolfram代码编译为C语言可执行代码,将在下一节中详细介绍Compile函数。当编译包含了已编译模块的代码时,需要配置编译选项为“CompilationOptions -> {"InlineExternalDefinitions" -> True}”,如图524所示。 图524Compile编译配置 在图524中“In[2]”所示的myf1函数中,使用选项“CompilationTarget->"C"”将其编译为C语言可执行代码; 在“In[3]”所示的myf2函数中,调用了myf1函数,故使用了选项“CompilationOptions -> {"InlineExternalDefinitions" -> True}”。在“In[4]”中调用“myf2[{2.3,1.6},{0.1,0.3}]”得到执行结果0.610677,如“Out[4]”所示。下一节将深入介绍Compile函数。 5.4.2Compile函数 Compile函数用于生成编译执行的代码模块,以提高代码的执行效率。目前Wolfram语言仅支持整型、实型、复数类型和逻辑型的数据处理方式的代码的编译。由于Wolfram语言是函数式的语言,而自定义函数必将调用内置函数和其他的自定义函数,所以,自定义函数的代码编译实质上是将其调用的内置函数和其他的自定义函数在可执行代码的级别上建立了接口。通过这种方式提高编译后的自定义函数(即Compile函数)的代码执行效率。 Compile函数的调用语法如下: (1) Compile[{x1,x2,…}, 语句组] 这里的“{x1,x2,…}”为Compile函数定义的函数的参数,默认为实数类型,“语句组”可以为Module函数或Block函数,也可以为由分号“;”分隔的大量语句组成,如图525所示。 图525Compile函数第一种语法实例 在图525的“In[2]”中,使用Compile函数定义了函数f1,具有x和y两个参数,默认为实型参数。Compile函数体由Block函数组成,这里Block函数实现了参数x与y相乘,积s作为返回值,也是Compile函数的返回值。而“CompilationTarget->"C"”为编译选项,如果省略该选项,相当于“CompilationTarget->"WVM"”,即编译为运行于Wolfram虚拟机上的机器代码。编译成功的函数如“Out[2]”所示,为“CompiledFunction”函数,显示的“Argument count: 2”表示具有2个参数,“Argument types: {_Real, _Real}”表示2个参数均为实数类型,单击“”将显示更多的编译后的函数信息。 (2) Compile[{{x1, 类型}, {x2, 类型}, …}, 语句组] 这里的“{{x1, 类型}, {x2, 类型}, …}”为Compile函数定义的函数的参数及其类型声明,注意: 这里的“类型”仅支持整型、实型、复数类型和逻辑型(True或False,逻辑型用“True|False”显式指定)。当这里的“类型”为整数、浮点数类型和复数类型时,分别用“_Integer”“_Real”和“_Complex”指定,表示为单个的数值参数。这里的“语句组”可以为Module函数或Block函数,也可以为由分号“;”分隔的任意多条语句。 在Compile函数中可以使用笔记本中定义的全局变量,但是一般情况下,不应在Compile函数中使用全局变量,而应全部使用局部变量,因此,Compile函数体常由Module或Block函数组成。Compile函数中也可以使用With函数,在使用With函数时需注意,With函数并不能定义局部变量,只能其将定义的局部“变量”替换为“常量”或“函数”符号进行后续计算,With函数用于Compile函数中,一般只用作函数替换使用。 Compile函数的这种语法实例如图526所示。 图526Compile函数的第二种语法实例 在图526中,“In[4]”使用Compile函数定义了一个函数f2,具有x、y和z三个参数,依次为整型、实型和复数类型的参数。在Module函数内部,使用了自定义函数f1,然后,添加了选项“CompilationOptions->{"InlineExternalDefinitions" -> True}”,表示本函数使用了编译了的自定义函数。自定义函数f2实现了x*y+|z|的计算。在“In[5]”中调用了“f2[3,5,2+3I]”,得到结果“18.6056”,如“Out[5]”所示。 (3) Compile[{{x1, 类型, 维度}, {x2, 类型, 维度}, …}, 语句组] Compile函数的这种语法比第(2)种语法多了一个参数的“维度”,如图527所示。 图527Compile函数的第三种用法 在图527中,“In[6]”使用Compile函数的第(3)种语法定义了函数f3,具有x、y和z三个参数,这里的参数x为一维列表,参数y为二维列表,参数z为三维列表。在Module函数内,“s=Total[x]”将列表x的元素之和赋给s,“t=Det[y]”将二维列表y(必须输入方阵)的行列式赋给t,“v=Flatten[z]”将三维列表z压平为一维列表,赋给v。“s*t+Total[v]”作为Module函数的返回结果,也是Compile函数的返回结果。在“In[7]”中,执行调用“f3[{1,2,3}, {{1,2},{3,4}}, {{{3},{4}},{{5},{6}}}]”,这里,第一个参数“{1,2,3}”为一维列表,第二个参数“{{1,2}, {3,4}}”为一个二维列表,第三个参数“{{{3},{4}},{{5},{6}}}”为一个三维列表,计算结果为“6.”,如“Out[7]”所示。 5.4.3Compile模块实例 本节将介绍两个使用Compile函数实现的实例,其一为生成Logistic混沌序列,其二为使用RC4进行数据流加密与解密。 实例一Logistic混沌序列发生器 这里,使用的离散Logistic映射的形式为xn+1=1-2,状态的取值范围为区间(-1,1),将状态序列转化为0~255的整数序列,程序如图528所示。 图528Logistic混沌序列发生器 在图528中,“In[2]”使用Compile函数定义了函数logistic,具有两个参数x0和n,分别表示Logistic映射的初值和生成的序列的长度,这里的“x0”和“n”在函数内部应视为常数,不能作为变量使用。在Module模块中,定义了局部变量x1(赋初值x0)、x2和dat,这里的“x1”和“x2”用于表示Logistic映射的两个状态,“dat”用于存储生成的伪随机序列。然后,“dat=ConstantArray[0,n]”将dat变量初始化为长度为n且元素值为0的一维列表。接着,在For循环中,循环变量i从1至n,循环执行语句组“x2=(1.0-2.0#2)&[x1]; dat[[i]]=Mod[Floor[x2*1010], 256]; x1=x2”n次,每次执行先由状态x1迭代得到状态x2,然后,由状态x2得到序列的第i个值dat[[i]],再将x2赋给x1进行下一次迭代。这里的Logistic映射使用了纯函数的形式“(1.0-2.0#2)&”。 在“In[3]”中,调用了logistic自定义函数“logistic[0.321,10]”,设置初始值参数为0321,序列长度参数为10,得到结果“{224,25,72,241,213,209,183,110,156,67}”,如“Out[3]”所示。 实例二RC4加密与解密数据流 RC4密码,全称为“Rivest Cipher 4”,是一种典型的分组密钥,习惯上称之为流密码,因为RC4可用于互联网中的实时数据流传输。RC4的密钥长度可为1~256B,建议实际保密通信应用中使用128字节以上的密钥。 这里,设p表示明文,k表示密钥,c表示密文,均为基于字节的向量。RC4加密过程如图529所示。 图529RC4加密过程 结合图529可知,对于RC4加密过程,输入为密钥k和长度为n个字节的明文p,输出为长度为n个字节的密文c。具体的加密步骤如下: 1. 密码流初始化 第1步: 将密钥k扩展为长度为256字节的key。设密钥k的长度为m个字节,则 key[i++]=k[(i++) mod m], i=0,1,2,…,255 第2步: 初始化长度为256字节的数组sbox,即sbox=[0,1,2,…,255]。 第3步: 循环变量i从0至255,循环执行以下两条语句: ① j=(j+sbox[i]+key[i]) mod 256; ② 互换sbox[i]与sbox[j]的值。 经过上述3步得到的sbox称为初始密码流。 2. 加密算法 已知明文p的长度为n。初始化变量i=0、j=0。变量u为0~n-1,循环执行以下语句: ① j=(j+sbox[i]) mod 256; ② 互换sbox[i]与sbox[j]的值; ③ t=(sbox[i]+sbox[j]) mod 256; ④ i=(i++) mod 256; ⑤ c[u]=sbox[t] 异或p[u]。 最后得到的c即为密文。 需要注意的是,RC4密码不是一次一密算法,使用RC4密码的通信双方在“密码流初始化”之后,将随着图529中循环变量u的增加持续加密过程。RC4可能的不安全性在于密码流的重复(或循环再现)。因此,RC4密码不宜长期使用,在使用一段时间(加密了足够长的数据)后,应借助于公钥技术替换RC4密码的密钥k。此外,RC4不宜加密大量的重复性内容,这种情况下即使密码流是变化的,仍然有信息泄露的危险。 RC4密码的解密过程与加密过程相似,但有两点不同: ①输入为密钥k和密文c,输出为还原后的明文p; ②图529中有灰色填充的方框中的内容由原来的“c[u]=sbox[t] 异或p[u]”变为“p[u]=sbox[t] 异或c[u]”。 借助于Compile函数实现的RC4密码(密钥为k、明文为p、密文为c)自定义函数如图530和图531所示。 图530密码流初始化 图531RC4密码 在图530中,“In[2]”用Compile函数定义了函数stream,该函数具有一个输入参数k,表示密钥,为在0~255取值的整数序列(一维列表),如“In[3]”所示。在“In[3]”使用了20个字节表示的整数序列作为密钥k,即密钥k的长度为160比特。 回到图530的“In[2]”,在Compile函数内部的Module函数中,定义了局部变量key保存密钥k扩展至256字节后的密钥,用语句“key=PadRight[k,256,k]”实现,这里的“PadRight”函数将密钥k向右填充为长度为256的列表,使用密钥k填充,相当于把密钥k循环扩展为长度为256的序列。Module函数还定义了局变量j,并初始化为0; 定义了局部变量sbox,在“sbox=Range[0,255]”,将sbox初始为列表{0, 1, …, 255}。在For循环“For[i=1, i<=256, i++, j=Mod[sbox[[i]] + key[[i]] + j, 256]; {sbox[[i]], sbox[[j+1]]}={sbox[[j+1]], sbox[[i]]};]”中,根据密钥key的值,打乱sbox中各个元素的顺序。然后,sbox作为stream函数的输出。 在图530中的“In[4]”中,输入密钥k,由“(sbox=stream[k])//Short”得到加密用的sbox,这里“Short”函数表示仅显示sbox的部分数据,sbox为一个长度为256的列表,元素为0~255共256个整数的特定的乱序排列。 现在,由图530过渡到图531。在图531中,由Compile函数定义了函数rc4,具有两个参数,一个为表示明文序列的p,另一个为表示密钥的k。在Module函数内部,定义了6个局部变量,其中,i和j均初始化为0; n存储明文序列的长度,即“n=Length[p]”; sbox由密钥k经图530的函数stream得到,即“sbox=stream[k]”; c保存密文,初始化为元素均为0的列表,即“c=ConstantArray[0,n]”,密文c的长度与明文p的长度相同。在Table函数中实现对明文序列p的加密,其中,u为循环变量(1~n),循环执行以下操作: (1) j=Mod[j+sbox[[i+1]], 256],根据i的值借助于sbox更新j的值; (2) {sbox[[i+1]], sbox[[j+1]]}={sbox[[j+1]], sbox[[i+1]]},交换sbox的第i+1个和第j+1个元素; (3) t=Mod[sbox[[i+1]]+sbox[[j+1]], 256],由sbox的第i+1个和第j+1个元素之和得到临时变量t的值; (4) i=Mod[i+1, 256],更新i的值; (5) c[[u]]=BitXor[sbox[[t+1]], p[[u]]],由sbox的第t+1个元素与明文p的第u个元素相异或得到密文c的第u个元素。这些操作对应着RC4密码的加密过程。 最后,Table函数的输出c为Compile函数的输出,c即为加密后的密文。 回到图530,密钥k为“{110, 110, 218, 136, 1, 54, 119, 10, 174, 198, 25, 79, 82, 226, 99, 3, 171, 173, 49, 147}”,在图531中,明文序列p为“{131, 116, 71, 92, 107, 199, 147, 167, 224, 173, 48, 221, 123, 202, 184, 200, 229, 117, 9, 237, 77, 225, 16, 144, 230, 170, 7, 247, 92, 117}”,如“In[6]”所示。在“In[7]”中调用“c=rc4[p,k]”,使用密钥k对明文p进行加密,得到密文c为“{50, 147, 164, 37, 121, 153, 154, 153, 248, 165, 94, 104, 125, 222, 76, 80, 244, 126, 22, 135, 19, 109, 240, 67, 45, 114, 89, 228, 9, 169}”,如“Out[7]”所示。 在图531的“In[8]”中,调用“r=rc4[c,k]”,使用密钥k对刚生成的密文序列c进行解密,解密后的文本r为“{131, 116, 71, 92, 107, 199, 147, 167, 224, 173, 48, 221, 123, 202, 184, 200, 229, 117, 9, 237, 77, 225, 16, 144, 230, 170, 7, 247, 92, 117}”,如“Out[8]”所示,与原始的明文序列p完全相同。这是因为,上述RC4密码的加密函数和解密函数是相同的。 由于在图531的rc4函数中调用了stream函数,故使用了选项“CompilationOptions-> {"InlineExternalDefinitions" -> True}”。 5.5并行编程 对于复杂的算法,设计其并行计算处理算法是件非常困难的事情,需要考虑算法执行过程中可能出现的诸多因素。在Wolfram语言中,并行计算由Wolfram内置函数Parallelize管理,这是一个为并行计算高效分配算法资源的函数。内置函数放在Parallelize函数中,Wolfram语言将自动进行并行化处理,以尽可能多的并行指令执行这些内置函数。 在Mathematica软件的笔记本中,选择菜单“Edit|Preferences…”,在弹出的窗口“Preferences”中选择选项卡“Parallel”,并在该页面选择“Local Kernels”页面,在这个页面可以设定执行并行计算的内核数。可设定的内核数受计算机的CPU内核总数以及Mathematica的版权限制,最多可支持16个并行内核。设定了并行内核个数后,下面介绍并行编程相关的内置函数,首先介绍并行计算函数Parallelize函数,然后介绍并行处理函数ParallelTable、ParallelMap、ParallelDo、ParallellSum、ParallelProduct和ParallelArray函数等。 5.5.1并行计算函数 Parallelize函数的语法为“Parallelize[语句组]”,自动使用并行计算方式计算“语句组”,Parallelize有一个Method选项,常用的选项为“Method->"FinestGrained"”和“Method-> "CoarsestGrained"”,分别表示将计算任务分成尽可能小的计算子单元(细粒度划分)和将计算任务按计算机并行内核的数量进行分隔(粗粒度划分)。 这里用Parallelize函数计算乘方xy,使用“平方—乘算法”。“平方—乘算法”实现乘方xy的方法为: 将y转化为二进制形式,例如y为10101110b,则xy为x10101110b。令s=1,从y的二进制序列的最左边为1的位开始向右遍历,当某位为1时,s的平方乘以x赋给s,即s=s2·x; 当某位为0时,s的平方赋给s,即s=s2。这一过程本质上由y的二进制序列构造y的十进制数的值,因此,“平方—乘算法”是正确的。 Parallelize函数的典型实例如图532所示。 图532Parallelize函数典型实例 图532中,“In[2]”为实现“平方—乘算法”的自定义函数sqmul,具有两个整型参数x和y,计算xy。在Module模块中,定义了两个局部变量h和s,h用于保存y的二进制数的各个数位,“h=IntegerDigits[y, 2]”; s用于保存最后的结果,初始化为1。Table函数实现“平方—乘算法”,“Table[s=s2; If[t==1, s=s*x], {t, h}]”,循环变量t在列表h中取值,当t为0时,s的平方赋给s; 当t为1时,s的平方乘以x赋给s。最后,s作为函数sqmul的返回值。在“In[4]”中将“{sqmul[2,20], sqmul[12,19]}”作为函数Parallelize的参数,执行结果如“Out[4]”所示。“In[5]”和“Out[5]”为Wolfram语言计算220和1219的结果,以供参考对比。 需要注意的是,并非使用了并行计算函数后,计算任务的执行时间一定会大幅度减少。并行计算需要将原来的计算任务划分为可并行执行的计算子单元,并需要考虑各个计算子单元间的数据通信,对于并不复杂的运算任务而言,这些并行预处理工作花费的时间可能比计算任务本身的执行时间更长。并行计算主要应用于非常复杂和耗时的计算任务中。 5.5.2并行处理函数 并行处理函数都有对应的单线程函数,例如,ParallelTable函数与单线程函数Table函数对应,其语法也类似; ParallelMap函数与单线程函数Map对应,其语法也类似。但是并行处理函数将其中的变量局部化,如果读取其中的变量,必须使用SetSharedVariable函数使这些变量全局化。 这里主要介绍常用的并行处理函数,其中,图533演示了ParallelTable和ParallelDo函数的典型用法。 图533并行处理函数ParallelTable和ParallelDo 在图533中,“In[2]”调用Table函数“Table[Pause[1];i^2,{i,4}]//AbsoluteTiming”并统计其执行的时间,这里“Pause[1]”为延时1秒,由于Table函数是单线程顺序执行语句组“Pause[1]; i2”4次(循环变量i为1~4),故运行时间至少为4秒,结果“Out[2]”显示运行时间为4.06139秒。而“In[4]”调用并行处理函数ParallelTable函数“ParallelTable[Pause[1]; i2, {i,4}]//AbsoluteTiming”,这里并行执行语句组“Pause[1]; i2”,执行时间如“Out[4]”所示,为1.02861秒。除了将变量局部化外,ParallelTable与Table的作用相同,在“In[5]”中,“ParallelTable[1/(i+j-1),{i,3},{j,3}]”计算了3阶Hilbert矩阵,如“Out[5]”所示。 在图533中,“In[6]”定义了全局变量m,为3×3的全0矩阵; 然后,调用函数“SetSharedVariable[m]”将m设为并行处理共享变量,接着,调用并行处理函数ParallelDo函数“ParallelDo[m[[i,j]]=1/(i+j-1),{i,3},{j,3}]”计算3阶Hilbert矩阵,最后,以矩阵形式输出m,如“Out[9]”所示。 使用并行处理函数时需要注意,如果程序本身是顺序执行的方式设计的,即前面语句组的执行结果将影响其后的语句组,此时使用并行处理函数可能得不到正确的结果,如图534所示。 图534并行处理函数的注意事项 在图534中,“In[6]”拟使用并行ParallelTable函数计算12+22+…+102,但是这个函数对于每个循环变量i的取值作并行处理,最后得到的结果s为i=10和s初始值为0的情况下的值,即100,如“Out[9]”所示。这不是设计的算法的正确结果。这时应使用“In[10]”的并行求和函数ParallelSum,即“ParallelSum[i2,{i,10}]”,计算结果为385,如“Out[10]”所示。此外,Wolfram语言还有并行乘法函数ParallelProduct,在“In[11]”中计算了10的阶乘,即“ParallelProduct[i,{i,10}]”,结果为3628800,如“Out[11]”所示。 最后需要介绍的三个常用并行处理函数为ParallelMap、ParallelArray和ParallelEvaluate,如图535所示。 图535并行处理函数ParallelMap、ParallelArray和ParallelEvaluate 在图535中,“In[3]”使用平行处理函数ParallelMap将纯函数“(#+1/#)&”作用于列表“Range[5]”(即“{1,2,3,4,5}”上,得到结果如“Out[3]”所示。ParallelMap函数是Map函数的并行版本,同样可以作用于列表的不同层上,例如“In[4]”中,ParallelMap函数将纯函数“(1+#2)&”作用于列表“{{1,2,3}}”第2层上的各个元素上,得到结果如“Out[4]”所示。 ParallelArray函数是Array函数的并行版本,图535中的“In[5]”使用ParallelArray函数将纯函数“(#+1/#)&”作用于{1,2,3,4,5}上,输出结果如“Out[5]”所示,与“In[3]”的计算结果相同。“In[6]”使用ParallelArray函数将纯函数“(#12+2#2)&”作用于{i,j}上,这里i={1,2},j={1,2,3},得到二维列表如“Out[6]”所示。 ParallelEvaluate函数是单线程Evaluate函数的并行版本,“ParallelEvaluate[表达式]”将使用所有的并行内核计算“表达式”的值,而且各个内核之间互不相关。在图535的“In[7]”中,“ParallelEvaluate[RandomInteger[{1,10}]]”使用所有的并行内核执行函数“RandomInteger[1, 10]”得到1~10的伪随机整数,如“Out[7]”所示,这里共16个内核,得到一个长度为16的伪随机整数列表。 本章小结 Wolfram语言有6000多个内置函数,涵盖了数学、物理、化学、生物和计算机科学等众多领域的常用计算方法,熟练掌握并灵活应用这些内置函数是精通Mathematica软件的关键,而有效地组织内置函数的方法是借助于模块编程技术。本章详细介绍了Wolfram语言的四种模块化编程方法,即Module模块、Block模块、With模块和Compile模块。Module模块可以定义局部变量并为局部变量赋初值,并可以组织任意多的内置函数共同完成特定的计算功能,是最常用的模块化编程手段。与Module模块类似,Block模块也可以定义局部变量并为局部变量赋初值,但是Block模块还可以将全局变量的值局部化,即给全局变量赋予新的值,并使这个值的作用域为整个Block模块,而Block模块外部,全局变量的值不受影响。一般地,认为Block模块可以替代Module模块,也就是说,所有的Module函数均可以直接变换为Block函数,其计算结果仍然正确且执行效率更高。Compile函数是针对数值计算的情况下,将Wolfram语言函数编译为机器代码以提高算法的执行效率。这种优化本质上为Compile函数中使用的内置函数提供了机器代码级别的接口,可以使用Wolfram语言自带的编译器生成Wolfram语言虚拟机上执行的机器代码,也可以借助于外部编译器(如Visual Studio或MinGW64)等对Compile函数进行编译优化。有时需要对比两个算法的运算速度,此时,借助于单线程的Compile模块可以相对公平地评测算法速度性能。