Undesired and inconsistent returning of NaN when instantiating a custom structure several times

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:

""" Auxilliary structure routines """
struct Poly1d{T}

	function Poly1d(c::Vector)
		deg = convert(Int64,length(c))-1;
		S = typeof(c[1]);

		return new{S}(deg,c)
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;

	return Poly1d(c)

""" 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];
		if k!=1
			Q[k] += Poly1d([-1.0,1.0])+P[k-1];
	return Q
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
				if ct==0
					@warn "At P[k+1]=P[$k+1] there were $ct many polys with nan's before continued" 
		return P[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!

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.

1 Like

It’d be better to write this like

function Poly1d(c::Vector{T}) where T
    new{T}(length(c)-1, c) # `length` already returns an `Int`

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:

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)

    return Poly1d(c)

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:

julia> NaN * 0

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> 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:

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)]

I did have a bad feeling about it, thanks! This should be better:

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)

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.


1 Like