Coverage Report

Created: 2018-06-24 14:39

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