keep track of number of Bernoulli sums
[barvinok.git] / zsolve / libzsolve.c
blob0b0ad3e139d0602393591f88288c7ab8e133bd6c
1 /*
2 4ti2 -- A software package for algebraic, geometric and combinatorial
3 problems on linear spaces.
5 Copyright (C) 2006 4ti2 team.
6 Main author(s): Matthias Walter.
8 This program is free software; you can redistribute it and/or
9 modify it under the terms of the GNU General Public License
10 as published by the Free Software Foundation; either version 2
11 of the License, or (at your option) any later version.
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 #include <assert.h>
24 #include <stdlib.h>
26 #include "libzsolve.h"
27 #include "valuetrees.h"
28 #include "lattice.h"
30 int chooseNext(ZSolveContext ctx, Vector possible, int value)
32 int i,j,zeros;
34 assert(possible);
36 if (value<0)
37 value = 0;
38 for (i=0; i<ctx->Variables; i++)
40 possible[i] = (possible[i]==value) ? 1 : 0;
41 zeros = 1;
42 for (j=0; j<ctx->Lattice->Size; j++)
43 if (ctx->Lattice->Data[j][i]==0)
44 zeros++;
45 possible[i] *= zeros;
48 j = -1;
49 for (i=0; i<ctx->Variables; i++)
50 if (possible[i]>0 && (j<0 || possible[i]>possible[j]))
51 j = i;
53 return j;
56 int nextVariable(ZSolveContext ctx)
58 IndexArray Zeros;
59 Vector gcds;
60 int i,j,k,col = -1;
61 int factor, value;
62 bool flag;
64 assert(ctx);
66 // Zeros is the list of vectors v with ||v|| = 0
67 Zeros = createIndexArray();
68 for (i=0; i<ctx->Lattice->Size; i++)
69 if (normVector(ctx->Lattice->Data[i], ctx->Current)==0)
70 appendToIndexArray(Zeros, i);
72 if (Zeros->Size)
74 gcds = createVector(ctx->Variables);
75 for (i=0; i<ctx->Variables; i++)
76 gcds[i] = -1;
77 value = -1;
78 for (i=ctx->Current; i<ctx->Variables; i++)
80 if (!ctx->Lattice->Properties[i].Free)
82 gcds[i] = abs(ctx->Lattice->Data[Zeros->Data[0]][i]);
83 for (j=1; j<Zeros->Size; j++)
84 gcds[i] = gcd(gcds[i], abs(ctx->Lattice->Data[Zeros->Data[j]][i]));
85 if ((value<0 || gcds[i]<value) && gcds[i]>0)
86 value = gcds[i];
90 col = chooseNext(ctx, gcds, value);
92 deleteVector(gcds);
93 if (col<0)
94 return -1;
95 // minimize Zeros
98 j = -1;
99 for (i=0; i<Zeros->Size; i++)
101 if (ctx->Lattice->Data[Zeros->Data[i]][col]!=0 && (j<0 || abs(ctx->Lattice->Data[Zeros->Data[i]][col])<abs(ctx->Lattice->Data[Zeros->Data[j]][col])))
102 j = i;
104 if (j<0)
105 break;
106 j = Zeros->Data[j];
108 printf("Current Norm[0]: ");
109 printVector(Zeros->Data, Zeros->Size, VECTOR_NEWLINE);
110 printf("ctx->Lattice:\n");
111 printVectorArray(ctx->Lattice, VECTORARRAY_INDENT);
112 printf("Reducing Slot: Row = %d, Col = %d\n\n", j, col);
114 flag = false;
115 for (i=0; i<ctx->Lattice->Size; i++)
117 if (ctx->Lattice->Data[i][col]!=0 && i!=j)
119 factor = ctx->Lattice->Data[i][col] / ctx->Lattice->Data[j][col];
120 for (k=0; k<ctx->Variables; k++)
121 if (ctx->Lattice->Data[i][k]+ctx->Lattice->Data[j][k]!=0)
122 break;
123 if (k==ctx->Variables)
124 continue;
125 if (factor!=0)
127 // printf("factor = %d\n", factor);
128 flag = true;
130 for (k=0; k<ctx->Variables; k++)
131 ctx->Lattice->Data[i][k] -= factor * ctx->Lattice->Data[j][k];
134 } while (flag);
136 else
138 gcds = createVector(ctx->Variables);
139 for (i=0; i<ctx->Variables; i++)
140 gcds[i] = (i<ctx->Current || ctx->Lattice->Properties[i].Free) ? 0 : 1;
141 col = chooseNext(ctx, gcds, 1);
142 deleteVector(gcds);
145 deleteIndexArray(Zeros);
147 return col;
150 void filterLimits(ZSolveContext ctx)
152 int i;
154 assert(ctx);
156 for (i=0; i<ctx->Lattice->Size; i++)
158 if (ctx->Lattice->Properties[ctx->Current].Lower>ctx->Lattice->Data[i][ctx->Current] || ctx->Lattice->Properties[ctx->Current].Upper<ctx->Lattice->Data[i][ctx->Current])
160 ctx->Lattice->Size--;
161 deleteVector(ctx->Lattice->Data[i]);
162 ctx->Lattice->Data[i] = ctx->Lattice->Data[ctx->Lattice->Size];
163 i--;
166 ctx->Lattice->Memory = ctx->Lattice->Size;
167 ctx->Lattice->Data = (Vector *)realloc(ctx->Lattice->Data, ctx->Lattice->Memory * sizeof(Vector));
170 void zsolveLogCallbackDefault(FILE *stream, int level, int type, int var, int sum, int norm, int vectors, CPUTime alltime, CPUTime steptime)
172 switch (level)
174 case 1:
175 if (type == ZSOLVE_LOG_VARIABLE_STARTED)
176 fprintf(stream, "Appending variable %d... ", var+1);
177 else if (type == ZSOLVE_LOG_VARIABLE_FINISHED)
179 fprintf(stream, "Solutions: %d, Step: ", vectors);
180 fprintCPUTime(stream, steptime);
181 fprintf(stream, ", Time: ");
182 fprintCPUTime(stream, alltime);
183 fprintf(stream, "\n");
185 break;
186 case 2:
187 if (type == ZSOLVE_LOG_VARIABLE_STARTED)
188 fprintf(stream, "%sAppending variable %d.\n\n", var ? "\n" : "", var+1);
189 else if (type == ZSOLVE_LOG_VARIABLE_FINISHED)
191 fprintf(stream, "\nFinished variable %d. Solutions: %d, Step: ", var+1, vectors);
192 fprintCPUTime(stream, steptime);
193 fprintf(stream, ", Time: ");
194 fprintCPUTime(stream, alltime);
195 fprintf(stream, "\n");
197 else if (type == ZSOLVE_LOG_SUM_STARTED)
198 fprintf(stream, " Var = %d, Sum = %d... ", var+1, sum);
199 else if (type == ZSOLVE_LOG_SUM_FINISHED)
201 fprintf(stream, "Solutions: %d, Step: ", vectors);
202 fprintCPUTime(stream, steptime);
203 fprintf(stream, ", Time: ");
204 fprintCPUTime(stream, alltime);
205 fprintf(stream, "\n");
207 break;
208 case 3:
209 if (type == ZSOLVE_LOG_VARIABLE_STARTED)
210 fprintf(stream, "%sAppending variable %d.\n\n", var ? "\n" : "", var+1);
211 else if (type == ZSOLVE_LOG_VARIABLE_FINISHED)
213 fprintf(stream, "\nFinished variable %d. Solutions: %d, Step: ", var+1, vectors);
214 fprintCPUTime(stream, steptime);
215 fprintf(stream, ", Time: ");
216 fprintCPUTime(stream, alltime);
217 fprintf(stream, "\n");
219 else if (type == ZSOLVE_LOG_SUM_STARTED)
220 fprintf(stream, "%s Var = %d, Sum = %d.\n\n", sum ? "\n" : "", var+1, sum);
221 else if (type == ZSOLVE_LOG_SUM_FINISHED)
223 fprintf(stream, "\n Finished sum %d. Solutions: %d, Step: ", sum, vectors);
224 fprintCPUTime(stream, steptime);
225 fprintf(stream, ", Time: ");
226 fprintCPUTime(stream, alltime);
227 fprintf(stream, "\n");
229 else if (type == ZSOLVE_LOG_NORM_STARTED)
230 fprintf(stream, " Var = %d, Norm = %d + %d... ", var+1, norm, sum-norm);
231 else if (type == ZSOLVE_LOG_NORM_FINISHED)
233 fprintf(stream, "Solutions: %d, Step: ", vectors);
234 fprintCPUTime(stream, steptime);
235 fprintf(stream, ", Time: ");
236 fprintCPUTime(stream, alltime);
237 fprintf(stream, "\n");
239 break;
241 if (level>0)
243 if (type == ZSOLVE_LOG_STARTED)
245 // need linearsystem - maybe with va
247 else if (type == ZSOLVE_LOG_RESUMED)
249 fprintf(stream, "Resuming from backup, starting at Var = %d, Norm = %d + %d... Solutions: %d, Time: ", var, norm, sum - norm, vectors);
250 fprintCPUTime(stream, getCPUTime() - alltime);
251 fprintf(stream, "\n\n");
253 else if (type == ZSOLVE_LOG_FINISHED)
255 // need homs / inhoms - maybe with va
258 if (stream)
259 fflush(stream);
263 void splitLog(ZSolveContext ctx, int type, CPUTime steptime)
265 CPUTime currenttime = getCPUTime();
267 if (ctx->LogCallback)
269 ctx->LogCallback(stdout, ctx->Verbosity, type, ctx->Current, ctx->SumNorm, ctx->FirstNorm, ctx->Lattice->Size, maxd(currenttime - ctx->AllTime, 0.0), maxd(steptime, 0.0));
270 if (ctx->LogFile)
271 ctx->LogCallback(ctx->LogFile, ctx->LogLevel, type, ctx->Current, ctx->SumNorm, ctx->FirstNorm, ctx->Lattice->Size, maxd(currenttime - ctx->AllTime, 0.0), maxd(steptime, 0.0));
275 void backupZSolveContext(FILE *stream, ZSolveContext ctx)
277 CPUTime currenttime = getCPUTime();
279 assert(stream);
280 assert(ctx);
282 fprintf(stream, "%d %d %d\n%d\n\n", ctx->Current, ctx->SumNorm, ctx->FirstNorm, ctx->Symmetric);
284 fprintf(stream, "%d %d %d\n\n", ctx->Verbosity, ctx->LogLevel, ctx->BackupTime);
286 fprintf(stream, "%.2f %.2f %.2f %.2f\n\n" , currenttime - ctx->AllTime, currenttime - ctx->VarTime, currenttime - ctx->SumTime, currenttime - ctx->NormTime);
288 fprintVectorArray(stream, ctx->Lattice, true);
291 int originalCompare(int a, int b)
293 if (a>=0 && b<0)
294 return -1;
295 if (a<0 && b>=0)
296 return 1;
297 return a-b;
300 ZSolveContext createZSolveContextFromSystem(LinearSystem initialsystem, FILE *logfile, int loglevel, int verbosity, ZSolveLogCallback logcallback, ZSolveBackupCallback backupcallback)
302 ZSolveContext ctx;
303 LinearSystem finalsystem;
305 ctx = (ZSolveContext)malloc(sizeof(zsolvecontext_t));
306 if (ctx==NULL)
308 fprintf(stderr, "Fatal Error (%s/%d): Could not allocate memory for ZSolveContext in createZSolveContextFromSystem!\n", __FILE__, __LINE__);
309 exit(1);
312 ctx->BackupTime = 0;
313 ctx->Verbosity = verbosity;
314 ctx->LogLevel = loglevel;
315 ctx->LogFile = logfile;
317 if (ctx->Verbosity>0)
319 printf("Linear integer system to solve:\n\n");
320 printLinearSystem(initialsystem);
322 if (ctx->LogLevel>0)
324 fprintf(ctx->LogFile, "Linear integer system to solve:\n\n");
325 fprintLinearSystem(ctx->LogFile, initialsystem);
328 finalsystem = homogenizeLinearSystem(initialsystem);
329 deleteLinearSystem(initialsystem);
331 if (ctx->Verbosity>0)
333 printf("\nLinear integer system of homogeneous equalities to solve:\n\n");
334 printLinearSystem(finalsystem);
335 printf("\n\n");
337 if (ctx->LogLevel>0)
339 fprintf(ctx->LogFile, "\nLinear integer system of homogeneous equalities to solve:\n\n");
340 fprintLinearSystem(ctx->LogFile, finalsystem);
341 fprintf(ctx->LogFile, "\n\n");
344 ctx->Lattice = generateLattice(finalsystem);
345 deleteLinearSystem(finalsystem);
347 // sort 0, 1, 2, ... , -2, -2, -2, -2, -1 (original order, then slackvars, then -b)
349 sortVectorArrayColumns(ctx->Lattice, originalCompare);
351 ctx->Variables = ctx->Lattice->Variables;
352 ctx->MaxNorm = 0;
353 ctx->Symmetric = true;
354 ctx->Current = 0;
355 ctx->SumNorm = 0;
356 ctx->FirstNorm = 0;
357 ctx->AllTime = getCPUTime();
358 ctx->BackupCallback = backupcallback;
359 ctx->LogCallback = logcallback;
360 ctx->Homs = NULL;
361 ctx->Inhoms = NULL;
362 ctx->Frees = NULL;
363 ctx->Graver = NULL;
365 return ctx;
368 ZSolveContext createZSolveContextFromLattice(VectorArray lattice, FILE *logfile, int loglevel, int verbosity, ZSolveLogCallback logcallback, ZSolveBackupCallback backupcallback)
370 ZSolveContext ctx;
372 assert(lattice);
374 ctx = (ZSolveContext)malloc(sizeof(zsolvecontext_t));
375 if (ctx==NULL)
377 fprintf(stderr, "Fatal Error (%s/%d): Could not allocate memory for ZSolveContext in createZSolveContextFromLattice!\n", __FILE__, __LINE__);
378 exit(1);
381 ctx->BackupTime = 0;
382 ctx->Verbosity = verbosity;
383 ctx->LogLevel = loglevel;
384 ctx->LogFile = logfile;
385 ctx->Lattice = lattice;
386 ctx->Variables = ctx->Lattice->Variables;
387 ctx->MaxNorm = 0;
388 ctx->Symmetric = true;
389 ctx->Current = 0;
390 ctx->SumNorm = 0;
391 ctx->FirstNorm = 0;
392 ctx->AllTime = getCPUTime();
393 ctx->BackupCallback = backupcallback;
394 ctx->LogCallback = logcallback;
395 ctx->Homs = NULL;
396 ctx->Inhoms = NULL;
397 ctx->Frees = NULL;
398 ctx->Graver = NULL;
400 return ctx;
403 ZSolveContext createZSolveContextFromBackup(FILE *stream, ZSolveLogCallback logcallback, ZSolveBackupCallback backupcallback)
405 ZSolveContext ctx;
406 CPUTime currenttime = getCPUTime();
408 assert(stream);
410 ctx = (ZSolveContext)malloc(sizeof(zsolvecontext_t));
411 if (ctx==NULL)
413 fprintf(stderr, "Fatal Error (%s/%d): Could not allocate memory for ZSolveContext in createZSolveContextFromBackup!\n", __FILE__, __LINE__);
414 exit(1);
417 fscanf(stream, "%d %d %d %d", &(ctx->Current), &(ctx->SumNorm), &(ctx->FirstNorm), &(ctx->Symmetric));
418 ctx->SecondNorm = ctx->SumNorm - ctx->FirstNorm;
420 fscanf(stream, "%d %d %d", &(ctx->Verbosity), &(ctx->LogLevel), &(ctx->BackupTime));
422 fscanf(stream, "%lf %lf %lf %lf", &(ctx->AllTime), &(ctx->VarTime), &(ctx->SumTime), &(ctx->NormTime));
423 ctx->AllTime = currenttime - ctx->AllTime;
424 ctx->VarTime = currenttime - ctx->VarTime;
425 ctx->SumTime = currenttime - ctx->SumTime;
426 ctx->NormTime = currenttime - ctx->NormTime;
428 ctx->Lattice = readVectorArray(stream, true);
430 ctx->Variables = ctx->Lattice->Variables;
432 ctx->BackupCallback = backupcallback;
433 ctx->LogCallback = logcallback;
434 ctx->Homs = NULL;
435 ctx->Inhoms = NULL;
436 ctx->Frees = NULL;
437 ctx->Graver = NULL;
439 return ctx;
442 void zsolveSystem(ZSolveContext ctx, bool appendnegatives)
444 int next;
445 int split;
446 int count;
447 int i,j,k;
448 bool is_hom, is_free, has_symmetric;
449 int lex_cmp;
450 Vector vector;
451 CPUTime lastbackup = getCPUTime();
453 if (appendnegatives)
454 appendVectorArrayNegatives(ctx->Lattice);
456 ctx->Sum = createVector(ctx->Variables);
458 if (ctx->SumNorm != 0)
459 createValueTrees(ctx);
461 while (ctx->Current < ctx->Lattice->Variables)
463 if (ctx->BackupTime > 0 && ctx->BackupCallback && (getCPUTime() - lastbackup >= (double)ctx->BackupTime))
465 ctx->BackupCallback(ctx);
466 lastbackup = getCPUTime();
469 if (ctx->SumNorm == 0) // start of var loop
471 ctx->VarTime = getCPUTime();
472 next = nextVariable(ctx);
473 // printf("next variable = %d\n", next);
474 if (next<0)
475 break;
476 splitLog(ctx, ZSOLVE_LOG_VARIABLE_STARTED, 0.0);
477 swapVectorArrayColumns(ctx->Lattice, next, ctx->Current);
478 // printf("after next_var_swap:\n");
479 // printVectorArray(ctx->Lattice, true);
480 createValueTrees(ctx);
483 if (ctx->FirstNorm == 0) // start of sum loop
485 ctx->SumTime = getCPUTime();
486 splitLog(ctx, ZSOLVE_LOG_SUM_STARTED, 0.0);
489 // start of norm loop
491 ctx->SecondNorm = ctx->SumNorm - ctx->FirstNorm;
492 ctx->NormTime = getCPUTime();
493 splitLog(ctx, ZSOLVE_LOG_NORM_STARTED, 0.0);
495 if (ctx->FirstNorm <= ctx->MaxNorm && ctx->SecondNorm <= ctx->MaxNorm && ctx->Norm[ctx->FirstNorm]!=NULL && ctx->Norm[ctx->SecondNorm]!=NULL)
497 // start completition procedure
498 completeValueTrees(ctx, ctx->FirstNorm, ctx->SumNorm - ctx->FirstNorm);
501 // end of norm loop
503 splitLog(ctx, ZSOLVE_LOG_NORM_FINISHED, ctx->NormTime);
504 ctx->FirstNorm++;
506 if (ctx->FirstNorm > ctx->SumNorm/2) // end of sum loop?
508 splitLog(ctx, ZSOLVE_LOG_SUM_FINISHED, ctx->SumTime);
509 ctx->SumNorm++;
511 if (ctx->SumNorm > 2*ctx->MaxNorm) // end of var loop?
513 deleteValueTrees(ctx);
514 if (ctx->Lattice->Properties[ctx->Current].Lower+ctx->Lattice->Properties[ctx->Current].Upper!=0)
515 ctx->Symmetric = false;
516 filterLimits(ctx);
517 splitLog(ctx, ZSOLVE_LOG_VARIABLE_FINISHED, ctx->VarTime);
518 ctx->Current++;
519 ctx->SumNorm = 0;
521 // printf("After COMPLETE and FILTER\n");
522 // printVectorArray(ctx->Lattice, true);
525 ctx->FirstNorm = 0;
530 // algorithm completed
532 split = -1;
533 count = 0;
535 // repermute
537 for (i = 0; i<ctx->Variables-1; i++) {
538 k = i;
539 for (j=i+1; j<ctx->Variables; j++) {
540 if (ctx->Lattice->Properties[k].Column < 0 || (ctx->Lattice->Properties[j].Column >=0 && ctx->Lattice->Properties[j].Column < ctx->Lattice->Properties[k].Column))
541 k = j;
543 swapVectorArrayColumns(ctx->Lattice, i, k);
546 // DEBUG
547 // printVectorArray(ctx->Lattice, true);
549 // split - variable for hom / inhom
550 // count - number of columns for output
552 for (i=0; i<ctx->Variables; i++)
554 if (ctx->Lattice->Properties[i].Column == -2)
555 split = i;
556 else if (ctx->Lattice->Properties[i].Column >=0 )
557 count = maxi(count, ctx->Lattice->Properties[i].Column+1);
560 ctx->Homs = createVectorArray(count);
561 ctx->Inhoms = createVectorArray(count);
562 ctx->Frees = createVectorArray(count);
563 ctx->Graver = createVectorArray(count);
565 // if no splitting, inhom only contains (0,...,0)
567 if (split<0)
568 appendToVectorArray(ctx->Inhoms, createZeroVector(count));
570 for (i=0; i<ctx->Lattice->Size; i++)
572 vector = ctx->Lattice->Data[i];
573 // hom depends on splitting variable
574 is_hom = split < 0 || vector[split] == 0;
575 // free, if all nonzero-vars are free
576 is_free = true;
577 for (j=0; j<ctx->Lattice->Variables; j++)
578 if (vector[j] != 0 && !ctx->Lattice->Properties[j].Free)
579 is_free = false;
580 // has_symmetric, if inverse fits the bounds
581 has_symmetric = true;
582 for (j=0; j<ctx->Lattice->Variables; j++)
583 if (!checkVariableBounds(ctx->Lattice->Properties, j, -vector[j]))
584 has_symmetric = false;
585 // lex compare of vector with its inverse
586 lex_cmp = lexCompareInverseVector(vector, ctx->Lattice->Variables);
588 assert(!is_free || has_symmetric);
590 if (is_free) {
591 if (!has_symmetric || lex_cmp>0) {
592 appendToVectorArray(ctx->Frees, copyVector(vector, count));
593 appendToVectorArray(ctx->Graver, copyVector(vector, count));
596 else {
597 if (is_hom) {
598 appendToVectorArray(ctx->Homs, copyVector(vector, count));
599 if (!has_symmetric || lex_cmp>0)
600 appendToVectorArray(ctx->Graver, copyVector(vector, count));
602 else {
603 appendToVectorArray(ctx->Inhoms, copyVector(vector, count));
608 if (ctx->Verbosity>0)
609 printf("\nFinal basis has %d inhomogeneous, %d homogeneous and %d free elements.\n", ctx->Inhoms->Size, ctx->Homs->Size, ctx->Frees->Size);
610 if (ctx->LogLevel>0)
611 fprintf(ctx->LogFile, "\nFinal basis has %d inhomogeneous, %d homogeneous and %d free elements.\n", ctx->Inhoms->Size, ctx->Homs->Size, ctx->Frees->Size);
614 void deleteZSolveContext(ZSolveContext ctx, bool deleteresult)
616 deleteVector(ctx->Sum);
617 deleteVectorArray(ctx->Lattice);
619 if (deleteresult)
621 deleteVectorArray(ctx->Homs);
622 deleteVectorArray(ctx->Inhoms);
623 deleteVectorArray(ctx->Frees);
624 deleteVectorArray(ctx->Graver);
627 free(ctx);