I have troubles in evaluating terms of a matrix nxn, where, starting with the first 3 rows, I compute the rest of the matrix with an algorithm which have a quite complicated loop. I'm trying to write the code for the generation of Number Walls, a pretty nice topic in mathematics Mathologer's youtube channel talked about. The code is wrote here:
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
I converted the most part of the calculations to fractions, so I avoid the problem with floating points. I think I summarized the totality of the cases I can potentially see, but my issue is: This code works for values of 'lato' below 20-21. After that, it starts to give a variety of errors like 'ERROR: InexactError: BigInt(-582176//73)'.
There are two reasons for me that should make this error impossible:
Is it possible that I forgot something in the code that prevent an eventual wrong evaluation of matrix cells?
Documentation on the formulas used: "The Number-Wall Algorithm: an LFSR Cookbook"
Your code is converting everything into Float32 values, and your code is demonstrating floating point math errors when lato becomes large. Declare data to be a Matrix{BigInt} or Matrix{Union{BigInt, Nothing}} instead, and adjust your code accordingly.
Float32 valuesYour first assumption, that your code makes it impossible not to have integer values inside the matrix, is incorrect. In fact, your data matrix cannot contain anything except Float32 values, because that is how you declared it:
data = zeros(Float32, 2, 2);
typeof(data) == Matrix{Float32} # true
eltype(data) == Float32 # true
You may think that by assigning values in the matrix to Int values will store Int values, but because you declared data to be a matrix of Float32 values, Julia will convert the values into the best possible Float32 representation it can:
data[1,:] = zeros(Int, 2);
typeof(data) == Matrix{Float32} # true
typeof(data[1,1]) == Float32 # true
data does not store Ints, it stores Float32s. It doesn't matter that you set some elements to Rational{BigInt} values: Julia will just convert those into Float32s as well:
data[2,1] = BigInt(1)//BigInt(10);
typeof(data) == Matrix{Float32} # true
typeof(data[2,1]) == Float32 # true
Now, hopefully, you can see how this becomes a problem in your code. You are summing up a lot of things that you may think are Rational{BigInt} values, but are actually Float32 values. This is a simple example of how things can go wrong:
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 differentlyIf you want data to be a matrix of BigInt values, just declare it that way:
data = zeros(BigInt, width, height)
eltype(data) == BigInt # true
(If you really want data to contain rational values, you can use Rational{BigInt} instead, but if you inspect your algorithm, I think you will find that you won't need them.)
Now, you may be falling into a python/pandas gotcha by assuming that, because you need a value to represent something undefined, you need data to be able to contain NaN and therefore require data to be a matrix of floating point values; however, Julia can handle unions of types, and a Union{T, Nothing} is a common way to handle this very case:
data = Matrix{Union{BigInt, Nothing}}(nothing, width, height)
This declares data as a matrix that can hold values of type BigInt or Nothing and initializes it to a width by height matrix of all nothing values. To check for uninitialized values, use data[i,j] == nothing or isnothing(data[i,j]) instead of checking against NaN.