计算矩阵数值项的问题。朱莉娅对结果的评估是错误的吗？

Question

我在评估矩阵 nxn 的项时遇到麻烦，其中从前 3 行开始，我使用具有相当复杂循环的算法计算矩阵的其余部分。我正在尝试编写用于生成“数字墙”的代码，这是数学家的 YouTube 频道中讨论的一个非常好的数学主题。代码写在这里：

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

我将大部分计算转换为分数，因此避免了浮点数的问题。我想我总结了我可能看到的所有情况，但我的问题是：此代码适用于低于 20-21 的“lato”值。之后，它开始给出各种错误，例如“ERROR: InexactError: BigInt(-582176//73)”。

对我来说有两个原因应该使这个错误不可能发生：

• 代码应该是我用来生成第三行值的序列，使得矩阵内不可能没有整数值。
• 朱莉娅表明的这些有问题的分数非常接近整数。

我是否可能忘记了代码中的某些内容，以防止矩阵单元最终出现错误评估？

所用公式的文档：“数墙算法：LFSR 食谱”

Answer 1

tl；博士

您的代码将所有内容转换为

Float32

值，并且当

lato

变大时，您的代码会出现浮点数学错误。将

data

声明为

Matrix{BigInt}

或

Matrix{Union{BigInt, Nothing}}

，并相应地调整您的代码。

您的代码正在将所有内容转换为

Float32

值

您的第一个假设，即您的代码使得矩阵中不可能不有整数值，这是不正确的。事实上，你的

data

矩阵不能包含任何 除了

Float32

值，因为这就是你声明它的方式：

data = zeros(Float32, 2, 2);
typeof(data) == Matrix{Float32}  # true
eltype(data) == Float32  # true

您可能认为通过将矩阵中的值分配给

Int

值将存储

Int

值，但由于您将

data

声明为

Float32

值的矩阵，Julia 会将这些值转换为最佳值

Float32

表示它可以:

data[1,:] = zeros(Int, 2);
typeof(data) == Matrix{Float32}  # true
typeof(data[1,1]) == Float32  # true

data

不存储

Int

，它存储

Float32

。将某些元素设置为

Rational{BigInt}

值并不重要：Julia 也会将它们转换为

Float32

：

data[2,1] = BigInt(1)//BigInt(10);
typeof(data) == Matrix{Float32}  # true
typeof(data[2,1]) == Float32  # true

您的代码正在演示浮点数学错误

现在，希望您能够看到这如何成为您的代码中的问题。你总结了很多你可能认为是

Rational{BigInt}

值的东西，但实际上是

Float32

值。这是一个简单的例子，说明事情可能会出错：

data[1,1] = BigInt(1)//BigInt(10)
data[1,2] = BigInt(0)

for _ in 1:10
    data[1,2] += data[1,1]
end

data[1,2] == BitInt(1)  # false, because Float32(1)/Float32(10) != 1//10

以不同方式声明

data

如果您希望

data

成为

BigInt

值的矩阵，只需这样声明即可：

data = zeros(BigInt, width, height)
eltype(data) == BigInt  # true

（如果你真的希望

data

包含有理值，你可以使用

Rational{BigInt}

代替，但是如果你检查你的算法，我想你会发现你不需要它们。）

现在，您可能会陷入 python/pandas 陷阱，因为您需要一个值来表示未定义的内容，因此需要

data

能够包含

NaN

，因此要求

data

是浮点值矩阵；然而，Julia 可以处理类型的联合，并且

Union{T, Nothing}

是处理这种情况的常见方法：

data = Matrix{Union{BigInt, Nothing}}(nothing, width, height)

This 将

data

声明为可以保存

BigInt

或

Nothing

类型值的矩阵，并将其初始化为所有

width

值的

height

by

nothing

矩阵。要检查未初始化的值，请使用

data[i,j] == nothing

或

isnothing(data[i,j])

而不是检查

NaN

。