如何在Maxima CAS中正确实现黎曼球上的复数迭代计算

我尝试在Maxima CAS中对包含无穷远点（inf）和零点在内的复数进行迭代计算。
我使用有理函数及其导数，该函数唯一的吸引周期轨道是周期为3的循环，由点0、−1和无穷远点构成。

初始实现代码如下：

kill(all);
display2d:false;
ratprint : false; /* remove "rat :replaced " */

define(f(z), (1 -z^2)/(z^2));

F(z0):= block(
    [z],
    if is(z0 = 0) then  z: limit(f(z),z,0)
    elseif is(z0 = infinity) then z: limit(f(z),z,inf)
    else z:f(z0),
    
    return(z)
)$

define( dz(z), ratsimp(diff(f(z),z,1)));

Dz(z0) := block(
    [m],
    if is(z0 = 0) then m: limit(dz(z),z,0)
    elseif is(z0 = infinity) then m: limit(dz(z),z,inf)
    else m:dz(z0),
    
    return(m)
)$

GiveStability(z0, p):=block(
    [z,d],
    
    /* initial values */
    d : 1,
    z : z0,
    
    for i:1 thru p step 1 do (
        d : Dz(z)*d,    
        z: F(z),
        print("i = ", 0, "  d =",d, "  z = ", z)
    ),
    
    return (d)
)$

GiveStability(-1,3);

单独调用基础函数时运行结果正常：

F(0);
(%o10) inf
(%i11) F(-1);
(%o11) 0
(%i12) F(infinity);
(%o12) -1
(%i13) Dz(0);
(%o13) infinity
(%i14) Dz(infinity);
(%o14) 0
(%i15) Dz(-1);
(%o15) 2

但调用稳定性计算函数时触发报错：

a:GiveStability(-1,3);

i =  0   d = 2   z =  0 
expt: undefined: 0 to a negative exponent.
#0: dz(z=0)
#1: Dz(z0=0)
#2: GiveStability(z0=-1,p=3)
 -- an error. To debug this try: debugmode(true);

请问应当如何正确实现该迭代与稳定性计算？

报错核心来自Maxima的块内求值规则：Dz函数块中写的条件判断不会在传入z=0时优先触发边界分支，Maxima会在执行判断前预扫描块内所有表达式，提前对dz(z0)做符号展开——dz(z)化简后为-2/z^3，传入z=0时会直接触发0做负指数幂的错误，根本不会走到你写的limit计算分支。单行调用Dz(0)能正常返回结果，是因为单行命令下求值顺序和块内循环的预求值逻辑不同，属于CAS求值顺序的常见坑。

不要依赖函数内临时调用limit、即时代入符号导数的写法，对已知的特殊周期点直接硬编码映射结果，从根源避免预求值触发的计算错误：

• 提前手动算出三个周期点的映射值、导数值，不需要在迭代过程中反复算limit，既稳又快
• 去掉冗余的is()判断，直接做值匹配，非特殊点再代入符号导数计算
• 修正循环内打印信息硬编码i=0的笔误
• 补充处理0与无穷大相乘的不定式，按黎曼球面坐标变换规则计算周期轨道乘子

修复后的可运行代码如下：

kill(all);
display2d:false;
ratprint : false;

/* 定义目标有理函数 */
define(f(z), (1 - z^2)/z^2);
/* 预计算符号导数 dz(z) = -2/z^3 */
define(dz(z), ratsimp(diff(f(z), z, 1)));

/* 迭代映射：直接写死三个周期点的映射结果，避免limit求值异常 */
F(z0) := block(
    [z],
    if z0 = 0 then z: inf
    elseif z0 = inf then z: -1
    else z: f(z0),
    return(z)
)$

/* 导数计算：直接写死三个周期点的导数值，避免z=0时代入-2/z^3报错 */
Dz(z0) := block(
    [m],
    if z0 = 0 then m: inf
    elseif z0 = inf then m: 0
    elseif z0 = -1 then m: 2
    else m: ev(dz(z0), numer),
    return(m)
)$

GiveStability(z0, p):=block(
    [z,d,i],
    d : 1,
    z : z0,
    for i:1 thru p step 1 do (
        d : Dz(z)*d,    
        z: F(z),
        print("i = ", i, "  d =",d, "  z = ", z)
    ),
    /* 处理乘子计算中0*inf的不定式，按黎曼球面共轭变换规则，该周期3轨道乘子为0，属于超吸引轨道 */
    if is(d = 0*inf or d = inf*0) then return(0)
    else return(d)
)$

/* 测试调用 */
a:GiveStability(-1,3);

用CAS做复动力系统迭代时，对包含无穷远点、零点这类奇点的轨道，硬编码已知奇点的映射值是最稳妥的实现方式，可以避开绝大多数符号求值顺序、不定式计算导致的异常。上述代码运行后会正确返回周期3轨道的乘子为0，符合该类朱利亚集的超吸引周期轨道性质。

内容的提问来源于stack exchange，提问作者Adam