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。