?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
iscc: support addition on a pair of a fold and a quasipolynomial.
[barvinok.git] / verify.c
blob108025972e1742f41e06eb2b0c0383906a3757f7
1 #include <assert.h>
2 #include <stdlib.h>
3 #include <barvinok/options.h>
4 #include <barvinok/util.h>
5 #include "verify.h"
6
7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
8
9 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
10 #define RANGE 50
11
12 /* SRANGE : small range for evalutations */
13 #define SRANGE 15
14
15 /* if dimension >= BIDDIM, use SRANGE */
16 #define BIGDIM 5
17
18 /* VSRANGE : very small range for evalutations */
19 #define VSRANGE 5
20
21 /* if dimension >= VBIDDIM, use VSRANGE */
22 #define VBIGDIM 8
23
24 static int set_m(void *opt, long val)
25 {
26 struct verify_options *options = (struct verify_options *)opt;
27 options->m = val;
28 options->verify = 1;
29 }
30
31 static int set_M(void *opt, long val)
32 {
33 struct verify_options *options = (struct verify_options *)opt;
34 options->M = val;
35 options->verify = 1;
36 }
37
38 static int set_r(void *opt, long val)
39 {
40 struct verify_options *options = (struct verify_options *)opt;
41 options->r = val;
42 options->m = -val;
43 options->M = val;
44 options->verify = 1;
45 }
46
47 struct isl_arg verify_options_arg[] = {
48 ISL_ARG_CHILD(struct verify_options, barvinok, NULL, barvinok_options_arg, NULL)
49 ISL_ARG_BOOL(struct verify_options, verify, 'T', "verify", 0, NULL)
50 ISL_ARG_BOOL(struct verify_options, exact, 'E', "exact", 0, NULL)
51 ISL_ARG_BOOL(struct verify_options, print_all, 'A', "print-all", 0, NULL)
52 ISL_ARG_BOOL(struct verify_options, continue_on_error, 'C',
53 "continue-on-error", 0, NULL)
54 ISL_ARG_USER_LONG(struct verify_options, m, 'm', "min", set_m, INT_MAX, NULL)
55 ISL_ARG_USER_LONG(struct verify_options, M, 'M', "max", set_M, INT_MIN, NULL)
56 ISL_ARG_USER_LONG(struct verify_options, r, 'r', "range", set_r, -1, NULL)
57 ISL_ARG_END
58 };
59
60 void verify_options_set_range(struct verify_options *options, int dim)
61 {
62 int r;
63
64 if (dim >= VBIGDIM)
65 r = VSRANGE;
66 else if (dim >= BIGDIM)
67 r = SRANGE;
68 else
69 r = RANGE;
70 /* If the user didn't set m or M, then we try to adjust the window
71 * to the context in check_poly_context_scan.
72 */
73 if (options->m == INT_MAX && options->M == INT_MIN)
74 options->r = r;
75 if (options->M == INT_MIN)
76 options->M = r;
77 if (options->m == INT_MAX)
78 options->m = -r;
79
80 if (options->verify && options->m > options->M) {
81 fprintf(stderr,"Nothing to do: min > max !\n");
82 exit(0);
83 }
84 }
85
86 static Polyhedron *project_on(Polyhedron *P, int i)
87 {
88 unsigned dim = P->Dimension;
89 Matrix *T;
90 Polyhedron *I;
91
92 if (dim == 1)
93 return Polyhedron_Copy(P);
94
95 T = Matrix_Alloc(2, dim+1);
96 value_set_si(T->p[0][i], 1);
97 value_set_si(T->p[1][dim], 1);
98 I = Polyhedron_Image(P, T, P->NbConstraints);
99 Matrix_Free(T);
100 return I;
101 }
102
103 static void set_bounds(Polyhedron *C, Value **rows, int i, unsigned nparam,
104 const struct verify_options *options)
105 {
106 Value min;
107 Value max;
108
109 value_init(min);
110 value_init(max);
111 value_set_si(min, options->m);
112 value_set_si(max, options->M);
113
114 if (options->r > 0) {
115 Polyhedron *I = project_on(C, i);
116 line_minmax(I, &min, &max);
117 if (value_cmp_si(min, options->M) >= 0)
118 value_add_int(max, min, options->r);
119 else if (value_cmp_si(max, options->m) <= 0)
120 value_sub_int(min, max, options->r);
121 else {
122 value_set_si(min, options->m);
123 value_set_si(max, options->M);
124 }
125 }
126
127 value_set_si(rows[0][0], 1);
128 value_set_si(rows[0][1+i], 1);
129 value_oppose(rows[0][1+nparam], min);
130 value_set_si(rows[1][0], 1);
131 value_set_si(rows[1][1+i], -1);
132 value_assign(rows[1][1+nparam], max);
133
134 value_clear(min);
135 value_clear(max);
136 }
137
138 Polyhedron *check_poly_context_scan(Polyhedron *P, Polyhedron **C,
139 unsigned nparam,
140 const struct verify_options *options)
141 {
142 int i;
143 Matrix *MM;
144 Polyhedron *CC, *CC2, *CS, *U;
145 unsigned MaxRays = options->barvinok->MaxRays;
146
147 if (nparam <= 0)
148 return NULL;
149
150 if (!P)
151 CC = *C;
152 else {
153 CC = Polyhedron_Project(P, nparam);
154 if (*C) {
155 CC2 = DomainIntersection(*C, CC, MaxRays);
156 Domain_Free(CC);
157 CC = CC2;
158 }
159 }
160
161 /* Intersect context with range */
162 MM = Matrix_Alloc(2*nparam, nparam+2);
163 for (i = 0; i < nparam; ++i)
164 set_bounds(CC, &MM->p[2*i], i, nparam, options);
165 CC2 = AddConstraints(MM->p[0], 2*nparam, CC, options->barvinok->MaxRays);
166 if (CC != *C)
167 Domain_Free(CC);
168 CC = CC2;
169 U = Universe_Polyhedron(0);
170 CS = Polyhedron_Scan(CC, U, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
171 Polyhedron_Free(U);
172 *C = CC;
173 Matrix_Free(MM);
174 return CS;
175 }
176
177 void check_poly_init(Polyhedron *C, struct verify_options *options)
178 {
179 int d, i;
180
181 if (options->print_all)
182 return;
183 if (C->Dimension <= 0)
184 return;
185
186 d = options->M - options->m;
187 if (d > 80)
188 options->st = 1+d/80;
189 else
190 options->st = 1;
191 for (i = options->m; i <= options->M; i += options->st)
192 printf(".");
193 printf( "\r" );
194 fflush(stdout);
195 }
196
197 static void print_rational(FILE *out, Value n, Value d)
198 {
199 value_print(out, VALUE_FMT, n);
200 if (value_notone_p(d)) {
201 fprintf(out, "/");
202 value_print(out, VALUE_FMT, d);
203 }
204 }
205
206 void check_poly_print(int ok, int nparam, Value *z,
207 Value want_n, Value want_d,
208 Value got_n, Value got_d,
209 const char *op, const char *check,
210 const char *long_op,
211 const struct verify_options *options)
212 {
213 int k;
214
215 if (options->print_all) {
216 printf("%s(", op);
217 value_print(stdout, VALUE_FMT, z[0]);
218 for (k = 1; k < nparam; ++k) {
219 printf(", ");
220 value_print(stdout, VALUE_FMT, z[k]);
221 }
222 printf(") = ");
223 print_rational(stdout, got_n, got_d);
224 printf(", %s = ", check);
225 print_rational(stdout, want_n, want_d);
226 printf(". ");
227 }
228
229 if (!ok) {
230 printf("\n");
231 fflush(stdout);
232 fprintf(stderr, "Error !\n");
233 fprintf(stderr, "%s(", op);
234 value_print(stderr, VALUE_FMT, z[0]);
235 for (k = 1; k < nparam; ++k) {
236 fprintf(stderr,", ");
237 value_print(stderr, VALUE_FMT, z[k]);
238 }
239 fprintf(stderr, ") should be ");
240 print_rational(stderr, want_n, want_d);
241 fprintf(stderr, ", while %s gives ", long_op);
242 print_rational(stderr, got_n, got_d);
243 fprintf(stderr, ".\n");
244 } else if (options->print_all)
245 printf("OK.\n");
246 }
247
248 /****************************************************/
249 /* function check_poly : */
250 /* scans the parameter space from Min to Max (all */
251 /* directions). Computes the number of points in */
252 /* the polytope using both methods, and compare them*/
253 /* returns 1 on success */
254 /****************************************************/
255
256 int check_poly(Polyhedron *CS, const struct check_poly_data *data,
257 int nparam, int pos, Value *z,
258 const struct verify_options *options)
259 {
260 if (pos == nparam) {
261 if (!data->check(data, nparam, z, options) && !options->continue_on_error)
262 return 0;
263 } else {
264 Value LB, UB;
265 int ok;
266 value_init(LB);
267 value_init(UB);
268 ok = !(lower_upper_bounds(1+pos, CS, z-1, &LB, &UB));
269 assert(ok);
270 for (; value_le(LB, UB); value_increment(LB, LB)) {
271 if (!options->print_all) {
272 int k = VALUE_TO_INT(LB);
273 if (!pos && !(k % options->st)) {
274 printf("o");
275 fflush(stdout);
276 }
277 }
278
279 value_assign(z[pos], LB);
280 if (!check_poly(CS->next, data, nparam, pos+1, z, options)) {
281 value_clear(LB);
282 value_clear(UB);
283 return 0;
284 }
285 }
286 value_set_si(z[pos], 0);
287 value_clear(LB);
288 value_clear(UB);
289 }
290 return 1;
291 } /* check_poly */
292
293 void check_EP_set_scan(struct check_EP_data *data, Polyhedron *C,
294 unsigned MaxRays)
295 {
296 const evalue *EP = data->EP;
297 int i;
298 int n_S = 0;
299
300 for (i = 0; i < EP->x.p->size/2; ++i) {
301 Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
302 for (; A; A = A->next)
303 ++n_S;
304 }
305
306 data->n_S = n_S;
307 data->S = ALLOCN(Polyhedron *, n_S);
308 n_S = 0;
309 for (i = 0; i < EP->x.p->size/2; ++i) {
310 Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
311 for (; A; A = A->next) {
312 Polyhedron *next = A->next;
313 A->next = NULL;
314 data->S[n_S++] = Polyhedron_Scan(A, C,
315 MaxRays & POL_NO_DUAL ? 0 : MaxRays);
316 A->next = next;
317 }
318 }
319 }
320
321 void check_EP_clear_scan(struct check_EP_data *data)
322 {
323 int i;
324
325 for (i = 0; i < data->n_S; ++i)
326 Domain_Free(data->S[i]);
327 free(data->S);
328 }
329
330 static int check_EP_on_poly(Polyhedron *P,
331 struct check_EP_data *data,
332 unsigned nvar, unsigned nparam,
333 struct verify_options *options)
334 {
335 Polyhedron *CS;
336 unsigned MaxRays = options->barvinok->MaxRays;
337 int ok = 1;
338 const evalue *EP = data->EP;
339
340 CS = check_poly_context_scan(NULL, &P, P->Dimension, options);
341
342 check_poly_init(P, options);
343
344 if (!(CS && emptyQ2(CS))) {
345 check_EP_set_scan(data, P, MaxRays);
346 ok = check_poly(CS, &data->cp, nparam, 0, data->cp.z+1+nvar, options);
347 check_EP_clear_scan(data);
348 }
349
350 if (!options->print_all)
351 printf("\n");
352
353 if (CS) {
354 Domain_Free(CS);
355 Domain_Free(P);
356 }
357
358 return ok;
359 }
360
361 /*
362 * Project on final dim dimensions
363 */
364 Polyhedron *DomainProject(Polyhedron *D, unsigned dim, unsigned MaxRays)
365 {
366 Polyhedron *P;
367 Polyhedron *R;
368
369 R = Polyhedron_Project(D, dim);
370 for (P = D->next; P; P = P->next) {
371 Polyhedron *R2 = Polyhedron_Project(P, dim);
372 Polyhedron *R3 = DomainUnion(R, R2, MaxRays);
373 Polyhedron_Free(R2);
374 Domain_Free(R);
375 R = R3;
376 }
377 return R;
378 }
379
380 static Polyhedron *evalue_parameter_domain(const evalue *e, unsigned nparam,
381 unsigned MaxRays)
382 {
383 int i;
384 Polyhedron *U = NULL;
385
386 if (EVALUE_IS_ZERO(*e))
387 return Universe_Polyhedron(0);
388
389 assert(value_zero_p(e->d));
390 assert(e->x.p->type == partition);
391 assert(e->x.p->size >= 2);
392
393 for (i = 0; i < e->x.p->size/2; ++i) {
394 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
395 Polyhedron *P = DomainProject(D, nparam, MaxRays);
396 if (!U) {
397 U = P;
398 } else {
399 Polyhedron *T = U;
400 U = DomainUnion(U, P, MaxRays);
401 Domain_Free(P);
402 Domain_Free(T);
403 }
404 }
405 return U;
406 }
407
408 int check_EP(struct check_EP_data *data, unsigned nvar, unsigned nparam,
409 struct verify_options *options)
410 {
411 Vector *p;
412 int i;
413 int ok = 1;
414 Polyhedron *D, *P;
415
416 p = Vector_Alloc(nvar+nparam+2);
417 value_set_si(p->p[nvar+nparam+1], 1);
418
419 data->cp.z = p->p;
420
421 D = evalue_parameter_domain(data->EP, nparam, options->barvinok->MaxRays);
422
423 for (P = D; P; P = P->next) {
424 ok = check_EP_on_poly(P, data, nvar, nparam, options);
425 if (!ok && !options->continue_on_error)
426 break;
427 }
428
429 Domain_Free(D);
430 Vector_Free(p);
431
432 return ok;
433 }
434
435 __isl_give isl_set *verify_context_set_bounds(__isl_take isl_set *set,
436 const struct verify_options *options)
437 {
438 int i;
439 unsigned nparam;
440 isl_point *pt, *pt2;
441 isl_set *box;
442
443 nparam = isl_set_dim(set, isl_dim_param);
444
445 if (options->r > 0) {
446 pt = isl_set_sample_point(isl_set_copy(set));
447 pt2 = isl_point_copy(pt);
448
449 for (i = 0; i < nparam; ++i) {
450 pt = isl_point_add_ui(pt, isl_dim_param, i, options->r);
451 pt2 = isl_point_sub_ui(pt2, isl_dim_param, i, options->r);
452 }
453 } else {
454 isl_int v;
455
456 isl_int_init(v);
457 pt = isl_point_zero(isl_set_get_dim(set));
458 isl_int_set_si(v, options->m);
459 for (i = 0; i < nparam; ++i)
460 pt = isl_point_set_coordinate(pt, isl_dim_param, i, v);
461
462 pt2 = isl_point_zero(isl_set_get_dim(set));
463 isl_int_set_si(v, options->M);
464 for (i = 0; i < nparam; ++i)
465 pt2 = isl_point_set_coordinate(pt2, isl_dim_param, i, v);
466
467 isl_int_clear(v);
468 }
469
470 box = isl_set_box_from_points(pt, pt2);
471
472 return isl_set_intersect(set, box);
473 }
474
475 int verify_point_data_init(struct verify_point_data *vpd,
476 __isl_keep isl_set *context)
477 {
478 isl_int v;
479 int i;
480 int r;
481
482 isl_int_init(v);
483 r = isl_set_count(context, &v);
484 vpd->n = isl_int_cmp_si(v, 200) < 0 ? isl_int_get_si(v) : 200;
485 isl_int_clear(v);
486
487 if (!vpd->options->print_all) {
488 vpd->s = vpd->n < 80 ? 1 : 1 + vpd->n/80;
489 for (i = 0; i < vpd->n; i += vpd->s)
490 printf(".");
491 printf("\r");
492 fflush(stdout);
493 }
494
495 vpd->error = r < 0 ? -1 : 0;
496
497 return r;
498 }
499
500 void verify_point_data_fini(struct verify_point_data *vpd)
501 {
502 if (!vpd->options->print_all)
503 printf("\n");
504
505 if (vpd->error)
506 fprintf(stderr, "Check failed !\n");
507 }
lattice point enumeration library