# Undesired and inconsistent returning of NaN when instantiating a custom structure several times
**URL:** https://discourse.julialang.org/t/undesired-and-inconsistent-returning-of-nan-when-instantiating-a-custom-structure-several-times/82529
**Category:** General Usage
**Tags:** structure
**Created:** [June 9, 2022, 9:31pm UTC](https://discourse.julialang.org/t/undesired-and-inconsistent-returning-of-nan-when-instantiating-a-custom-structure-several-times/82529 "2022-06-09T21:31:34Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@klaus.68](https://discourse.julialang.org/u/klaus.68)
#### Post date: [June 9, 2022, 9:31pm UTC](https://discourse.julialang.org/t/undesired-and-inconsistent-returning-of-nan-when-instantiating-a-custom-structure-several-times/82529/1 "2022-06-09T21:31:34Z")
Hi,
This is my first post to the boards about a behavior I haven’t seen yet in a forum (appreciate any direction to a post if it exists!). I was attempting to build a b-spline structure atop a polynomial structure, and to do this I have to define a family of polynomials recursively. What I’m finding is that when I call the code from REPL, eg regbn(7), I get seemingly random occurences of NaN in the polynomial coefficients. Here is my attempt at a much simplified MWE demonstrating the issue:
```julia
""" Auxilliary structure routines """
struct Poly1d{T}
deg::Int64
c::Vector{T}
function Poly1d(c::Vector)
deg = convert(Int64,length(c))-1;
S = typeof(c[1]);
return new{S}(deg,c)
end
end
function Base.:+(p₁::Poly1d,p₂::Poly1d)
deg = p₁.deg<=p₂.deg ? p₂.deg : p₁.deg;
T = typeof(p₁.c[1]+p₂.c[1]);
c = Vector{T}(undef,deg+1);
@inbounds for i=1:(deg+1)
c[i] *= 0;
c[i] += (p₁.deg+1)<=i ? p₁.c[i] : 0;
c[i] += (p₂.deg+1)<=i ? p₂.c[i] : 0;
end
return Poly1d(c)
end
""" Recursion relation """
function regbⁿ(P::Vector{Poly1d{Float64}})
n = length(P);
Q = fill(Poly1d([0.0]),n+1);
@inbounds for k=1:n+1
if k!=n+1
Q[k] += Poly1d([0.0,1.0])+P[k];
end
if k!=1
Q[k] += Poly1d([-1.0,1.0])+P[k-1];
end
end
return Q
end
function regbⁿ(deg::Integer)
P = Vector{Vector{Poly1d{Float64}}}(undef,deg+1);
P[1] = [Poly1d([1.0])];
if deg==0
return P[1]
elseif deg>0
@inbounds for k=1:deg
flagnan = true
while flagnan
P[k+1] = regbⁿ(P[k])
# Catch the erroneous NaN's
ct = 0
for i=1:length(P[k+1])
flag = isnan.(P[k+1][i].c).==true
if sum(flag)!=0
ct+=1;
end
end
if ct==0
flagnan=false
else
@warn "At P[k+1]=P[$k+1] there were $ct many polys with nan's before continued"
end
end
end
return P[end]
end
end
```
You’ll notice I’ve wrapped the part producing random NaN’s in a while loop that only terminates when no NaN’s are produced and warns the user how many times NaN’s occurred. I’ve also attached a screenshot of my REPL showing how the output changes despite being given the exact same input. I believe I have avoided any accidental argument mutation. The coefficients that do appear seem to be correct, and sometimes the code does run without returning any NaN’s. On the other hand, at other times calls will produce many more NaN’s than what the screenshot shows.
In the screenshot you’ll see in the first and third calls, the NaN’s appeared in different polynomials (k=4 vs k=7) and in the second call no NaN was produced.
The while loop-fix may be a workable solution, but I wish I understood what was happening better to guard against in the future.
Appreciate any insights!
---
### Author:  [@jondea](https://discourse.julialang.org/u/jondea)
#### Post date: [June 9, 2022, 9:37pm UTC](https://discourse.julialang.org/t/undesired-and-inconsistent-returning-of-nan-when-instantiating-a-custom-structure-several-times/82529/2 "2022-06-09T21:37:59Z")
The vector `c` is initialised as undefined so it may indeed sometimes be NaN. Should you be setting it to zero instead? Either each entry in the loop, or using `zeros` when the vector is initialised.
---
### Author:  [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [June 9, 2022, 9:40pm UTC](https://discourse.julialang.org/t/undesired-and-inconsistent-returning-of-nan-when-instantiating-a-custom-structure-several-times/82529/3 "2022-06-09T21:40:08Z")
> [@klaus.68](#):
>
> ```julia
> function Poly1d(c::Vector)
> deg = convert(Int64,length(c))-1;
> S = typeof(c[1]);
>
> return new{S}(deg,c)
> end
>
> ```
It’d be better to write this like
```julia
function Poly1d(c::Vector{T}) where T
new{T}(length(c)-1, c) # `length` already returns an `Int`
end
```
This gives you a `Poly1d{T}` with the same type parameter as the element type of the vector. The same goes for your `Base.:+` - you can use the type parameters here again:
```julia
function Base.:+(p1::Poly1d{T1}, p2::Poly1d{T2}) where {T1, T2}
deg = max(p1.deg, p2.deg)
rettype = promote_type(T1, T2)
c = zeros(rettype, deg+1)
for idx in eachindex(c)
c[idx] += (p1.deg + 1) <= idx ? p1.c[idx] : zero(rettype)
c[idx] += (p2.deg + 1) <= idx ? p2.c[idx] : zero(rettype)
end
return Poly1d(c)
end
```
Notice that I replaced the `Vector{T}(undef, deg+1)` with `zeros(rettype, deg+1)` - this is to prevent getting `NaN` from uninitialized memory and instead get a zero initialized array. You could also replace the `c[i] *= 0` with just `c[i] = zero(rettype)`, since `NaN*0 == NaN`, following [IEEE floating point NaN semantics](https://en.wikipedia.org/wiki/NaN#Operations_generating_NaN):
```julia
julia> NaN * 0
NaN
```
---
### Author:  [@mcabbott](https://discourse.julialang.org/u/mcabbott)
#### Post date: [June 9, 2022, 9:47pm UTC](https://discourse.julialang.org/t/undesired-and-inconsistent-returning-of-nan-when-instantiating-a-custom-structure-several-times/82529/4 "2022-06-09T21:47:30Z")
I think this still reads out of bounds. `(p₁.deg+1)<=i ? p₁.c[i] : 0` tests that `i` is the last element of `p₁.c` or later, and then chooses to read from it. Deleting `@inbounds` (or using @Sukera’s version):
```julia
julia> Poly1d([1,2,3]) + Poly1d([4,5,6,7,8])
ERROR: BoundsError: attempt to access 3-element Vector{Int64} at index [4]
```
That may give you NaN (or segfaults) in addition to `isnan(NaN * 0)` as mentioned.
Edit: If the intention is to pad with zeros, then one tidy way is:
```julia
Base.:+(p1::Poly1d, p2::Poly1d) = Poly1d(padadd(p1.c, p2.c))
function padadd(x::AbstractVector, y::AbstractVector)
length(x) >= length(y) || return padadd(y, x)
[x[i] + get(y, i, zero(eltype(y))) for i in eachindex(x)]
end
```
---
### Author:  [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [June 9, 2022, 9:51pm UTC](https://discourse.julialang.org/t/undesired-and-inconsistent-returning-of-nan-when-instantiating-a-custom-structure-several-times/82529/5 "2022-06-09T21:51:42Z")
I did have a bad feeling about it, thanks! This should be better:
```julia
for idx in eachindex(c)
c[idx] += checkbounds(Bool, p1.c, idx) ? p1.c[idx] : zero(rettype)
c[idx] += checkbounds(Bool, p2.c, idx) ? p2.c[idx] : zero(rettype)
end
```
---
### Author:  [@klaus.68](https://discourse.julialang.org/u/klaus.68)
#### Post date: [June 10, 2022, 3:31am UTC](https://discourse.julialang.org/t/undesired-and-inconsistent-returning-of-nan-when-instantiating-a-custom-structure-several-times/82529/6 "2022-06-10T03:31:16Z")
Thanks everyone! Much clearer now. I hadn’t realized undef initialization can also generate NaN’s. And also @Sukera thanks for the style tips. I’m still building intuition for where-keyword uses as well as (custom) type promotions. I hadn’t directly used zeros (when I probably should have) because I didn’t know the cleanest way of initializing a zero element of a custom type. The comments help with that too.
Cheers