注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐

SylvesterMatrix的Mathematica实现  

2012-04-05 14:49:33|  分类: 学习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

即以此功德,庄严佛净土。上报四重恩,下济三途苦。惟愿见闻者,悉发菩提心。在世富贵全,往生极乐国。

缘起:读懂Sylvester Matrix的Mathematica实现,验证SylvesterMatrix的行列式就是两个多项式的Resultant。

1) SylvesterMatrix 的实现

SylvesterMatrix1[poly1_, poly2_, var_] :=
 Function[{coeffs1, coeffs2},
   With[{l1 = Length[coeffs1], l2 = Length[coeffs2]},
    Join[NestList[RotateRight, PadRight[coeffs1, l1 + l2 - 2],
      l2 - 2],
     NestList[RotateRight, PadRight[coeffs2, l1 + l2 - 2], l1 - 2]]]][
  Reverse[CoefficientList[poly1, var]],
  Reverse[CoefficientList[poly2, var]]]

执行

SylvesterMatrix1[a0 + a1 x + a2 x^2 + a3 x^3,
  b0 + b1 x + b2 x^2 + b3 x^3, x] // MatrixForm

返回

{{a3, a2, a1, a0, 0, 0}, {0, a3, a2, a1, a0, 0}, {0, 0, a3, a2, a1,
  a0}, {b3, b2, b1, b0, 0, 0}, {0, b3, b2, b1, b0, 0}, {0, 0, b3, b2,
  b1, b0}}

2)验证Sylvester Matrix的行列式就是Resultant。依次执行如下语句,返回 True

SylvesterMatrix1[poly1_, poly2_, var_] :=
 Function[{coeffs1, coeffs2},  (*以前没看出来这两个参数就是形参*)
   With[{l1 = Length[coeffs1], l2 = Length[coeffs2]},
    Join[NestList[RotateRight, PadRight[coeffs1, l1 + l2 - 2],
      l2 - 2],  (*以前没看出来这里加粗部分是纯函数主体*)
     NestList[RotateRight, PadRight[coeffs2, l1 + l2 - 2], l1 - 2]]]
][
  Reverse[CoefficientList[poly1, var]],  (* 以前没看出来这两个参数就是实参 *)
  Reverse[CoefficientList[poly2, var]]]

SylvesterMatrix1[a0 + a1 x + a2 x^2 + a3 x^3,
  b0 + b1 x + b2 x^2 + b3 x^3, x] // MatrixForm

Det[%]

Resultant[a0 + a1*x + a2*x^2 + a3*x^3, b0 + b1*x + b2*x^2 + b3*x^3,
  x] === %3

3)SylvesterMatrix1函数的编写困扰了我整整1天半,原因是Pure Function的概念没有彻底搞明白。原来Function[{u, v}, u^2 + v^4][x, y]中{u,v}是形参,而[x, y]就是实参。再如g[#, #^2] & /@ {x, y, z} 返回 {g[x, x^2], g[y, y^2], g[z, z^2]},我不知道/@后面就是实参列表。

  1. 执行 SylvesterMatrix1[a0 + a1 x + a2 x^2 + a3 x^3, b0 + b1 x + b2 x^2 + b3 x^3, x]
    两个实参Reverse[CoefficientList[poly1, var]], Reverse[CoefficientList[poly2, var]]分别传递给Pure Function的形参{coeffs1, coeffs2}
  2. With[ ]中给l1和l2赋值,PadRight[coeffs1, l1 + l2 - 2]的意思是将coeffs1这个list扩充到长度为l1+l2-2的list,其右边用0填充。第一个多项式有l1个系数,为l1-1次多项式,第二个多项式有l2个系数,为l2-1次多项式。SylvesterMatrix是一个l1-1+l2-1的方阵。
  3. 对于NestList[RotateRight, PadRight[coeffs1, l1 + l2 - 2], l2 - 2]就是将PadRight[coeffs1, l1 + l2 - 2]的结果RotateRight执行l2-2次,因为coeffs1的长度本身就是l1,注意NestList的第一个结果是移动0次(即没有执行操作)的结果;NestList[RotateRight, PadRight[coeffs2, l1 + l2 - 2], l1 - 2]]就是将{b3, b2, b1, b0, 0, 0}循环右移l1-2次直到b0碰到矩阵的右边界
  4. Join[] 的结果为{{a3,a2,a1,a0,0,0},{0,a3,a2,a1,a0,0},{0,0,a3,a2,a1,a0},{b3,b2,b1,b0,0,0},{0,b3,b2,b1,b0,0},{0,0,b3,b2,b1,b0}}

