multable( p, P, g ) =
{
P = Mod(1,p) *  P;
g = Mod( Mod(1,p) *  g, P );
d = poldegree( P );
n = p^d;

/*
    This version makes B[2] = Mod( 1, P ) which is wrong
    B = concat( [0], vector( n-1, i, g^(i-1) ) );

    The following line corrects the mistake
*/

B = concat( [0], vector( n-1, i, Mod( Mod(1,p), P ) * g^(i-1) ) );
R = vector( n );
for( i = 1, n,
	R[ 1 + subst( lift(lift( B[i] )), x, p )] = i-1;
);

for( a = 1, p^d,
 	ae = B[a];
	for ( b = 1, p^d,
		be = B[b];
		re = ae * be;
		print1( R[ 1 + subst( lift(lift( re )), x, p )], "\t" );
	);
	print;
);

}