Problem with `cat`

General Usage
1

I have a problem with cat. I am using polynomials (my own package but this is not the problem)

julia> q=Pol(:q)
Pol{Int64}: q

julia> (q+1)^3
Pol{Int64}: q³+3q²+3q+1

and I am all the time using cat to construct block-diagonal matrices:

julia> m=[q 0; 1 q+1]
2×2 Array{Pol{Int64},2}:
 q  0  
 1  q+1

julia> cat(m,m;dims=(1,2))
4×4 Array{Pol{Int64},2}:
    q       0    #undef  #undef
    1       q+1  #undef  #undef
 #undef  #undef     q       0  
 #undef  #undef     1       q+1

here is the problem: why is cat putting #undef in my matrix instead of polynomial 0?

julia> zero(q)
Pol{Int64}: 0

I tracked down the problem to lines 1443–1445 of abstractarray.jl:

  if T <: Number && count(!iszero, catdims) > 1
        fill!(A, zero(T))
    end

Is not the test T<:Number wrong? Should not the right test be that T has the method zero?

Thank you for help on what to do…

2

Can’t you simply declare your polynomial type to be <: Number, or is that abstraction wrong?

3

Why not simply roll your own function to construct a diagonal matrix?

function diagm(m, n)
     l = size(m)[1]
     nd = l * n
     d = fill(0, nd, nd);
     for j = 1 : n
           idx1 = 1 + (j-1) * l
           idx2 = idx1 + (l-1)
           d[idx1:idx2, idx1:idx2] = m
      end
      return d 
end

where the 0 in the fill would have to be replaced with the zero for your polynomial.

4

I would say it is wrong. I distinguish polynomials from their coefficients by saying that the coefficients are <:Number

5

That reads like an implementation detail, not like a hard requirement of Poly{T} where T.

In any case, you can also of course overload cat for your type.

EDIT:

I agree with you in that this probably should not just check for <: Number, but making your type <: Number is a valid fix until this can be fixed in Base. Moreover, polynomials do very much feel like numbers to me, so not having that parent type and all the associated methods available seems weird to me.

6

The question is not whether I can write myself functions from julia Base (I can) but whether cat is properly designed. By the way your code is wrong

  • It assumes m is square
  • You cannot initialize d to 0. One should rather use zero(promote_type(eltype(m),eltype(n)))
7

Well, if the cat code uses the zero function, then it makes sense to check for Number subtypes, because that’s how zero is defined (in addition to some other types).

8

I tried to make a system where I do not distinguish polynomials from their coefficients, to get free
multi-variate polynomials (polynomials with coefficients another polynomial) but this did not work (too long to explain here why). In my system I have many numbers: integers, rationals, BigInts, Cyclotomics, field extensions; I think my type system is ok.

9

I don’t know how your package is structured, but are you sure things break as soon as you add <: Number to the type definition? Because other than that, I don’t see how this could be resolved, if the <: Number check for zero is intended.

10

In my mind Number means “some subfield of the complex numbers” or perhaps also “element of a finite field”. If you think it should mean “any element of a ring” then that’s a completely different viewpoint.
Anyway, thinking in terms of “duck typing”, I think the logical test should be the existence or not of a zero method. I would like some authoritative answer (like a sound argument why it should not be so).