Català-Valencià – Catalan中文 – Chinese (Simplified)中文 – Chinese (Traditional)Česky – CzechDansk – DanishNederlands – DutchEnglish – EnglishSuomi – FinnishFrançais – FrenchDeutsch – Germanעברית – Hebrewहिंदी – HindiMagyar – HungarianBahasa Indonesia – IndonesianItaliano – Italian日本語 – Japanese한국어 – KoreanМакедонски – Macedonianमराठी – MarathiNorsk – NorwegianPolski – PolishPortuguês – PortuguesePortuguês – Portuguese (Brazil)Русский – RussianSlovenčina – SlovakSlovenščina – SlovenianEspañol – SpanishSvenska – SwedishTürkçe – TurkishУкраїнська – UkrainianOëzbekcha – Uzbek

nexmon – Blame information for rev 1

Rev	Author	Line No.	Line
1	office	1	/*
		2	* Copyright 2010 INRIA Saclay
		3	*
		4	* Use of this software is governed by the GNU LGPLv2.1 license
		5	*
		6	* Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
		7	* Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
		8	* 91893 Orsay, France
		9	*/
		10
		11	#include <stdlib.h>
		12	#define ISL_DIM_H
		13	#include <isl_ctx_private.h>
		14	#include <isl_map_private.h>
		15	#include <isl_factorization.h>
		16	#include <isl/lp.h>
		17	#include <isl/seq.h>
		18	#include <isl_union_map_private.h>
		19	#include <isl_constraint_private.h>
		20	#include <isl_polynomial_private.h>
		21	#include <isl_point_private.h>
		22	#include <isl_space_private.h>
		23	#include <isl_mat_private.h>
		24	#include <isl_range.h>
		25	#include <isl_local_space_private.h>
		26	#include <isl_aff_private.h>
		27	#include <isl_config.h>
		28
		29	static unsigned pos(__isl_keep isl_space *dim, enum isl_dim_type type)
		30	{
		31	switch (type) {
		32	case isl_dim_param: return 0;
		33	case isl_dim_in: return dim->nparam;
		34	case isl_dim_out: return dim->nparam + dim->n_in;
		35	default: return 0;
		36	}
		37	}
		38
		39	int isl_upoly_is_cst(__isl_keep struct isl_upoly *up)
		40	{
		41	if (!up)
		42	return -1;
		43
		44	return up->var < 0;
		45	}
		46
		47	__isl_keep struct isl_upoly_cst *isl_upoly_as_cst(__isl_keep struct isl_upoly *up)
		48	{
		49	if (!up)
		50	return NULL;
		51
		52	isl_assert(up->ctx, up->var < 0, return NULL);
		53
		54	return (struct isl_upoly_cst *)up;
		55	}
		56
		57	__isl_keep struct isl_upoly_rec *isl_upoly_as_rec(__isl_keep struct isl_upoly *up)
		58	{
		59	if (!up)
		60	return NULL;
		61
		62	isl_assert(up->ctx, up->var >= 0, return NULL);
		63
		64	return (struct isl_upoly_rec *)up;
		65	}
		66
		67	int isl_upoly_is_equal(__isl_keep struct isl_upoly *up1,
		68	__isl_keep struct isl_upoly *up2)
		69	{
		70	int i;
		71	struct isl_upoly_rec *rec1, *rec2;
		72
		73	if (!up1 || !up2)
		74	return -1;
		75	if (up1 == up2)
		76	return 1;
		77	if (up1->var != up2->var)
		78	return 0;
		79	if (isl_upoly_is_cst(up1)) {
		80	struct isl_upoly_cst *cst1, *cst2;
		81	cst1 = isl_upoly_as_cst(up1);
		82	cst2 = isl_upoly_as_cst(up2);
		83	if (!cst1 || !cst2)
		84	return -1;
		85	return isl_int_eq(cst1->n, cst2->n) &&
		86	isl_int_eq(cst1->d, cst2->d);
		87	}
		88
		89	rec1 = isl_upoly_as_rec(up1);
		90	rec2 = isl_upoly_as_rec(up2);
		91	if (!rec1 || !rec2)
		92	return -1;
		93
		94	if (rec1->n != rec2->n)
		95	return 0;
		96
		97	for (i = 0; i < rec1->n; ++i) {
		98	int eq = isl_upoly_is_equal(rec1->p[i], rec2->p[i]);
		99	if (eq < 0 || !eq)
		100	return eq;
		101	}
		102
		103	return 1;
		104	}
		105
		106	int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
		107	{
		108	struct isl_upoly_cst *cst;
		109
		110	if (!up)
		111	return -1;
		112	if (!isl_upoly_is_cst(up))
		113	return 0;
		114
		115	cst = isl_upoly_as_cst(up);
		116	if (!cst)
		117	return -1;
		118
		119	return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
		120	}
		121
		122	int isl_upoly_sgn(__isl_keep struct isl_upoly *up)
		123	{
		124	struct isl_upoly_cst *cst;
		125
		126	if (!up)
		127	return 0;
		128	if (!isl_upoly_is_cst(up))
		129	return 0;
		130
		131	cst = isl_upoly_as_cst(up);
		132	if (!cst)
		133	return 0;
		134
		135	return isl_int_sgn(cst->n);
		136	}
		137
		138	int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
		139	{
		140	struct isl_upoly_cst *cst;
		141
		142	if (!up)
		143	return -1;
		144	if (!isl_upoly_is_cst(up))
		145	return 0;
		146
		147	cst = isl_upoly_as_cst(up);
		148	if (!cst)
		149	return -1;
		150
		151	return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
		152	}
		153
		154	int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
		155	{
		156	struct isl_upoly_cst *cst;
		157
		158	if (!up)
		159	return -1;
		160	if (!isl_upoly_is_cst(up))
		161	return 0;
		162
		163	cst = isl_upoly_as_cst(up);
		164	if (!cst)
		165	return -1;
		166
		167	return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
		168	}
		169
		170	int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
		171	{
		172	struct isl_upoly_cst *cst;
		173
		174	if (!up)
		175	return -1;
		176	if (!isl_upoly_is_cst(up))
		177	return 0;
		178
		179	cst = isl_upoly_as_cst(up);
		180	if (!cst)
		181	return -1;
		182
		183	return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
		184	}
		185
		186	int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
		187	{
		188	struct isl_upoly_cst *cst;
		189
		190	if (!up)
		191	return -1;
		192	if (!isl_upoly_is_cst(up))
		193	return 0;
		194
		195	cst = isl_upoly_as_cst(up);
		196	if (!cst)
		197	return -1;
		198
		199	return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
		200	}
		201
		202	int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
		203	{
		204	struct isl_upoly_cst *cst;
		205
		206	if (!up)
		207	return -1;
		208	if (!isl_upoly_is_cst(up))
		209	return 0;
		210
		211	cst = isl_upoly_as_cst(up);
		212	if (!cst)
		213	return -1;
		214
		215	return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
		216	}
		217
		218	__isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
		219	{
		220	struct isl_upoly_cst *cst;
		221
		222	cst = isl_alloc_type(ctx, struct isl_upoly_cst);
		223	if (!cst)
		224	return NULL;
		225
		226	cst->up.ref = 1;
		227	cst->up.ctx = ctx;
		228	isl_ctx_ref(ctx);
		229	cst->up.var = -1;
		230
		231	isl_int_init(cst->n);
		232	isl_int_init(cst->d);
		233
		234	return cst;
		235	}
		236
		237	__isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
		238	{
		239	struct isl_upoly_cst *cst;
		240
		241	cst = isl_upoly_cst_alloc(ctx);
		242	if (!cst)
		243	return NULL;
		244
		245	isl_int_set_si(cst->n, 0);
		246	isl_int_set_si(cst->d, 1);
		247
		248	return &cst->up;
		249	}
		250
		251	__isl_give struct isl_upoly *isl_upoly_one(struct isl_ctx *ctx)
		252	{
		253	struct isl_upoly_cst *cst;
		254
		255	cst = isl_upoly_cst_alloc(ctx);
		256	if (!cst)
		257	return NULL;
		258
		259	isl_int_set_si(cst->n, 1);
		260	isl_int_set_si(cst->d, 1);
		261
		262	return &cst->up;
		263	}
		264
		265	__isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
		266	{
		267	struct isl_upoly_cst *cst;
		268
		269	cst = isl_upoly_cst_alloc(ctx);
		270	if (!cst)
		271	return NULL;
		272
		273	isl_int_set_si(cst->n, 1);
		274	isl_int_set_si(cst->d, 0);
		275
		276	return &cst->up;
		277	}
		278
		279	__isl_give struct isl_upoly *isl_upoly_neginfty(struct isl_ctx *ctx)
		280	{
		281	struct isl_upoly_cst *cst;
		282
		283	cst = isl_upoly_cst_alloc(ctx);
		284	if (!cst)
		285	return NULL;
		286
		287	isl_int_set_si(cst->n, -1);
		288	isl_int_set_si(cst->d, 0);
		289
		290	return &cst->up;
		291	}
		292
		293	__isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
		294	{
		295	struct isl_upoly_cst *cst;
		296
		297	cst = isl_upoly_cst_alloc(ctx);
		298	if (!cst)
		299	return NULL;
		300
		301	isl_int_set_si(cst->n, 0);
		302	isl_int_set_si(cst->d, 0);
		303
		304	return &cst->up;
		305	}
		306
		307	__isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
		308	isl_int n, isl_int d)
		309	{
		310	struct isl_upoly_cst *cst;
		311
		312	cst = isl_upoly_cst_alloc(ctx);
		313	if (!cst)
		314	return NULL;
		315
		316	isl_int_set(cst->n, n);
		317	isl_int_set(cst->d, d);
		318
		319	return &cst->up;
		320	}
		321
		322	__isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
		323	int var, int size)
		324	{
		325	struct isl_upoly_rec *rec;
		326
		327	isl_assert(ctx, var >= 0, return NULL);
		328	isl_assert(ctx, size >= 0, return NULL);
		329	rec = isl_calloc(ctx, struct isl_upoly_rec,
		330	sizeof(struct isl_upoly_rec) +
		331	size * sizeof(struct isl_upoly *));
		332	if (!rec)
		333	return NULL;
		334
		335	rec->up.ref = 1;
		336	rec->up.ctx = ctx;
		337	isl_ctx_ref(ctx);
		338	rec->up.var = var;
		339
		340	rec->n = 0;
		341	rec->size = size;
		342
		343	return rec;
		344	}
		345
		346	__isl_give isl_qpolynomial *isl_qpolynomial_reset_domain_space(
		347	__isl_take isl_qpolynomial *qp, __isl_take isl_space *dim)
		348	{
		349	qp = isl_qpolynomial_cow(qp);
		350	if (!qp || !dim)
		351	goto error;
		352
		353	isl_space_free(qp->dim);
		354	qp->dim = dim;
		355
		356	return qp;
		357	error:
		358	isl_qpolynomial_free(qp);
		359	isl_space_free(dim);
		360	return NULL;
		361	}
		362
		363	/* Reset the space of "qp". This function is called from isl_pw_templ.c
		364	* and doesn't know if the space of an element object is represented
		365	* directly or through its domain. It therefore passes along both.
		366	*/
		367	__isl_give isl_qpolynomial *isl_qpolynomial_reset_space_and_domain(
		368	__isl_take isl_qpolynomial *qp, __isl_take isl_space *space,
		369	__isl_take isl_space *domain)
		370	{
		371	isl_space_free(space);
		372	return isl_qpolynomial_reset_domain_space(qp, domain);
		373	}
		374
		375	isl_ctx *isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial *qp)
		376	{
		377	return qp ? qp->dim->ctx : NULL;
		378	}
		379
		380	__isl_give isl_space *isl_qpolynomial_get_domain_space(
		381	__isl_keep isl_qpolynomial *qp)
		382	{
		383	return qp ? isl_space_copy(qp->dim) : NULL;
		384	}
		385
		386	__isl_give isl_space *isl_qpolynomial_get_space(__isl_keep isl_qpolynomial *qp)
		387	{
		388	isl_space *space;
		389	if (!qp)
		390	return NULL;
		391	space = isl_space_copy(qp->dim);
		392	space = isl_space_from_domain(space);
		393	space = isl_space_add_dims(space, isl_dim_out, 1);
		394	return space;
		395	}
		396
		397	/* Externally, an isl_qpolynomial has a map space, but internally, the
		398	* ls field corresponds to the domain of that space.
		399	*/
		400	unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial *qp,
		401	enum isl_dim_type type)
		402	{
		403	if (!qp)
		404	return 0;
		405	if (type == isl_dim_out)
		406	return 1;
		407	if (type == isl_dim_in)
		408	type = isl_dim_set;
		409	return isl_space_dim(qp->dim, type);
		410	}
		411
		412	int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
		413	{
		414	return qp ? isl_upoly_is_zero(qp->upoly) : -1;
		415	}
		416
		417	int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
		418	{
		419	return qp ? isl_upoly_is_one(qp->upoly) : -1;
		420	}
		421
		422	int isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp)
		423	{
		424	return qp ? isl_upoly_is_nan(qp->upoly) : -1;
		425	}
		426
		427	int isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp)
		428	{
		429	return qp ? isl_upoly_is_infty(qp->upoly) : -1;
		430	}
		431
		432	int isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp)
		433	{
		434	return qp ? isl_upoly_is_neginfty(qp->upoly) : -1;
		435	}
		436
		437	int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp)
		438	{
		439	return qp ? isl_upoly_sgn(qp->upoly) : 0;
		440	}
		441
		442	static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
		443	{
		444	isl_int_clear(cst->n);
		445	isl_int_clear(cst->d);
		446	}
		447
		448	static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
		449	{
		450	int i;
		451
		452	for (i = 0; i < rec->n; ++i)
		453	isl_upoly_free(rec->p[i]);
		454	}
		455
		456	__isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
		457	{
		458	if (!up)
		459	return NULL;
		460
		461	up->ref++;
		462	return up;
		463	}
		464
		465	__isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
		466	{
		467	struct isl_upoly_cst *cst;
		468	struct isl_upoly_cst *dup;
		469
		470	cst = isl_upoly_as_cst(up);
		471	if (!cst)
		472	return NULL;
		473
		474	dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
		475	if (!dup)
		476	return NULL;
		477	isl_int_set(dup->n, cst->n);
		478	isl_int_set(dup->d, cst->d);
		479
		480	return &dup->up;
		481	}
		482
		483	__isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
		484	{
		485	int i;
		486	struct isl_upoly_rec *rec;
		487	struct isl_upoly_rec *dup;
		488
		489	rec = isl_upoly_as_rec(up);
		490	if (!rec)
		491	return NULL;
		492
		493	dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
		494	if (!dup)
		495	return NULL;
		496
		497	for (i = 0; i < rec->n; ++i) {
		498	dup->p[i] = isl_upoly_copy(rec->p[i]);
		499	if (!dup->p[i])
		500	goto error;
		501	dup->n++;
		502	}
		503
		504	return &dup->up;
		505	error:
		506	isl_upoly_free(&dup->up);
		507	return NULL;
		508	}
		509
		510	__isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
		511	{
		512	if (!up)
		513	return NULL;
		514
		515	if (isl_upoly_is_cst(up))
		516	return isl_upoly_dup_cst(up);
		517	else
		518	return isl_upoly_dup_rec(up);
		519	}
		520
		521	__isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
		522	{
		523	if (!up)
		524	return NULL;
		525
		526	if (up->ref == 1)
		527	return up;
		528	up->ref--;
		529	return isl_upoly_dup(up);
		530	}
		531
		532	void isl_upoly_free(__isl_take struct isl_upoly *up)
		533	{
		534	if (!up)
		535	return;
		536
		537	if (--up->ref > 0)
		538	return;
		539
		540	if (up->var < 0)
		541	upoly_free_cst((struct isl_upoly_cst *)up);
		542	else
		543	upoly_free_rec((struct isl_upoly_rec *)up);
		544
		545	isl_ctx_deref(up->ctx);
		546	free(up);
		547	}
		548
		549	static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
		550	{
		551	isl_int gcd;
		552
		553	isl_int_init(gcd);
		554	isl_int_gcd(gcd, cst->n, cst->d);
		555	if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
		556	isl_int_divexact(cst->n, cst->n, gcd);
		557	isl_int_divexact(cst->d, cst->d, gcd);
		558	}
		559	isl_int_clear(gcd);
		560	}
		561
		562	__isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
		563	__isl_take struct isl_upoly *up2)
		564	{
		565	struct isl_upoly_cst *cst1;
		566	struct isl_upoly_cst *cst2;
		567
		568	up1 = isl_upoly_cow(up1);
		569	if (!up1 || !up2)
		570	goto error;
		571
		572	cst1 = isl_upoly_as_cst(up1);
		573	cst2 = isl_upoly_as_cst(up2);
		574
		575	if (isl_int_eq(cst1->d, cst2->d))
		576	isl_int_add(cst1->n, cst1->n, cst2->n);
		577	else {
		578	isl_int_mul(cst1->n, cst1->n, cst2->d);
		579	isl_int_addmul(cst1->n, cst2->n, cst1->d);
		580	isl_int_mul(cst1->d, cst1->d, cst2->d);
		581	}
		582
		583	isl_upoly_cst_reduce(cst1);
		584
		585	isl_upoly_free(up2);
		586	return up1;
		587	error:
		588	isl_upoly_free(up1);
		589	isl_upoly_free(up2);
		590	return NULL;
		591	}
		592
		593	static __isl_give struct isl_upoly *replace_by_zero(
		594	__isl_take struct isl_upoly *up)
		595	{
		596	struct isl_ctx *ctx;
		597
		598	if (!up)
		599	return NULL;
		600	ctx = up->ctx;
		601	isl_upoly_free(up);
		602	return isl_upoly_zero(ctx);
		603	}
		604
		605	static __isl_give struct isl_upoly *replace_by_constant_term(
		606	__isl_take struct isl_upoly *up)
		607	{
		608	struct isl_upoly_rec *rec;
		609	struct isl_upoly *cst;
		610
		611	if (!up)
		612	return NULL;
		613
		614	rec = isl_upoly_as_rec(up);
		615	if (!rec)
		616	goto error;
		617	cst = isl_upoly_copy(rec->p[0]);
		618	isl_upoly_free(up);
		619	return cst;
		620	error:
		621	isl_upoly_free(up);
		622	return NULL;
		623	}
		624
		625	__isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
		626	__isl_take struct isl_upoly *up2)
		627	{
		628	int i;
		629	struct isl_upoly_rec *rec1, *rec2;
		630
		631	if (!up1 || !up2)
		632	goto error;
		633
		634	if (isl_upoly_is_nan(up1)) {
		635	isl_upoly_free(up2);
		636	return up1;
		637	}
		638
		639	if (isl_upoly_is_nan(up2)) {
		640	isl_upoly_free(up1);
		641	return up2;
		642	}
		643
		644	if (isl_upoly_is_zero(up1)) {
		645	isl_upoly_free(up1);
		646	return up2;
		647	}
		648
		649	if (isl_upoly_is_zero(up2)) {
		650	isl_upoly_free(up2);
		651	return up1;
		652	}
		653
		654	if (up1->var < up2->var)
		655	return isl_upoly_sum(up2, up1);
		656
		657	if (up2->var < up1->var) {
		658	struct isl_upoly_rec *rec;
		659	if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
		660	isl_upoly_free(up1);
		661	return up2;
		662	}
		663	up1 = isl_upoly_cow(up1);
		664	rec = isl_upoly_as_rec(up1);
		665	if (!rec)
		666	goto error;
		667	rec->p[0] = isl_upoly_sum(rec->p[0], up2);
		668	if (rec->n == 1)
		669	up1 = replace_by_constant_term(up1);
		670	return up1;
		671	}
		672
		673	if (isl_upoly_is_cst(up1))
		674	return isl_upoly_sum_cst(up1, up2);
		675
		676	rec1 = isl_upoly_as_rec(up1);
		677	rec2 = isl_upoly_as_rec(up2);
		678	if (!rec1 || !rec2)
		679	goto error;
		680
		681	if (rec1->n < rec2->n)
		682	return isl_upoly_sum(up2, up1);
		683
		684	up1 = isl_upoly_cow(up1);
		685	rec1 = isl_upoly_as_rec(up1);
		686	if (!rec1)
		687	goto error;
		688
		689	for (i = rec2->n - 1; i >= 0; --i) {
		690	rec1->p[i] = isl_upoly_sum(rec1->p[i],
		691	isl_upoly_copy(rec2->p[i]));
		692	if (!rec1->p[i])
		693	goto error;
		694	if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
		695	isl_upoly_free(rec1->p[i]);
		696	rec1->n--;
		697	}
		698	}
		699
		700	if (rec1->n == 0)
		701	up1 = replace_by_zero(up1);
		702	else if (rec1->n == 1)
		703	up1 = replace_by_constant_term(up1);
		704
		705	isl_upoly_free(up2);
		706
		707	return up1;
		708	error:
		709	isl_upoly_free(up1);
		710	isl_upoly_free(up2);
		711	return NULL;
		712	}
		713
		714	__isl_give struct isl_upoly *isl_upoly_cst_add_isl_int(
		715	__isl_take struct isl_upoly *up, isl_int v)
		716	{
		717	struct isl_upoly_cst *cst;
		718
		719	up = isl_upoly_cow(up);
		720	if (!up)
		721	return NULL;
		722
		723	cst = isl_upoly_as_cst(up);
		724
		725	isl_int_addmul(cst->n, cst->d, v);
		726
		727	return up;
		728	}
		729
		730	__isl_give struct isl_upoly *isl_upoly_add_isl_int(
		731	__isl_take struct isl_upoly *up, isl_int v)
		732	{
		733	struct isl_upoly_rec *rec;
		734
		735	if (!up)
		736	return NULL;
		737
		738	if (isl_upoly_is_cst(up))
		739	return isl_upoly_cst_add_isl_int(up, v);
		740
		741	up = isl_upoly_cow(up);
		742	rec = isl_upoly_as_rec(up);
		743	if (!rec)
		744	goto error;
		745
		746	rec->p[0] = isl_upoly_add_isl_int(rec->p[0], v);
		747	if (!rec->p[0])
		748	goto error;
		749
		750	return up;
		751	error:
		752	isl_upoly_free(up);
		753	return NULL;
		754	}
		755
		756	__isl_give struct isl_upoly *isl_upoly_cst_mul_isl_int(
		757	__isl_take struct isl_upoly *up, isl_int v)
		758	{
		759	struct isl_upoly_cst *cst;
		760
		761	if (isl_upoly_is_zero(up))
		762	return up;
		763
		764	up = isl_upoly_cow(up);
		765	if (!up)
		766	return NULL;
		767
		768	cst = isl_upoly_as_cst(up);
		769
		770	isl_int_mul(cst->n, cst->n, v);
		771
		772	return up;
		773	}
		774
		775	__isl_give struct isl_upoly *isl_upoly_mul_isl_int(
		776	__isl_take struct isl_upoly *up, isl_int v)
		777	{
		778	int i;
		779	struct isl_upoly_rec *rec;
		780
		781	if (!up)
		782	return NULL;
		783
		784	if (isl_upoly_is_cst(up))
		785	return isl_upoly_cst_mul_isl_int(up, v);
		786
		787	up = isl_upoly_cow(up);
		788	rec = isl_upoly_as_rec(up);
		789	if (!rec)
		790	goto error;
		791
		792	for (i = 0; i < rec->n; ++i) {
		793	rec->p[i] = isl_upoly_mul_isl_int(rec->p[i], v);
		794	if (!rec->p[i])
		795	goto error;
		796	}
		797
		798	return up;
		799	error:
		800	isl_upoly_free(up);
		801	return NULL;
		802	}
		803
		804	__isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
		805	__isl_take struct isl_upoly *up2)
		806	{
		807	struct isl_upoly_cst *cst1;
		808	struct isl_upoly_cst *cst2;
		809
		810	up1 = isl_upoly_cow(up1);
		811	if (!up1 || !up2)
		812	goto error;
		813
		814	cst1 = isl_upoly_as_cst(up1);
		815	cst2 = isl_upoly_as_cst(up2);
		816
		817	isl_int_mul(cst1->n, cst1->n, cst2->n);
		818	isl_int_mul(cst1->d, cst1->d, cst2->d);
		819
		820	isl_upoly_cst_reduce(cst1);
		821
		822	isl_upoly_free(up2);
		823	return up1;
		824	error:
		825	isl_upoly_free(up1);
		826	isl_upoly_free(up2);
		827	return NULL;
		828	}
		829
		830	__isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
		831	__isl_take struct isl_upoly *up2)
		832	{
		833	struct isl_upoly_rec *rec1;
		834	struct isl_upoly_rec *rec2;
		835	struct isl_upoly_rec *res = NULL;
		836	int i, j;
		837	int size;
		838
		839	rec1 = isl_upoly_as_rec(up1);
		840	rec2 = isl_upoly_as_rec(up2);
		841	if (!rec1 || !rec2)
		842	goto error;
		843	size = rec1->n + rec2->n - 1;
		844	res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
		845	if (!res)
		846	goto error;
		847
		848	for (i = 0; i < rec1->n; ++i) {
		849	res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
		850	isl_upoly_copy(rec1->p[i]));
		851	if (!res->p[i])
		852	goto error;
		853	res->n++;
		854	}
		855	for (; i < size; ++i) {
		856	res->p[i] = isl_upoly_zero(up1->ctx);
		857	if (!res->p[i])
		858	goto error;
		859	res->n++;
		860	}
		861	for (i = 0; i < rec1->n; ++i) {
		862	for (j = 1; j < rec2->n; ++j) {
		863	struct isl_upoly *up;
		864	up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
		865	isl_upoly_copy(rec1->p[i]));
		866	res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
		867	if (!res->p[i + j])
		868	goto error;
		869	}
		870	}
		871
		872	isl_upoly_free(up1);
		873	isl_upoly_free(up2);
		874
		875	return &res->up;
		876	error:
		877	isl_upoly_free(up1);
		878	isl_upoly_free(up2);
		879	isl_upoly_free(&res->up);
		880	return NULL;
		881	}
		882
		883	__isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
		884	__isl_take struct isl_upoly *up2)
		885	{
		886	if (!up1 || !up2)
		887	goto error;
		888
		889	if (isl_upoly_is_nan(up1)) {
		890	isl_upoly_free(up2);
		891	return up1;
		892	}
		893
		894	if (isl_upoly_is_nan(up2)) {
		895	isl_upoly_free(up1);
		896	return up2;
		897	}
		898
		899	if (isl_upoly_is_zero(up1)) {
		900	isl_upoly_free(up2);
		901	return up1;
		902	}
		903
		904	if (isl_upoly_is_zero(up2)) {
		905	isl_upoly_free(up1);
		906	return up2;
		907	}
		908
		909	if (isl_upoly_is_one(up1)) {
		910	isl_upoly_free(up1);
		911	return up2;
		912	}
		913
		914	if (isl_upoly_is_one(up2)) {
		915	isl_upoly_free(up2);
		916	return up1;
		917	}
		918
		919	if (up1->var < up2->var)
		920	return isl_upoly_mul(up2, up1);
		921
		922	if (up2->var < up1->var) {
		923	int i;
		924	struct isl_upoly_rec *rec;
		925	if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
		926	isl_ctx *ctx = up1->ctx;
		927	isl_upoly_free(up1);
		928	isl_upoly_free(up2);
		929	return isl_upoly_nan(ctx);
		930	}
		931	up1 = isl_upoly_cow(up1);
		932	rec = isl_upoly_as_rec(up1);
		933	if (!rec)
		934	goto error;
		935
		936	for (i = 0; i < rec->n; ++i) {
		937	rec->p[i] = isl_upoly_mul(rec->p[i],
		938	isl_upoly_copy(up2));
		939	if (!rec->p[i])
		940	goto error;
		941	}
		942	isl_upoly_free(up2);
		943	return up1;
		944	}
		945
		946	if (isl_upoly_is_cst(up1))
		947	return isl_upoly_mul_cst(up1, up2);
		948
		949	return isl_upoly_mul_rec(up1, up2);
		950	error:
		951	isl_upoly_free(up1);
		952	isl_upoly_free(up2);
		953	return NULL;
		954	}
		955
		956	__isl_give struct isl_upoly *isl_upoly_pow(__isl_take struct isl_upoly *up,
		957	unsigned power)
		958	{
		959	struct isl_upoly *res;
		960
		961	if (!up)
		962	return NULL;
		963	if (power == 1)
		964	return up;
		965
		966	if (power % 2)
		967	res = isl_upoly_copy(up);
		968	else
		969	res = isl_upoly_one(up->ctx);
		970
		971	while (power >>= 1) {
		972	up = isl_upoly_mul(up, isl_upoly_copy(up));
		973	if (power % 2)
		974	res = isl_upoly_mul(res, isl_upoly_copy(up));
		975	}
		976
		977	isl_upoly_free(up);
		978	return res;
		979	}
		980
		981	__isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_space *dim,
		982	unsigned n_div, __isl_take struct isl_upoly *up)
		983	{
		984	struct isl_qpolynomial *qp = NULL;
		985	unsigned total;
		986
		987	if (!dim || !up)
		988	goto error;
		989
		990	if (!isl_space_is_set(dim))
		991	isl_die(isl_space_get_ctx(dim), isl_error_invalid,
		992	"domain of polynomial should be a set", goto error);
		993
		994	total = isl_space_dim(dim, isl_dim_all);
		995
		996	qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
		997	if (!qp)
		998	goto error;
		999
		1000	qp->ref = 1;
		1001	qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
		1002	if (!qp->div)
		1003	goto error;
		1004
		1005	qp->dim = dim;
		1006	qp->upoly = up;
		1007
		1008	return qp;
		1009	error:
		1010	isl_space_free(dim);
		1011	isl_upoly_free(up);
		1012	isl_qpolynomial_free(qp);
		1013	return NULL;
		1014	}
		1015
		1016	__isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
		1017	{
		1018	if (!qp)
		1019	return NULL;
		1020
		1021	qp->ref++;
		1022	return qp;
		1023	}
		1024
		1025	__isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
		1026	{
		1027	struct isl_qpolynomial *dup;
		1028
		1029	if (!qp)
		1030	return NULL;
		1031
		1032	dup = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row,
		1033	isl_upoly_copy(qp->upoly));
		1034	if (!dup)
		1035	return NULL;
		1036	isl_mat_free(dup->div);
		1037	dup->div = isl_mat_copy(qp->div);
		1038	if (!dup->div)
		1039	goto error;
		1040
		1041	return dup;
		1042	error:
		1043	isl_qpolynomial_free(dup);
		1044	return NULL;
		1045	}
		1046
		1047	__isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
		1048	{
		1049	if (!qp)
		1050	return NULL;
		1051
		1052	if (qp->ref == 1)
		1053	return qp;
		1054	qp->ref--;
		1055	return isl_qpolynomial_dup(qp);
		1056	}
		1057
		1058	void *isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
		1059	{
		1060	if (!qp)
		1061	return NULL;
		1062
		1063	if (--qp->ref > 0)
		1064	return NULL;
		1065
		1066	isl_space_free(qp->dim);
		1067	isl_mat_free(qp->div);
		1068	isl_upoly_free(qp->upoly);
		1069
		1070	free(qp);
		1071	return NULL;
		1072	}
		1073
		1074	__isl_give struct isl_upoly *isl_upoly_var_pow(isl_ctx *ctx, int pos, int power)
		1075	{
		1076	int i;
		1077	struct isl_upoly_rec *rec;
		1078	struct isl_upoly_cst *cst;
		1079
		1080	rec = isl_upoly_alloc_rec(ctx, pos, 1 + power);
		1081	if (!rec)
		1082	return NULL;
		1083	for (i = 0; i < 1 + power; ++i) {
		1084	rec->p[i] = isl_upoly_zero(ctx);
		1085	if (!rec->p[i])
		1086	goto error;
		1087	rec->n++;
		1088	}
		1089	cst = isl_upoly_as_cst(rec->p[power]);
		1090	isl_int_set_si(cst->n, 1);
		1091
		1092	return &rec->up;
		1093	error:
		1094	isl_upoly_free(&rec->up);
		1095	return NULL;
		1096	}
		1097
		1098	/* r array maps original positions to new positions.
		1099	*/
		1100	static __isl_give struct isl_upoly *reorder(__isl_take struct isl_upoly *up,
		1101	int *r)
		1102	{
		1103	int i;
		1104	struct isl_upoly_rec *rec;
		1105	struct isl_upoly *base;
		1106	struct isl_upoly *res;
		1107
		1108	if (isl_upoly_is_cst(up))
		1109	return up;
		1110
		1111	rec = isl_upoly_as_rec(up);
		1112	if (!rec)
		1113	goto error;
		1114
		1115	isl_assert(up->ctx, rec->n >= 1, goto error);
		1116
		1117	base = isl_upoly_var_pow(up->ctx, r[up->var], 1);
		1118	res = reorder(isl_upoly_copy(rec->p[rec->n - 1]), r);
		1119
		1120	for (i = rec->n - 2; i >= 0; --i) {
		1121	res = isl_upoly_mul(res, isl_upoly_copy(base));
		1122	res = isl_upoly_sum(res, reorder(isl_upoly_copy(rec->p[i]), r));
		1123	}
		1124
		1125	isl_upoly_free(base);
		1126	isl_upoly_free(up);
		1127
		1128	return res;
		1129	error:
		1130	isl_upoly_free(up);
		1131	return NULL;
		1132	}
		1133
		1134	static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
		1135	{
		1136	int n_row, n_col;
		1137	int equal;
		1138
		1139	isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
		1140	div1->n_col >= div2->n_col, return -1);
		1141
		1142	if (div1->n_row == div2->n_row)
		1143	return isl_mat_is_equal(div1, div2);
		1144
		1145	n_row = div1->n_row;
		1146	n_col = div1->n_col;
		1147	div1->n_row = div2->n_row;
		1148	div1->n_col = div2->n_col;
		1149
		1150	equal = isl_mat_is_equal(div1, div2);
		1151
		1152	div1->n_row = n_row;
		1153	div1->n_col = n_col;
		1154
		1155	return equal;
		1156	}
		1157
		1158	static int cmp_row(__isl_keep isl_mat *div, int i, int j)
		1159	{
		1160	int li, lj;
		1161
		1162	li = isl_seq_last_non_zero(div->row[i], div->n_col);
		1163	lj = isl_seq_last_non_zero(div->row[j], div->n_col);
		1164
		1165	if (li != lj)
		1166	return li - lj;
		1167
		1168	return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
		1169	}
		1170
		1171	struct isl_div_sort_info {
		1172	isl_mat *div;
		1173	int row;
		1174	};
		1175
		1176	static int div_sort_cmp(const void *p1, const void *p2)
		1177	{
		1178	const struct isl_div_sort_info *i1, *i2;
		1179	i1 = (const struct isl_div_sort_info *) p1;
		1180	i2 = (const struct isl_div_sort_info *) p2;
		1181
		1182	return cmp_row(i1->div, i1->row, i2->row);
		1183	}
		1184
		1185	/* Sort divs and remove duplicates.
		1186	*/
		1187	static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
		1188	{
		1189	int i;
		1190	int skip;
		1191	int len;
		1192	struct isl_div_sort_info *array = NULL;
		1193	int *pos = NULL, *at = NULL;
		1194	int *reordering = NULL;
		1195	unsigned div_pos;
		1196
		1197	if (!qp)
		1198	return NULL;
		1199	if (qp->div->n_row <= 1)
		1200	return qp;
		1201
		1202	div_pos = isl_space_dim(qp->dim, isl_dim_all);
		1203
		1204	array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info,
		1205	qp->div->n_row);
		1206	pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
		1207	at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
		1208	len = qp->div->n_col - 2;
		1209	reordering = isl_alloc_array(qp->div->ctx, int, len);
		1210	if (!array || !pos || !at || !reordering)
		1211	goto error;
		1212
		1213	for (i = 0; i < qp->div->n_row; ++i) {
		1214	array[i].div = qp->div;
		1215	array[i].row = i;
		1216	pos[i] = i;
		1217	at[i] = i;
		1218	}
		1219
		1220	qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
		1221	div_sort_cmp);
		1222
		1223	for (i = 0; i < div_pos; ++i)
		1224	reordering[i] = i;
		1225
		1226	for (i = 0; i < qp->div->n_row; ++i) {
		1227	if (pos[array[i].row] == i)
		1228	continue;
		1229	qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
		1230	pos[at[i]] = pos[array[i].row];
		1231	at[pos[array[i].row]] = at[i];
		1232	at[i] = array[i].row;
		1233	pos[array[i].row] = i;
		1234	}
		1235
		1236	skip = 0;
		1237	for (i = 0; i < len - div_pos; ++i) {
		1238	if (i > 0 &&
		1239	isl_seq_eq(qp->div->row[i - skip - 1],
		1240	qp->div->row[i - skip], qp->div->n_col)) {
		1241	qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
		1242	isl_mat_col_add(qp->div, 2 + div_pos + i - skip - 1,
		1243	2 + div_pos + i - skip);
		1244	qp->div = isl_mat_drop_cols(qp->div,
		1245	2 + div_pos + i - skip, 1);
		1246	skip++;
		1247	}
		1248	reordering[div_pos + array[i].row] = div_pos + i - skip;
		1249	}
		1250
		1251	qp->upoly = reorder(qp->upoly, reordering);
		1252
		1253	if (!qp->upoly || !qp->div)
		1254	goto error;
		1255
		1256	free(at);
		1257	free(pos);
		1258	free(array);
		1259	free(reordering);
		1260
		1261	return qp;
		1262	error:
		1263	free(at);
		1264	free(pos);
		1265	free(array);
		1266	free(reordering);
		1267	isl_qpolynomial_free(qp);
		1268	return NULL;
		1269	}
		1270
		1271	static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
		1272	int *exp, int first)
		1273	{
		1274	int i;
		1275	struct isl_upoly_rec *rec;
		1276
		1277	if (isl_upoly_is_cst(up))
		1278	return up;
		1279
		1280	if (up->var < first)
		1281	return up;
		1282
		1283	if (exp[up->var - first] == up->var - first)
		1284	return up;
		1285
		1286	up = isl_upoly_cow(up);
		1287	if (!up)
		1288	goto error;
		1289
		1290	up->var = exp[up->var - first] + first;
		1291
		1292	rec = isl_upoly_as_rec(up);
		1293	if (!rec)
		1294	goto error;
		1295
		1296	for (i = 0; i < rec->n; ++i) {
		1297	rec->p[i] = expand(rec->p[i], exp, first);
		1298	if (!rec->p[i])
		1299	goto error;
		1300	}
		1301
		1302	return up;
		1303	error:
		1304	isl_upoly_free(up);
		1305	return NULL;
		1306	}
		1307
		1308	static __isl_give isl_qpolynomial *with_merged_divs(
		1309	__isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
		1310	__isl_take isl_qpolynomial *qp2),
		1311	__isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
		1312	{
		1313	int *exp1 = NULL;
		1314	int *exp2 = NULL;
		1315	isl_mat *div = NULL;
		1316
		1317	qp1 = isl_qpolynomial_cow(qp1);
		1318	qp2 = isl_qpolynomial_cow(qp2);
		1319
		1320	if (!qp1 || !qp2)
		1321	goto error;
		1322
		1323	isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
		1324	qp1->div->n_col >= qp2->div->n_col, goto error);
		1325
		1326	exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
		1327	exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
		1328	if (!exp1 || !exp2)
		1329	goto error;
		1330
		1331	div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2);
		1332	if (!div)
		1333	goto error;
		1334
		1335	isl_mat_free(qp1->div);
		1336	qp1->div = isl_mat_copy(div);
		1337	isl_mat_free(qp2->div);
		1338	qp2->div = isl_mat_copy(div);
		1339
		1340	qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
		1341	qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
		1342
		1343	if (!qp1->upoly || !qp2->upoly)
		1344	goto error;
		1345
		1346	isl_mat_free(div);
		1347	free(exp1);
		1348	free(exp2);
		1349
		1350	return fn(qp1, qp2);
		1351	error:
		1352	isl_mat_free(div);
		1353	free(exp1);
		1354	free(exp2);
		1355	isl_qpolynomial_free(qp1);
		1356	isl_qpolynomial_free(qp2);
		1357	return NULL;
		1358	}
		1359
		1360	__isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
		1361	__isl_take isl_qpolynomial *qp2)
		1362	{
		1363	qp1 = isl_qpolynomial_cow(qp1);
		1364
		1365	if (!qp1 || !qp2)
		1366	goto error;
		1367
		1368	if (qp1->div->n_row < qp2->div->n_row)
		1369	return isl_qpolynomial_add(qp2, qp1);
		1370
		1371	isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
		1372	if (!compatible_divs(qp1->div, qp2->div))
		1373	return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
		1374
		1375	qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
		1376	if (!qp1->upoly)
		1377	goto error;
		1378
		1379	isl_qpolynomial_free(qp2);
		1380
		1381	return qp1;
		1382	error:
		1383	isl_qpolynomial_free(qp1);
		1384	isl_qpolynomial_free(qp2);
		1385	return NULL;
		1386	}
		1387
		1388	__isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
		1389	__isl_keep isl_set *dom,
		1390	__isl_take isl_qpolynomial *qp1,
		1391	__isl_take isl_qpolynomial *qp2)
		1392	{
		1393	qp1 = isl_qpolynomial_add(qp1, qp2);
		1394	qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom));
		1395	return qp1;
		1396	}
		1397
		1398	__isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
		1399	__isl_take isl_qpolynomial *qp2)
		1400	{
		1401	return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
		1402	}
		1403
		1404	__isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int(
		1405	__isl_take isl_qpolynomial *qp, isl_int v)
		1406	{
		1407	if (isl_int_is_zero(v))
		1408	return qp;
		1409
		1410	qp = isl_qpolynomial_cow(qp);
		1411	if (!qp)
		1412	return NULL;
		1413
		1414	qp->upoly = isl_upoly_add_isl_int(qp->upoly, v);
		1415	if (!qp->upoly)
		1416	goto error;
		1417
		1418	return qp;
		1419	error:
		1420	isl_qpolynomial_free(qp);
		1421	return NULL;
		1422
		1423	}
		1424
		1425	__isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
		1426	{
		1427	if (!qp)
		1428	return NULL;
		1429
		1430	return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone);
		1431	}
		1432
		1433	__isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int(
		1434	__isl_take isl_qpolynomial *qp, isl_int v)
		1435	{
		1436	if (isl_int_is_one(v))
		1437	return qp;
		1438
		1439	if (qp && isl_int_is_zero(v)) {
		1440	isl_qpolynomial *zero;
		1441	zero = isl_qpolynomial_zero_on_domain(isl_space_copy(qp->dim));
		1442	isl_qpolynomial_free(qp);
		1443	return zero;
		1444	}
		1445
		1446	qp = isl_qpolynomial_cow(qp);
		1447	if (!qp)
		1448	return NULL;
		1449
		1450	qp->upoly = isl_upoly_mul_isl_int(qp->upoly, v);
		1451	if (!qp->upoly)
		1452	goto error;
		1453
		1454	return qp;
		1455	error:
		1456	isl_qpolynomial_free(qp);
		1457	return NULL;
		1458	}
		1459
		1460	__isl_give isl_qpolynomial *isl_qpolynomial_scale(
		1461	__isl_take isl_qpolynomial *qp, isl_int v)
		1462	{
		1463	return isl_qpolynomial_mul_isl_int(qp, v);
		1464	}
		1465
		1466	__isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
		1467	__isl_take isl_qpolynomial *qp2)
		1468	{
		1469	qp1 = isl_qpolynomial_cow(qp1);
		1470
		1471	if (!qp1 || !qp2)
		1472	goto error;
		1473
		1474	if (qp1->div->n_row < qp2->div->n_row)
		1475	return isl_qpolynomial_mul(qp2, qp1);
		1476
		1477	isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
		1478	if (!compatible_divs(qp1->div, qp2->div))
		1479	return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
		1480
		1481	qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
		1482	if (!qp1->upoly)
		1483	goto error;
		1484
		1485	isl_qpolynomial_free(qp2);
		1486
		1487	return qp1;
		1488	error:
		1489	isl_qpolynomial_free(qp1);
		1490	isl_qpolynomial_free(qp2);
		1491	return NULL;
		1492	}
		1493
		1494	__isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp,
		1495	unsigned power)
		1496	{
		1497	qp = isl_qpolynomial_cow(qp);
		1498
		1499	if (!qp)
		1500	return NULL;
		1501
		1502	qp->upoly = isl_upoly_pow(qp->upoly, power);
		1503	if (!qp->upoly)
		1504	goto error;
		1505
		1506	return qp;
		1507	error:
		1508	isl_qpolynomial_free(qp);
		1509	return NULL;
		1510	}
		1511
		1512	__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_pow(
		1513	__isl_take isl_pw_qpolynomial *pwqp, unsigned power)
		1514	{
		1515	int i;
		1516
		1517	if (power == 1)
		1518	return pwqp;
		1519
		1520	pwqp = isl_pw_qpolynomial_cow(pwqp);
		1521	if (!pwqp)
		1522	return NULL;
		1523
		1524	for (i = 0; i < pwqp->n; ++i) {
		1525	pwqp->p[i].qp = isl_qpolynomial_pow(pwqp->p[i].qp, power);
		1526	if (!pwqp->p[i].qp)
		1527	return isl_pw_qpolynomial_free(pwqp);
		1528	}
		1529
		1530	return pwqp;
		1531	}
		1532
		1533	__isl_give isl_qpolynomial *isl_qpolynomial_zero_on_domain(
		1534	__isl_take isl_space *dim)
		1535	{
		1536	if (!dim)
		1537	return NULL;
		1538	return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
		1539	}
		1540
		1541	__isl_give isl_qpolynomial *isl_qpolynomial_one_on_domain(
		1542	__isl_take isl_space *dim)
		1543	{
		1544	if (!dim)
		1545	return NULL;
		1546	return isl_qpolynomial_alloc(dim, 0, isl_upoly_one(dim->ctx));
		1547	}
		1548
		1549	__isl_give isl_qpolynomial *isl_qpolynomial_infty_on_domain(
		1550	__isl_take isl_space *dim)
		1551	{
		1552	if (!dim)
		1553	return NULL;
		1554	return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
		1555	}
		1556
		1557	__isl_give isl_qpolynomial *isl_qpolynomial_neginfty_on_domain(
		1558	__isl_take isl_space *dim)
		1559	{
		1560	if (!dim)
		1561	return NULL;
		1562	return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
		1563	}
		1564
		1565	__isl_give isl_qpolynomial *isl_qpolynomial_nan_on_domain(
		1566	__isl_take isl_space *dim)
		1567	{
		1568	if (!dim)
		1569	return NULL;
		1570	return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
		1571	}
		1572
		1573	__isl_give isl_qpolynomial *isl_qpolynomial_cst_on_domain(
		1574	__isl_take isl_space *dim,
		1575	isl_int v)
		1576	{
		1577	struct isl_qpolynomial *qp;
		1578	struct isl_upoly_cst *cst;
		1579
		1580	if (!dim)
		1581	return NULL;
		1582
		1583	qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
		1584	if (!qp)
		1585	return NULL;
		1586
		1587	cst = isl_upoly_as_cst(qp->upoly);
		1588	isl_int_set(cst->n, v);
		1589
		1590	return qp;
		1591	}
		1592
		1593	int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
		1594	isl_int *n, isl_int *d)
		1595	{
		1596	struct isl_upoly_cst *cst;
		1597
		1598	if (!qp)
		1599	return -1;
		1600
		1601	if (!isl_upoly_is_cst(qp->upoly))
		1602	return 0;
		1603
		1604	cst = isl_upoly_as_cst(qp->upoly);
		1605	if (!cst)
		1606	return -1;
		1607
		1608	if (n)
		1609	isl_int_set(*n, cst->n);
		1610	if (d)
		1611	isl_int_set(*d, cst->d);
		1612
		1613	return 1;
		1614	}
		1615
		1616	int isl_upoly_is_affine(__isl_keep struct isl_upoly *up)
		1617	{
		1618	int is_cst;
		1619	struct isl_upoly_rec *rec;
		1620
		1621	if (!up)
		1622	return -1;
		1623
		1624	if (up->var < 0)
		1625	return 1;
		1626
		1627	rec = isl_upoly_as_rec(up);
		1628	if (!rec)
		1629	return -1;
		1630
		1631	if (rec->n > 2)
		1632	return 0;
		1633
		1634	isl_assert(up->ctx, rec->n > 1, return -1);
		1635
		1636	is_cst = isl_upoly_is_cst(rec->p[1]);
		1637	if (is_cst < 0)
		1638	return -1;
		1639	if (!is_cst)
		1640	return 0;
		1641
		1642	return isl_upoly_is_affine(rec->p[0]);
		1643	}
		1644
		1645	int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
		1646	{
		1647	if (!qp)
		1648	return -1;
		1649
		1650	if (qp->div->n_row > 0)
		1651	return 0;
		1652
		1653	return isl_upoly_is_affine(qp->upoly);
		1654	}
		1655
		1656	static void update_coeff(__isl_keep isl_vec *aff,
		1657	__isl_keep struct isl_upoly_cst *cst, int pos)
		1658	{
		1659	isl_int gcd;
		1660	isl_int f;
		1661
		1662	if (isl_int_is_zero(cst->n))
		1663	return;
		1664
		1665	isl_int_init(gcd);
		1666	isl_int_init(f);
		1667	isl_int_gcd(gcd, cst->d, aff->el[0]);
		1668	isl_int_divexact(f, cst->d, gcd);
		1669	isl_int_divexact(gcd, aff->el[0], gcd);
		1670	isl_seq_scale(aff->el, aff->el, f, aff->size);
		1671	isl_int_mul(aff->el[1 + pos], gcd, cst->n);
		1672	isl_int_clear(gcd);
		1673	isl_int_clear(f);
		1674	}
		1675
		1676	int isl_upoly_update_affine(__isl_keep struct isl_upoly *up,
		1677	__isl_keep isl_vec *aff)
		1678	{
		1679	struct isl_upoly_cst *cst;
		1680	struct isl_upoly_rec *rec;
		1681
		1682	if (!up || !aff)
		1683	return -1;
		1684
		1685	if (up->var < 0) {
		1686	struct isl_upoly_cst *cst;
		1687
		1688	cst = isl_upoly_as_cst(up);
		1689	if (!cst)
		1690	return -1;
		1691	update_coeff(aff, cst, 0);
		1692	return 0;
		1693	}
		1694
		1695	rec = isl_upoly_as_rec(up);
		1696	if (!rec)
		1697	return -1;
		1698	isl_assert(up->ctx, rec->n == 2, return -1);
		1699
		1700	cst = isl_upoly_as_cst(rec->p[1]);
		1701	if (!cst)
		1702	return -1;
		1703	update_coeff(aff, cst, 1 + up->var);
		1704
		1705	return isl_upoly_update_affine(rec->p[0], aff);
		1706	}
		1707
		1708	__isl_give isl_vec *isl_qpolynomial_extract_affine(
		1709	__isl_keep isl_qpolynomial *qp)
		1710	{
		1711	isl_vec *aff;
		1712	unsigned d;
		1713
		1714	if (!qp)
		1715	return NULL;
		1716
		1717	d = isl_space_dim(qp->dim, isl_dim_all);
		1718	aff = isl_vec_alloc(qp->div->ctx, 2 + d + qp->div->n_row);
		1719	if (!aff)
		1720	return NULL;
		1721
		1722	isl_seq_clr(aff->el + 1, 1 + d + qp->div->n_row);
		1723	isl_int_set_si(aff->el[0], 1);
		1724
		1725	if (isl_upoly_update_affine(qp->upoly, aff) < 0)
		1726	goto error;
		1727
		1728	return aff;
		1729	error:
		1730	isl_vec_free(aff);
		1731	return NULL;
		1732	}
		1733
		1734	int isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1,
		1735	__isl_keep isl_qpolynomial *qp2)
		1736	{
		1737	int equal;
		1738
		1739	if (!qp1 || !qp2)
		1740	return -1;
		1741
		1742	equal = isl_space_is_equal(qp1->dim, qp2->dim);
		1743	if (equal < 0 || !equal)
		1744	return equal;
		1745
		1746	equal = isl_mat_is_equal(qp1->div, qp2->div);
		1747	if (equal < 0 || !equal)
		1748	return equal;
		1749
		1750	return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
		1751	}
		1752
		1753	static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
		1754	{
		1755	int i;
		1756	struct isl_upoly_rec *rec;
		1757
		1758	if (isl_upoly_is_cst(up)) {
		1759	struct isl_upoly_cst *cst;
		1760	cst = isl_upoly_as_cst(up);
		1761	if (!cst)
		1762	return;
		1763	isl_int_lcm(*d, *d, cst->d);
		1764	return;
		1765	}
		1766
		1767	rec = isl_upoly_as_rec(up);
		1768	if (!rec)
		1769	return;
		1770
		1771	for (i = 0; i < rec->n; ++i)
		1772	upoly_update_den(rec->p[i], d);
		1773	}
		1774
		1775	void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
		1776	{
		1777	isl_int_set_si(*d, 1);
		1778	if (!qp)
		1779	return;
		1780	upoly_update_den(qp->upoly, d);
		1781	}
		1782
		1783	__isl_give isl_qpolynomial *isl_qpolynomial_var_pow_on_domain(
		1784	__isl_take isl_space *dim, int pos, int power)
		1785	{
		1786	struct isl_ctx *ctx;
		1787
		1788	if (!dim)
		1789	return NULL;
		1790
		1791	ctx = dim->ctx;
		1792
		1793	return isl_qpolynomial_alloc(dim, 0, isl_upoly_var_pow(ctx, pos, power));
		1794	}
		1795
		1796	__isl_give isl_qpolynomial *isl_qpolynomial_var_on_domain(__isl_take isl_space *dim,
		1797	enum isl_dim_type type, unsigned pos)
		1798	{
		1799	if (!dim)
		1800	return NULL;
		1801
		1802	isl_assert(dim->ctx, isl_space_dim(dim, isl_dim_in) == 0, goto error);
		1803	isl_assert(dim->ctx, pos < isl_space_dim(dim, type), goto error);
		1804
		1805	if (type == isl_dim_set)
		1806	pos += isl_space_dim(dim, isl_dim_param);
		1807
		1808	return isl_qpolynomial_var_pow_on_domain(dim, pos, 1);
		1809	error:
		1810	isl_space_free(dim);
		1811	return NULL;
		1812	}
		1813
		1814	__isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
		1815	unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
		1816	{
		1817	int i;
		1818	struct isl_upoly_rec *rec;
		1819	struct isl_upoly *base, *res;
		1820
		1821	if (!up)
		1822	return NULL;
		1823
		1824	if (isl_upoly_is_cst(up))
		1825	return up;
		1826
		1827	if (up->var < first)
		1828	return up;
		1829
		1830	rec = isl_upoly_as_rec(up);
		1831	if (!rec)
		1832	goto error;
		1833
		1834	isl_assert(up->ctx, rec->n >= 1, goto error);
		1835
		1836	if (up->var >= first + n)
		1837	base = isl_upoly_var_pow(up->ctx, up->var, 1);
		1838	else
		1839	base = isl_upoly_copy(subs[up->var - first]);
		1840
		1841	res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
		1842	for (i = rec->n - 2; i >= 0; --i) {
		1843	struct isl_upoly *t;
		1844	t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
		1845	res = isl_upoly_mul(res, isl_upoly_copy(base));
		1846	res = isl_upoly_sum(res, t);
		1847	}
		1848
		1849	isl_upoly_free(base);
		1850	isl_upoly_free(up);
		1851
		1852	return res;
		1853	error:
		1854	isl_upoly_free(up);
		1855	return NULL;
		1856	}
		1857
		1858	__isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
		1859	isl_int denom, unsigned len)
		1860	{
		1861	int i;
		1862	struct isl_upoly *up;
		1863
		1864	isl_assert(ctx, len >= 1, return NULL);
		1865
		1866	up = isl_upoly_rat_cst(ctx, f[0], denom);
		1867	for (i = 0; i < len - 1; ++i) {
		1868	struct isl_upoly *t;
		1869	struct isl_upoly *c;
		1870
		1871	if (isl_int_is_zero(f[1 + i]))
		1872	continue;
		1873
		1874	c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
		1875	t = isl_upoly_var_pow(ctx, i, 1);
		1876	t = isl_upoly_mul(c, t);
		1877	up = isl_upoly_sum(up, t);
		1878	}
		1879
		1880	return up;
		1881	}
		1882
		1883	/* Remove common factor of non-constant terms and denominator.
		1884	*/
		1885	static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
		1886	{
		1887	isl_ctx *ctx = qp->div->ctx;
		1888	unsigned total = qp->div->n_col - 2;
		1889
		1890	isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd);
		1891	isl_int_gcd(ctx->normalize_gcd,
		1892	ctx->normalize_gcd, qp->div->row[div][0]);
		1893	if (isl_int_is_one(ctx->normalize_gcd))
		1894	return;
		1895
		1896	isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2,
		1897	ctx->normalize_gcd, total);
		1898	isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0],
		1899	ctx->normalize_gcd);
		1900	isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1],
		1901	ctx->normalize_gcd);
		1902	}
		1903
		1904	/* Replace the integer division identified by "div" by the polynomial "s".
		1905	* The integer division is assumed not to appear in the definition
		1906	* of any other integer divisions.
		1907	*/
		1908	static __isl_give isl_qpolynomial *substitute_div(
		1909	__isl_take isl_qpolynomial *qp,
		1910	int div, __isl_take struct isl_upoly *s)
		1911	{
		1912	int i;
		1913	int total;
		1914	int *reordering;
		1915
		1916	if (!qp || !s)
		1917	goto error;
		1918
		1919	qp = isl_qpolynomial_cow(qp);
		1920	if (!qp)
		1921	goto error;
		1922
		1923	total = isl_space_dim(qp->dim, isl_dim_all);
		1924	qp->upoly = isl_upoly_subs(qp->upoly, total + div, 1, &s);
		1925	if (!qp->upoly)
		1926	goto error;
		1927
		1928	reordering = isl_alloc_array(qp->dim->ctx, int, total + qp->div->n_row);
		1929	if (!reordering)
		1930	goto error;
		1931	for (i = 0; i < total + div; ++i)
		1932	reordering[i] = i;
		1933	for (i = total + div + 1; i < total + qp->div->n_row; ++i)
		1934	reordering[i] = i - 1;
		1935	qp->div = isl_mat_drop_rows(qp->div, div, 1);
		1936	qp->div = isl_mat_drop_cols(qp->div, 2 + total + div, 1);
		1937	qp->upoly = reorder(qp->upoly, reordering);
		1938	free(reordering);
		1939
		1940	if (!qp->upoly || !qp->div)
		1941	goto error;
		1942
		1943	isl_upoly_free(s);
		1944	return qp;
		1945	error:
		1946	isl_qpolynomial_free(qp);
		1947	isl_upoly_free(s);
		1948	return NULL;
		1949	}
		1950
		1951	/* Replace all integer divisions [e/d] that turn out to not actually be integer
		1952	* divisions because d is equal to 1 by their definition, i.e., e.
		1953	*/
		1954	static __isl_give isl_qpolynomial *substitute_non_divs(
		1955	__isl_take isl_qpolynomial *qp)
		1956	{
		1957	int i, j;
		1958	int total;
		1959	struct isl_upoly *s;
		1960
		1961	if (!qp)
		1962	return NULL;
		1963
		1964	total = isl_space_dim(qp->dim, isl_dim_all);
		1965	for (i = 0; qp && i < qp->div->n_row; ++i) {
		1966	if (!isl_int_is_one(qp->div->row[i][0]))
		1967	continue;
		1968	for (j = i + 1; j < qp->div->n_row; ++j) {
		1969	if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
		1970	continue;
		1971	isl_seq_combine(qp->div->row[j] + 1,
		1972	qp->div->ctx->one, qp->div->row[j] + 1,
		1973	qp->div->row[j][2 + total + i],
		1974	qp->div->row[i] + 1, 1 + total + i);
		1975	isl_int_set_si(qp->div->row[j][2 + total + i], 0);
		1976	normalize_div(qp, j);
		1977	}
		1978	s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
		1979	qp->div->row[i][0], qp->div->n_col - 1);
		1980	qp = substitute_div(qp, i, s);
		1981	--i;
		1982	}
		1983
		1984	return qp;
		1985	}
		1986
		1987	/* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
		1988	* with d the denominator. When replacing the coefficient e of x by
		1989	* d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
		1990	* inside the division, so we need to add floor(e/d) * x outside.
		1991	* That is, we replace q by q' + floor(e/d) * x and we therefore need
		1992	* to adjust the coefficient of x in each later div that depends on the
		1993	* current div "div" and also in the affine expression "aff"
		1994	* (if it too depends on "div").
		1995	*/
		1996	static void reduce_div(__isl_keep isl_qpolynomial *qp, int div,
		1997	__isl_keep isl_vec *aff)
		1998	{
		1999	int i, j;
		2000	isl_int v;
		2001	unsigned total = qp->div->n_col - qp->div->n_row - 2;
		2002
		2003	isl_int_init(v);
		2004	for (i = 0; i < 1 + total + div; ++i) {
		2005	if (isl_int_is_nonneg(qp->div->row[div][1 + i]) &&
		2006	isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0]))
		2007	continue;
		2008	isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]);
		2009	isl_int_fdiv_r(qp->div->row[div][1 + i],
		2010	qp->div->row[div][1 + i], qp->div->row[div][0]);
		2011	if (!isl_int_is_zero(aff->el[1 + total + div]))
		2012	isl_int_addmul(aff->el[i], v, aff->el[1 + total + div]);
		2013	for (j = div + 1; j < qp->div->n_row; ++j) {
		2014	if (isl_int_is_zero(qp->div->row[j][2 + total + div]))
		2015	continue;
		2016	isl_int_addmul(qp->div->row[j][1 + i],
		2017	v, qp->div->row[j][2 + total + div]);
		2018	}
		2019	}
		2020	isl_int_clear(v);
		2021	}
		2022
		2023	/* Check if the last non-zero coefficient is bigger that half of the
		2024	* denominator. If so, we will invert the div to further reduce the number
		2025	* of distinct divs that may appear.
		2026	* If the last non-zero coefficient is exactly half the denominator,
		2027	* then we continue looking for earlier coefficients that are bigger
		2028	* than half the denominator.
		2029	*/
		2030	static int needs_invert(__isl_keep isl_mat *div, int row)
		2031	{
		2032	int i;
		2033	int cmp;
		2034
		2035	for (i = div->n_col - 1; i >= 1; --i) {
		2036	if (isl_int_is_zero(div->row[row][i]))
		2037	continue;
		2038	isl_int_mul_ui(div->row[row][i], div->row[row][i], 2);
		2039	cmp = isl_int_cmp(div->row[row][i], div->row[row][0]);
		2040	isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2);
		2041	if (cmp)
		2042	return cmp > 0;
		2043	if (i == 1)
		2044	return 1;
		2045	}
		2046
		2047	return 0;
		2048	}
		2049
		2050	/* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
		2051	* We only invert the coefficients of e (and the coefficient of q in
		2052	* later divs and in "aff"). After calling this function, the
		2053	* coefficients of e should be reduced again.
		2054	*/
		2055	static void invert_div(__isl_keep isl_qpolynomial *qp, int div,
		2056	__isl_keep isl_vec *aff)
		2057	{
		2058	unsigned total = qp->div->n_col - qp->div->n_row - 2;
		2059
		2060	isl_seq_neg(qp->div->row[div] + 1,
		2061	qp->div->row[div] + 1, qp->div->n_col - 1);
		2062	isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1);
		2063	isl_int_add(qp->div->row[div][1],
		2064	qp->div->row[div][1], qp->div->row[div][0]);
		2065	if (!isl_int_is_zero(aff->el[1 + total + div]))
		2066	isl_int_neg(aff->el[1 + total + div], aff->el[1 + total + div]);
		2067	isl_mat_col_mul(qp->div, 2 + total + div,
		2068	qp->div->ctx->negone, 2 + total + div);
		2069	}
		2070
		2071	/* Assuming "qp" is a monomial, reduce all its divs to have coefficients
		2072	* in the interval [0, d-1], with d the denominator and such that the
		2073	* last non-zero coefficient that is not equal to d/2 is smaller than d/2.
		2074	*
		2075	* After the reduction, some divs may have become redundant or identical,
		2076	* so we call substitute_non_divs and sort_divs. If these functions
		2077	* eliminate divs or merge two or more divs into one, the coefficients
		2078	* of the enclosing divs may have to be reduced again, so we call
		2079	* ourselves recursively if the number of divs decreases.
		2080	*/
		2081	static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp)
		2082	{
		2083	int i;
		2084	isl_vec *aff = NULL;
		2085	struct isl_upoly *s;
		2086	unsigned n_div;
		2087
		2088	if (!qp)
		2089	return NULL;
		2090
		2091	aff = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
		2092	aff = isl_vec_clr(aff);
		2093	if (!aff)
		2094	goto error;
		2095
		2096	isl_int_set_si(aff->el[1 + qp->upoly->var], 1);
		2097
		2098	for (i = 0; i < qp->div->n_row; ++i) {
		2099	normalize_div(qp, i);
		2100	reduce_div(qp, i, aff);
		2101	if (needs_invert(qp->div, i)) {
		2102	invert_div(qp, i, aff);
		2103	reduce_div(qp, i, aff);
		2104	}
		2105	}
		2106
		2107	s = isl_upoly_from_affine(qp->div->ctx, aff->el,
		2108	qp->div->ctx->one, aff->size);
		2109	qp->upoly = isl_upoly_subs(qp->upoly, qp->upoly->var, 1, &s);
		2110	isl_upoly_free(s);
		2111	if (!qp->upoly)
		2112	goto error;
		2113
		2114	isl_vec_free(aff);
		2115
		2116	n_div = qp->div->n_row;
		2117	qp = substitute_non_divs(qp);
		2118	qp = sort_divs(qp);
		2119	if (qp && qp->div->n_row < n_div)
		2120	return reduce_divs(qp);
		2121
		2122	return qp;
		2123	error:
		2124	isl_qpolynomial_free(qp);
		2125	isl_vec_free(aff);
		2126	return NULL;
		2127	}
		2128
		2129	__isl_give isl_qpolynomial *isl_qpolynomial_rat_cst_on_domain(
		2130	__isl_take isl_space *dim, const isl_int n, const isl_int d)
		2131	{
		2132	struct isl_qpolynomial *qp;
		2133	struct isl_upoly_cst *cst;
		2134
		2135	if (!dim)
		2136	return NULL;
		2137
		2138	qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
		2139	if (!qp)
		2140	return NULL;
		2141
		2142	cst = isl_upoly_as_cst(qp->upoly);
		2143	isl_int_set(cst->n, n);
		2144	isl_int_set(cst->d, d);
		2145
		2146	return qp;
		2147	}
		2148
		2149	static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
		2150	{
		2151	struct isl_upoly_rec *rec;
		2152	int i;
		2153
		2154	if (!up)
		2155	return -1;
		2156
		2157	if (isl_upoly_is_cst(up))
		2158	return 0;
		2159
		2160	if (up->var < d)
		2161	active[up->var] = 1;
		2162
		2163	rec = isl_upoly_as_rec(up);
		2164	for (i = 0; i < rec->n; ++i)
		2165	if (up_set_active(rec->p[i], active, d) < 0)
		2166	return -1;
		2167
		2168	return 0;
		2169	}
		2170
		2171	static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
		2172	{
		2173	int i, j;
		2174	int d = isl_space_dim(qp->dim, isl_dim_all);
		2175
		2176	if (!qp || !active)
		2177	return -1;
		2178
		2179	for (i = 0; i < d; ++i)
		2180	for (j = 0; j < qp->div->n_row; ++j) {
		2181	if (isl_int_is_zero(qp->div->row[j][2 + i]))
		2182	continue;
		2183	active[i] = 1;
		2184	break;
		2185	}
		2186
		2187	return up_set_active(qp->upoly, active, d);
		2188	}
		2189
		2190	int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
		2191	enum isl_dim_type type, unsigned first, unsigned n)
		2192	{
		2193	int i;
		2194	int *active = NULL;
		2195	int involves = 0;
		2196
		2197	if (!qp)
		2198	return -1;
		2199	if (n == 0)
		2200	return 0;
		2201
		2202	isl_assert(qp->dim->ctx,
		2203	first + n <= isl_qpolynomial_dim(qp, type), return -1);
		2204	isl_assert(qp->dim->ctx, type == isl_dim_param ||
		2205	type == isl_dim_in, return -1);
		2206
		2207	active = isl_calloc_array(qp->dim->ctx, int,
		2208	isl_space_dim(qp->dim, isl_dim_all));
		2209	if (set_active(qp, active) < 0)
		2210	goto error;
		2211
		2212	if (type == isl_dim_in)
		2213	first += isl_space_dim(qp->dim, isl_dim_param);
		2214	for (i = 0; i < n; ++i)
		2215	if (active[first + i]) {
		2216	involves = 1;
		2217	break;
		2218	}
		2219
		2220	free(active);
		2221
		2222	return involves;
		2223	error:
		2224	free(active);
		2225	return -1;
		2226	}
		2227
		2228	/* Remove divs that do not appear in the quasi-polynomial, nor in any
		2229	* of the divs that do appear in the quasi-polynomial.
		2230	*/
		2231	static __isl_give isl_qpolynomial *remove_redundant_divs(
		2232	__isl_take isl_qpolynomial *qp)
		2233	{
		2234	int i, j;
		2235	int d;
		2236	int len;
		2237	int skip;
		2238	int *active = NULL;
		2239	int *reordering = NULL;
		2240	int redundant = 0;
		2241	int n_div;
		2242	isl_ctx *ctx;
		2243
		2244	if (!qp)
		2245	return NULL;
		2246	if (qp->div->n_row == 0)
		2247	return qp;
		2248
		2249	d = isl_space_dim(qp->dim, isl_dim_all);
		2250	len = qp->div->n_col - 2;
		2251	ctx = isl_qpolynomial_get_ctx(qp);
		2252	active = isl_calloc_array(ctx, int, len);
		2253	if (!active)
		2254	goto error;
		2255
		2256	if (up_set_active(qp->upoly, active, len) < 0)
		2257	goto error;
		2258
		2259	for (i = qp->div->n_row - 1; i >= 0; --i) {
		2260	if (!active[d + i]) {
		2261	redundant = 1;
		2262	continue;
		2263	}
		2264	for (j = 0; j < i; ++j) {
		2265	if (isl_int_is_zero(qp->div->row[i][2 + d + j]))
		2266	continue;
		2267	active[d + j] = 1;
		2268	break;
		2269	}
		2270	}
		2271
		2272	if (!redundant) {
		2273	free(active);
		2274	return qp;
		2275	}
		2276
		2277	reordering = isl_alloc_array(qp->div->ctx, int, len);
		2278	if (!reordering)
		2279	goto error;
		2280
		2281	for (i = 0; i < d; ++i)
		2282	reordering[i] = i;
		2283
		2284	skip = 0;
		2285	n_div = qp->div->n_row;
		2286	for (i = 0; i < n_div; ++i) {
		2287	if (!active[d + i]) {
		2288	qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
		2289	qp->div = isl_mat_drop_cols(qp->div,
		2290	2 + d + i - skip, 1);
		2291	skip++;
		2292	}
		2293	reordering[d + i] = d + i - skip;
		2294	}
		2295
		2296	qp->upoly = reorder(qp->upoly, reordering);
		2297
		2298	if (!qp->upoly || !qp->div)
		2299	goto error;
		2300
		2301	free(active);
		2302	free(reordering);
		2303
		2304	return qp;
		2305	error:
		2306	free(active);
		2307	free(reordering);
		2308	isl_qpolynomial_free(qp);
		2309	return NULL;
		2310	}
		2311
		2312	__isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
		2313	unsigned first, unsigned n)
		2314	{
		2315	int i;
		2316	struct isl_upoly_rec *rec;
		2317
		2318	if (!up)
		2319	return NULL;
		2320	if (n == 0 || up->var < 0 || up->var < first)
		2321	return up;
		2322	if (up->var < first + n) {
		2323	up = replace_by_constant_term(up);
		2324	return isl_upoly_drop(up, first, n);
		2325	}
		2326	up = isl_upoly_cow(up);
		2327	if (!up)
		2328	return NULL;
		2329	up->var -= n;
		2330	rec = isl_upoly_as_rec(up);
		2331	if (!rec)
		2332	goto error;
		2333
		2334	for (i = 0; i < rec->n; ++i) {
		2335	rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
		2336	if (!rec->p[i])
		2337	goto error;
		2338	}
		2339
		2340	return up;
		2341	error:
		2342	isl_upoly_free(up);
		2343	return NULL;
		2344	}
		2345
		2346	__isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
		2347	__isl_take isl_qpolynomial *qp,
		2348	enum isl_dim_type type, unsigned pos, const char *s)
		2349	{
		2350	qp = isl_qpolynomial_cow(qp);
		2351	if (!qp)
		2352	return NULL;
		2353	qp->dim = isl_space_set_dim_name(qp->dim, type, pos, s);
		2354	if (!qp->dim)
		2355	goto error;
		2356	return qp;
		2357	error:
		2358	isl_qpolynomial_free(qp);
		2359	return NULL;
		2360	}
		2361
		2362	__isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
		2363	__isl_take isl_qpolynomial *qp,
		2364	enum isl_dim_type type, unsigned first, unsigned n)
		2365	{
		2366	if (!qp)
		2367	return NULL;
		2368	if (type == isl_dim_out)
		2369	isl_die(qp->dim->ctx, isl_error_invalid,
		2370	"cannot drop output/set dimension",
		2371	goto error);
		2372	if (type == isl_dim_in)
		2373	type = isl_dim_set;
		2374	if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
		2375	return qp;
		2376
		2377	qp = isl_qpolynomial_cow(qp);
		2378	if (!qp)
		2379	return NULL;
		2380
		2381	isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
		2382	goto error);
		2383	isl_assert(qp->dim->ctx, type == isl_dim_param ||
		2384	type == isl_dim_set, goto error);
		2385
		2386	qp->dim = isl_space_drop_dims(qp->dim, type, first, n);
		2387	if (!qp->dim)
		2388	goto error;
		2389
		2390	if (type == isl_dim_set)
		2391	first += isl_space_dim(qp->dim, isl_dim_param);
		2392
		2393	qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
		2394	if (!qp->div)
		2395	goto error;
		2396
		2397	qp->upoly = isl_upoly_drop(qp->upoly, first, n);
		2398	if (!qp->upoly)
		2399	goto error;
		2400
		2401	return qp;
		2402	error:
		2403	isl_qpolynomial_free(qp);
		2404	return NULL;
		2405	}
		2406
		2407	/* Project the domain of the quasi-polynomial onto its parameter space.
		2408	* The quasi-polynomial may not involve any of the domain dimensions.
		2409	*/
		2410	__isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params(
		2411	__isl_take isl_qpolynomial *qp)
		2412	{
		2413	isl_space *space;
		2414	unsigned n;
		2415	int involves;
		2416
		2417	n = isl_qpolynomial_dim(qp, isl_dim_in);
		2418	involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n);
		2419	if (involves < 0)
		2420	return isl_qpolynomial_free(qp);
		2421	if (involves)
		2422	isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
		2423	"polynomial involves some of the domain dimensions",
		2424	return isl_qpolynomial_free(qp));
		2425	qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n);
		2426	space = isl_qpolynomial_get_domain_space(qp);
		2427	space = isl_space_params(space);
		2428	qp = isl_qpolynomial_reset_domain_space(qp, space);
		2429	return qp;
		2430	}
		2431
		2432	static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
		2433	__isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
		2434	{
		2435	int i, j, k;
		2436	isl_int denom;
		2437	unsigned total;
		2438	unsigned n_div;
		2439	struct isl_upoly *up;
		2440
		2441	if (!eq)
		2442	goto error;
		2443	if (eq->n_eq == 0) {
		2444	isl_basic_set_free(eq);
		2445	return qp;
		2446	}
		2447
		2448	qp = isl_qpolynomial_cow(qp);
		2449	if (!qp)
		2450	goto error;
		2451	qp->div = isl_mat_cow(qp->div);
		2452	if (!qp->div)
		2453	goto error;
		2454
		2455	total = 1 + isl_space_dim(eq->dim, isl_dim_all);
		2456	n_div = eq->n_div;
		2457	isl_int_init(denom);
		2458	for (i = 0; i < eq->n_eq; ++i) {
		2459	j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
		2460	if (j < 0 || j == 0 || j >= total)
		2461	continue;
		2462
		2463	for (k = 0; k < qp->div->n_row; ++k) {
		2464	if (isl_int_is_zero(qp->div->row[k][1 + j]))
		2465	continue;
		2466	isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
		2467	&qp->div->row[k][0]);
		2468	normalize_div(qp, k);
		2469	}
		2470
		2471	if (isl_int_is_pos(eq->eq[i][j]))
		2472	isl_seq_neg(eq->eq[i], eq->eq[i], total);
		2473	isl_int_abs(denom, eq->eq[i][j]);
		2474	isl_int_set_si(eq->eq[i][j], 0);
		2475
		2476	up = isl_upoly_from_affine(qp->dim->ctx,
		2477	eq->eq[i], denom, total);
		2478	qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
		2479	isl_upoly_free(up);
		2480	}
		2481	isl_int_clear(denom);
		2482
		2483	if (!qp->upoly)
		2484	goto error;
		2485
		2486	isl_basic_set_free(eq);
		2487
		2488	qp = substitute_non_divs(qp);
		2489	qp = sort_divs(qp);
		2490
		2491	return qp;
		2492	error:
		2493	isl_basic_set_free(eq);
		2494	isl_qpolynomial_free(qp);
		2495	return NULL;
		2496	}
		2497
		2498	/* Exploit the equalities in "eq" to simplify the quasi-polynomial.
		2499	*/
		2500	__isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
		2501	__isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
		2502	{
		2503	if (!qp || !eq)
		2504	goto error;
		2505	if (qp->div->n_row > 0)
		2506	eq = isl_basic_set_add(eq, isl_dim_set, qp->div->n_row);
		2507	return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
		2508	error:
		2509	isl_basic_set_free(eq);
		2510	isl_qpolynomial_free(qp);
		2511	return NULL;
		2512	}
		2513
		2514	static __isl_give isl_basic_set *add_div_constraints(
		2515	__isl_take isl_basic_set *bset, __isl_take isl_mat *div)
		2516	{
		2517	int i;
		2518	unsigned total;
		2519
		2520	if (!bset || !div)
		2521	goto error;
		2522
		2523	bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
		2524	if (!bset)
		2525	goto error;
		2526	total = isl_basic_set_total_dim(bset);
		2527	for (i = 0; i < div->n_row; ++i)
		2528	if (isl_basic_set_add_div_constraints_var(bset,
		2529	total - div->n_row + i, div->row[i]) < 0)
		2530	goto error;
		2531
		2532	isl_mat_free(div);
		2533	return bset;
		2534	error:
		2535	isl_mat_free(div);
		2536	isl_basic_set_free(bset);
		2537	return NULL;
		2538	}
		2539
		2540	/* Look for equalities among the variables shared by context and qp
		2541	* and the integer divisions of qp, if any.
		2542	* The equalities are then used to eliminate variables and/or integer
		2543	* divisions from qp.
		2544	*/
		2545	__isl_give isl_qpolynomial *isl_qpolynomial_gist(
		2546	__isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
		2547	{
		2548	isl_basic_set *aff;
		2549
		2550	if (!qp)
		2551	goto error;
		2552	if (qp->div->n_row > 0) {
		2553	isl_basic_set *bset;
		2554	context = isl_set_add_dims(context, isl_dim_set,
		2555	qp->div->n_row);
		2556	bset = isl_basic_set_universe(isl_set_get_space(context));
		2557	bset = add_div_constraints(bset, isl_mat_copy(qp->div));
		2558	context = isl_set_intersect(context,
		2559	isl_set_from_basic_set(bset));
		2560	}
		2561
		2562	aff = isl_set_affine_hull(context);
		2563	return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
		2564	error:
		2565	isl_qpolynomial_free(qp);
		2566	isl_set_free(context);
		2567	return NULL;
		2568	}
		2569
		2570	__isl_give isl_qpolynomial *isl_qpolynomial_gist_params(
		2571	__isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
		2572	{
		2573	isl_space *space = isl_qpolynomial_get_domain_space(qp);
		2574	isl_set *dom_context = isl_set_universe(space);
		2575	dom_context = isl_set_intersect_params(dom_context, context);
		2576	return isl_qpolynomial_gist(qp, dom_context);
		2577	}
		2578
		2579	__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_qpolynomial(
		2580	__isl_take isl_qpolynomial *qp)
		2581	{
		2582	isl_set *dom;
		2583
		2584	if (!qp)
		2585	return NULL;
		2586	if (isl_qpolynomial_is_zero(qp)) {
		2587	isl_space *dim = isl_qpolynomial_get_space(qp);
		2588	isl_qpolynomial_free(qp);
		2589	return isl_pw_qpolynomial_zero(dim);
		2590	}
		2591
		2592	dom = isl_set_universe(isl_qpolynomial_get_domain_space(qp));
		2593	return isl_pw_qpolynomial_alloc(dom, qp);
		2594	}
		2595
		2596	#undef PW
		2597	#define PW isl_pw_qpolynomial
		2598	#undef EL
		2599	#define EL isl_qpolynomial
		2600	#undef EL_IS_ZERO
		2601	#define EL_IS_ZERO is_zero
		2602	#undef ZERO
		2603	#define ZERO zero
		2604	#undef IS_ZERO
		2605	#define IS_ZERO is_zero
		2606	#undef FIELD
		2607	#define FIELD qp
		2608	#undef DEFAULT_IS_ZERO
		2609	#define DEFAULT_IS_ZERO 1
		2610
		2611	#include <isl_pw_templ.c>
		2612
		2613	#undef UNION
		2614	#define UNION isl_union_pw_qpolynomial
		2615	#undef PART
		2616	#define PART isl_pw_qpolynomial
		2617	#undef PARTS
		2618	#define PARTS pw_qpolynomial
		2619	#define ALIGN_DOMAIN
		2620
		2621	#include <isl_union_templ.c>
		2622
		2623	int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
		2624	{
		2625	if (!pwqp)
		2626	return -1;
		2627
		2628	if (pwqp->n != -1)
		2629	return 0;
		2630
		2631	if (!isl_set_plain_is_universe(pwqp->p[0].set))
		2632	return 0;
		2633
		2634	return isl_qpolynomial_is_one(pwqp->p[0].qp);
		2635	}
		2636
		2637	__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
		2638	__isl_take isl_pw_qpolynomial *pwqp1,
		2639	__isl_take isl_pw_qpolynomial *pwqp2)
		2640	{
		2641	return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
		2642	}
		2643
		2644	__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
		2645	__isl_take isl_pw_qpolynomial *pwqp1,
		2646	__isl_take isl_pw_qpolynomial *pwqp2)
		2647	{
		2648	int i, j, n;
		2649	struct isl_pw_qpolynomial *res;
		2650
		2651	if (!pwqp1 || !pwqp2)
		2652	goto error;
		2653
		2654	isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
		2655	goto error);
		2656
		2657	if (isl_pw_qpolynomial_is_zero(pwqp1)) {
		2658	isl_pw_qpolynomial_free(pwqp2);
		2659	return pwqp1;
		2660	}
		2661
		2662	if (isl_pw_qpolynomial_is_zero(pwqp2)) {
		2663	isl_pw_qpolynomial_free(pwqp1);
		2664	return pwqp2;
		2665	}
		2666
		2667	if (isl_pw_qpolynomial_is_one(pwqp1)) {
		2668	isl_pw_qpolynomial_free(pwqp1);
		2669	return pwqp2;
		2670	}
		2671
		2672	if (isl_pw_qpolynomial_is_one(pwqp2)) {
		2673	isl_pw_qpolynomial_free(pwqp2);
		2674	return pwqp1;
		2675	}
		2676
		2677	n = pwqp1->n * pwqp2->n;
		2678	res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
		2679
		2680	for (i = 0; i < pwqp1->n; ++i) {
		2681	for (j = 0; j < pwqp2->n; ++j) {
		2682	struct isl_set *common;
		2683	struct isl_qpolynomial *prod;
		2684	common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
		2685	isl_set_copy(pwqp2->p[j].set));
		2686	if (isl_set_plain_is_empty(common)) {
		2687	isl_set_free(common);
		2688	continue;
		2689	}
		2690
		2691	prod = isl_qpolynomial_mul(
		2692	isl_qpolynomial_copy(pwqp1->p[i].qp),
		2693	isl_qpolynomial_copy(pwqp2->p[j].qp));
		2694
		2695	res = isl_pw_qpolynomial_add_piece(res, common, prod);
		2696	}
		2697	}
		2698
		2699	isl_pw_qpolynomial_free(pwqp1);
		2700	isl_pw_qpolynomial_free(pwqp2);
		2701
		2702	return res;
		2703	error:
		2704	isl_pw_qpolynomial_free(pwqp1);
		2705	isl_pw_qpolynomial_free(pwqp2);
		2706	return NULL;
		2707	}
		2708
		2709	__isl_give struct isl_upoly *isl_upoly_eval(
		2710	__isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
		2711	{
		2712	int i;
		2713	struct isl_upoly_rec *rec;
		2714	struct isl_upoly *res;
		2715	struct isl_upoly *base;
		2716
		2717	if (isl_upoly_is_cst(up)) {
		2718	isl_vec_free(vec);
		2719	return up;
		2720	}
		2721
		2722	rec = isl_upoly_as_rec(up);
		2723	if (!rec)
		2724	goto error;
		2725
		2726	isl_assert(up->ctx, rec->n >= 1, goto error);
		2727
		2728	base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
		2729
		2730	res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
		2731	isl_vec_copy(vec));
		2732
		2733	for (i = rec->n - 2; i >= 0; --i) {
		2734	res = isl_upoly_mul(res, isl_upoly_copy(base));
		2735	res = isl_upoly_sum(res,
		2736	isl_upoly_eval(isl_upoly_copy(rec->p[i]),
		2737	isl_vec_copy(vec)));
		2738	}
		2739
		2740	isl_upoly_free(base);
		2741	isl_upoly_free(up);
		2742	isl_vec_free(vec);
		2743	return res;
		2744	error:
		2745	isl_upoly_free(up);
		2746	isl_vec_free(vec);
		2747	return NULL;
		2748	}
		2749
		2750	__isl_give isl_qpolynomial *isl_qpolynomial_eval(
		2751	__isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
		2752	{
		2753	isl_vec *ext;
		2754	struct isl_upoly *up;
		2755	isl_space *dim;
		2756
		2757	if (!qp || !pnt)
		2758	goto error;
		2759	isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
		2760
		2761	if (qp->div->n_row == 0)
		2762	ext = isl_vec_copy(pnt->vec);
		2763	else {
		2764	int i;
		2765	unsigned dim = isl_space_dim(qp->dim, isl_dim_all);
		2766	ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
		2767	if (!ext)
		2768	goto error;
		2769
		2770	isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
		2771	for (i = 0; i < qp->div->n_row; ++i) {
		2772	isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
		2773	1 + dim + i, &ext->el[1+dim+i]);
		2774	isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
		2775	qp->div->row[i][0]);
		2776	}
		2777	}
		2778
		2779	up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
		2780	if (!up)
		2781	goto error;
		2782
		2783	dim = isl_space_copy(qp->dim);
		2784	isl_qpolynomial_free(qp);
		2785	isl_point_free(pnt);
		2786
		2787	return isl_qpolynomial_alloc(dim, 0, up);
		2788	error:
		2789	isl_qpolynomial_free(qp);
		2790	isl_point_free(pnt);
		2791	return NULL;
		2792	}
		2793
		2794	int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
		2795	__isl_keep struct isl_upoly_cst *cst2)
		2796	{
		2797	int cmp;
		2798	isl_int t;
		2799	isl_int_init(t);
		2800	isl_int_mul(t, cst1->n, cst2->d);
		2801	isl_int_submul(t, cst2->n, cst1->d);
		2802	cmp = isl_int_sgn(t);
		2803	isl_int_clear(t);
		2804	return cmp;
		2805	}
		2806
		2807	int isl_qpolynomial_le_cst(__isl_keep isl_qpolynomial *qp1,
		2808	__isl_keep isl_qpolynomial *qp2)
		2809	{
		2810	struct isl_upoly_cst *cst1, *cst2;
		2811
		2812	if (!qp1 || !qp2)
		2813	return -1;
		2814	isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), return -1);
		2815	isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), return -1);
		2816	if (isl_qpolynomial_is_nan(qp1))
		2817	return -1;
		2818	if (isl_qpolynomial_is_nan(qp2))
		2819	return -1;
		2820	cst1 = isl_upoly_as_cst(qp1->upoly);
		2821	cst2 = isl_upoly_as_cst(qp2->upoly);
		2822
		2823	return isl_upoly_cmp(cst1, cst2) <= 0;
		2824	}
		2825
		2826	__isl_give isl_qpolynomial *isl_qpolynomial_min_cst(
		2827	__isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
		2828	{
		2829	struct isl_upoly_cst *cst1, *cst2;
		2830	int cmp;
		2831
		2832	if (!qp1 || !qp2)
		2833	goto error;
		2834	isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
		2835	isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
		2836	cst1 = isl_upoly_as_cst(qp1->upoly);
		2837	cst2 = isl_upoly_as_cst(qp2->upoly);
		2838	cmp = isl_upoly_cmp(cst1, cst2);
		2839
		2840	if (cmp <= 0) {
		2841	isl_qpolynomial_free(qp2);
		2842	} else {
		2843	isl_qpolynomial_free(qp1);
		2844	qp1 = qp2;
		2845	}
		2846	return qp1;
		2847	error:
		2848	isl_qpolynomial_free(qp1);
		2849	isl_qpolynomial_free(qp2);
		2850	return NULL;
		2851	}
		2852
		2853	__isl_give isl_qpolynomial *isl_qpolynomial_max_cst(
		2854	__isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
		2855	{
		2856	struct isl_upoly_cst *cst1, *cst2;
		2857	int cmp;
		2858
		2859	if (!qp1 || !qp2)
		2860	goto error;
		2861	isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
		2862	isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
		2863	cst1 = isl_upoly_as_cst(qp1->upoly);
		2864	cst2 = isl_upoly_as_cst(qp2->upoly);
		2865	cmp = isl_upoly_cmp(cst1, cst2);
		2866
		2867	if (cmp >= 0) {
		2868	isl_qpolynomial_free(qp2);
		2869	} else {
		2870	isl_qpolynomial_free(qp1);
		2871	qp1 = qp2;
		2872	}
		2873	return qp1;
		2874	error:
		2875	isl_qpolynomial_free(qp1);
		2876	isl_qpolynomial_free(qp2);
		2877	return NULL;
		2878	}
		2879
		2880	__isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
		2881	__isl_take isl_qpolynomial *qp, enum isl_dim_type type,
		2882	unsigned first, unsigned n)
		2883	{
		2884	unsigned total;
		2885	unsigned g_pos;
		2886	int *exp;
		2887
		2888	if (!qp)
		2889	return NULL;
		2890	if (type == isl_dim_out)
		2891	isl_die(qp->div->ctx, isl_error_invalid,
		2892	"cannot insert output/set dimensions",
		2893	goto error);
		2894	if (type == isl_dim_in)
		2895	type = isl_dim_set;
		2896	if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
		2897	return qp;
		2898
		2899	qp = isl_qpolynomial_cow(qp);
		2900	if (!qp)
		2901	return NULL;
		2902
		2903	isl_assert(qp->div->ctx, first <= isl_space_dim(qp->dim, type),
		2904	goto error);
		2905
		2906	g_pos = pos(qp->dim, type) + first;
		2907
		2908	qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
		2909	if (!qp->div)
		2910	goto error;
		2911
		2912	total = qp->div->n_col - 2;
		2913	if (total > g_pos) {
		2914	int i;
		2915	exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
		2916	if (!exp)
		2917	goto error;
		2918	for (i = 0; i < total - g_pos; ++i)
		2919	exp[i] = i + n;
		2920	qp->upoly = expand(qp->upoly, exp, g_pos);
		2921	free(exp);
		2922	if (!qp->upoly)
		2923	goto error;
		2924	}
		2925
		2926	qp->dim = isl_space_insert_dims(qp->dim, type, first, n);
		2927	if (!qp->dim)
		2928	goto error;
		2929
		2930	return qp;
		2931	error:
		2932	isl_qpolynomial_free(qp);
		2933	return NULL;
		2934	}
		2935
		2936	__isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
		2937	__isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
		2938	{
		2939	unsigned pos;
		2940
		2941	pos = isl_qpolynomial_dim(qp, type);
		2942
		2943	return isl_qpolynomial_insert_dims(qp, type, pos, n);
		2944	}
		2945
		2946	__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
		2947	__isl_take isl_pw_qpolynomial *pwqp,
		2948	enum isl_dim_type type, unsigned n)
		2949	{
		2950	unsigned pos;
		2951
		2952	pos = isl_pw_qpolynomial_dim(pwqp, type);
		2953
		2954	return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
		2955	}
		2956
		2957	static int *reordering_move(isl_ctx *ctx,
		2958	unsigned len, unsigned dst, unsigned src, unsigned n)
		2959	{
		2960	int i;
		2961	int *reordering;
		2962
		2963	reordering = isl_alloc_array(ctx, int, len);
		2964	if (!reordering)
		2965	return NULL;
		2966
		2967	if (dst <= src) {
		2968	for (i = 0; i < dst; ++i)
		2969	reordering[i] = i;
		2970	for (i = 0; i < n; ++i)
		2971	reordering[src + i] = dst + i;
		2972	for (i = 0; i < src - dst; ++i)
		2973	reordering[dst + i] = dst + n + i;
		2974	for (i = 0; i < len - src - n; ++i)
		2975	reordering[src + n + i] = src + n + i;
		2976	} else {
		2977	for (i = 0; i < src; ++i)
		2978	reordering[i] = i;
		2979	for (i = 0; i < n; ++i)
		2980	reordering[src + i] = dst + i;
		2981	for (i = 0; i < dst - src; ++i)
		2982	reordering[src + n + i] = src + i;
		2983	for (i = 0; i < len - dst - n; ++i)
		2984	reordering[dst + n + i] = dst + n + i;
		2985	}
		2986
		2987	return reordering;
		2988	}
		2989
		2990	__isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
		2991	__isl_take isl_qpolynomial *qp,
		2992	enum isl_dim_type dst_type, unsigned dst_pos,
		2993	enum isl_dim_type src_type, unsigned src_pos, unsigned n)
		2994	{
		2995	unsigned g_dst_pos;
		2996	unsigned g_src_pos;
		2997	int *reordering;
		2998
		2999	qp = isl_qpolynomial_cow(qp);
		3000	if (!qp)
		3001	return NULL;
		3002
		3003	if (dst_type == isl_dim_out || src_type == isl_dim_out)
		3004	isl_die(qp->dim->ctx, isl_error_invalid,
		3005	"cannot move output/set dimension",
		3006	goto error);
		3007	if (dst_type == isl_dim_in)
		3008	dst_type = isl_dim_set;
		3009	if (src_type == isl_dim_in)
		3010	src_type = isl_dim_set;
		3011
		3012	isl_assert(qp->dim->ctx, src_pos + n <= isl_space_dim(qp->dim, src_type),
		3013	goto error);
		3014
		3015	g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
		3016	g_src_pos = pos(qp->dim, src_type) + src_pos;
		3017	if (dst_type > src_type)
		3018	g_dst_pos -= n;
		3019
		3020	qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
		3021	if (!qp->div)
		3022	goto error;
		3023	qp = sort_divs(qp);
		3024	if (!qp)
		3025	goto error;
		3026
		3027	reordering = reordering_move(qp->dim->ctx,
		3028	qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
		3029	if (!reordering)
		3030	goto error;
		3031
		3032	qp->upoly = reorder(qp->upoly, reordering);
		3033	free(reordering);
		3034	if (!qp->upoly)
		3035	goto error;
		3036
		3037	qp->dim = isl_space_move_dims(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
		3038	if (!qp->dim)
		3039	goto error;
		3040
		3041	return qp;
		3042	error:
		3043	isl_qpolynomial_free(qp);
		3044	return NULL;
		3045	}
		3046
		3047	__isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_space *dim,
		3048	isl_int *f, isl_int denom)
		3049	{
		3050	struct isl_upoly *up;
		3051
		3052	dim = isl_space_domain(dim);
		3053	if (!dim)
		3054	return NULL;
		3055
		3056	up = isl_upoly_from_affine(dim->ctx, f, denom,
		3057	1 + isl_space_dim(dim, isl_dim_all));
		3058
		3059	return isl_qpolynomial_alloc(dim, 0, up);
		3060	}
		3061
		3062	__isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
		3063	{
		3064	isl_ctx *ctx;
		3065	struct isl_upoly *up;
		3066	isl_qpolynomial *qp;
		3067
		3068	if (!aff)
		3069	return NULL;
		3070
		3071	ctx = isl_aff_get_ctx(aff);
		3072	up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
		3073	aff->v->size - 1);
		3074
		3075	qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
		3076	aff->ls->div->n_row, up);
		3077	if (!qp)
		3078	goto error;
		3079
		3080	isl_mat_free(qp->div);
		3081	qp->div = isl_mat_copy(aff->ls->div);
		3082	qp->div = isl_mat_cow(qp->div);
		3083	if (!qp->div)
		3084	goto error;
		3085
		3086	isl_aff_free(aff);
		3087	qp = reduce_divs(qp);
		3088	qp = remove_redundant_divs(qp);
		3089	return qp;
		3090	error:
		3091	isl_aff_free(aff);
		3092	return NULL;
		3093	}
		3094
		3095	__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
		3096	__isl_take isl_pw_aff *pwaff)
		3097	{
		3098	int i;
		3099	isl_pw_qpolynomial *pwqp;
		3100
		3101	if (!pwaff)
		3102	return NULL;
		3103
		3104	pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
		3105	pwaff->n);
		3106
		3107	for (i = 0; i < pwaff->n; ++i) {
		3108	isl_set *dom;
		3109	isl_qpolynomial *qp;
		3110
		3111	dom = isl_set_copy(pwaff->p[i].set);
		3112	qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
		3113	pwqp = isl_pw_qpolynomial_add_piece(pwqp, dom, qp);
		3114	}
		3115
		3116	isl_pw_aff_free(pwaff);
		3117	return pwqp;
		3118	}
		3119
		3120	__isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
		3121	__isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
		3122	{
		3123	isl_aff *aff;
		3124
		3125	aff = isl_constraint_get_bound(c, type, pos);
		3126	isl_constraint_free(c);
		3127	return isl_qpolynomial_from_aff(aff);
		3128	}
		3129
		3130	/* For each 0 <= i < "n", replace variable "first" + i of type "type"
		3131	* in "qp" by subs[i].
		3132	*/
		3133	__isl_give isl_qpolynomial *isl_qpolynomial_substitute(
		3134	__isl_take isl_qpolynomial *qp,
		3135	enum isl_dim_type type, unsigned first, unsigned n,
		3136	__isl_keep isl_qpolynomial **subs)
		3137	{
		3138	int i;
		3139	struct isl_upoly **ups;
		3140
		3141	if (n == 0)
		3142	return qp;
		3143
		3144	qp = isl_qpolynomial_cow(qp);
		3145	if (!qp)
		3146	return NULL;
		3147
		3148	if (type == isl_dim_out)
		3149	isl_die(qp->dim->ctx, isl_error_invalid,
		3150	"cannot substitute output/set dimension",
		3151	goto error);
		3152	if (type == isl_dim_in)
		3153	type = isl_dim_set;
		3154
		3155	for (i = 0; i < n; ++i)
		3156	if (!subs[i])
		3157	goto error;
		3158
		3159	isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
		3160	goto error);
		3161
		3162	for (i = 0; i < n; ++i)
		3163	isl_assert(qp->dim->ctx, isl_space_is_equal(qp->dim, subs[i]->dim),
		3164	goto error);
		3165
		3166	isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
		3167	for (i = 0; i < n; ++i)
		3168	isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
		3169
		3170	first += pos(qp->dim, type);
		3171
		3172	ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
		3173	if (!ups)
		3174	goto error;
		3175	for (i = 0; i < n; ++i)
		3176	ups[i] = subs[i]->upoly;
		3177
		3178	qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
		3179
		3180	free(ups);
		3181
		3182	if (!qp->upoly)
		3183	goto error;
		3184
		3185	return qp;
		3186	error:
		3187	isl_qpolynomial_free(qp);
		3188	return NULL;
		3189	}
		3190
		3191	/* Extend "bset" with extra set dimensions for each integer division
		3192	* in "qp" and then call "fn" with the extended bset and the polynomial
		3193	* that results from replacing each of the integer divisions by the
		3194	* corresponding extra set dimension.
		3195	*/
		3196	int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
		3197	__isl_keep isl_basic_set *bset,
		3198	int (*fn)(__isl_take isl_basic_set *bset,
		3199	__isl_take isl_qpolynomial *poly, void *user), void *user)
		3200	{
		3201	isl_space *dim;
		3202	isl_mat *div;
		3203	isl_qpolynomial *poly;
		3204
		3205	if (!qp || !bset)
		3206	goto error;
		3207	if (qp->div->n_row == 0)
		3208	return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
		3209	user);
		3210
		3211	div = isl_mat_copy(qp->div);
		3212	dim = isl_space_copy(qp->dim);
		3213	dim = isl_space_add_dims(dim, isl_dim_set, qp->div->n_row);
		3214	poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
		3215	bset = isl_basic_set_copy(bset);
		3216	bset = isl_basic_set_add(bset, isl_dim_set, qp->div->n_row);
		3217	bset = add_div_constraints(bset, div);
		3218
		3219	return fn(bset, poly, user);
		3220	error:
		3221	return -1;
		3222	}
		3223
		3224	/* Return total degree in variables first (inclusive) up to last (exclusive).
		3225	*/
		3226	int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
		3227	{
		3228	int deg = -1;
		3229	int i;
		3230	struct isl_upoly_rec *rec;
		3231
		3232	if (!up)
		3233	return -2;
		3234	if (isl_upoly_is_zero(up))
		3235	return -1;
		3236	if (isl_upoly_is_cst(up) || up->var < first)
		3237	return 0;
		3238
		3239	rec = isl_upoly_as_rec(up);
		3240	if (!rec)
		3241	return -2;
		3242
		3243	for (i = 0; i < rec->n; ++i) {
		3244	int d;
		3245
		3246	if (isl_upoly_is_zero(rec->p[i]))
		3247	continue;
		3248	d = isl_upoly_degree(rec->p[i], first, last);
		3249	if (up->var < last)
		3250	d += i;
		3251	if (d > deg)
		3252	deg = d;
		3253	}
		3254
		3255	return deg;
		3256	}
		3257
		3258	/* Return total degree in set variables.
		3259	*/
		3260	int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
		3261	{
		3262	unsigned ovar;
		3263	unsigned nvar;
		3264
		3265	if (!poly)
		3266	return -2;
		3267
		3268	ovar = isl_space_offset(poly->dim, isl_dim_set);
		3269	nvar = isl_space_dim(poly->dim, isl_dim_set);
		3270	return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
		3271	}
		3272
		3273	__isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
		3274	unsigned pos, int deg)
		3275	{
		3276	int i;
		3277	struct isl_upoly_rec *rec;
		3278
		3279	if (!up)
		3280	return NULL;
		3281
		3282	if (isl_upoly_is_cst(up) || up->var < pos) {
		3283	if (deg == 0)
		3284	return isl_upoly_copy(up);
		3285	else
		3286	return isl_upoly_zero(up->ctx);
		3287	}
		3288
		3289	rec = isl_upoly_as_rec(up);
		3290	if (!rec)
		3291	return NULL;
		3292
		3293	if (up->var == pos) {
		3294	if (deg < rec->n)
		3295	return isl_upoly_copy(rec->p[deg]);
		3296	else
		3297	return isl_upoly_zero(up->ctx);
		3298	}
		3299
		3300	up = isl_upoly_copy(up);
		3301	up = isl_upoly_cow(up);
		3302	rec = isl_upoly_as_rec(up);
		3303	if (!rec)
		3304	goto error;
		3305
		3306	for (i = 0; i < rec->n; ++i) {
		3307	struct isl_upoly *t;
		3308	t = isl_upoly_coeff(rec->p[i], pos, deg);
		3309	if (!t)
		3310	goto error;
		3311	isl_upoly_free(rec->p[i]);
		3312	rec->p[i] = t;
		3313	}
		3314
		3315	return up;
		3316	error:
		3317	isl_upoly_free(up);
		3318	return NULL;
		3319	}
		3320
		3321	/* Return coefficient of power "deg" of variable "t_pos" of type "type".
		3322	*/
		3323	__isl_give isl_qpolynomial *isl_qpolynomial_coeff(
		3324	__isl_keep isl_qpolynomial *qp,
		3325	enum isl_dim_type type, unsigned t_pos, int deg)
		3326	{
		3327	unsigned g_pos;
		3328	struct isl_upoly *up;
		3329	isl_qpolynomial *c;
		3330
		3331	if (!qp)
		3332	return NULL;
		3333
		3334	if (type == isl_dim_out)
		3335	isl_die(qp->div->ctx, isl_error_invalid,
		3336	"output/set dimension does not have a coefficient",
		3337	return NULL);
		3338	if (type == isl_dim_in)
		3339	type = isl_dim_set;
		3340
		3341	isl_assert(qp->div->ctx, t_pos < isl_space_dim(qp->dim, type),
		3342	return NULL);
		3343
		3344	g_pos = pos(qp->dim, type) + t_pos;
		3345	up = isl_upoly_coeff(qp->upoly, g_pos, deg);
		3346
		3347	c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, up);
		3348	if (!c)
		3349	return NULL;
		3350	isl_mat_free(c->div);
		3351	c->div = isl_mat_copy(qp->div);
		3352	if (!c->div)
		3353	goto error;
		3354	return c;
		3355	error:
		3356	isl_qpolynomial_free(c);
		3357	return NULL;
		3358	}
		3359
		3360	/* Homogenize the polynomial in the variables first (inclusive) up to
		3361	* last (exclusive) by inserting powers of variable first.
		3362	* Variable first is assumed not to appear in the input.
		3363	*/
		3364	__isl_give struct isl_upoly *isl_upoly_homogenize(
		3365	__isl_take struct isl_upoly *up, int deg, int target,
		3366	int first, int last)
		3367	{
		3368	int i;
		3369	struct isl_upoly_rec *rec;
		3370
		3371	if (!up)
		3372	return NULL;
		3373	if (isl_upoly_is_zero(up))
		3374	return up;
		3375	if (deg == target)
		3376	return up;
		3377	if (isl_upoly_is_cst(up) || up->var < first) {
		3378	struct isl_upoly *hom;
		3379
		3380	hom = isl_upoly_var_pow(up->ctx, first, target - deg);
		3381	if (!hom)
		3382	goto error;
		3383	rec = isl_upoly_as_rec(hom);
		3384	rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
		3385
		3386	return hom;
		3387	}
		3388
		3389	up = isl_upoly_cow(up);
		3390	rec = isl_upoly_as_rec(up);
		3391	if (!rec)
		3392	goto error;
		3393
		3394	for (i = 0; i < rec->n; ++i) {
		3395	if (isl_upoly_is_zero(rec->p[i]))
		3396	continue;
		3397	rec->p[i] = isl_upoly_homogenize(rec->p[i],
		3398	up->var < last ? deg + i : i, target,
		3399	first, last);
		3400	if (!rec->p[i])
		3401	goto error;
		3402	}
		3403
		3404	return up;
		3405	error:
		3406	isl_upoly_free(up);
		3407	return NULL;
		3408	}
		3409
		3410	/* Homogenize the polynomial in the set variables by introducing
		3411	* powers of an extra set variable at position 0.
		3412	*/
		3413	__isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
		3414	__isl_take isl_qpolynomial *poly)
		3415	{
		3416	unsigned ovar;
		3417	unsigned nvar;
		3418	int deg = isl_qpolynomial_degree(poly);
		3419
		3420	if (deg < -1)
		3421	goto error;
		3422
		3423	poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
		3424	poly = isl_qpolynomial_cow(poly);
		3425	if (!poly)
		3426	goto error;
		3427
		3428	ovar = isl_space_offset(poly->dim, isl_dim_set);
		3429	nvar = isl_space_dim(poly->dim, isl_dim_set);
		3430	poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
		3431	ovar, ovar + nvar);
		3432	if (!poly->upoly)
		3433	goto error;
		3434
		3435	return poly;
		3436	error:
		3437	isl_qpolynomial_free(poly);
		3438	return NULL;
		3439	}
		3440
		3441	__isl_give isl_term *isl_term_alloc(__isl_take isl_space *dim,
		3442	__isl_take isl_mat *div)
		3443	{
		3444	isl_term *term;
		3445	int n;
		3446
		3447	if (!dim || !div)
		3448	goto error;
		3449
		3450	n = isl_space_dim(dim, isl_dim_all) + div->n_row;
		3451
		3452	term = isl_calloc(dim->ctx, struct isl_term,
		3453	sizeof(struct isl_term) + (n - 1) * sizeof(int));
		3454	if (!term)
		3455	goto error;
		3456
		3457	term->ref = 1;
		3458	term->dim = dim;
		3459	term->div = div;
		3460	isl_int_init(term->n);
		3461	isl_int_init(term->d);
		3462
		3463	return term;
		3464	error:
		3465	isl_space_free(dim);
		3466	isl_mat_free(div);
		3467	return NULL;
		3468	}
		3469
		3470	__isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
		3471	{
		3472	if (!term)
		3473	return NULL;
		3474
		3475	term->ref++;
		3476	return term;
		3477	}
		3478
		3479	__isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
		3480	{
		3481	int i;
		3482	isl_term *dup;
		3483	unsigned total;
		3484
		3485	if (!term)
		3486	return NULL;
		3487
		3488	total = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
		3489
		3490	dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
		3491	if (!dup)
		3492	return NULL;
		3493
		3494	isl_int_set(dup->n, term->n);
		3495	isl_int_set(dup->d, term->d);
		3496
		3497	for (i = 0; i < total; ++i)
		3498	dup->pow[i] = term->pow[i];
		3499
		3500	return dup;
		3501	}
		3502
		3503	__isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
		3504	{
		3505	if (!term)
		3506	return NULL;
		3507
		3508	if (term->ref == 1)
		3509	return term;
		3510	term->ref--;
		3511	return isl_term_dup(term);
		3512	}
		3513
		3514	void isl_term_free(__isl_take isl_term *term)
		3515	{
		3516	if (!term)
		3517	return;
		3518
		3519	if (--term->ref > 0)
		3520	return;
		3521
		3522	isl_space_free(term->dim);
		3523	isl_mat_free(term->div);
		3524	isl_int_clear(term->n);
		3525	isl_int_clear(term->d);
		3526	free(term);
		3527	}
		3528
		3529	unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
		3530	{
		3531	if (!term)
		3532	return 0;
		3533
		3534	switch (type) {
		3535	case isl_dim_param:
		3536	case isl_dim_in:
		3537	case isl_dim_out: return isl_space_dim(term->dim, type);
		3538	case isl_dim_div: return term->div->n_row;
		3539	case isl_dim_all: return isl_space_dim(term->dim, isl_dim_all) +
		3540	term->div->n_row;
		3541	default: return 0;
		3542	}
		3543	}
		3544
		3545	isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
		3546	{
		3547	return term ? term->dim->ctx : NULL;
		3548	}
		3549
		3550	void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
		3551	{
		3552	if (!term)
		3553	return;
		3554	isl_int_set(*n, term->n);
		3555	}
		3556
		3557	void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
		3558	{
		3559	if (!term)
		3560	return;
		3561	isl_int_set(*d, term->d);
		3562	}
		3563
		3564	int isl_term_get_exp(__isl_keep isl_term *term,
		3565	enum isl_dim_type type, unsigned pos)
		3566	{
		3567	if (!term)
		3568	return -1;
		3569
		3570	isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
		3571
		3572	if (type >= isl_dim_set)
		3573	pos += isl_space_dim(term->dim, isl_dim_param);
		3574	if (type >= isl_dim_div)
		3575	pos += isl_space_dim(term->dim, isl_dim_set);
		3576
		3577	return term->pow[pos];
		3578	}
		3579
		3580	__isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
		3581	{
		3582	isl_local_space *ls;
		3583	isl_aff *aff;
		3584
		3585	if (!term)
		3586	return NULL;
		3587
		3588	isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
		3589	return NULL);
		3590
		3591	ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
		3592	isl_mat_copy(term->div));
		3593	aff = isl_aff_alloc(ls);
		3594	if (!aff)
		3595	return NULL;
		3596
		3597	isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
		3598
		3599	aff = isl_aff_normalize(aff);
		3600
		3601	return aff;
		3602	}
		3603
		3604	__isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
		3605	int (*fn)(__isl_take isl_term *term, void *user),
		3606	__isl_take isl_term *term, void *user)
		3607	{
		3608	int i;
		3609	struct isl_upoly_rec *rec;
		3610
		3611	if (!up || !term)
		3612	goto error;
		3613
		3614	if (isl_upoly_is_zero(up))
		3615	return term;
		3616
		3617	isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
		3618	isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
		3619	isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
		3620
		3621	if (isl_upoly_is_cst(up)) {
		3622	struct isl_upoly_cst *cst;
		3623	cst = isl_upoly_as_cst(up);
		3624	if (!cst)
		3625	goto error;
		3626	term = isl_term_cow(term);
		3627	if (!term)
		3628	goto error;
		3629	isl_int_set(term->n, cst->n);
		3630	isl_int_set(term->d, cst->d);
		3631	if (fn(isl_term_copy(term), user) < 0)
		3632	goto error;
		3633	return term;
		3634	}
		3635
		3636	rec = isl_upoly_as_rec(up);
		3637	if (!rec)
		3638	goto error;
		3639
		3640	for (i = 0; i < rec->n; ++i) {
		3641	term = isl_term_cow(term);
		3642	if (!term)
		3643	goto error;
		3644	term->pow[up->var] = i;
		3645	term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
		3646	if (!term)
		3647	goto error;
		3648	}
		3649	term->pow[up->var] = 0;
		3650
		3651	return term;
		3652	error:
		3653	isl_term_free(term);
		3654	return NULL;
		3655	}
		3656
		3657	int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
		3658	int (*fn)(__isl_take isl_term *term, void *user), void *user)
		3659	{
		3660	isl_term *term;
		3661
		3662	if (!qp)
		3663	return -1;
		3664
		3665	term = isl_term_alloc(isl_space_copy(qp->dim), isl_mat_copy(qp->div));
		3666	if (!term)
		3667	return -1;
		3668
		3669	term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
		3670
		3671	isl_term_free(term);
		3672
		3673	return term ? 0 : -1;
		3674	}
		3675
		3676	__isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
		3677	{
		3678	struct isl_upoly *up;
		3679	isl_qpolynomial *qp;
		3680	int i, n;
		3681
		3682	if (!term)
		3683	return NULL;
		3684
		3685	n = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
		3686
		3687	up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
		3688	for (i = 0; i < n; ++i) {
		3689	if (!term->pow[i])
		3690	continue;
		3691	up = isl_upoly_mul(up,
		3692	isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
		3693	}
		3694
		3695	qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), term->div->n_row, up);
		3696	if (!qp)
		3697	goto error;
		3698	isl_mat_free(qp->div);
		3699	qp->div = isl_mat_copy(term->div);
		3700	if (!qp->div)
		3701	goto error;
		3702
		3703	isl_term_free(term);
		3704	return qp;
		3705	error:
		3706	isl_qpolynomial_free(qp);
		3707	isl_term_free(term);
		3708	return NULL;
		3709	}
		3710
		3711	__isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
		3712	__isl_take isl_space *dim)
		3713	{
		3714	int i;
		3715	int extra;
		3716	unsigned total;
		3717
		3718	if (!qp || !dim)
		3719	goto error;
		3720
		3721	if (isl_space_is_equal(qp->dim, dim)) {
		3722	isl_space_free(dim);
		3723	return qp;
		3724	}
		3725
		3726	qp = isl_qpolynomial_cow(qp);
		3727	if (!qp)
		3728	goto error;
		3729
		3730	extra = isl_space_dim(dim, isl_dim_set) -
		3731	isl_space_dim(qp->dim, isl_dim_set);
		3732	total = isl_space_dim(qp->dim, isl_dim_all);
		3733	if (qp->div->n_row) {
		3734	int *exp;
		3735
		3736	exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
		3737	if (!exp)
		3738	goto error;
		3739	for (i = 0; i < qp->div->n_row; ++i)
		3740	exp[i] = extra + i;
		3741	qp->upoly = expand(qp->upoly, exp, total);
		3742	free(exp);
		3743	if (!qp->upoly)
		3744	goto error;
		3745	}
		3746	qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
		3747	if (!qp->div)
		3748	goto error;
		3749	for (i = 0; i < qp->div->n_row; ++i)
		3750	isl_seq_clr(qp->div->row[i] + 2 + total, extra);
		3751
		3752	isl_space_free(qp->dim);
		3753	qp->dim = dim;
		3754
		3755	return qp;
		3756	error:
		3757	isl_space_free(dim);
		3758	isl_qpolynomial_free(qp);
		3759	return NULL;
		3760	}
		3761
		3762	/* For each parameter or variable that does not appear in qp,
		3763	* first eliminate the variable from all constraints and then set it to zero.
		3764	*/
		3765	static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
		3766	__isl_keep isl_qpolynomial *qp)
		3767	{
		3768	int *active = NULL;
		3769	int i;
		3770	int d;
		3771	unsigned nparam;
		3772	unsigned nvar;
		3773
		3774	if (!set || !qp)
		3775	goto error;
		3776
		3777	d = isl_space_dim(set->dim, isl_dim_all);
		3778	active = isl_calloc_array(set->ctx, int, d);
		3779	if (set_active(qp, active) < 0)
		3780	goto error;
		3781
		3782	for (i = 0; i < d; ++i)
		3783	if (!active[i])
		3784	break;
		3785
		3786	if (i == d) {
		3787	free(active);
		3788	return set;
		3789	}
		3790
		3791	nparam = isl_space_dim(set->dim, isl_dim_param);
		3792	nvar = isl_space_dim(set->dim, isl_dim_set);
		3793	for (i = 0; i < nparam; ++i) {
		3794	if (active[i])
		3795	continue;
		3796	set = isl_set_eliminate(set, isl_dim_param, i, 1);
		3797	set = isl_set_fix_si(set, isl_dim_param, i, 0);
		3798	}
		3799	for (i = 0; i < nvar; ++i) {
		3800	if (active[nparam + i])
		3801	continue;
		3802	set = isl_set_eliminate(set, isl_dim_set, i, 1);
		3803	set = isl_set_fix_si(set, isl_dim_set, i, 0);
		3804	}
		3805
		3806	free(active);
		3807
		3808	return set;
		3809	error:
		3810	free(active);
		3811	isl_set_free(set);
		3812	return NULL;
		3813	}
		3814
		3815	struct isl_opt_data {
		3816	isl_qpolynomial *qp;
		3817	int first;
		3818	isl_qpolynomial *opt;
		3819	int max;
		3820	};
		3821
		3822	static int opt_fn(__isl_take isl_point *pnt, void *user)
		3823	{
		3824	struct isl_opt_data *data = (struct isl_opt_data *)user;
		3825	isl_qpolynomial *val;
		3826
		3827	val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
		3828	if (data->first) {
		3829	data->first = 0;
		3830	data->opt = val;
		3831	} else if (data->max) {
		3832	data->opt = isl_qpolynomial_max_cst(data->opt, val);
		3833	} else {
		3834	data->opt = isl_qpolynomial_min_cst(data->opt, val);
		3835	}
		3836
		3837	return 0;
		3838	}
		3839
		3840	__isl_give isl_qpolynomial *isl_qpolynomial_opt_on_domain(
		3841	__isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
		3842	{
		3843	struct isl_opt_data data = { NULL, 1, NULL, max };
		3844
		3845	if (!set || !qp)
		3846	goto error;
		3847
		3848	if (isl_upoly_is_cst(qp->upoly)) {
		3849	isl_set_free(set);
		3850	return qp;
		3851	}
		3852
		3853	set = fix_inactive(set, qp);
		3854
		3855	data.qp = qp;
		3856	if (isl_set_foreach_point(set, opt_fn, &data) < 0)
		3857	goto error;
		3858
		3859	if (data.first) {
		3860	isl_space *space = isl_qpolynomial_get_domain_space(qp);
		3861	data.opt = isl_qpolynomial_zero_on_domain(space);
		3862	}
		3863
		3864	isl_set_free(set);
		3865	isl_qpolynomial_free(qp);
		3866	return data.opt;
		3867	error:
		3868	isl_set_free(set);
		3869	isl_qpolynomial_free(qp);
		3870	isl_qpolynomial_free(data.opt);
		3871	return NULL;
		3872	}
		3873
		3874	__isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
		3875	__isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
		3876	{
		3877	int i;
		3878	int n_sub;
		3879	isl_ctx *ctx;
		3880	struct isl_upoly **subs;
		3881	isl_mat *mat, *diag;
		3882
		3883	qp = isl_qpolynomial_cow(qp);
		3884	if (!qp || !morph)
		3885	goto error;
		3886
		3887	ctx = qp->dim->ctx;
		3888	isl_assert(ctx, isl_space_is_equal(qp->dim, morph->dom->dim), goto error);
		3889
		3890	n_sub = morph->inv->n_row - 1;
		3891	if (morph->inv->n_row != morph->inv->n_col)
		3892	n_sub += qp->div->n_row;
		3893	subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
		3894	if (!subs)
		3895	goto error;
		3896
		3897	for (i = 0; 1 + i < morph->inv->n_row; ++i)
		3898	subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
		3899	morph->inv->row[0][0], morph->inv->n_col);
		3900	if (morph->inv->n_row != morph->inv->n_col)
		3901	for (i = 0; i < qp->div->n_row; ++i)
		3902	subs[morph->inv->n_row - 1 + i] =
		3903	isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
		3904
		3905	qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
		3906
		3907	for (i = 0; i < n_sub; ++i)
		3908	isl_upoly_free(subs[i]);
		3909	free(subs);
		3910
		3911	diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
		3912	mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
		3913	diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
		3914	mat = isl_mat_diagonal(mat, diag);
		3915	qp->div = isl_mat_product(qp->div, mat);
		3916	isl_space_free(qp->dim);
		3917	qp->dim = isl_space_copy(morph->ran->dim);
		3918
		3919	if (!qp->upoly || !qp->div || !qp->dim)
		3920	goto error;
		3921
		3922	isl_morph_free(morph);
		3923
		3924	return qp;
		3925	error:
		3926	isl_qpolynomial_free(qp);
		3927	isl_morph_free(morph);
		3928	return NULL;
		3929	}
		3930
		3931	static int neg_entry(void **entry, void *user)
		3932	{
		3933	isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
		3934
		3935	*pwqp = isl_pw_qpolynomial_neg(*pwqp);
		3936
		3937	return *pwqp ? 0 : -1;
		3938	}
		3939
		3940	__isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_neg(
		3941	__isl_take isl_union_pw_qpolynomial *upwqp)
		3942	{
		3943	upwqp = isl_union_pw_qpolynomial_cow(upwqp);
		3944	if (!upwqp)
		3945	return NULL;
		3946
		3947	if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
		3948	&neg_entry, NULL) < 0)
		3949	goto error;
		3950
		3951	return upwqp;
		3952	error:
		3953	isl_union_pw_qpolynomial_free(upwqp);
		3954	return NULL;
		3955	}
		3956
		3957	__isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_sub(
		3958	__isl_take isl_union_pw_qpolynomial *upwqp1,
		3959	__isl_take isl_union_pw_qpolynomial *upwqp2)
		3960	{
		3961	return isl_union_pw_qpolynomial_add(upwqp1,
		3962	isl_union_pw_qpolynomial_neg(upwqp2));
		3963	}
		3964
		3965	__isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
		3966	__isl_take isl_union_pw_qpolynomial *upwqp1,
		3967	__isl_take isl_union_pw_qpolynomial *upwqp2)
		3968	{
		3969	return match_bin_op(upwqp1, upwqp2, &isl_pw_qpolynomial_mul);
		3970	}
		3971
		3972	/* Reorder the columns of the given div definitions according to the
		3973	* given reordering.
		3974	*/
		3975	static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
		3976	__isl_take isl_reordering *r)
		3977	{
		3978	int i, j;
		3979	isl_mat *mat;
		3980	int extra;
		3981
		3982	if (!div || !r)
		3983	goto error;
		3984
		3985	extra = isl_space_dim(r->dim, isl_dim_all) + div->n_row - r->len;
		3986	mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
		3987	if (!mat)
		3988	goto error;
		3989
		3990	for (i = 0; i < div->n_row; ++i) {
		3991	isl_seq_cpy(mat->row[i], div->row[i], 2);
		3992	isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
		3993	for (j = 0; j < r->len; ++j)
		3994	isl_int_set(mat->row[i][2 + r->pos[j]],
		3995	div->row[i][2 + j]);
		3996	}
		3997
		3998	isl_reordering_free(r);
		3999	isl_mat_free(div);
		4000	return mat;
		4001	error:
		4002	isl_reordering_free(r);
		4003	isl_mat_free(div);
		4004	return NULL;
		4005	}
		4006
		4007	/* Reorder the dimension of "qp" according to the given reordering.
		4008	*/
		4009	__isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
		4010	__isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
		4011	{
		4012	qp = isl_qpolynomial_cow(qp);
		4013	if (!qp)
		4014	goto error;
		4015
		4016	r = isl_reordering_extend(r, qp->div->n_row);
		4017	if (!r)
		4018	goto error;
		4019
		4020	qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
		4021	if (!qp->div)
		4022	goto error;
		4023
		4024	qp->upoly = reorder(qp->upoly, r->pos);
		4025	if (!qp->upoly)
		4026	goto error;
		4027
		4028	qp = isl_qpolynomial_reset_domain_space(qp, isl_space_copy(r->dim));
		4029
		4030	isl_reordering_free(r);
		4031	return qp;
		4032	error:
		4033	isl_qpolynomial_free(qp);
		4034	isl_reordering_free(r);
		4035	return NULL;
		4036	}
		4037
		4038	__isl_give isl_qpolynomial *isl_qpolynomial_align_params(
		4039	__isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
		4040	{
		4041	if (!qp || !model)
		4042	goto error;
		4043
		4044	if (!isl_space_match(qp->dim, isl_dim_param, model, isl_dim_param)) {
		4045	isl_reordering *exp;
		4046
		4047	model = isl_space_drop_dims(model, isl_dim_in,
		4048	0, isl_space_dim(model, isl_dim_in));
		4049	model = isl_space_drop_dims(model, isl_dim_out,
		4050	0, isl_space_dim(model, isl_dim_out));
		4051	exp = isl_parameter_alignment_reordering(qp->dim, model);
		4052	exp = isl_reordering_extend_space(exp,
		4053	isl_qpolynomial_get_domain_space(qp));
		4054	qp = isl_qpolynomial_realign_domain(qp, exp);
		4055	}
		4056
		4057	isl_space_free(model);
		4058	return qp;
		4059	error:
		4060	isl_space_free(model);
		4061	isl_qpolynomial_free(qp);
		4062	return NULL;
		4063	}
		4064
		4065	struct isl_split_periods_data {
		4066	int max_periods;
		4067	isl_pw_qpolynomial *res;
		4068	};
		4069
		4070	/* Create a slice where the integer division "div" has the fixed value "v".
		4071	* In particular, if "div" refers to floor(f/m), then create a slice
		4072	*
		4073	* m v <= f <= m v + (m - 1)
		4074	*
		4075	* or
		4076	*
		4077	* f - m v >= 0
		4078	* -f + m v + (m - 1) >= 0
		4079	*/
		4080	static __isl_give isl_set *set_div_slice(__isl_take isl_space *dim,
		4081	__isl_keep isl_qpolynomial *qp, int div, isl_int v)
		4082	{
		4083	int total;
		4084	isl_basic_set *bset = NULL;
		4085	int k;
		4086
		4087	if (!dim || !qp)
		4088	goto error;
		4089
		4090	total = isl_space_dim(dim, isl_dim_all);
		4091	bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0, 0, 2);
		4092
		4093	k = isl_basic_set_alloc_inequality(bset);
		4094	if (k < 0)
		4095	goto error;
		4096	isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
		4097	isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
		4098
		4099	k = isl_basic_set_alloc_inequality(bset);
		4100	if (k < 0)
		4101	goto error;
		4102	isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
		4103	isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
		4104	isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
		4105	isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
		4106
		4107	isl_space_free(dim);
		4108	return isl_set_from_basic_set(bset);
		4109	error:
		4110	isl_basic_set_free(bset);
		4111	isl_space_free(dim);
		4112	return NULL;
		4113	}
		4114
		4115	static int split_periods(__isl_take isl_set *set,
		4116	__isl_take isl_qpolynomial *qp, void *user);
		4117
		4118	/* Create a slice of the domain "set" such that integer division "div"
		4119	* has the fixed value "v" and add the results to data->res,
		4120	* replacing the integer division by "v" in "qp".
		4121	*/
		4122	static int set_div(__isl_take isl_set *set,
		4123	__isl_take isl_qpolynomial *qp, int div, isl_int v,
		4124	struct isl_split_periods_data *data)
		4125	{
		4126	int i;
		4127	int total;
		4128	isl_set *slice;
		4129	struct isl_upoly *cst;
		4130
		4131	slice = set_div_slice(isl_set_get_space(set), qp, div, v);
		4132	set = isl_set_intersect(set, slice);
		4133
		4134	if (!qp)
		4135	goto error;
		4136
		4137	total = isl_space_dim(qp->dim, isl_dim_all);
		4138
		4139	for (i = div + 1; i < qp->div->n_row; ++i) {
		4140	if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
		4141	continue;
		4142	isl_int_addmul(qp->div->row[i][1],
		4143	qp->div->row[i][2 + total + div], v);
		4144	isl_int_set_si(qp->div->row[i][2 + total + div], 0);
		4145	}
		4146
		4147	cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
		4148	qp = substitute_div(qp, div, cst);
		4149
		4150	return split_periods(set, qp, data);
		4151	error:
		4152	isl_set_free(set);
		4153	isl_qpolynomial_free(qp);
		4154	return -1;
		4155	}
		4156
		4157	/* Split the domain "set" such that integer division "div"
		4158	* has a fixed value (ranging from "min" to "max") on each slice
		4159	* and add the results to data->res.
		4160	*/
		4161	static int split_div(__isl_take isl_set *set,
		4162	__isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
		4163	struct isl_split_periods_data *data)
		4164	{
		4165	for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
		4166	isl_set *set_i = isl_set_copy(set);
		4167	isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
		4168
		4169	if (set_div(set_i, qp_i, div, min, data) < 0)
		4170	goto error;
		4171	}
		4172	isl_set_free(set);
		4173	isl_qpolynomial_free(qp);
		4174	return 0;
		4175	error:
		4176	isl_set_free(set);
		4177	isl_qpolynomial_free(qp);
		4178	return -1;
		4179	}
		4180
		4181	/* If "qp" refers to any integer division
		4182	* that can only attain "max_periods" distinct values on "set"
		4183	* then split the domain along those distinct values.
		4184	* Add the results (or the original if no splitting occurs)
		4185	* to data->res.
		4186	*/
		4187	static int split_periods(__isl_take isl_set *set,
		4188	__isl_take isl_qpolynomial *qp, void *user)
		4189	{
		4190	int i;
		4191	isl_pw_qpolynomial *pwqp;
		4192	struct isl_split_periods_data *data;
		4193	isl_int min, max;
		4194	int total;
		4195	int r = 0;
		4196
		4197	data = (struct isl_split_periods_data *)user;
		4198
		4199	if (!set || !qp)
		4200	goto error;
		4201
		4202	if (qp->div->n_row == 0) {
		4203	pwqp = isl_pw_qpolynomial_alloc(set, qp);
		4204	data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
		4205	return 0;
		4206	}
		4207
		4208	isl_int_init(min);
		4209	isl_int_init(max);
		4210	total = isl_space_dim(qp->dim, isl_dim_all);
		4211	for (i = 0; i < qp->div->n_row; ++i) {
		4212	enum isl_lp_result lp_res;
		4213
		4214	if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
		4215	qp->div->n_row) != -1)
		4216	continue;
		4217
		4218	lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
		4219	set->ctx->one, &min, NULL, NULL);
		4220	if (lp_res == isl_lp_error)
		4221	goto error2;
		4222	if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
		4223	continue;
		4224	isl_int_fdiv_q(min, min, qp->div->row[i][0]);
		4225
		4226	lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
		4227	set->ctx->one, &max, NULL, NULL);
		4228	if (lp_res == isl_lp_error)
		4229	goto error2;
		4230	if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
		4231	continue;
		4232	isl_int_fdiv_q(max, max, qp->div->row[i][0]);
		4233
		4234	isl_int_sub(max, max, min);
		4235	if (isl_int_cmp_si(max, data->max_periods) < 0) {
		4236	isl_int_add(max, max, min);
		4237	break;
		4238	}
		4239	}
		4240
		4241	if (i < qp->div->n_row) {
		4242	r = split_div(set, qp, i, min, max, data);
		4243	} else {
		4244	pwqp = isl_pw_qpolynomial_alloc(set, qp);
		4245	data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
		4246	}
		4247
		4248	isl_int_clear(max);
		4249	isl_int_clear(min);
		4250
		4251	return r;
		4252	error2:
		4253	isl_int_clear(max);
		4254	isl_int_clear(min);
		4255	error:
		4256	isl_set_free(set);
		4257	isl_qpolynomial_free(qp);
		4258	return -1;
		4259	}
		4260
		4261	/* If any quasi-polynomial in pwqp refers to any integer division
		4262	* that can only attain "max_periods" distinct values on its domain
		4263	* then split the domain along those distinct values.
		4264	*/
		4265	__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
		4266	__isl_take isl_pw_qpolynomial *pwqp, int max_periods)
		4267	{
		4268	struct isl_split_periods_data data;
		4269
		4270	data.max_periods = max_periods;
		4271	data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
		4272
		4273	if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
		4274	goto error;
		4275
		4276	isl_pw_qpolynomial_free(pwqp);
		4277
		4278	return data.res;
		4279	error:
		4280	isl_pw_qpolynomial_free(data.res);
		4281	isl_pw_qpolynomial_free(pwqp);
		4282	return NULL;
		4283	}
		4284
		4285	/* Construct a piecewise quasipolynomial that is constant on the given
		4286	* domain. In particular, it is
		4287	* 0 if cst == 0
		4288	* 1 if cst == 1
		4289	* infinity if cst == -1
		4290	*/
		4291	static __isl_give isl_pw_qpolynomial *constant_on_domain(
		4292	__isl_take isl_basic_set *bset, int cst)
		4293	{
		4294	isl_space *dim;
		4295	isl_qpolynomial *qp;
		4296
		4297	if (!bset)
		4298	return NULL;
		4299
		4300	bset = isl_basic_set_params(bset);
		4301	dim = isl_basic_set_get_space(bset);
		4302	if (cst < 0)
		4303	qp = isl_qpolynomial_infty_on_domain(dim);
		4304	else if (cst == 0)
		4305	qp = isl_qpolynomial_zero_on_domain(dim);
		4306	else
		4307	qp = isl_qpolynomial_one_on_domain(dim);
		4308	return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
		4309	}
		4310
		4311	/* Factor bset, call fn on each of the factors and return the product.
		4312	*
		4313	* If no factors can be found, simply call fn on the input.
		4314	* Otherwise, construct the factors based on the factorizer,
		4315	* call fn on each factor and compute the product.
		4316	*/
		4317	static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
		4318	__isl_take isl_basic_set *bset,
		4319	__isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
		4320	{
		4321	int i, n;
		4322	isl_space *dim;
		4323	isl_set *set;
		4324	isl_factorizer *f;
		4325	isl_qpolynomial *qp;
		4326	isl_pw_qpolynomial *pwqp;
		4327	unsigned nparam;
		4328	unsigned nvar;
		4329
		4330	f = isl_basic_set_factorizer(bset);
		4331	if (!f)
		4332	goto error;
		4333	if (f->n_group == 0) {
		4334	isl_factorizer_free(f);
		4335	return fn(bset);
		4336	}
		4337
		4338	nparam = isl_basic_set_dim(bset, isl_dim_param);
		4339	nvar = isl_basic_set_dim(bset, isl_dim_set);
		4340
		4341	dim = isl_basic_set_get_space(bset);
		4342	dim = isl_space_domain(dim);
		4343	set = isl_set_universe(isl_space_copy(dim));
		4344	qp = isl_qpolynomial_one_on_domain(dim);
		4345	pwqp = isl_pw_qpolynomial_alloc(set, qp);
		4346
		4347	bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
		4348
		4349	for (i = 0, n = 0; i < f->n_group; ++i) {
		4350	isl_basic_set *bset_i;
		4351	isl_pw_qpolynomial *pwqp_i;
		4352
		4353	bset_i = isl_basic_set_copy(bset);
		4354	bset_i = isl_basic_set_drop_constraints_involving(bset_i,
		4355	nparam + n + f->len[i], nvar - n - f->len[i]);
		4356	bset_i = isl_basic_set_drop_constraints_involving(bset_i,
		4357	nparam, n);
		4358	bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
		4359	n + f->len[i], nvar - n - f->len[i]);
		4360	bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
		4361
		4362	pwqp_i = fn(bset_i);
		4363	pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
		4364
		4365	n += f->len[i];
		4366	}
		4367
		4368	isl_basic_set_free(bset);
		4369	isl_factorizer_free(f);
		4370
		4371	return pwqp;
		4372	error:
		4373	isl_basic_set_free(bset);
		4374	return NULL;
		4375	}
		4376
		4377	/* Factor bset, call fn on each of the factors and return the product.
		4378	* The function is assumed to evaluate to zero on empty domains,
		4379	* to one on zero-dimensional domains and to infinity on unbounded domains
		4380	* and will not be called explicitly on zero-dimensional or unbounded domains.
		4381	*
		4382	* We first check for some special cases and remove all equalities.
		4383	* Then we hand over control to compressed_multiplicative_call.
		4384	*/
		4385	__isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
		4386	__isl_take isl_basic_set *bset,
		4387	__isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
		4388	{
		4389	int bounded;
		4390	isl_morph *morph;
		4391	isl_pw_qpolynomial *pwqp;
		4392
		4393	if (!bset)
		4394	return NULL;
		4395
		4396	if (isl_basic_set_plain_is_empty(bset))
		4397	return constant_on_domain(bset, 0);
		4398
		4399	if (isl_basic_set_dim(bset, isl_dim_set) == 0)
		4400	return constant_on_domain(bset, 1);
		4401
		4402	bounded = isl_basic_set_is_bounded(bset);
		4403	if (bounded < 0)
		4404	goto error;
		4405	if (!bounded)
		4406	return constant_on_domain(bset, -1);
		4407
		4408	if (bset->n_eq == 0)
		4409	return compressed_multiplicative_call(bset, fn);
		4410
		4411	morph = isl_basic_set_full_compression(bset);
		4412	bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
		4413
		4414	pwqp = compressed_multiplicative_call(bset, fn);
		4415
		4416	morph = isl_morph_dom_params(morph);
		4417	morph = isl_morph_ran_params(morph);
		4418	morph = isl_morph_inverse(morph);
		4419
		4420	pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
		4421
		4422	return pwqp;
		4423	error:
		4424	isl_basic_set_free(bset);
		4425	return NULL;
		4426	}
		4427
		4428	/* Drop all floors in "qp", turning each integer division [a/m] into
		4429	* a rational division a/m. If "down" is set, then the integer division
		4430	* is replaces by (a-(m-1))/m instead.
		4431	*/
		4432	static __isl_give isl_qpolynomial *qp_drop_floors(
		4433	__isl_take isl_qpolynomial *qp, int down)
		4434	{
		4435	int i;
		4436	struct isl_upoly *s;
		4437
		4438	if (!qp)
		4439	return NULL;
		4440	if (qp->div->n_row == 0)
		4441	return qp;
		4442
		4443	qp = isl_qpolynomial_cow(qp);
		4444	if (!qp)
		4445	return NULL;
		4446
		4447	for (i = qp->div->n_row - 1; i >= 0; --i) {
		4448	if (down) {
		4449	isl_int_sub(qp->div->row[i][1],
		4450	qp->div->row[i][1], qp->div->row[i][0]);
		4451	isl_int_add_ui(qp->div->row[i][1],
		4452	qp->div->row[i][1], 1);
		4453	}
		4454	s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
		4455	qp->div->row[i][0], qp->div->n_col - 1);
		4456	qp = substitute_div(qp, i, s);
		4457	if (!qp)
		4458	return NULL;
		4459	}
		4460
		4461	return qp;
		4462	}
		4463
		4464	/* Drop all floors in "pwqp", turning each integer division [a/m] into
		4465	* a rational division a/m.
		4466	*/
		4467	static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
		4468	__isl_take isl_pw_qpolynomial *pwqp)
		4469	{
		4470	int i;
		4471
		4472	if (!pwqp)
		4473	return NULL;
		4474
		4475	if (isl_pw_qpolynomial_is_zero(pwqp))
		4476	return pwqp;
		4477
		4478	pwqp = isl_pw_qpolynomial_cow(pwqp);
		4479	if (!pwqp)
		4480	return NULL;
		4481
		4482	for (i = 0; i < pwqp->n; ++i) {
		4483	pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
		4484	if (!pwqp->p[i].qp)
		4485	goto error;
		4486	}
		4487
		4488	return pwqp;
		4489	error:
		4490	isl_pw_qpolynomial_free(pwqp);
		4491	return NULL;
		4492	}
		4493
		4494	/* Adjust all the integer divisions in "qp" such that they are at least
		4495	* one over the given orthant (identified by "signs"). This ensures
		4496	* that they will still be non-negative even after subtracting (m-1)/m.
		4497	*
		4498	* In particular, f is replaced by f' + v, changing f = [a/m]
		4499	* to f' = [(a - m v)/m].
		4500	* If the constant term k in a is smaller than m,
		4501	* the constant term of v is set to floor(k/m) - 1.
		4502	* For any other term, if the coefficient c and the variable x have
		4503	* the same sign, then no changes are needed.
		4504	* Otherwise, if the variable is positive (and c is negative),
		4505	* then the coefficient of x in v is set to floor(c/m).
		4506	* If the variable is negative (and c is positive),
		4507	* then the coefficient of x in v is set to ceil(c/m).
		4508	*/
		4509	static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
		4510	int *signs)
		4511	{
		4512	int i, j;
		4513	int total;
		4514	isl_vec *v = NULL;
		4515	struct isl_upoly *s;
		4516
		4517	qp = isl_qpolynomial_cow(qp);
		4518	if (!qp)
		4519	return NULL;
		4520	qp->div = isl_mat_cow(qp->div);
		4521	if (!qp->div)
		4522	goto error;
		4523
		4524	total = isl_space_dim(qp->dim, isl_dim_all);
		4525	v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
		4526
		4527	for (i = 0; i < qp->div->n_row; ++i) {
		4528	isl_int *row = qp->div->row[i];
		4529	v = isl_vec_clr(v);
		4530	if (!v)
		4531	goto error;
		4532	if (isl_int_lt(row[1], row[0])) {
		4533	isl_int_fdiv_q(v->el[0], row[1], row[0]);
		4534	isl_int_sub_ui(v->el[0], v->el[0], 1);
		4535	isl_int_submul(row[1], row[0], v->el[0]);
		4536	}
		4537	for (j = 0; j < total; ++j) {
		4538	if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
		4539	continue;
		4540	if (signs[j] < 0)
		4541	isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
		4542	else
		4543	isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
		4544	isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
		4545	}
		4546	for (j = 0; j < i; ++j) {
		4547	if (isl_int_sgn(row[2 + total + j]) >= 0)
		4548	continue;
		4549	isl_int_fdiv_q(v->el[1 + total + j],
		4550	row[2 + total + j], row[0]);
		4551	isl_int_submul(row[2 + total + j],
		4552	row[0], v->el[1 + total + j]);
		4553	}
		4554	for (j = i + 1; j < qp->div->n_row; ++j) {
		4555	if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
		4556	continue;
		4557	isl_seq_combine(qp->div->row[j] + 1,
		4558	qp->div->ctx->one, qp->div->row[j] + 1,
		4559	qp->div->row[j][2 + total + i], v->el, v->size);
		4560	}
		4561	isl_int_set_si(v->el[1 + total + i], 1);
		4562	s = isl_upoly_from_affine(qp->dim->ctx, v->el,
		4563	qp->div->ctx->one, v->size);
		4564	qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
		4565	isl_upoly_free(s);
		4566	if (!qp->upoly)
		4567	goto error;
		4568	}
		4569
		4570	isl_vec_free(v);
		4571	return qp;
		4572	error:
		4573	isl_vec_free(v);
		4574	isl_qpolynomial_free(qp);
		4575	return NULL;
		4576	}
		4577
		4578	struct isl_to_poly_data {
		4579	int sign;
		4580	isl_pw_qpolynomial *res;
		4581	isl_qpolynomial *qp;
		4582	};
		4583
		4584	/* Appoximate data->qp by a polynomial on the orthant identified by "signs".
		4585	* We first make all integer divisions positive and then split the
		4586	* quasipolynomials into terms with sign data->sign (the direction
		4587	* of the requested approximation) and terms with the opposite sign.
		4588	* In the first set of terms, each integer division [a/m] is
		4589	* overapproximated by a/m, while in the second it is underapproximated
		4590	* by (a-(m-1))/m.
		4591	*/
		4592	static int to_polynomial_on_orthant(__isl_take isl_set *orthant, int *signs,
		4593	void *user)
		4594	{
		4595	struct isl_to_poly_data *data = user;
		4596	isl_pw_qpolynomial *t;
		4597	isl_qpolynomial *qp, *up, *down;
		4598
		4599	qp = isl_qpolynomial_copy(data->qp);
		4600	qp = make_divs_pos(qp, signs);
		4601
		4602	up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
		4603	up = qp_drop_floors(up, 0);
		4604	down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
		4605	down = qp_drop_floors(down, 1);
		4606
		4607	isl_qpolynomial_free(qp);
		4608	qp = isl_qpolynomial_add(up, down);
		4609
		4610	t = isl_pw_qpolynomial_alloc(orthant, qp);
		4611	data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
		4612
		4613	return 0;
		4614	}
		4615
		4616	/* Approximate each quasipolynomial by a polynomial. If "sign" is positive,
		4617	* the polynomial will be an overapproximation. If "sign" is negative,
		4618	* it will be an underapproximation. If "sign" is zero, the approximation
		4619	* will lie somewhere in between.
		4620	*
		4621	* In particular, is sign == 0, we simply drop the floors, turning
		4622	* the integer divisions into rational divisions.
		4623	* Otherwise, we split the domains into orthants, make all integer divisions
		4624	* positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
		4625	* depending on the requested sign and the sign of the term in which
		4626	* the integer division appears.
		4627	*/
		4628	__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
		4629	__isl_take isl_pw_qpolynomial *pwqp, int sign)
		4630	{
		4631	int i;
		4632	struct isl_to_poly_data data;
		4633
		4634	if (sign == 0)
		4635	return pwqp_drop_floors(pwqp);
		4636
		4637	if (!pwqp)
		4638	return NULL;
		4639
		4640	data.sign = sign;
		4641	data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
		4642
		4643	for (i = 0; i < pwqp->n; ++i) {
		4644	if (pwqp->p[i].qp->div->n_row == 0) {
		4645	isl_pw_qpolynomial *t;
		4646	t = isl_pw_qpolynomial_alloc(
		4647	isl_set_copy(pwqp->p[i].set),
		4648	isl_qpolynomial_copy(pwqp->p[i].qp));
		4649	data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
		4650	continue;
		4651	}
		4652	data.qp = pwqp->p[i].qp;
		4653	if (isl_set_foreach_orthant(pwqp->p[i].set,
		4654	&to_polynomial_on_orthant, &data) < 0)
		4655	goto error;
		4656	}
		4657
		4658	isl_pw_qpolynomial_free(pwqp);
		4659
		4660	return data.res;
		4661	error:
		4662	isl_pw_qpolynomial_free(pwqp);
		4663	isl_pw_qpolynomial_free(data.res);
		4664	return NULL;
		4665	}
		4666
		4667	static int poly_entry(void **entry, void *user)
		4668	{
		4669	int *sign = user;
		4670	isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
		4671
		4672	*pwqp = isl_pw_qpolynomial_to_polynomial(*pwqp, *sign);
		4673
		4674	return *pwqp ? 0 : -1;
		4675	}
		4676
		4677	__isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
		4678	__isl_take isl_union_pw_qpolynomial *upwqp, int sign)
		4679	{
		4680	upwqp = isl_union_pw_qpolynomial_cow(upwqp);
		4681	if (!upwqp)
		4682	return NULL;
		4683
		4684	if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
		4685	&poly_entry, &sign) < 0)
		4686	goto error;
		4687
		4688	return upwqp;
		4689	error:
		4690	isl_union_pw_qpolynomial_free(upwqp);
		4691	return NULL;
		4692	}
		4693
		4694	__isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
		4695	__isl_take isl_qpolynomial *qp)
		4696	{
		4697	int i, k;
		4698	isl_space *dim;
		4699	isl_vec *aff = NULL;
		4700	isl_basic_map *bmap = NULL;
		4701	unsigned pos;
		4702	unsigned n_div;
		4703
		4704	if (!qp)
		4705	return NULL;
		4706	if (!isl_upoly_is_affine(qp->upoly))
		4707	isl_die(qp->dim->ctx, isl_error_invalid,
		4708	"input quasi-polynomial not affine", goto error);
		4709	aff = isl_qpolynomial_extract_affine(qp);
		4710	if (!aff)
		4711	goto error;
		4712	dim = isl_qpolynomial_get_space(qp);
		4713	pos = 1 + isl_space_offset(dim, isl_dim_out);
		4714	n_div = qp->div->n_row;
		4715	bmap = isl_basic_map_alloc_space(dim, n_div, 1, 2 * n_div);
		4716
		4717	for (i = 0; i < n_div; ++i) {
		4718	k = isl_basic_map_alloc_div(bmap);
		4719	if (k < 0)
		4720	goto error;
		4721	isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
		4722	isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
		4723	if (isl_basic_map_add_div_constraints(bmap, k) < 0)
		4724	goto error;
		4725	}
		4726	k = isl_basic_map_alloc_equality(bmap);
		4727	if (k < 0)
		4728	goto error;
		4729	isl_int_neg(bmap->eq[k][pos], aff->el[0]);
		4730	isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
		4731	isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
		4732
		4733	isl_vec_free(aff);
		4734	isl_qpolynomial_free(qp);
		4735	bmap = isl_basic_map_finalize(bmap);
		4736	return bmap;
		4737	error:
		4738	isl_vec_free(aff);
		4739	isl_qpolynomial_free(qp);
		4740	isl_basic_map_free(bmap);
		4741	return NULL;
		4742	}