def kron(b)
      b = NArray.cast(b)
      nda = ndim
      ndb = b.ndim
      shpa = shape
      shpb = b.shape
      adim = [:new]*(2*[ndb-nda,0].max) + [true,:new]*nda
      bdim = [:new]*(2*[nda-ndb,0].max) + [:new,true]*ndb
      shpr = (-[nda,ndb].max..-1).map{|i| (shpa[i]||1) * (shpb[i]||1)}
      (self[*adim] * b[*bdim]).reshape(*shpr)
    end
def kron(b)
      b = NArray.cast(b)
      nda = ndim
      ndb = b.ndim
      shpa = shape
      shpb = b.shape
      adim = [:new]*(2*[ndb-nda,0].max) + [true,:new]*nda
      bdim = [:new]*(2*[nda-ndb,0].max) + [:new,true]*ndb
      shpr = (-[nda,ndb].max..-1).map{|i| (shpa[i]||1) * (shpb[i]||1)}
      (self[*adim] * b[*bdim]).reshape(*shpr)
    end
Kronecker product of two arrays.

     kron(a,b)[k_0, k_1, ...] = a[i_0, i_1, ...] * b[j_0, j_1, ...]
        where:  k_n = i_n * b.shape[n] + j_n

 @param b [Numo::NArray]
 @return [Numo::NArray]  return Kronecker product
 @example
   Numo::DFloat[1,10,100].kron([5,6,7])
   => Numo::DFloat#shape=[9]
   [5, 6, 7, 50, 60, 70, 500, 600, 700]
   Numo::DFloat[5,6,7].kron([1,10,100])
   => Numo::DFloat#shape=[9]
   [5, 50, 500, 6, 60, 600, 7, 70, 700]
   Numo::DFloat.eye(2).kron(Numo::DFloat.ones(2,2))
   => Numo::DFloat#shape=[4,4]
   [[1, 1, 0, 0],
    [1, 1, 0, 0],
    [0, 0, 1, 1],
    [0, 0, 1, 1]]