brintos / llvm-project-archived public Read only

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