4)执行Trace[SylvesterMatrix1[a0 + a1 x + a2 x^2 + a3 x^3,
   b0 + b1 x + b2 x^2 + b3 x^3, x]]

查看函数执行过程,返回

{SylvesterMatrix1[a0+a1 x+a2 x^2+a3 x^3,b0+b1 x+b2 x^2+b3 x^3,x],Function[{coeffs1,coeffs2},With[{l1=Length[coeffs1],l2=Length[coeffs2]},Join[NestList[RotateRight,PadRight[coeffs1,l1+l2-2],l2-2],NestList[RotateRight,PadRight[coeffs2,l1+l2-2],l1-2]]]][Reverse[CoefficientList[a0+a1 x+a2 x^2+a3 x^3,x]],Reverse[CoefficientList[b0+b1 x+b2 x^2+b3 x^3,x]]],{{CoefficientList[a0+a1 x+a2 x^2+a3 x^3,x],{a0,a1,a2,a3}},Reverse[{a0,a1,a2,a3}],{a3,a2,a1,a0}},{{CoefficientList[b0+b1 x+b2 x^2+b3 x^3,x],{b0,b1,b2,b3}},Reverse[{b0,b1,b2,b3}],{b3,b2,b1,b0}},Function[{coeffs1,coeffs2},With[{l1=Length[coeffs1],l2=Length[coeffs2]},Join[NestList[RotateRight,PadRight[coeffs1,l1+l2-2],l2-2],NestList[RotateRight,PadRight[coeffs2,l1+l2-2],l1-2]]]][{a3,a2,a1,a0},{b3,b2,b1,b0}],With[{l1$=Length[{a3,a2,a1,a0}],l2$=Length[{b3,b2,b1,b0}]},Join[NestList[RotateRight,PadRight[{a3,a2,a1,a0},l1$+l2$-2],l2$-2],NestList[RotateRight,PadRight[{b3,b2,b1,b0},l1$+l2$-2],l1$-2]]],{Length[{a3,a2,a1,a0}],4},{Length[{b3,b2,b1,b0}],4},Join[NestList[RotateRight,PadRight[{a3,a2,a1,a0},4+4-2],4-2],NestList[RotateRight,PadRight[{b3,b2,b1,b0},4+4-2],4-2]],{{{4+4-2,6},PadRight[{a3,a2,a1,a0},6],{a3,a2,a1,a0,0,0}},{4-2,2},NestList[RotateRight,{a3,a2,a1,a0,0,0},2],{RotateRight[{a3,a2,a1,a0,0,0}],{0,a3,a2,a1,a0,0}},{RotateRight[{0,a3,a2,a1,a0,0}],{0,0,a3,a2,a1,a0}},{{a3,a2,a1,a0,0,0},{0,a3,a2,a1,a0,0},{0,0,a3,a2,a1,a0}}},{{{4+4-2,6},PadRight[{b3,b2,b1,b0},6],{b3,b2,b1,b0,0,0}},{4-2,2},NestList[RotateRight,{b3,b2,b1,b0,0,0},2],{RotateRight[{b3,b2,b1,b0,0,0}],{0,b3,b2,b1,b0,0}},{RotateRight[{0,b3,b2,b1,b0,0}],{0,0,b3,b2,b1,b0}},{{b3,b2,b1,b0,0,0},{0,b3,b2,b1,b0,0},{0,0,b3,b2,b1,b0}}},Join[{{a3,a2,a1,a0,0,0},{0,a3,a2,a1,a0,0},{0,0,a3,a2,a1,a0}},{{b3,b2,b1,b0,0,0},{0,b3,b2,b1,b0,0},{0,0,b3,b2,b1,b0}}],{{a3,a2,a1,a0,0,0},{0,a3,a2,a1,a0,0},{0,0,a3,a2,a1,a0},{b3,b2,b1,b0,0,0},{0,b3,b2,b1,b0,0},{0,0,b3,b2,b1,b0}}}


 5)另外一个版本

