Input and output not allowed to be same instance in composing?

[kaelingre]

Are the first two polynomials in the function fmpz_mpoly_compose_fmpz_mpoly allowed to be the same? I haven't found anything in the documentation that would tell me otherwise.
If so, the following (C++) code leads to random (incorrect) output:

fmpz_mpoly_ctx_struct ctx;
fmpz_mpoly_ctx_init(&ctx, 2, ORD_LEX);

fmpz_mpoly_struct n1;
fmpz_mpoly_init(&n1, &ctx);
fmpz_mpoly_gen(&n1, 0, &ctx);
fmpz_mpoly_struct n2;
fmpz_mpoly_init(&n2, &ctx);
fmpz_mpoly_gen(&n2, 1, &ctx);

fmpz_mpoly_struct p;
fmpz_mpoly_init(&p, &ctx);

fmpz_mpoly_add(&p, &n1, &n2, &ctx);
std::cout << "p = ";
fmpz_mpoly_print_pretty(&p, nullptr, &ctx);
std::cout << std::endl;

std::vector<fmpz_mpoly_struct> replData(2);
std::vector<fmpz_mpoly_struct*> repl(2);
for (std::size_t i = 0; i < 2; ++i) {
	auto& v = replData[i];
	repl[i] = &v;
	fmpz_mpoly_init(&v, &ctx);
	fmpz_mpoly_gen(&v, i, &ctx);
	//fmpz_mpoly_scalar_mul_si(&v, &v, 2, &ctx);
}

fmpz_mpoly_struct p2;
fmpz_mpoly_init(&p2, &ctx);

std::cout << "success: " << fmpz_mpoly_compose_fmpz_mpoly(&p, &p, repl.data(), &ctx, &ctx) << std::endl;

std::cout << "p = ";
fmpz_mpoly_print_pretty(&p, nullptr, &ctx);
std::cout << std::endl;

Leads for example to an output of the form:

p = x1+x2
success: 1
p = x1^1969438400*x2^1928772650+x1

Also note that uncommenting the line in the for loop will produce a correct output.
If one is not intended to have input=output it would be useful to know, for which functions this is forbidden.

[albinahlback]

Could you provide a minimal working example written in C?

[kaelingre]

Here's the simplest C example where this happened to me:

#include <flint/fmpz_mpoly.h>

int main() {
	fmpz_mpoly_ctx_struct ctx;
	fmpz_mpoly_ctx_init(&ctx, 2, ORD_LEX);

	fmpz_mpoly_struct x1;
	fmpz_mpoly_init(&x1, &ctx);
	fmpz_mpoly_gen(&x1, 0, &ctx);
	fmpz_mpoly_struct x2;
	fmpz_mpoly_init(&x2, &ctx);
	fmpz_mpoly_gen(&x2, 1, &ctx);

	//p = x1 + x2
	fmpz_mpoly_struct p;
	fmpz_mpoly_init(&p, &ctx);
	fmpz_mpoly_add(&p, &x1, &x2, &ctx);
	printf("p = ");
	fmpz_mpoly_print_pretty(&p, NULL, &ctx);
	printf("\n");

	//make the trivial replacement x1->x1, x2->x2
	fmpz_mpoly_struct replData[2];
	fmpz_mpoly_struct* repl[2];
	repl[0] = &replData[0];
	repl[1] = &replData[1];
	fmpz_mpoly_init(repl[0], &ctx);
	fmpz_mpoly_gen(repl[0], 0, &ctx);
	fmpz_mpoly_init(repl[1], &ctx);
	fmpz_mpoly_gen(repl[1], 0, &ctx);
		
	printf("success: %d\n", fmpz_mpoly_compose_fmpz_mpoly(&p, &p, repl, &ctx, &ctx));

	printf("p = ");
	fmpz_mpoly_print_pretty(&p, NULL, &ctx);
	printf("\n");

	return 0;
}

which outputs

p = x1+x2
success: 1
p = x1^3305670466+x1

If this is indeed not the intended behavior I'd also be happy to track down the issue.

[fredrik-johansson]

Yeah, it looks like there should be an aliasing check here.

[kaelingre]

Alright, the problem is that the function _fmpz_mpoly_compose_mat(A,B,...) requires A!=B. There's an assert which would have caught the above problem, had I turned on asserts.