1	/*
2	* Copyright 2010 INRIA Saclay
3	*
4	* Use of this software is governed by the MIT license
5	*
6	* Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7	* Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
8	* 91893 Orsay, France
9	*/
10
11	#include <stdlib.h>
12	#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_BASE
32	#define EL_BASE qpolynomial
33
34	#include <isl_list_templ.c>
35
36	#undef EL_BASE
37	#define EL_BASE pw_qpolynomial
38
39	#include <isl_list_templ.c>
40
41	static 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
51	isl_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(b: 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	*/
84	static 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(poly: 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(poly: poly1);
106	cst2 = isl_poly_as_cst(poly: 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(poly: poly1);
116	rec2 = isl_poly_as_rec(poly: 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(poly1: rec1->p[i], poly2: rec2->p[i]);
125	if (cmp != 0)
126	return cmp;
127	}
128
129	return 0;
130	}
131
132	isl_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(poly: 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(poly: poly1);
150	cst2 = isl_poly_as_cst(poly: 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(b: r);
156	}
157
158	rec1 = isl_poly_as_rec(poly: poly1);
159	rec2 = isl_poly_as_rec(poly: 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(poly1: rec1->p[i], poly2: rec2->p[i]);
168	if (eq < 0 || !eq)
169	return eq;
170	}
171
172	return isl_bool_true;
173	}
174
175	isl_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
191	int 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
207	isl_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
223	isl_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
239	isl_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
255	isl_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(b: r);
271	}
272
273	isl_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(space: qp->dim);
423	qp->dim = space;
424
425	return qp;
426	error:
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.c
433	* and doesn't know if the space of an element object is represented
434	* 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, space: domain);
442	}
443
444	isl_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	*/
451	static __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(space: isl_qpolynomial_peek_domain_space(qp));
463	}
464
465	#undef TYPE
466	#define TYPE isl_qpolynomial
467	#undef PEEK_SPACE
468	#define PEEK_SPACE peek_domain_space
469
470	static
471	#include "isl_type_has_equal_space_bin_templ.c"
472	static
473	#include "isl_type_check_equal_space_templ.c"
474
475	#undef PEEK_SPACE
476
477	/* Return a copy of the local space on which "qp" is defined.
478	*/
479	static __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, div: isl_mat_copy(mat: 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(space: qp->dim);
497	space = isl_space_from_domain(space);
498	space = isl_space_add_dims(space, type: isl_dim_out, n: 1);
499	return space;
500	}
501
502	/* Return the number of variables of the given type in the domain of "qp".
503	*/
504	isl_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, type: 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 or
533	* isl_dim_param.
534	*/
535	static 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, the
541	* ls field corresponds to the domain of that space.
542	*/
543	isl_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" within
555	* the variables of the domain of "qp".
556	*/
557	static 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, type: 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" in
578	* the domain of "qp".
579	*/
580	unsigned 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
595	isl_bool isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
596	{
597	return qp ? isl_poly_is_zero(poly: qp->poly) : isl_bool_error;
598	}
599
600	isl_bool isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
601	{
602	return qp ? isl_poly_is_one(poly: qp->poly) : isl_bool_error;
603	}
604
605	isl_bool isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp)
606	{
607	return qp ? isl_poly_is_nan(poly: qp->poly) : isl_bool_error;
608	}
609
610	isl_bool isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp)
611	{
612	return qp ? isl_poly_is_infty(poly: qp->poly) : isl_bool_error;
613	}
614
615	isl_bool isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp)
616	{
617	return qp ? isl_poly_is_neginfty(poly: qp->poly) : isl_bool_error;
618	}
619
620	int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp)
621	{
622	return qp ? isl_poly_sgn(poly: qp->poly) : 0;
623	}
624
625	static 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
631	static 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(poly: 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(poly: isl_poly_zero(ctx: 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(ctx: poly->ctx, var: poly->var, size: rec->n);
677	if (!dup)
678	return NULL;
679
680	for (i = 0; i < rec->n; ++i) {
681	dup->p[i] = isl_poly_copy(poly: rec->p[i]);
682	if (!dup->p[i])
683	goto error;
684	dup->n++;
685	}
686
687	return &dup->poly;
688	error:
689	isl_poly_free(poly: &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	else
703	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(cst: (isl_poly_cst *) poly);
727	else
728	poly_free_rec(rec: (isl_poly_rec *) poly);
729
730	isl_ctx_deref(ctx: poly->ctx);
731	free(ptr: poly);
732	return NULL;
733	}
734
735	static 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(poly: poly1);
755	if (!poly1 || !poly2)
756	goto error;
757
758	cst1 = isl_poly_as_cst(poly: poly1);
759	cst2 = isl_poly_as_cst(poly: 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(cst: cst1);
770
771	isl_poly_free(poly: poly2);
772	return poly1;
773	error:
774	isl_poly_free(poly: poly1);
775	isl_poly_free(poly: poly2);
776	return NULL;
777	}
778
779	static __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
790	static __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(poly: rec->p[0]);
802	isl_poly_free(poly);
803	return cst;
804	error:
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(poly: poly1);
820	if (is_nan < 0)
821	goto error;
822	if (is_nan) {
823	isl_poly_free(poly: poly2);
824	return poly1;
825	}
826
827	is_nan = isl_poly_is_nan(poly: poly2);
828	if (is_nan < 0)
829	goto error;
830	if (is_nan) {
831	isl_poly_free(poly: poly1);
832	return poly2;
833	}
834
835	is_zero = isl_poly_is_zero(poly: poly1);
836	if (is_zero < 0)
837	goto error;
838	if (is_zero) {
839	isl_poly_free(poly: poly1);
840	return poly2;
841	}
842
843	is_zero = isl_poly_is_zero(poly: poly2);
844	if (is_zero < 0)
845	goto error;
846	if (is_zero) {
847	isl_poly_free(poly: poly2);
848	return poly1;
849	}
850
851	if (poly1->var < poly2->var)
852	return isl_poly_sum(poly1: poly2, 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(poly: poly2);
859	if (is_infty >= 0 && !is_infty)
860	is_infty = isl_poly_is_neginfty(poly: poly2);
861	if (is_infty < 0)
862	goto error;
863	if (is_infty) {
864	isl_poly_free(poly: poly1);
865	return poly2;
866	}
867	poly1 = isl_poly_cow(poly: poly1);
868	rec = isl_poly_as_rec(poly: poly1);
869	if (!rec)
870	goto error;
871	rec->p[0] = isl_poly_sum(poly1: rec->p[0], poly2);
872	if (rec->n == 1)
873	poly1 = replace_by_constant_term(poly: poly1);
874	return poly1;
875	}
876
877	is_cst = isl_poly_is_cst(poly: 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(poly: poly1);
884	rec2 = isl_poly_as_rec(poly: poly2);
885	if (!rec1 || !rec2)
886	goto error;
887
888	if (rec1->n < rec2->n)
889	return isl_poly_sum(poly1: poly2, poly2: poly1);
890
891	poly1 = isl_poly_cow(poly: poly1);
892	rec1 = isl_poly_as_rec(poly: 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(poly1: rec1->p[i],
900	poly2: isl_poly_copy(poly: 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(poly: rec1->p[i]);
906	if (is_zero < 0)
907	goto error;
908	if (is_zero) {
909	isl_poly_free(poly: rec1->p[i]);
910	rec1->n--;
911	}
912	}
913
914	if (rec1->n == 0)
915	poly1 = replace_by_zero(poly: poly1);
916	else if (rec1->n == 1)
917	poly1 = replace_by_constant_term(poly: poly1);
918
919	isl_poly_free(poly: poly2);
920
921	return poly1;
922	error:
923	isl_poly_free(poly: poly1);
924	isl_poly_free(poly: 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(poly: rec->p[0], v);
961	if (!rec->p[0])
962	goto error;
963
964	return poly;
965	error:
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(poly: rec->p[i], v);
1012	if (!rec->p[i])
1013	goto error;
1014	}
1015
1016	return poly;
1017	error:
1018	isl_poly_free(poly);
1019	return NULL;
1020	}
1021
1022	/* Multiply the constant polynomial "poly" by "v".
1023	*/
1024	static __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	*/
1051	static __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(poly: rec->p[i], v);
1071	if (!rec->p[i])
1072	goto error;
1073	}
1074
1075	return poly;
1076	error:
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(poly: poly1);
1088	if (!poly1 || !poly2)
1089	goto error;
1090
1091	cst1 = isl_poly_as_cst(poly: poly1);
1092	cst2 = isl_poly_as_cst(poly: 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(cst: cst1);
1098
1099	isl_poly_free(poly: poly2);
1100	return poly1;
1101	error:
1102	isl_poly_free(poly: poly1);
1103	isl_poly_free(poly: 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(poly: poly1);
1117	rec2 = isl_poly_as_rec(poly: poly2);
1118	if (!rec1 || !rec2)
1119	goto error;
1120	size = rec1->n + rec2->n - 1;
1121	res = isl_poly_alloc_rec(ctx: poly1->ctx, var: 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(poly1: isl_poly_copy(poly: rec2->p[0]),
1127	poly2: isl_poly_copy(poly: 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(ctx: 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(poly1: isl_poly_copy(poly: rec2->p[j]),
1142	poly2: isl_poly_copy(poly: rec1->p[i]));
1143	res->p[i + j] = isl_poly_sum(poly1: res->p[i + j], poly2: poly);
1144	if (!res->p[i + j])
1145	goto error;
1146	}
1147	}
1148
1149	isl_poly_free(poly: poly1);
1150	isl_poly_free(poly: poly2);
1151
1152	return &res->poly;
1153	error:
1154	isl_poly_free(poly: poly1);
1155	isl_poly_free(poly: poly2);
1156	isl_poly_free(poly: &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(poly: poly1);
1169	if (is_nan < 0)
1170	goto error;
1171	if (is_nan) {
1172	isl_poly_free(poly: poly2);
1173	return poly1;
1174	}
1175
1176	is_nan = isl_poly_is_nan(poly: poly2);
1177	if (is_nan < 0)
1178	goto error;
1179	if (is_nan) {
1180	isl_poly_free(poly: poly1);
1181	return poly2;
1182	}
1183
1184	is_zero = isl_poly_is_zero(poly: poly1);
1185	if (is_zero < 0)
1186	goto error;
1187	if (is_zero) {
1188	isl_poly_free(poly: poly2);
1189	return poly1;
1190	}
1191
1192	is_zero = isl_poly_is_zero(poly: poly2);
1193	if (is_zero < 0)
1194	goto error;
1195	if (is_zero) {
1196	isl_poly_free(poly: poly1);
1197	return poly2;
1198	}
1199
1200	is_one = isl_poly_is_one(poly: poly1);
1201	if (is_one < 0)
1202	goto error;
1203	if (is_one) {
1204	isl_poly_free(poly: poly1);
1205	return poly2;
1206	}
1207
1208	is_one = isl_poly_is_one(poly: poly2);
1209	if (is_one < 0)
1210	goto error;
1211	if (is_one) {
1212	isl_poly_free(poly: poly2);
1213	return poly1;
1214	}
1215
1216	if (poly1->var < poly2->var)
1217	return isl_poly_mul(poly1: poly2, 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(poly: poly2);
1225	if (is_infty >= 0 && !is_infty)
1226	is_infty = isl_poly_is_neginfty(poly: poly2);
1227	if (is_infty < 0)
1228	goto error;
1229	if (is_infty) {
1230	isl_ctx *ctx = poly1->ctx;
1231	isl_poly_free(poly: poly1);
1232	isl_poly_free(poly: poly2);
1233	return isl_poly_nan(ctx);
1234	}
1235	poly1 = isl_poly_cow(poly: poly1);
1236	rec = isl_poly_as_rec(poly: poly1);
1237	if (!rec)
1238	goto error;
1239
1240	for (i = 0; i < rec->n; ++i) {
1241	rec->p[i] = isl_poly_mul(poly1: rec->p[i],
1242	poly2: isl_poly_copy(poly: poly2));
1243	if (!rec->p[i])
1244	goto error;
1245	}
1246	isl_poly_free(poly: poly2);
1247	return poly1;
1248	}
1249
1250	is_cst = isl_poly_is_cst(poly: 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);
1257	error:
1258	isl_poly_free(poly: poly1);
1259	isl_poly_free(poly: 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	else
1275	res = isl_poly_one(ctx: poly->ctx);
1276
1277	while (power >>= 1) {
1278	poly = isl_poly_mul(poly1: poly, poly2: isl_poly_copy(poly));
1279	if (power % 2)
1280	res = isl_poly_mul(poly1: res, poly2: 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, type: 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(ctx: space->ctx, n_row: n_div, n_col: 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;
1314	error:
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(space: isl_space_copy(space: qp->dim), n_div: qp->div->n_row,
1338	poly: isl_poly_copy(poly: qp->poly));
1339	if (!dup)
1340	return NULL;
1341	isl_mat_free(mat: dup->div);
1342	dup->div = isl_mat_copy(mat: qp->div);
1343	if (!dup->div)
1344	goto error;
1345
1346	return dup;
1347	error:
1348	isl_qpolynomial_free(qp: 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(space: qp->dim);
1373	isl_mat_free(mat: qp->div);
1374	isl_poly_free(poly: qp->poly);
1375
1376	free(ptr: 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, var: pos, size: 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(poly: rec->p[power]);
1396	isl_int_set_si(cst->n, 1);
1397
1398	return &rec->poly;
1399	error:
1400	isl_poly_free(poly: &rec->poly);
1401	return NULL;
1402	}
1403
1404	/* r array maps original positions to new positions.
1405	*/
1406	static __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(ctx: poly->ctx, pos: r[poly->var], power: 1);
1427	res = reorder(poly: isl_poly_copy(poly: rec->p[rec->n - 1]), r);
1428
1429	for (i = rec->n - 2; i >= 0; --i) {
1430	res = isl_poly_mul(poly1: res, poly2: isl_poly_copy(poly: base));
1431	res = isl_poly_sum(poly1: res, poly2: reorder(poly: isl_poly_copy(poly: rec->p[i]), r));
1432	}
1433
1434	isl_poly_free(poly: base);
1435	isl_poly_free(poly);
1436
1437	return res;
1438	error:
1439	isl_poly_free(poly);
1440	return NULL;
1441	}
1442
1443	static 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(mat1: div1, mat2: 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(mat1: div1, mat2: div2);
1462
1463	div1->n_row = n_row;
1464	div1->n_col = n_col;
1465
1466	return equal;
1467	}
1468
1469	static 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(p: div->row[i], len: div->n_col);
1474	lj = isl_seq_last_non_zero(p: div->row[j], len: div->n_col);
1475
1476	if (li != lj)
1477	return li - lj;
1478
1479	return isl_seq_cmp(p1: div->row[i], p2: div->row[j], len: div->n_col);
1480	}
1481
1482	struct isl_div_sort_info {
1483	isl_mat *div;
1484	int row;
1485	};
1486
1487	static 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(div: i1->div, i: i1->row, j: i2->row);
1494	}
1495
1496	/* Sort divs and remove duplicates.
1497	*/
1498	static __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, type: 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(base: array, nmemb: qp->div->n_row, size: sizeof(struct isl_div_sort_info),
1534	compar: 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(mat: qp->div, i, j: 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(p1: qp->div->row[i - skip - 1],
1553	p2: qp->div->row[i - skip], len: qp->div->n_col)) {
1554	qp->div = isl_mat_drop_rows(mat: qp->div, row: i - skip, n: 1);
1555	isl_mat_col_add(mat: qp->div, dst_col: 2 + div_pos + i - skip - 1,
1556	src_col: 2 + div_pos + i - skip);
1557	qp->div = isl_mat_drop_cols(mat: qp->div,
1558	col: 2 + div_pos + i - skip, n: 1);
1559	skip++;
1560	}
1561	reordering[div_pos + array[i].row] = div_pos + i - skip;
1562	}
1563
1564	qp->poly = reorder(poly: qp->poly, r: reordering);
1565
1566	if (!qp->poly || !qp->div)
1567	goto error;
1568
1569	free(ptr: at);
1570	free(ptr: pos);
1571	free(ptr: array);
1572	free(ptr: reordering);
1573
1574	return qp;
1575	error:
1576	free(ptr: at);
1577	free(ptr: pos);
1578	free(ptr: array);
1579	free(ptr: reordering);
1580	isl_qpolynomial_free(qp);
1581	return NULL;
1582	}
1583
1584	static __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(poly: rec->p[i], exp, first);
1615	if (!rec->p[i])
1616	goto error;
1617	}
1618
1619	return poly;
1620	error:
1621	isl_poly_free(poly);
1622	return NULL;
1623	}
1624
1625	static __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(qp: qp1);
1636	qp2 = isl_qpolynomial_cow(qp: 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(div1: qp1->div, div2: qp2->div, exp1, exp2);
1652	if (!div)
1653	goto error;
1654
1655	isl_mat_free(mat: qp1->div);
1656	qp1->div = isl_mat_copy(mat: div);
1657	isl_mat_free(mat: qp2->div);
1658	qp2->div = isl_mat_copy(mat: div);
1659
1660	qp1->poly = expand(poly: qp1->poly, exp: exp1, first: div->n_col - div->n_row - 2);
1661	qp2->poly = expand(poly: qp2->poly, exp: exp2, first: div->n_col - div->n_row - 2);
1662
1663	if (!qp1->poly || !qp2->poly)
1664	goto error;
1665
1666	isl_mat_free(mat: div);
1667	free(ptr: exp1);
1668	free(ptr: exp2);
1669
1670	return fn(qp1, qp2);
1671	error:
1672	isl_mat_free(mat: div);
1673	free(ptr: exp1);
1674	free(ptr: exp2);
1675	isl_qpolynomial_free(qp: qp1);
1676	isl_qpolynomial_free(qp: 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(qp: qp1);
1686
1687	if (isl_qpolynomial_check_equal_space(obj1: qp1, obj2: qp2) < 0)
1688	goto error;
1689
1690	if (qp1->div->n_row < qp2->div->n_row)
1691	return isl_qpolynomial_add(qp1: qp2, qp2: qp1);
1692
1693	compatible = compatible_divs(div1: qp1->div, div2: qp2->div);
1694	if (compatible < 0)
1695	goto error;
1696	if (!compatible)
1697	return with_merged_divs(fn: isl_qpolynomial_add, qp1, qp2);
1698
1699	qp1->poly = isl_poly_sum(poly1: qp1->poly, poly2: isl_poly_copy(poly: qp2->poly));
1700	if (!qp1->poly)
1701	goto error;
1702
1703	isl_qpolynomial_free(qp: qp2);
1704
1705	return qp1;
1706	error:
1707	isl_qpolynomial_free(qp: qp1);
1708	isl_qpolynomial_free(qp: 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(qp: qp1, context: isl_set_copy(set: 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, qp2: isl_qpolynomial_neg(qp: 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(poly: qp->poly, v);
1739	if (!qp->poly)
1740	goto error;
1741
1742	return qp;
1743	error:
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, v: 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(domain: isl_space_copy(space: 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(poly: qp->poly, v);
1775	if (!qp->poly)
1776	goto error;
1777
1778	return qp;
1779	error:
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(domain: space);
1814	}
1815
1816	qp = isl_qpolynomial_cow(qp);
1817	if (!qp)
1818	goto error;
1819
1820	qp->poly = isl_poly_scale_val(poly: qp->poly, v);
1821	if (!qp->poly)
1822	qp = isl_qpolynomial_free(qp);
1823
1824	isl_val_free(v);
1825	return qp;
1826	error:
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, v: isl_val_inv(v));
1848	error:
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(qp: qp1);
1860
1861	if (isl_qpolynomial_check_equal_space(obj1: qp1, obj2: qp2) < 0)
1862	goto error;
1863
1864	if (qp1->div->n_row < qp2->div->n_row)
1865	return isl_qpolynomial_mul(qp1: qp2, qp2: qp1);
1866
1867	compatible = compatible_divs(div1: qp1->div, div2: qp2->div);
1868	if (compatible < 0)
1869	goto error;
1870	if (!compatible)
1871	return with_merged_divs(fn: isl_qpolynomial_mul, qp1, qp2);
1872
1873	qp1->poly = isl_poly_mul(poly1: qp1->poly, poly2: isl_poly_copy(poly: qp2->poly));
1874	if (!qp1->poly)
1875	goto error;
1876
1877	isl_qpolynomial_free(qp: qp2);
1878
1879	return qp1;
1880	error:
1881	isl_qpolynomial_free(qp: qp1);
1882	isl_qpolynomial_free(qp: 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(poly: qp->poly, power);
1895	if (!qp->poly)
1896	goto error;
1897
1898	return qp;
1899	error:
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(qp: 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(space: domain, n_div: 0, poly: isl_poly_zero(ctx: 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(space: domain, n_div: 0, poly: isl_poly_one(ctx: 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(space: domain, n_div: 0, poly: isl_poly_infty(ctx: 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(space: domain, n_div: 0, poly: isl_poly_neginfty(ctx: 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(space: domain, n_div: 0, poly: isl_poly_nan(ctx: 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(poly: qp->poly);
1977	isl_int_set(cst->n, v);
1978
1979	return qp;
1980	}
1981
1982	isl_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(poly: qp->poly);
1992	if (is_cst < 0 || !is_cst)
1993	return is_cst;
1994
1995	cst = isl_poly_as_cst(poly: 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	*/
2009	static __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(ctx: cst->poly.ctx, n: cst->n, d: 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(poly: qp->poly);
2043	}
2044
2045	isl_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(poly: rec->p[1]);
2066	if (is_cst < 0 || !is_cst)
2067	return is_cst;
2068
2069	return isl_poly_is_affine(poly: rec->p[0]);
2070	}
2071
2072	isl_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(poly: qp->poly);
2081	}
2082
2083	static 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(dst: aff->el, src: aff->el, f, len: 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
2103	int 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, pos: 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(poly: rec->p[1]);
2127	if (!cst)
2128	return -1;
2129	update_coeff(aff, cst, pos: 1 + poly->var);
2130
2131	return isl_poly_update_affine(poly: 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, type: isl_dim_all);
2141	if (d < 0)
2142	return NULL;
2143
2144	aff = isl_vec_alloc(ctx: qp->div->ctx, size: 2 + d);
2145	if (!aff)
2146	return NULL;
2147
2148	isl_seq_clr(p: aff->el + 1, len: 1 + d);
2149	isl_int_set_si(aff->el[0], 1);
2150
2151	if (isl_poly_update_affine(poly: qp->poly, aff) < 0)
2152	goto error;
2153
2154	return aff;
2155	error:
2156	isl_vec_free(vec: 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	*/
2165	int 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(space1: qp1->dim, space2: qp2->dim);
2178	if (cmp != 0)
2179	return cmp;
2180
2181	cmp = isl_local_cmp(local1: qp1->div, local2: qp2->div);
2182	if (cmp != 0)
2183	return cmp;
2184
2185	return isl_poly_plain_cmp(poly1: qp1->poly, poly2: 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	*/
2192	isl_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(qp: qp1) || isl_qpolynomial_is_nan(qp: qp2))
2201	return isl_bool_false;
2202
2203	equal = isl_space_is_equal(space1: qp1->dim, space2: qp2->dim);
2204	if (equal < 0 || !equal)
2205	return equal;
2206
2207	equal = isl_mat_is_equal(mat1: qp1->div, mat2: qp2->div);
2208	if (equal < 0 || !equal)
2209	return equal;
2210
2211	return isl_poly_is_equal(poly1: qp1->poly, poly2: qp2->poly);
2212	}
2213
2214	static 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(poly: 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(ctx: isl_qpolynomial_get_ctx(qp));
2249	if (!d)
2250	return NULL;
2251	if (poly_update_den(poly: qp->poly, d: &d->n) < 0)
2252	return isl_val_free(v: 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(space: domain, n_div: 0,
2267	poly: 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(space: domain ) < 0)
2274	goto error;
2275	if (isl_space_check_range(space: domain, type, first: pos, n: 1) < 0)
2276	goto error;
2277
2278	pos += isl_space_offset(space: domain, type);
2279
2280	return isl_qpolynomial_var_pow_on_domain(domain, pos, power: 1);
2281	error:
2282	isl_space_free(space: 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(ctx: poly->ctx, pos: poly->var, power: 1);
2311	else
2312	base = isl_poly_copy(poly: subs[poly->var - first]);
2313
2314	res = isl_poly_subs(poly: isl_poly_copy(poly: 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(poly: isl_poly_copy(poly: rec->p[i]), first, n, subs);
2318	res = isl_poly_mul(poly1: res, poly2: isl_poly_copy(poly: base));
2319	res = isl_poly_sum(poly1: res, poly2: t);
2320	}
2321
2322	isl_poly_free(poly: base);
2323	isl_poly_free(poly);
2324
2325	return res;
2326	error:
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, n: f[0], d: 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, n: f[1 + i], d: denom);
2348	t = isl_poly_var_pow(ctx, pos: i, power: 1);
2349	t = isl_poly_mul(poly1: c, poly2: t);
2350	poly = isl_poly_sum(poly1: poly, poly2: t);
2351	}
2352
2353	return poly;
2354	}
2355
2356	/* Remove common factor of non-constant terms and denominator.
2357	*/
2358	static 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(p: qp->div->row[div] + 2, len: total, gcd: &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(dst: qp->div->row[div] + 2, src: qp->div->row[div] + 2,
2370	f: ctx->normalize_gcd, len: 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 definition
2379	* of any other integer divisions.
2380	*/
2381	static __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, type: isl_dim_div);
2397	if (div_pos < 0)
2398	goto error;
2399	qp->poly = isl_poly_subs(poly: qp->poly, first: div_pos + div, n: 1, subs: &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(mat: qp->div, row: div, n: 1);
2412	qp->div = isl_mat_drop_cols(mat: qp->div, col: 2 + div_pos + div, n: 1);
2413	qp->poly = reorder(poly: qp->poly, r: reordering);
2414	free(ptr: reordering);
2415
2416	if (!qp->poly || !qp->div)
2417	goto error;
2418
2419	isl_poly_free(poly: s);
2420	return qp;
2421	error:
2422	isl_qpolynomial_free(qp);
2423	isl_poly_free(poly: s);
2424	return NULL;
2425	}
2426
2427	/* Replace all integer divisions [e/d] that turn out to not actually be integer
2428	* divisions because d is equal to 1 by their definition, i.e., e.
2429	*/
2430	static __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, type: 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(dst: qp->div->row[j] + 1,
2448	m1: qp->div->ctx->one, src1: qp->div->row[j] + 1,
2449	m2: qp->div->row[j][2 + div_pos + i],
2450	src2: qp->div->row[i] + 1, len: 1 + div_pos + i);
2451	isl_int_set_si(qp->div->row[j][2 + div_pos + i], 0);
2452	normalize_div(qp, div: j);
2453	}
2454	s = isl_poly_from_affine(ctx: qp->dim->ctx, f: qp->div->row[i] + 1,
2455	denom: qp->div->row[i][0], len: qp->div->n_col - 1);
2456	qp = substitute_div(qp, div: 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 by
2465	* d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
2466	* 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 need
2468	* to adjust the coefficient of x in each later div that depends on the
2469	* current div "div" and also in the affine expressions in the rows of "mat"
2470	* (if they too depend on "div").
2471	*/
2472	static 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: *mat, dst_col: i, f: v, src_col: 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 the
2499	* denominator. If so, we will invert the div to further reduce the number
2500	* 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 bigger
2503	* than half the denominator.
2504	*/
2505	static 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 in
2527	* later divs and in the rows of "mat"). After calling this function, the
2528	* coefficients of e should be reduced again.
2529	*/
2530	static 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(dst: qp->div->row[div] + 1,
2536	src: qp->div->row[div] + 1, len: 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: *mat, col: 1 + total + div);
2541	isl_mat_col_mul(mat: qp->div, dst_col: 2 + total + div,
2542	f: qp->div->ctx->negone, src_col: 2 + total + div);
2543	}
2544
2545	/* Reduce all divs of "qp" to have coefficients
2546	* in the interval [0, d-1], with d the denominator and such that the
2547	* 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 reflected
2549	* in the factors of the polynomial that refer to the original
2550	* integer divisions. To this end, the modifications are collected
2551	* 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 functions
2555	* eliminate divs or merge two or more divs into one, the coefficients
2556	* of the enclosing divs may have to be reduced again, so we call
2557	* ourselves recursively if the number of divs decreases.
2558	*/
2559	static __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, type: isl_dim_all);
2569	n_div = isl_qpolynomial_domain_dim(qp, type: isl_dim_div);
2570	o_div = isl_qpolynomial_domain_offset(qp, type: 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_row: n_div, n_col: 1 + total);
2575
2576	for (i = 0; i < n_div; ++i)
2577	mat = isl_mat_set_element_si(mat, row: i, col: o_div + i, v: 1);
2578
2579	for (i = 0; i < qp->div->n_row; ++i) {
2580	normalize_div(qp, div: i);
2581	reduce_div(qp, div: i, mat: &mat);
2582	if (needs_invert(div: qp->div, row: i)) {
2583	invert_div(qp, div: i, mat: &mat);
2584	reduce_div(qp, div: i, mat: &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, f: mat->row[i], denom: ctx->one,
2595	len: 1 + total);
2596	qp->poly = isl_poly_subs(poly: qp->poly, first: o_div - 1, n: n_div, subs: s);
2597	for (i = 0; i < n_div; ++i)
2598	isl_poly_free(poly: s[i]);
2599	free(ptr: 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, type: 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;
2614	error:
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(poly: 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(poly: qp->poly);
2650	isl_int_set(cst->n, val->n);
2651	isl_int_set(cst->d, val->d);
2652
2653	isl_val_free(v: val);
2654	return qp;
2655	error:
2656	isl_val_free(v: val);
2657	isl_qpolynomial_free(qp);
2658	return NULL;
2659	}
2660
2661	static 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(poly: rec->p[i], active, d) < 0)
2679	return isl_stat_error;
2680
2681	return isl_stat_ok;
2682	}
2683
2684	static 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, type: 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(poly: qp->poly, active, d);
2704	}
2705
2706	#undef TYPE
2707	#define TYPE isl_qpolynomial
2708	static
2709	#include "check_type_range_templ.c"
2710
2711	isl_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(obj: 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, type: 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, type: 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(ptr: active);
2750
2751	return involves;
2752	error:
2753	free(ptr: active);
2754	return isl_bool_error;
2755	}
2756
2757	/* Remove divs that do not appear in the quasi-polynomial, nor in any
2758	* of the divs that do appear in the quasi-polynomial.
2759	*/
2760	static __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, type: 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(poly: qp->poly, active, d: 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(ptr: 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(mat: qp->div, row: i - skip, n: 1);
2820	qp->div = isl_mat_drop_cols(mat: qp->div,
2821	col: 2 + div_pos + i - skip, n: 1);
2822	skip++;
2823	}
2824	reordering[div_pos + i] = div_pos + i - skip;
2825	}
2826
2827	qp->poly = reorder(poly: qp->poly, r: reordering);
2828
2829	if (!qp->poly || !qp->div)
2830	goto error;
2831
2832	free(ptr: active);
2833	free(ptr: reordering);
2834
2835	return qp;
2836	error:
2837	free(ptr: active);
2838	free(ptr: 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(poly: rec->p[i], first, n);
2867	if (!rec->p[i])
2868	goto error;
2869	}
2870
2871	return poly;
2872	error:
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(space: qp->dim, type, pos, name: s);
2890	if (!qp->dim)
2891	goto error;
2892	return qp;
2893	error:
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(obj: 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(space: 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(space: qp->dim, type, first, num: 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(mat: qp->div, col: 2 + first, n);
2933	if (!qp->div)
2934	goto error;
2935
2936	qp->poly = isl_poly_drop(poly: qp->poly, first, n);
2937	if (!qp->poly)
2938	goto error;
2939
2940	return qp;
2941	error:
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, type: isl_dim_in);
2957	if (n < 0)
2958	return isl_qpolynomial_free(qp);
2959	involves = isl_qpolynomial_involves_dims(qp, type: isl_dim_in, first: 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, type: isl_dim_in, first: 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
2973	static __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(bset: eq);
2986	return qp;
2987	}
2988
2989	qp = isl_qpolynomial_cow(qp);
2990	if (!qp)
2991	goto error;
2992	qp->div = isl_mat_cow(mat: qp->div);
2993	if (!qp->div)
2994	goto error;
2995
2996	total = isl_basic_set_offset(bset: eq, type: 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(p: eq->eq[i], len: 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(dst: qp->div->row[k] + 1, src: eq->eq[i], pos: j, len: total,
3008	m: &qp->div->row[k][0]);
3009	normalize_div(qp, div: k);
3010	}
3011
3012	if (isl_int_is_pos(eq->eq[i][j]))
3013	isl_seq_neg(dst: eq->eq[i], src: eq->eq[i], len: 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(ctx: qp->dim->ctx,
3018	f: eq->eq[i], denom, len: total);
3019	qp->poly = isl_poly_subs(poly: qp->poly, first: j - 1, n: 1, subs: &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(bset: eq);
3028
3029	qp = substitute_non_divs(qp);
3030	qp = sort_divs(qp);
3031
3032	return qp;
3033	error:
3034	isl_basic_set_free(bset: 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(bset: eq, type: isl_dim_set, n: qp->div->n_row);
3048	return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
3049	error:
3050	isl_basic_set_free(bset: eq);
3051	isl_qpolynomial_free(qp);
3052	return NULL;
3053	}
3054
3055	/* Look for equalities among the variables shared by context and qp
3056	* and the integer divisions of qp, if any.
3057	* The equalities are then used to eliminate variables and/or integer
3058	* 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, set: context);
3068
3069	aff = isl_set_affine_hull(set: context);
3070	return isl_qpolynomial_substitute_equalities_lifted(qp, eq: 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(set: dom_context, params: context);
3079	return isl_qpolynomial_gist(qp, context: 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 uniform
3085	* interface over all piecewise types.
3086	*/
3087	static __isl_give isl_qpolynomial *isl_qpolynomial_zero_in_space(
3088	__isl_take isl_space *space)
3089	{
3090	return isl_qpolynomial_zero_on_domain(domain: isl_space_domain(space));
3091	}
3092
3093	#define isl_qpolynomial_involves_nan isl_qpolynomial_is_nan
3094
3095	#undef PW
3096	#define PW isl_pw_qpolynomial
3097	#undef BASE
3098	#define BASE qpolynomial
3099	#undef EL_IS_ZERO
3100	#define EL_IS_ZERO is_zero
3101	#undef ZERO
3102	#define ZERO zero
3103	#undef IS_ZERO
3104	#define IS_ZERO is_zero
3105	#undef FIELD
3106	#define FIELD qp
3107	#undef DEFAULT_IS_ZERO
3108	#define DEFAULT_IS_ZERO 1
3109
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 BASE
3126	#define BASE pw_qpolynomial
3127
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
3133	int 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(set: pwqp->p[0].set))
3142	return 0;
3143
3144	return isl_qpolynomial_is_one(qp: 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_(pw1: pwqp1, pw2: 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(pw: pwqp1)) {
3168	isl_pw_qpolynomial_free(pw: pwqp2);
3169	return pwqp1;
3170	}
3171
3172	if (isl_pw_qpolynomial_is_zero(pw: pwqp2)) {
3173	isl_pw_qpolynomial_free(pw: pwqp1);
3174	return pwqp2;
3175	}
3176
3177	if (isl_pw_qpolynomial_is_one(pwqp: pwqp1)) {
3178	isl_pw_qpolynomial_free(pw: pwqp1);
3179	return pwqp2;
3180	}
3181
3182	if (isl_pw_qpolynomial_is_one(pwqp: pwqp2)) {
3183	isl_pw_qpolynomial_free(pw: pwqp2);
3184	return pwqp1;
3185	}
3186
3187	n = pwqp1->n * pwqp2->n;
3188	res = isl_pw_qpolynomial_alloc_size(space: isl_space_copy(space: 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(set1: isl_set_copy(set: pwqp1->p[i].set),
3195	set2: isl_set_copy(set: pwqp2->p[j].set));
3196	if (isl_set_plain_is_empty(set: common)) {
3197	isl_set_free(set: common);
3198	continue;
3199	}
3200
3201	prod = isl_qpolynomial_mul(
3202	qp1: isl_qpolynomial_copy(qp: pwqp1->p[i].qp),
3203	qp2: isl_qpolynomial_copy(qp: pwqp2->p[j].qp));
3204
3205	res = isl_pw_qpolynomial_add_piece(pw: res, set: common, el: prod);
3206	}
3207	}
3208
3209	isl_pw_qpolynomial_free(pw: pwqp1);
3210	isl_pw_qpolynomial_free(pw: pwqp2);
3211
3212	return res;
3213	error:
3214	isl_pw_qpolynomial_free(pw: pwqp1);
3215	isl_pw_qpolynomial_free(pw: 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(ctx: poly->ctx,
3245	n: vec->el[1 + poly->var], d: vec->el[0]);
3246
3247	res = isl_poly_eval(poly: isl_poly_copy(poly: rec->p[rec->n - 1]),
3248	vec: isl_vec_copy(vec));
3249
3250	for (i = rec->n - 2; i >= 0; --i) {
3251	res = isl_val_mul(v1: res, v2: isl_val_copy(v: base));
3252	res = isl_val_add(v1: res, v2: isl_poly_eval(poly: isl_poly_copy(poly: rec->p[i]),
3253	vec: isl_vec_copy(vec)));
3254	}
3255
3256	isl_val_free(v: base);
3257	isl_poly_free(poly);
3258	isl_vec_free(vec);
3259	return res;
3260	error:
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	*/
3269	static __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(local: qp->div, v: isl_vec_copy(vec: pnt->vec));
3297
3298	v = isl_poly_eval(poly: isl_poly_copy(poly: qp->poly), vec: ext);
3299
3300	isl_qpolynomial_free(qp);
3301	isl_point_free(pnt);
3302
3303	return v;
3304	error:
3305	isl_qpolynomial_free(qp);
3306	isl_point_free(pnt);
3307	return NULL;
3308	}
3309
3310	int 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(obj: qp, type, first, n: 0) < 0)
3337	return isl_qpolynomial_free(qp);
3338	type = domain_type(type);
3339	if (n == 0 && !isl_space_is_named_or_nested(space: qp->dim, type))
3340	return qp;
3341
3342	qp = isl_qpolynomial_cow(qp);
3343	if (!qp)
3344	return NULL;
3345
3346	g_pos = pos(space: qp->dim, type) + first;
3347
3348	qp->div = isl_mat_insert_zero_cols(mat: qp->div, first: 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(poly: qp->poly, exp, first: g_pos);
3361	free(ptr: exp);
3362	if (!qp->poly)
3363	goto error;
3364	}
3365
3366	qp->dim = isl_space_insert_dims(space: qp->dim, type, pos: first, n);
3367	if (!qp->dim)
3368	goto error;
3369
3370	return qp;
3371	error:
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, first: pos, n);
3386	}
3387
3388	static 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(obj: qp, type: src_type, first: 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(space: qp->dim, type: src_type) &&
3446	!isl_space_is_named_or_nested(space: qp->dim, type: dst_type))
3447	return qp;
3448
3449	qp = isl_qpolynomial_cow(qp);
3450	if (!qp)
3451	return NULL;
3452
3453	g_dst_pos = pos(space: qp->dim, type: dst_type) + dst_pos;
3454	g_src_pos = pos(space: qp->dim, type: src_type) + src_pos;
3455	if (dst_type > src_type)
3456	g_dst_pos -= n;
3457
3458	qp->div = isl_mat_move_cols(mat: qp->div, dst_col: 2 + g_dst_pos, src_col: 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(ctx: qp->dim->ctx,
3466	len: qp->div->n_col - 2, dst: g_dst_pos, src: g_src_pos, n);
3467	if (!reordering)
3468	goto error;
3469
3470	qp->poly = reorder(poly: qp->poly, r: reordering);
3471	free(ptr: reordering);
3472	if (!qp->poly)
3473	goto error;
3474
3475	qp->dim = isl_space_move_dims(space: qp->dim, dst_type, dst_pos, src_type, src_pos, n);
3476	if (!qp->dim)
3477	goto error;
3478
3479	return qp;
3480	error:
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, type: isl_dim_all);
3496	poly = d < 0 ? NULL : isl_poly_from_affine(ctx: space->ctx, f, denom, len: 1 + d);
3497
3498	return isl_qpolynomial_alloc(space, n_div: 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, f: aff->v->el + 1, denom: aff->v->el[0],
3512	len: aff->v->size - 1);
3513
3514	qp = isl_qpolynomial_alloc(space: isl_aff_get_domain_space(aff),
3515	n_div: aff->ls->div->n_row, poly);
3516	if (!qp)
3517	goto error;
3518
3519	isl_mat_free(mat: qp->div);
3520	qp->div = isl_mat_copy(mat: aff->ls->div);
3521	qp->div = isl_mat_cow(mat: 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;
3529	error:
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(space: isl_pw_aff_get_space(pwaff),
3544	n: 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(set: pwaff->p[i].set);
3551	qp = isl_qpolynomial_from_aff(aff: isl_aff_copy(aff: pwaff->p[i].aff));
3552	pwqp = isl_pw_qpolynomial_add_piece(pw: pwqp, set: dom, el: 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(constraint: 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(obj: 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(obj1: qp, obj2: 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(space: 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(poly: qp->poly, first, n, subs: polys);
3616
3617	free(ptr: polys);
3618
3619	if (!qp->poly)
3620	goto error;
3621
3622	return qp;
3623	error:
3624	isl_qpolynomial_free(qp);
3625	return NULL;
3626	}
3627
3628	/* Extend "bset" with extra set dimensions for each integer division
3629	* in "qp" and then call "fn" with the extended bset and the polynomial
3630	* that results from replacing each of the integer divisions by the
3631	* corresponding extra set dimension.
3632	*/
3633	isl_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(space: qp->dim);
3649	space = isl_space_add_dims(space, type: isl_dim_set, n: qp->div->n_row);
3650	poly = isl_qpolynomial_alloc(space, n_div: 0, poly: isl_poly_copy(poly: 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	*/
3660	int 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(poly: rec->p[i]);
3686	if (is_zero < 0)
3687	return -2;
3688	if (is_zero)
3689	continue;
3690	d = isl_poly_degree(poly: 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	*/
3702	int 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(space: poly->dim, type: isl_dim_set);
3711	nvar = isl_space_dim(space: poly->dim, type: isl_dim_set);
3712	if (nvar < 0)
3713	return -2;
3714	return isl_poly_degree(poly: poly->poly, first: ovar, last: 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	else
3731	return isl_poly_zero(ctx: 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(poly: rec->p[deg]);
3741	else
3742	return isl_poly_zero(ctx: 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(poly: rec->p[i], pos, deg);
3754	if (!t)
3755	goto error;
3756	isl_poly_free(poly: rec->p[i]);
3757	rec->p[i] = t;
3758	}
3759
3760	return poly;
3761	error:
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(obj: qp, type, first: t_pos, n: 1) < 0)
3784	return NULL;
3785	type = domain_type(type);
3786
3787	g_pos = pos(space: qp->dim, type) + t_pos;
3788	poly = isl_poly_coeff(poly: qp->poly, pos: g_pos, deg);
3789
3790	c = isl_qpolynomial_alloc(space: isl_space_copy(space: qp->dim),
3791	n_div: qp->div->n_row, poly);
3792	if (!c)
3793	return NULL;
3794	isl_mat_free(mat: c->div);
3795	c->div = isl_mat_copy(mat: qp->div);
3796	if (!c->div)
3797	goto error;
3798	return c;
3799	error:
3800	isl_qpolynomial_free(qp: c);
3801	return NULL;
3802	}
3803
3804	/* Homogenize the polynomial in the variables first (inclusive) up to
3805	* 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(ctx: poly->ctx, pos: first, power: target - deg);
3829	if (!hom)
3830	goto error;
3831	rec = isl_poly_as_rec(poly: hom);
3832	rec->p[target - deg] = isl_poly_mul(poly1: rec->p[target - deg], poly2: 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(poly: 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(poly: rec->p[i],
3849	deg: poly->var < last ? deg + i : i, target,
3850	first, last);
3851	if (!rec->p[i])
3852	goto error;
3853	}
3854
3855	return poly;
3856	error:
3857	isl_poly_free(poly);
3858	return NULL;
3859	}
3860
3861	/* Homogenize the polynomial in the set variables by introducing
3862	* 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(qp: poly, type: isl_dim_in, first: 0, n: 1);
3875	poly = isl_qpolynomial_cow(qp: poly);
3876	if (!poly)
3877	goto error;
3878
3879	ovar = isl_space_offset(space: poly->dim, type: isl_dim_set);
3880	nvar = isl_space_dim(space: poly->dim, type: isl_dim_set);
3881	if (nvar < 0)
3882	return isl_qpolynomial_free(qp: poly);
3883	poly->poly = isl_poly_homogenize(poly: poly->poly, deg: 0, target: deg, first: ovar, last: ovar + nvar);
3884	if (!poly->poly)
3885	goto error;
3886
3887	return poly;
3888	error:
3889	isl_qpolynomial_free(qp: 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, type: 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;
3918	error:
3919	isl_space_free(space);
3920	isl_mat_free(mat: 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, type: isl_dim_all);
3940	if (total < 0)
3941	return NULL;
3942
3943	dup = isl_term_alloc(space: isl_space_copy(space: term->dim), div: isl_mat_copy(mat: 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(space: term->dim);
3976	isl_mat_free(mat: term->div);
3977	isl_int_clear(term->n);
3978	isl_int_clear(term->d);
3979	free(ptr: term);
3980
3981	return NULL;
3982	}
3983
3984	isl_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(space: term->dim, type);
3995	case isl_dim_div: return term->div->n_row;
3996	case isl_dim_all: dim = isl_space_dim(space: term->dim, type: 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	*/
4006	static __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" within
4012	* the variables of "term".
4013	*/
4014	static 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, type: 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
4033	isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
4034	{
4035	return term ? term->dim->ctx : NULL;
4036	}
4037
4038	void 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(ctx: isl_term_get_ctx(term),
4053	n: term->n, d: term->d);
4054	}
4055
4056	#undef TYPE
4057	#define TYPE isl_term
4058	static
4059	#include "check_type_range_templ.c"
4060
4061	isl_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(obj: term, type, first: pos, n: 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(obj: term, type: isl_dim_div, first: pos, n: 1) < 0)
4081	return NULL;
4082
4083	ls = isl_local_space_alloc_div(space: isl_space_copy(space: term->dim),
4084	div: isl_mat_copy(mat: term->div));
4085	aff = isl_aff_alloc(ls);
4086	if (!aff)
4087	return NULL;
4088
4089	isl_seq_cpy(dst: aff->v->el, src: term->div->row[pos], len: 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(poly: 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;
4158	error:
4159	isl_term_free(term);
4160	return NULL;
4161	}
4162
4163	isl_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(space: isl_space_copy(space: qp->dim), div: isl_mat_copy(mat: qp->div));
4172	if (!term)
4173	return isl_stat_error;
4174
4175	term = isl_poly_foreach_term(poly: 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, type: 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(ctx: term->dim->ctx, n: term->n, d: term->d);
4196	for (i = 0; i < n; ++i) {
4197	if (!term->pow[i])
4198	continue;
4199	poly = isl_poly_mul(poly1: poly,
4200	poly2: isl_poly_var_pow(ctx: term->dim->ctx, pos: i, power: term->pow[i]));
4201	}
4202
4203	qp = isl_qpolynomial_alloc(space: isl_space_copy(space: term->dim),
4204	n_div: term->div->n_row, poly);
4205	if (!qp)
4206	goto error;
4207	isl_mat_free(mat: qp->div);
4208	qp->div = isl_mat_copy(mat: term->div);
4209	if (!qp->div)
4210	goto error;
4211
4212	isl_term_free(term);
4213	return qp;
4214	error:
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(space1: qp->dim, space2: 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, type: isl_dim_set);
4240	d_qp = isl_qpolynomial_domain_dim(qp, type: isl_dim_set);
4241	extra = d_set - d_qp;
4242	total = isl_space_dim(space: qp->dim, type: 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(poly: qp->poly, exp, first: total);
4254	free(ptr: exp);
4255	if (!qp->poly)
4256	goto error;
4257	}
4258	qp->div = isl_mat_insert_cols(mat: qp->div, col: 2 + total, n: extra);
4259	if (!qp->div)
4260	goto error;
4261	for (i = 0; i < qp->div->n_row; ++i)
4262	isl_seq_clr(p: qp->div->row[i] + 2 + total, len: extra);
4263
4264	isl_space_free(space: qp->dim);
4265	qp->dim = space;
4266
4267	return qp;
4268	error:
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	*/
4277	static __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, type: 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(ptr: active);
4300	return set;
4301	}
4302
4303	nparam = isl_set_dim(set, type: isl_dim_param);
4304	nvar = isl_set_dim(set, type: 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, type: isl_dim_param, first: i, n: 1);
4311	set = isl_set_fix_si(set, type: isl_dim_param, pos: i, value: 0);
4312	}
4313	for (i = 0; i < nvar; ++i) {
4314	if (active[nparam + i])
4315	continue;
4316	set = isl_set_eliminate(set, type: isl_dim_set, first: i, n: 1);
4317	set = isl_set_fix_si(set, type: isl_dim_set, pos: i, value: 0);
4318	}
4319
4320	free(ptr: active);
4321
4322	return set;
4323	error:
4324	free(ptr: active);
4325	isl_set_free(set);
4326	return NULL;
4327	}
4328
4329	struct isl_opt_data {
4330	isl_qpolynomial *qp;
4331	int first;
4332	isl_val *opt;
4333	int max;
4334	};
4335
4336	static 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(qp: isl_qpolynomial_copy(qp: 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(v1: data->opt, v2: val);
4347	} else {
4348	data->opt = isl_val_min(v1: data->opt, v2: 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(poly: 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, fn: opt_fn, user: &data) < 0)
4377	goto error;
4378
4379	if (data.first)
4380	data.opt = isl_val_zero(ctx: isl_set_get_ctx(set));
4381
4382	isl_set_free(set);
4383	isl_qpolynomial_free(qp);
4384	return data.opt;
4385	error:
4386	isl_set_free(set);
4387	isl_qpolynomial_free(qp);
4388	isl_val_free(v: 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, f: morph->inv->row[1 + i],
4418	denom: morph->inv->row[0][0], len: 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, pos: morph->inv->n_col - 1 + i, power: 1);
4423
4424	qp->poly = isl_poly_subs(poly: qp->poly, first: 0, n: n_sub, subs);
4425
4426	for (i = 0; i < n_sub; ++i)
4427	isl_poly_free(poly: subs[i]);
4428	free(ptr: subs);
4429
4430	diag = isl_mat_diag(ctx, n_row: 1, d: morph->inv->row[0][0]);
4431	mat = isl_mat_diagonal(mat1: diag, mat2: isl_mat_copy(mat: morph->inv));
4432	diag = isl_mat_diag(ctx, n_row: qp->div->n_row, d: morph->inv->row[0][0]);
4433	mat = isl_mat_diagonal(mat1: mat, mat2: diag);
4434	qp->div = isl_mat_product(left: qp->div, right: mat);
4435	isl_space_free(space: qp->dim);
4436	qp->dim = isl_space_copy(space: 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;
4444	error:
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(u1: upwqp1, u2: upwqp2,
4455	fn: &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(exp: r, extra: qp->div->n_row);
4470	if (!r)
4471	goto error;
4472
4473	qp->div = isl_local_reorder(local: qp->div, r: isl_reordering_copy(exp: r));
4474	if (!qp->div)
4475	goto error;
4476
4477	qp->poly = reorder(poly: qp->poly, r: 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(exp: r);
4485	return qp;
4486	error:
4487	isl_qpolynomial_free(qp);
4488	isl_reordering_free(exp: 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(space1: domain_space, space2: model);
4500	if (equal_params < 0)
4501	goto error;
4502	if (!equal_params) {
4503	isl_reordering *exp;
4504
4505	exp = isl_parameter_alignment_reordering(alignee: domain_space, aligner: model);
4506	qp = isl_qpolynomial_realign_domain(qp, r: exp);
4507	}
4508
4509	isl_space_free(space: model);
4510	return qp;
4511	error:
4512	isl_space_free(space: model);
4513	isl_qpolynomial_free(qp);
4514	return NULL;
4515	}
4516
4517	struct 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 slice
4524	*
4525	* m v <= f <= m v + (m - 1)
4526	*
4527	* or
4528	*
4529	* f - m v >= 0
4530	* -f + m v + (m - 1) >= 0
4531	*/
4532	static __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, type: isl_dim_all);
4540	if (total < 0 || !qp)
4541	goto error;
4542
4543	bset = isl_basic_set_alloc_space(space: isl_space_copy(space), extra: 0, n_eq: 0, n_ineq: 2);
4544
4545	k = isl_basic_set_alloc_inequality(bset);
4546	if (k < 0)
4547	goto error;
4548	isl_seq_cpy(dst: bset->ineq[k], src: qp->div->row[div] + 1, len: 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(dst: bset->ineq[k], src: qp->div->row[div] + 1, len: 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);
4561	error:
4562	isl_basic_set_free(bset);
4563	isl_space_free(space);
4564	return NULL;
4565	}
4566
4567	static 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	*/
4574	static 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(space: isl_set_get_space(set), qp, div, v);
4584	set = isl_set_intersect(set1: set, set2: slice);
4585
4586	div_pos = isl_qpolynomial_domain_var_offset(qp, type: 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(ctx: qp->dim->ctx, n: v, d: qp->dim->ctx->one);
4599	qp = substitute_div(qp, div, s: cst);
4600
4601	return split_periods(set, qp, user: data);
4602	error:
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 slice
4610	* and add the results to data->res.
4611	*/
4612	static 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: set_i, qp: qp_i, div, v: min, data) < 0)
4621	goto error;
4622	}
4623	isl_set_free(set);
4624	isl_qpolynomial_free(qp);
4625	return isl_stat_ok;
4626	error:
4627	isl_set_free(set);
4628	isl_qpolynomial_free(qp);
4629	return isl_stat_error;
4630	}
4631
4632	/* If "qp" refers to any integer division
4633	* 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	*/
4638	static 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, el: qp);
4655	data->res = isl_pw_qpolynomial_add_disjoint(pw1: data->res, pw2: pwqp);
4656	return isl_stat_ok;
4657	}
4658
4659	div_pos = isl_qpolynomial_domain_var_offset(qp, type: 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(p: qp->div->row[i] + 2 + div_pos,
4669	len: qp->div->n_row) != -1)
4670	continue;
4671
4672	lp_res = isl_set_solve_lp(set, max: 0, f: qp->div->row[i] + 1,
4673	denom: set->ctx->one, opt: &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, max: 1, f: qp->div->row[i] + 1,
4681	denom: set->ctx->one, opt: &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, div: i, min, max, data);
4697	} else {
4698	pwqp = isl_pw_qpolynomial_alloc(set, el: qp);
4699	data->res = isl_pw_qpolynomial_add_disjoint(pw1: data->res, pw2: pwqp);
4700	}
4701
4702	isl_int_clear(max);
4703	isl_int_clear(min);
4704
4705	return r;
4706	error2:
4707	isl_int_clear(max);
4708	isl_int_clear(min);
4709	error:
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 division
4716	* that can only attain "max_periods" distinct values on its domain
4717	* 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(space: isl_pw_qpolynomial_get_space(pw: pwqp));
4726
4727	if (isl_pw_qpolynomial_foreach_piece(pw: pwqp, fn: &split_periods, user: &data) < 0)
4728	goto error;
4729
4730	isl_pw_qpolynomial_free(pw: pwqp);
4731
4732	return data.res;
4733	error:
4734	isl_pw_qpolynomial_free(pw: data.res);
4735	isl_pw_qpolynomial_free(pw: pwqp);
4736	return NULL;
4737	}
4738
4739	/* Construct a piecewise quasipolynomial that is constant on the given
4740	* domain. In particular, it is
4741	* 0 if cst == 0
4742	* 1 if cst == 1
4743	* infinity if cst == -1
4744	*
4745	* If cst == -1, then explicitly check whether the domain is empty and,
4746	* if so, return 0 instead.
4747	*/
4748	static __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(domain: space);
4763	else if (cst == 0)
4764	qp = isl_qpolynomial_zero_on_domain(domain: space);
4765	else
4766	qp = isl_qpolynomial_one_on_domain(domain: space);
4767	return isl_pw_qpolynomial_alloc(set: isl_set_from_basic_set(bset), el: 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	*/
4774	struct 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 or
4781	* implicit equality constraints.
4782	* If so, there may be further opportunities for detecting factors or
4783	* removing equality constraints, so recursively call
4784	* the top-level isl_basic_set_multiplicative_call.
4785	*/
4786	static __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 applies
4807	* data->fn to the factor "bset" and multiplies in the result
4808	* in data->pwqp.
4809	*/
4810	static 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, fn: data->fn);
4818	data->pwqp = isl_pw_qpolynomial_mul(pwqp1: data->pwqp, pwqp2: 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	*/
4831	static __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(space: isl_space_copy(space));
4853	qp = isl_qpolynomial_one_on_domain(domain: space);
4854	data.pwqp = isl_pw_qpolynomial_alloc(set, el: qp);
4855
4856	every = isl_factorizer_every_factor_basic_set(f,
4857	test: &multiplicative_call_factor_pw_qpolynomial, user: &data);
4858	if (every < 0)
4859	data.pwqp = isl_pw_qpolynomial_free(pw: data.pwqp);
4860
4861	isl_basic_set_free(bset);
4862	isl_factorizer_free(f);
4863
4864	return data.pwqp;
4865	error:
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 domains
4873	* 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, cst: 0);
4892
4893	dim = isl_basic_set_dim(bset, type: isl_dim_set);
4894	if (dim < 0)
4895	goto error;
4896	if (dim == 0)
4897	return constant_on_domain(bset, cst: 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, cst: -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(morph: 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(pw: pwqp, morph);
4918
4919	return pwqp;
4920	error:
4921	isl_basic_set_free(bset);
4922	return NULL;
4923	}
4924
4925	/* Drop all floors in "qp", turning each integer division [a/m] into
4926	* a rational division a/m. If "down" is set, then the integer division
4927	* is replaced by (a-(m-1))/m instead.
4928	*/
4929	static __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(ctx: qp->dim->ctx, f: qp->div->row[i] + 1,
4952	denom: qp->div->row[i][0], len: qp->div->n_col - 1);
4953	qp = substitute_div(qp, div: 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] into
4962	* a rational division a/m.
4963	*/
4964	static __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(pw: pwqp))
4973	return pwqp;
4974
4975	pwqp = isl_pw_qpolynomial_cow(pw: pwqp);
4976	if (!pwqp)
4977	return NULL;
4978
4979	for (i = 0; i < pwqp->n; ++i) {
4980	pwqp->p[i].qp = qp_drop_floors(qp: pwqp->p[i].qp, down: 0);
4981	if (!pwqp->p[i].qp)
4982	goto error;
4983	}
4984
4985	return pwqp;
4986	error:
4987	isl_pw_qpolynomial_free(pw: pwqp);
4988	return NULL;
4989	}
4990
4991	/* Adjust all the integer divisions in "qp" such that they are at least
4992	* one over the given orthant (identified by "signs"). This ensures
4993	* 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 have
5000	* 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	*/
5006	static __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, type: isl_dim_div);
5016	if (div_pos < 0)
5017	return isl_qpolynomial_free(qp);
5018	qp->div = isl_mat_cow(mat: qp->div);
5019	if (!qp->div)
5020	goto error;
5021
5022	v = isl_vec_alloc(ctx: qp->div->ctx, size: 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(vec: 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	else
5040	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(dst: qp->div->row[j] + 1,
5055	m1: qp->div->ctx->one, src1: qp->div->row[j] + 1,
5056	m2: qp->div->row[j][2 + div_pos + i], src2: v->el,
5057	len: v->size);
5058	}
5059	isl_int_set_si(v->el[1 + div_pos + i], 1);
5060	s = isl_poly_from_affine(ctx: qp->dim->ctx, f: v->el,
5061	denom: qp->div->ctx->one, len: v->size);
5062	qp->poly = isl_poly_subs(poly: qp->poly, first: div_pos + i, n: 1, subs: &s);
5063	isl_poly_free(poly: s);
5064	if (!qp->poly)
5065	goto error;
5066	}
5067
5068	isl_vec_free(vec: v);
5069	return qp;
5070	error:
5071	isl_vec_free(vec: v);
5072	isl_qpolynomial_free(qp);
5073	return NULL;
5074	}
5075
5076	struct 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 the
5084	* quasipolynomials into terms with sign data->sign (the direction
5085	* of the requested approximation) and terms with the opposite sign.
5086	* In the first set of terms, each integer division [a/m] is
5087	* overapproximated by a/m, while in the second it is underapproximated
5088	* by (a-(m-1))/m.
5089	*/
5090	static 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(qp: data->qp);
5098	qp = make_divs_pos(qp, signs);
5099
5100	up = isl_qpolynomial_terms_of_sign(poly: qp, signs, sign: data->sign);
5101	up = qp_drop_floors(qp: up, down: 0);
5102	down = isl_qpolynomial_terms_of_sign(poly: qp, signs, sign: -data->sign);
5103	down = qp_drop_floors(qp: down, down: 1);
5104
5105	isl_qpolynomial_free(qp);
5106	qp = isl_qpolynomial_add(qp1: up, qp2: down);
5107
5108	t = isl_pw_qpolynomial_alloc(set: orthant, el: qp);
5109	data->res = isl_pw_qpolynomial_add_disjoint(pw1: data->res, pw2: 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 approximation
5117	* will lie somewhere in between.
5118	*
5119	* In particular, is sign == 0, we simply drop the floors, turning
5120	* the integer divisions into rational divisions.
5121	* Otherwise, we split the domains into orthants, make all integer divisions
5122	* 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 which
5124	* 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(space: isl_pw_qpolynomial_get_space(pw: 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	set: isl_set_copy(set: pwqp->p[i].set),
5146	el: isl_qpolynomial_copy(qp: pwqp->p[i].qp));
5147	data.res = isl_pw_qpolynomial_add_disjoint(pw1: data.res, pw2: t);
5148	continue;
5149	}
5150	data.qp = pwqp->p[i].qp;
5151	if (isl_set_foreach_orthant(set: pwqp->p[i].set,
5152	fn: &to_polynomial_on_orthant, user: &data) < 0)
5153	goto error;
5154	}
5155
5156	isl_pw_qpolynomial_free(pw: pwqp);
5157
5158	return data.res;
5159	error:
5160	isl_pw_qpolynomial_free(pw: pwqp);
5161	isl_pw_qpolynomial_free(pw: data.res);
5162	return NULL;
5163	}
5164
5165	static __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: *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(u: upwqp,
5177	fn: &poly_entry, user: &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(poly: 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, type: isl_dim_out);
5204	n_div = qp->div->n_row;
5205	bmap = isl_basic_map_alloc_space(space, extra: n_div, n_eq: 1, n_ineq: 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(dst: bmap->div[k], src: qp->div->row[i], len: 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, div: 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(dst: bmap->eq[k], src: aff->el + 1, len: pos);
5220	isl_seq_cpy(dst: bmap->eq[k] + pos + 1, src: aff->el + 1 + pos, len: n_div);
5221
5222	isl_vec_free(vec: aff);
5223	isl_qpolynomial_free(qp);
5224	bmap = isl_basic_map_finalize(bmap);
5225	return bmap;
5226	error:
5227	isl_vec_free(vec: aff);
5228	isl_qpolynomial_free(qp);
5229	isl_basic_map_free(bmap);
5230	return NULL;
5231	}
5232