矩阵数值计算问题：Julia生成Number Walls时的整数转换错误

问题描述 #

我在Julia中实现了数学领域的Number Walls（数字墙）生成算法，需计算n×n矩阵元素：前三行初始化后，通过复杂循环生成剩余元素。为避免浮点精度问题，已将多数计算转为分数形式。但当参数lato小于20-21时代码正常运行，超过该值后会出现如ERROR: InexactError: BigInt(-582176//73)的错误。根据算法逻辑与初始序列，矩阵元素应均为整数，且报错的分数与整数极为接近，想确认代码是否存在未考虑的疏漏。

实现代码 #

using Plotly, PlotlyJS, Plots
using Statistics, ProgressMeter, ImageShow

riempi_matrice_3 = function(lato)
    
    width = lato
    height = width

    data = zeros(Float32, width, height)

    data[data .== 0] .= NaN

    data[1,:] = zeros(Int, width)
    data[2,:] = ones(Int, width)

    seqq = [0, 1, 1, 1, 0]
    w = 0
    while w <13
        seqq2 = .-1 .* seqq .+ 1
        seqq = vcat(seqq,seqq2)
        w += 1
    end

    data[3,:]= seqq[1:width]

    total_steps = (width - 3) * (width - 2)
    p = Progress(total_steps)

    for i in 3 : (height-1)
        for j in 2:(width-1)
            if !isnan(data[i, j-1]) && !isnan(data[i, j+1]) 
                if data[i,j] != 0 && data[i-1, j] != 0 

                    data[i+1, j] = BigInt(BigInt(data[i, j]^2 - data[i, j-1] * data[i, j+1])//BigInt(data[i-1, j]))  # = F1

                elseif data[i,j] != 0 && data[i-1, j] == 0 

                    l_lato_finestra = 0
                    posiz_sul_lato_dal_fondo = 1
                    posiz_sul_lato_dalla_cima = 1

                    for n in 1:(i)
                        if data[i-n,j] == 0
                            l_lato_finestra +=1
                            n += 1
                        else
                            break
                        end
                    end

                    for n in 1:(i)
                        if data[i-1,j+n] == 0 
                            posiz_sul_lato_dal_fondo +=1
                            n += 1
                        else
                            break
                        end
                    end

                    for n in 1:(i)
                        if data[i-1,j-n] == 0
                            posiz_sul_lato_dalla_cima += 1
                            n += 1
                        else
                            break
                        end
                    end

                    lato_sinistro = BigInt(j-posiz_sul_lato_dalla_cima)

                    rapp_geom_up = BigInt(data[i-l_lato_finestra-1,j-posiz_sul_lato_dalla_cima+1])//BigInt(data[i-l_lato_finestra-1,j-posiz_sul_lato_dalla_cima])                        
                    rapp_geom_left = BigInt(data[i-l_lato_finestra, lato_sinistro])//BigInt(data[i-l_lato_finestra-1, lato_sinistro]) #cambiato il valore 1
                    rapp_geom_right = BigInt(data[i-l_lato_finestra-1, j + posiz_sul_lato_dal_fondo])//BigInt(data[i-l_lato_finestra, j + posiz_sul_lato_dal_fondo]) #cambiato il valore 1
                    rapp_geom_down = BigInt(data[i,j-1])//BigInt(data[i,j])

                    A = 0
                    B = 0
                    C = 0
                    if data[i-l_lato_finestra-1,j+posiz_sul_lato_dal_fondo-posiz_sul_lato_dalla_cima] ==0 

                        return["Il denominatore del rapporto superiore è zero", data]
                    else
                        A=rapp_geom_left*BigInt(data[i-l_lato_finestra-2,j+posiz_sul_lato_dal_fondo-posiz_sul_lato_dalla_cima])//BigInt(data[i-l_lato_finestra-1,j+posiz_sul_lato_dal_fondo-posiz_sul_lato_dalla_cima])
                    end

                    if data[i-posiz_sul_lato_dalla_cima,lato_sinistro] == 0 

                        return["Il denominatore del rapporto a sinistra è zero", data]
                    else
                        B=rapp_geom_up*BigInt(data[i-posiz_sul_lato_dalla_cima,lato_sinistro-1])//BigInt(data[i-posiz_sul_lato_dalla_cima,lato_sinistro]*(-1)^posiz_sul_lato_dal_fondo)
                    end

                    if data[i-posiz_sul_lato_dal_fondo,j+posiz_sul_lato_dal_fondo] ==0 

                        retunr("Il denominatore del rapporto a destra è zero")
                    else
                        C=rapp_geom_down*BigInt(data[i-posiz_sul_lato_dal_fondo,j+posiz_sul_lato_dal_fondo+1])//BigInt(data[i-posiz_sul_lato_dal_fondo,j+posiz_sul_lato_dal_fondo]*(-1)^posiz_sul_lato_dal_fondo)
                    end 

                    D=(BigInt(data[i,j])//rapp_geom_right)

                    # A +- B = x/D +- C 

                    data[i+1,j] = BigInt((A+B-C)*D)

                    # F4

                elseif data[i,j] == 0 

                    posiz_dalla_prima_riga = 1

                    for n in 1:(i)
                        if data[i-n,j] == 0
                            posiz_dalla_prima_riga += 1
                            n += 1
                        else
                            break
                        end
                    end

                    posiz_sul_lato_dal_fondo = 1
                    posiz_sul_lato_dalla_cima = 1



                    for n in 1:(i)

                        if j+n < width 
                            if data[i-posiz_dalla_prima_riga+1,j+n] == 0
                                posiz_sul_lato_dal_fondo += 1
                                n += 1
                            else
                                break
                            end
                        end
                    end

                    for n in 1:(i)
                        if data[i-posiz_dalla_prima_riga+1,j-n] == 0
                            posiz_sul_lato_dalla_cima += 1
                            n += 1
                        else
                            break
                        end
                    end

                    l_lato_finestra = BigInt(posiz_sul_lato_dal_fondo + posiz_sul_lato_dalla_cima - 1)
          
                    if posiz_dalla_prima_riga == l_lato_finestra

                        lato_sinistro = BigInt(j+posiz_sul_lato_dalla_cima)  

                        up = BigInt(data[i-l_lato_finestra,j-posiz_sul_lato_dalla_cima+posiz_sul_lato_dal_fondo])
                        left = BigInt(data[i-posiz_sul_lato_dalla_cima+1,j-posiz_sul_lato_dalla_cima])
                        right = BigInt(data[i-l_lato_finestra+posiz_sul_lato_dalla_cima,j+posiz_sul_lato_dal_fondo])

                        data[i+1,j] = BigInt((-1)^(l_lato_finestra*posiz_sul_lato_dal_fondo)*left*right//up)

                    else

                        data[i+1,j] = Int(0)

                    end
                                               
                end
            end
            next!(p)
        end
    end

    color_matrix = [data[Int(i),Int(j)] == 0 ? RGB(1,1,1) : (isnan(data[Int(i),Int(j)]) ? RGB(1,0,0) : RGB(0, 0, 0)) for i in 1:height/2+2, j in 1:width]
end

问题定位与修复方案 #

核心问题 #

1. 浮点类型的精度污染
   初始化data为Float32矩阵，即便后续用BigInt转换，大整数存入浮点类型时会丢失精度，迭代中误差放大，导致分数计算无法整除。
2. 强制类型转换的风险
   直接将Rational类型强制转为BigInt，若因精度问题分数非整数，就触发InexactError。
3. 循环索引逻辑错误
   循环中手动执行n +=1会跳过元素，导致索引判断错误，引发错误的分数计算。
4. 代码笔误
   存在retunr拼写错误，导致该错误分支逻辑失效。

修复步骤 #

1. 改用整数类型存储矩阵
   避免浮点精度损失，初始化时用支持缺失值的整数矩阵：

   data = Matrix{Union{BigInt, Missing}}(missing, width, height)
   data[1,:] .= BigInt(0)
   data[2,:] .= BigInt(1)
   data[3,:] .= BigInt.(seqq[1:width])
   

   同时将所有isnan判断改为ismissing，data[i,j] !=0改为!ismissing(data[i,j]) && data[i,j] != BigInt(0)。

2. 验证分数整除性后转换
   在转为BigInt前先检查分数是否为整数：

   frac = (BigInt(data[i, j]^2 - data[i, j-1] * data[i, j+1]) // BigInt(data[i-1, j]))
   if denominator(frac) != 1
       # 可添加调试逻辑定位问题
       error("非整数结果：$frac at ($i, $j)")
   end
   data[i+1, j] = numerator(frac)
   

3. 修正循环索引逻辑
   替换手动递增的for循环为while循环，避免跳过元素：

   l_lato_finestra = 0
   n = 1
   while n <= i && !ismissing(data[i-n,j]) && data[i-n,j] == BigInt(0)
       l_lato_finestra +=1
       n +=1
   end
   

4. 修复拼写错误
   将retunr改为return，确保错误分支正常执行。

参考文档翻译 #

《数字墙算法：LFSR实用指南》

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