MatrixSylvester[poly1_, poly2_, var_] :=
Function[{coeffs1, coeffs2},
With[{l1 = Length[coeffs1], l2 = Length[coeffs2]},
(* form the matrix *)
Join @@ (Table[Join[Table[0, {j}], #1,   (* Join括号中的两个Table都是前后补零用的 *)
Table[0, {#2 - 2 - j}]], {j, 0, #2 - 2}]& @@@
{ {coeffs1, l1}, {coeffs2, l2} })]][
(* the coefficients of the two polynomials *)
Sequence @@ (Reverse[CoefficientList[#, var]]& /@ {poly1, poly2})]

MatrixSylvester[a0 + a1 x + a2 x^2 + a3 x^3,
  b0 + b1 x + b2 x^2 + b3 x^3, x] // MatrixForm

  1. 首先明白 f @@@ {{a, b, c}, {d, e}}的结果为{f[a, b, c], f[d, e]},对比f @@ {{a, b, c}, {d, e}},得到 f[{a, b, c}, {d, e}],明白为什么用@@@而不是@@。在Table[Join[Table[0, {j}], #1, Table[0, {#2 - 2 - j}]], {j,
         0, #2 - 2}] & @@@ {{coeffs1, l1}, {coeffs2, l2}}中@@@的作用就是依次将coeffs1作为#1,l1作为#2,然后coeffs2作为#1,l2作为#2。
  2. Sequence @@ (Reverse[CoefficientList[#, var]]& /@ {poly1, poly2})中#代表#1,列表{poly1, poly2}中的项依次代入#,然后Sequence @@的作用就是将Reverse[CoefficientList[#, var]]& /@ {poly1, poly2}的依次执行的两个结果一次性的作为两个实参传递给Pure Function的两个形参{coeffs1, coeffs2},这个以后编程时模仿着做就是了,没有必要搞的太明白,只要It Works就OK
  3. Join的作用就是将最外边的那层大括号在两两之间缝合起来!Join[{{a, b, c}}, {x, y}, {u, v, w}]返回{{a, b, c}, x, y, u, v, w}
  4. 执行 Trace[MatrixSylvester[a0 + a1 x + a2 x^2 + a3 x^3,
       b0 + b1 x + b2 x^2 + b3 x^3, x]]
    返回 {MatrixSylvester[a0+a1 x+a2 x^2+a3 x^3,b0+b1 x+b2 x^2+b3 x^3,x],Function[{coeffs1,coeffs2},With[{l1=Length[coeffs1],l2=Length[coeffs2]},Join@@Apply[Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&,{{coeffs1,l1},{coeffs2,l2}},{1}]]][Sequence@@(Reverse[CoefficientList[#1,x]]&)/@{a0+a1 x+a2 x^2+a3 x^3,b0+b1 x+b2 x^2+b3 x^3}],{{(Reverse[CoefficientList[#1,x]]&)/@{a0+a1 x+a2 x^2+a3 x^3,b0+b1 x+b2 x^2+b3 x^3},{(Reverse[CoefficientList[#1,x]]&)[a0+a1 x+a2 x^2+a3 x^3],(Reverse[CoefficientList[#1,x]]&)[b0+b1 x+b2 x^2+b3 x^3]},{(Reverse[CoefficientList[#1,x]]&)[a0+a1 x+a2 x^2+a3 x^3],Reverse[CoefficientList[a0+a1 x+a2 x^2+a3 x^3,x]],{CoefficientList[a0+a1 x+a2 x^2+a3 x^3,x],{a0,a1,a2,a3}},Reverse[{a0,a1,a2,a3}],{a3,a2,a1,a0}},{(Reverse[CoefficientList[#1,x]]&)[b0+b1 x+b2 x^2+b3 x^3],Reverse[CoefficientList[b0+b1 x+b2 x^2+b3 x^3,x]],{CoefficientList[b0+b1 x+b2 x^2+b3 x^3,x],{b0,b1,b2,b3}},Reverse[{b0,b1,b2,b3}],{b3,b2,b1,b0}},{{a3,a2,a1,a0},{b3,b2,b1,b0}}},Sequence@@{{a3,a2,a1,a0},{b3,b2,b1,b0}},HoldForm[{a3,a2,a1,a0},{b3,b2,b1,b0}]},Function[{coeffs1,coeffs2},With[{l1=Length[coeffs1],l2=Length[coeffs2]},Join@@Apply[Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&,{{coeffs1,l1},{coeffs2,l2}},{1}]]][Sequence[{a3,a2,a1,a0},{b3,b2,b1,b0}]],Function[{coeffs1,coeffs2},With[{l1=Length[coeffs1],l2=Length[coeffs2]},Join@@Apply[Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&,{{coeffs1,l1},{coeffs2,l2}},{1}]]][{a3,a2,a1,a0},{b3,b2,b1,b0}],With[{l1$=Length[{a3,a2,a1,a0}],l2$=Length[{b3,b2,b1,b0}]},Join@@Apply[Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&,{{{a3,a2,a1,a0},l1$},{{b3,b2,b1,b0},l2$}},{1}]],{Length[{a3,a2,a1,a0}],4},{Length[{b3,b2,b1,b0}],4},Join@@Apply[Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&,{{{a3,a2,a1,a0},4},{{b3,b2,b1,b0},4}},{1}],{Apply[Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&,{{{a3,a2,a1,a0},4},{{b3,b2,b1,b0},4}},{1}],{(Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&)[{a3,a2,a1,a0},4],(Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&)[{b3,b2,b1,b0},4]},{(Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&)[{a3,a2,a1,a0},4],Table[Join[Table[0,{j}],{a3,a2,a1,a0},Table[0,{4-2-j}]],{j,0,4-2}],{{Table[0,{j}],{}},{Table[0,{4-2-j}],{0,0}},Join[{},{a3,a2,a1,a0},{0,0}],{a3,a2,a1,a0,0,0}},{{Table[0,{j}],{0}},{Table[0,{4-2-j}],{0}},Join[{0},{a3,a2,a1,a0},{0}],{0,a3,a2,a1,a0,0}},{{Table[0,{j}],{0,0}},{Table[0,{4-2-j}],{}},Join[{0,0},{a3,a2,a1,a0},{}],{0,0,a3,a2,a1,a0}},{{a3,a2,a1,a0,0,0},{0,a3,a2,a1,a0,0},{0,0,a3,a2,a1,a0}}},{(Table[Join[Table[0,{j}],#1,Table[0,{#2-2-j}]],{j,0,#2-2}]&)[{b3,b2,b1,b0},4],Table[Join[Table[0,{j}],{b3,b2,b1,b0},Table[0,{4-2-j}]],{j,0,4-2}],{{Table[0,{j}],{}},{Table[0,{4-2-j}],{0,0}},Join[{},{b3,b2,b1,b0},{0,0}],{b3,b2,b1,b0,0,0}},{{Table[0,{j}],{0}},{Table[0,{4-2-j}],{0}},Join[{0},{b3,b2,b1,b0},{0}],{0,b3,b2,b1,b0,0}},{{Table[0,{j}],{0,0}},{Table[0,{4-2-j}],{}},Join[{0,0},{b3,b2,b1,b0},{}],{0,0,b3,b2,b1,b0}},{{b3,b2,b1,b0,0,0},{0,b3,b2,b1,b0,0},{0,0,b3,b2,b1,b0}}},{{{a3,a2,a1,a0,0,0},{0,a3,a2,a1,a0,0},{0,0,a3,a2,a1,a0}},{{b3,b2,b1,b0,0,0},{0,b3,b2,b1,b0,0},{0,0,b3,b2,b1,b0}}}},Join@@{{{a3,a2,a1,a0,0,0},{0,a3,a2,a1,a0,0},{0,0,a3,a2,a1,a0}},{{b3,b2,b1,b0,0,0},{0,b3,b2,b1,b0,0},{0,0,b3,b2,b1,b0}}},Join[{{a3,a2,a1,a0,0,0},{0,a3,a2,a1,a0,0},{0,0,a3,a2,a1,a0}},{{b3,b2,b1,b0,0,0},{0,b3,b2,b1,b0,0},{0,0,b3,b2,b1,b0}}],{{a3,a2,a1,a0,0,0},{0,a3,a2,a1,a0,0},{0,0,a3,a2,a1,a0},{b3,b2,b1,b0,0,0},{0,b3,b2,b1,b0,0},{0,0,b3,b2,b1,b0}}}
  评论这张
 
阅读(591)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016