Coverage Report

Created: 2017-03-28 09:59

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