Coverage Report

Created: 2018-02-20 12:54

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