Coverage Report

Created: 2018-04-23 18:20

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