公告:本站正式转型为非交互式静态网站!
转型:本站将通过笔记和博客的形式继续为大家服务,关于 Mathematica 问答服务请移步至QQ群:365716997
联系:如有问题请联系QQ群管理员,或发送邮件至:lixuan.xyz@qq.com。
感谢:最后非常感谢大家多年来的支持与帮助!
参考《互联网跟帖评论服务管理规定》《中华人民共和国网络安全法》《网络信息内容生态治理规定》《互联网用户账号信息管理规定》

—— 2022-11-27

欢迎来到 Mathematica 问答社区

提问时请贴上文本代码

语法高亮:在编辑器中点击

被禁止的话题:广告破解

请阅读:《提问的智慧》

备用域名:mma.ooo

支持LaTex数学公式
行内公式标识符:\$ 或“$\backslash ($”+“$\backslash )$”,
行间公式标识符:\$\$ 或 “$\backslash [$”+“$\backslash ]$”

社区建议QQ群:365716997

分类

0 投票
3.2k 浏览
In[52]:= CompileSuperLog = 
  Compile[{{x, _Complex}, {z, _Real}}, 
   Block[{dec, mtx, y1, y2, terms1, terms2},
    mtx = Table[If[k == 12, If[j == 0, 1, 0],
       k^j/Gamma[k + 1] - If[j == k, Log[x]^-k, 0]], {k, 1, 12}, {j, 
       0, 10}]; mtx[[Range[10]]] = RowReduce[mtx[[Range[10]]]];
    terms1 =  Most[mtx][[All, 12]];
    y1 =  Which[
      z < 0,
        Sum[
        terms1[[k]]*x^(z*k)/Gamma[k + 1], {k, 1, Length[terms1]}] - 2,
      z == 0, -1,
      0 < z < 1,
        Sum[terms1[[k]]*z^k/Gamma[k + 1], {k, 1, Length[terms1]}] - 1,
      z == 1, 0,
      z > 1,
        Block[{i = -1}, 
            i + Sum[terms1[[k]]*
                (NestWhile[Log[x, #] &, z, (i++; # > 1) &])^k/
            Gamma[k + 1], {k, 1, Length[terms1]}] - 1],
      True, 0.];
    terms2 = RowReduce[mtx][[All, 12]];
    y2 =  Which[
      z < 0,
        Sum[
        terms2[[k]]*x^(z*k)/Gamma[k + 1], {k, 1, Length[terms2]}] - 2,
      z == 0, -1,
      0 < z < 1,
        Sum[terms2[[k]]*z^k/Gamma[k + 1], {k, 1, Length[terms2]}] - 1,
      z == 1, 0,
      z > 1,
        Block[{j = -1}, 
            j + Sum[terms2[[k]]*
                (NestWhile[Log[x, #] &, z, (j++; # > 1) &])^k/
            Gamma[k + 1], {k, 1, Length[terms2]}] - 1],
      True, 0.];
    If[y1 == y2, y2,
     dec = Floor[-Log[10, Abs[y2 - y1]]];
     (Floor@Re[y2*10^dec] + Floor@Im[y2*10^dec]*I)/10^dec]]];
CompilePrint@CompileSuperLog

Out[53]= "
        2 arguments
        6 Boolean registers
        16 Integer registers
        8 Real registers
        11 Complex registers
        5 Tensor registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking off
        RuntimeAttributes -> {Listable}

(*略过非重点*)

137	C7 = Part[ T(C1)3, I9]
138	R1I9T(C1)3I4T(C2)0 = MainEvaluate[ Function[{x, z, kCompile$42, \
terms1Compile$43, iCompile$44, mtxCompile$45}, Block[{k = \
kCompile$42, terms1 = terms1Compile$43, i = iCompile$44, mtx = \
mtxCompile$45}, {NestWhile[Log[x, #1] & , z, (i++; #1 > 1) & ], k, \
terms1, i, mtx}]][ C0, R0, I9, T(C1)3, I4, T(C2)0]]
139	R5 = Power[ R1, I9]

(*......*)

238	C7 = Part[ T(C1)4, I9]
239	R5C6I9T(C1)3T(C1)4I4T(C2)0 = MainEvaluate[ Function[{x, z, \
y1Compile$46, kCompile$47, terms1Compile$48, terms2Compile$49, \
jCompile$50, mtxCompile$51}, Block[{y1 = y1Compile$46, k = \
kCompile$47, terms1 = terms1Compile$48, terms2 = terms2Compile$49, j \
= jCompile$50, mtx = mtxCompile$51}, {NestWhile[Log[x, #1] & , z, \
(j++; #1 > 1) & ], y1, k, terms1, terms2, j, mtx}]][ C0, R0, C6, I9, \
T(C1)3, T(C1)4, I4, T(C2)0]]
240	R4 = Power[ R5, I9]

可以看到会出现无法编译的NestWhile。但是明明矩阵是标准的,所有的函数都是能被编译的并且没有忘记声明的中间项啊。。。难道说NestWhile不可被编译?

联动:http://mmaqa.com/index.php/734

================================================================

重构NestWhile:

CompileSuperLog2 = 
  Compile[{{x, _Complex}, {z, _Real}}, 
   Block[{dec, mtx, y1, y2, terms1, terms2},
    mtx = Table[If[k == 12, If[j == 0, 1, 0],
       k^j/Gamma[k + 1] - If[j == k, Log[x]^-k, 0]], {k, 1, 12}, {j, 
       0, 10}]; mtx[[Range[10]]] = RowReduce[mtx[[Range[10]]]];
    terms1 =  Most[mtx][[All, 12]];
    y1 =  Which[
      z < 0,
      	Sum[
        terms1[[k]]*x^(z*k)/Gamma[k + 1], {k, 1, Length[terms1]}] - 2,
      z == 0, -1,
      0 < z < 1,
      	Sum[terms1[[k]]*z^k/Gamma[k + 1], {k, 1, Length[terms1]}] - 1,
      z == 1, 0,
      z > 1,
      	Block[{i = -1, l = z}, 
       		i + Sum[terms1[[k]]*
          		(While[i++; If[l > 1, False, True, True], l = Log[x, l]]; 
             l)^k/
           	Gamma[k + 1], {k, 1, Length[terms1]}] - 1],
      True, 0.];
    terms2 = RowReduce[mtx][[All, 12]];
    y2 =  Which[
      z < 0,
      	Sum[
        terms2[[k]]*x^(z*k)/Gamma[k + 1], {k, 1, Length[terms2]}] - 2,
      z == 0, -1,
      0 < z < 1,
      	Sum[terms2[[k]]*z^k/Gamma[k + 1], {k, 1, Length[terms2]}] - 1,
      z == 1, 0,
      z > 1,
      	Block[{j = -1, m = z}, 
       		j + Sum[terms2[[k]]*
          		(While[i++; If[m > 1, False, True, True], m = Log[x, m]]; 
             m)^k/
           	Gamma[k + 1], {k, 1, Length[terms2]}] - 1],
      True, 0.];
    If[y1 == y2, y2,
     dec = Floor[-Log[10, Abs[y2 - y1]]];
     (Floor@Re[y2*10^dec] + Floor@Im[y2*10^dec]*I)/10^dec]]];
CompilePrint@CompileSuperLog2

关键更改:

(* 更改前 *)
	Block[{i = -1}, 
 		i + Sum[terms$$$[[k]]*
    		(NestWhile[Log[x, #] &, z, (i++; # > 1) &])^k/
     	Gamma[k + 1], {k, 1, Length[terms$$$]}] - 1]
(* 更改后 *)
Block[{i = -1, l = z}, 
 		i + Sum[terms$$$[[k]]*
    		(While[i++; If[l > 1, False, True, True], l = Log[x, l]]; l)^k/
     	Gamma[k + 1], {k, 1, Length[terms$$$]}] - 1]

结果:

(* 省略非重点 *)
133	R4I4T(C1)3R2T(C2)0 = MainEvaluate[                                \
                                                                      \
                                                                      \
                                                                      \
  k
                                                                      \
                                                                      \
                                           (While[i++; If[l > 1, \
False, True, True], l = Log[x, l]]; l)
Function[{x, z, iCompile$21, terms1Compile$22, lCompile$23, \
mtxCompile$24}, Block[{i = iCompile$21, terms1 = terms1Compile$22, l \
= lCompile$23, mtx = mtxCompile$24}, {Sum[terms1[[k]] \
-------------------------------------------------------------, {k, 1, \
Length[terms1]}], i, terms1, l, mtx}]]
                                                                      \
                                                                      \
                                                                   Gam\
ma[k + 1][ C0, R0, I4, T(C1)3, R2, T(C2)0]]

(* ...... *)

216	R4T(C1)4T(C1)3C2R2I4T(C2)0 = MainEvaluate[                        \
                                                                      \
                                                                      \
                                                                      \
                                                                      \
                  k
                                                                      \
                                                                      \
                                                                      \
                                                   (While[i++; If[m > \
1, False, True, True], m = Log[x, m]]; m)
Function[{x, z, terms2Compile$25, terms1Compile$26, y1Compile$27, \
mCompile$28, jCompile$29, mtxCompile$30}, Block[{terms2 = \
terms2Compile$25, terms1 = terms1Compile$26, y1 = y1Compile$27, m = \
mCompile$28, j = jCompile$29, mtx = mtxCompile$30}, {Sum[terms2[[k]] -------------------------------------------------------------\
, {k, 1, Length[terms2]}], terms2, terms1, y1, m, j, mtx}]]
                                                                      \
                                                                      \
                                                                      \
                                                                      \
     Gamma[k + 1][ C0, R0, T(C1)4, T(C1)3, C2, R2, I4, T(C2)0]]

并不成功。

===============================================================

没错,# > 1并不正确。应改为Abs@# > 1。

分类:程序包 | 用户: EmberEdison (806 分)
修改于 用户:EmberEdison
NestWhile[Log[x, #] &, z, (Abs@# > 1) &]], x != 1&& z>1.至少我x从-5到5,z的1.1到3,还试了几个复数都没有锁死了啊....
之前看错了,不是“随便写个z和a”就死循环,是有时候会死循环,例如Abs@NestList[Log[2 + I, #] &, 20., 50],死循环。
这个完全应该是判断的问题啊...2.79491 - 1.61032 I的Abs已经大于1了...奇怪的bug真多
2.79491 - 1.61032 I // Abs是3.2啊。。Abs作用在复数上是取模呀。
我有点错乱了,我明天检查下,但是限定为实数应该是没问题的

登录 或者 注册 后回答这个问题。

...