A = [ 1 2 ]
    [ 3 4 ]

B = [ 5 6 ]
    [ 7 8 ]

A x B = [ 1 x [ 5 6 ]   2 x [ 5 6 ] ]
        [     [ 7 8 ]       [ 7 8 ] ]
        [ 3 x [ 5 6 ]   4 x [ 5 6 ] ]
        [     [ 7 8 ]       [ 7 8 ] ]

      = [ 1x5 1x6 2x5 2x6 ]
        [ 1x7 1x8 2x7 2x8 ]
        [ 3x5 3x6 4x5 4x6 ]
        [ 3x7 3x8 4x7 4x8 ]

      = [  5  6 10 12 ]
        [  7  8 14 16 ]
        [ 15 18 20 24 ]
        [ 21 24 28 32 ]

##</code><code>##

#!/usr/bin/perl
use strict;
use warnings;
use PDL;
use PDL::NiceSlice;
 
sub kron{
	my $A = shift;
	my $B = shift;
	my ($r0, $c0) = $A->dims;
	my ($r1, $c1) = $B->dims;
	my $kron = zeroes($r0 * $r1, $c0 * $c1);
	for(my $i = 0; $i < $r0; ++$i){
		for(my $j = 0; $j < $c0; ++$j){
			$kron(
				($i * $r1) : (($i + 1) * $r1 - 1),
				($j * $c1) : (($j + 1) * $c1 - 1)
			) .= $A($i,$j) * $B;
		}
	}
	return $kron;
}
##</code><code>##

[ 1 2 ]
[ 3 4 ]

##</code><code>##

[ 1 1 ]   [ 2 2 ]
[ 1 1 ]   [ 2 2 ]

[ 3 3 ]   [ 4 4 ]
[ 3 3 ]   [ 4 4 ]

##</code><code>##

[  5  6 ] [ 10 12 ]
[  7  8 ] [ 14 16 ]

[ 15 18 ] [ 20 24 ]
[ 21 24 ] [ 28 32 ]

##</code><code>##

sub kronecker_product {
    my ($x, $y) = @_;
    my ($x0, $x1) = $x->dims;
    my ($y0, $y1) = $y->dims;

    return (
        $y * $x->dummy(0, $y0)->dummy(1, $y1)
    )->xchg(1, 2)->reshape($x0 * $y0, $x1 * $y1)
}