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.
26 #include "libzsolve.h"
27 #include "valuetrees.h"
30 int chooseNext(ZSolveContext ctx
, Vector possible
, int value
)
38 for (i
=0; i
<ctx
->Variables
; i
++)
40 possible
[i
] = (possible
[i
]==value
) ? 1 : 0;
42 for (j
=0; j
<ctx
->Lattice
->Size
; j
++)
43 if (ctx
->Lattice
->Data
[j
][i
]==0)
49 for (i
=0; i
<ctx
->Variables
; i
++)
50 if (possible
[i
]>0 && (j
<0 || possible
[i
]>possible
[j
]))
56 int nextVariable(ZSolveContext 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
);
74 gcds
= createVector(ctx
->Variables
);
75 for (i
=0; i
<ctx
->Variables
; i
++)
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)
90 col
= chooseNext(ctx
, gcds
, value
);
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
])))
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);
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)
123 if (k
==ctx
->Variables
)
127 // printf("factor = %d\n", factor);
130 for (k
=0; k
<ctx
->Variables
; k
++)
131 ctx
->Lattice
->Data
[i
][k
] -= factor
* ctx
->Lattice
->Data
[j
][k
];
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);
145 deleteIndexArray(Zeros
);
150 void filterLimits(ZSolveContext 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
];
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
)
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");
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");
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");
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
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));
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();
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
)
300 ZSolveContext
createZSolveContextFromSystem(LinearSystem initialsystem
, FILE *logfile
, int loglevel
, int verbosity
, ZSolveLogCallback logcallback
, ZSolveBackupCallback backupcallback
)
303 LinearSystem finalsystem
;
305 ctx
= (ZSolveContext
)malloc(sizeof(zsolvecontext_t
));
308 fprintf(stderr
, "Fatal Error (%s/%d): Could not allocate memory for ZSolveContext in createZSolveContextFromSystem!\n", __FILE__
, __LINE__
);
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
);
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
);
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
;
353 ctx
->Symmetric
= true;
357 ctx
->AllTime
= getCPUTime();
358 ctx
->BackupCallback
= backupcallback
;
359 ctx
->LogCallback
= logcallback
;
368 ZSolveContext
createZSolveContextFromLattice(VectorArray lattice
, FILE *logfile
, int loglevel
, int verbosity
, ZSolveLogCallback logcallback
, ZSolveBackupCallback backupcallback
)
374 ctx
= (ZSolveContext
)malloc(sizeof(zsolvecontext_t
));
377 fprintf(stderr
, "Fatal Error (%s/%d): Could not allocate memory for ZSolveContext in createZSolveContextFromLattice!\n", __FILE__
, __LINE__
);
382 ctx
->Verbosity
= verbosity
;
383 ctx
->LogLevel
= loglevel
;
384 ctx
->LogFile
= logfile
;
385 ctx
->Lattice
= lattice
;
386 ctx
->Variables
= ctx
->Lattice
->Variables
;
388 ctx
->Symmetric
= true;
392 ctx
->AllTime
= getCPUTime();
393 ctx
->BackupCallback
= backupcallback
;
394 ctx
->LogCallback
= logcallback
;
403 ZSolveContext
createZSolveContextFromBackup(FILE *stream
, ZSolveLogCallback logcallback
, ZSolveBackupCallback backupcallback
)
406 CPUTime currenttime
= getCPUTime();
410 ctx
= (ZSolveContext
)malloc(sizeof(zsolvecontext_t
));
413 fprintf(stderr
, "Fatal Error (%s/%d): Could not allocate memory for ZSolveContext in createZSolveContextFromBackup!\n", __FILE__
, __LINE__
);
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
;
442 void zsolveSystem(ZSolveContext ctx
, bool appendnegatives
)
448 bool is_hom
, is_free
, has_symmetric
;
451 CPUTime lastbackup
= getCPUTime();
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);
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
);
503 splitLog(ctx
, ZSOLVE_LOG_NORM_FINISHED
, ctx
->NormTime
);
506 if (ctx
->FirstNorm
> ctx
->SumNorm
/2) // end of sum loop?
508 splitLog(ctx
, ZSOLVE_LOG_SUM_FINISHED
, ctx
->SumTime
);
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;
517 splitLog(ctx
, ZSOLVE_LOG_VARIABLE_FINISHED
, ctx
->VarTime
);
521 // printf("After COMPLETE and FILTER\n");
522 // printVectorArray(ctx->Lattice, true);
530 // algorithm completed
537 for (i
= 0; i
<ctx
->Variables
-1; 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
))
543 swapVectorArrayColumns(ctx
->Lattice
, i
, k
);
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)
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)
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
577 for (j
=0; j
<ctx
->Lattice
->Variables
; j
++)
578 if (vector
[j
] != 0 && !ctx
->Lattice
->Properties
[j
].Free
)
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
);
591 if (!has_symmetric
|| lex_cmp
>0) {
592 appendToVectorArray(ctx
->Frees
, copyVector(vector
, count
));
593 appendToVectorArray(ctx
->Graver
, copyVector(vector
, count
));
598 appendToVectorArray(ctx
->Homs
, copyVector(vector
, count
));
599 if (!has_symmetric
|| lex_cmp
>0)
600 appendToVectorArray(ctx
->Graver
, copyVector(vector
, count
));
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
);
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
);
621 deleteVectorArray(ctx
->Homs
);
622 deleteVectorArray(ctx
->Inhoms
);
623 deleteVectorArray(ctx
->Frees
);
624 deleteVectorArray(ctx
->Graver
);