2 * Copyright (c) 2016-2019, S. Gilles <sgilles@math.umd.edu>
4 * Permission to use, copy, modify, and/or distribute this software
5 * for any purpose with or without fee is hereby granted, provided
6 * that the above copyright notice and this permission notice appear
9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
10 * WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
11 * WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
12 * AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR
13 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
14 * OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
15 * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
16 * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
28 #define MAX_SUPPORTED_EDGE_NUM (((size_t) -1) >> 2)
29 #define MAX_SUPPORTED_VERTEX_NUM (((size_t) -1) >> 2)
31 /* Primes for reducing fractions */
32 static uint32_t primes
[] = {
34 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
35 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
36 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223,
37 227, 229, 233, 239, 241, 251
42 gcd(uint_fast32_t x
, uint_fast32_t y
)
60 /* Attempt to store a fraction in out, reducing if possible */
62 reduce_fraction(int_fast64_t n
, uint_fast64_t d
, struct rational
*out
, const
67 *out
= (struct rational
) { .p
= 0, .q
= 1 };
72 /* We special-case 1-3 since, in those cases, this actually needs to go fast */
93 for (size_t j
= 0; j
< ((sizeof primes
) / (sizeof primes
[0])); ++j
) {
94 uint32_t p
= primes
[j
];
110 if (n
> INT_FAST32_MAX
||
111 n
< INT_FAST32_MIN
||
113 d
> UINT_FAST32_MAX
) {
114 IF_NZ_SET(out_errstr
, L("unrepresentable fraction"));
119 *out
= (struct rational
) { .p
= n
, .q
= d
};
126 add_to_fraction(int_fast32_t a
, uint_fast32_t b
, struct rational
*out
, const
129 int_fast32_t n
= a
* out
->q
+ out
->p
* b
;
130 uint_fast32_t d
= b
* out
->q
;
132 return reduce_fraction(n
, d
, out
, out_errstr
);
135 /* Add an arrow of weight a/b from i -> j, affecting e(i,j) and e(j,i) */
137 quiver_add_to_edge(struct quiver
*q
, size_t i
, size_t j
, int_fast32_t a
,
138 uint_fast32_t b
, const char **out_errstr
)
143 IF_NZ_SET(out_errstr
, L("nonexistant quiver"));
150 IF_NZ_SET(out_errstr
, L("edge includes nonexistant vertex"));
156 IF_NZ_SET(out_errstr
, L("loops are not allowed"));
161 struct rational
*eij
= &(q
->e
[i
* q
->v_len
+ j
]);
162 struct rational
*eji
= &(q
->e
[j
* q
->v_len
+ i
]);
163 uint_fast32_t d
= gcd(q
->v
[i
].fatness
, q
->v
[j
].fatness
);
165 if ((ret
= add_to_fraction(a
* q
->v
[j
].fatness
, b
* d
, eij
,
170 if ((ret
= add_to_fraction(-1 * a
* q
->v
[i
].fatness
, b
* d
, eji
,
178 /* Add a vertex with a name and weight */
180 quiver_add_vertex(struct quiver
*q
, size_t *out_i
, const char *name
,
181 uint_fast16_t fatness
, int x
, int y
, uint32_t argb
, const
189 IF_NZ_SET(out_errstr
, L("invalid quiver"));
195 len
= snprintf(0, 0, "%zu", q
->v_num
);
198 IF_NZ_SET(out_errstr
, L("overflow"));
203 if (!(newname
= malloc(len
+ 1))) {
206 IF_NZ_SET(out_errstr
, L("malloc"));
211 sprintf(newname
, "%zu", q
->v_num
);
214 * Should be no need to check overflow, either name
215 * was malloc'd, or came from a file which was read
216 * in with a length cap.
218 if (!(newname
= malloc(strlen(name
) + 1))) {
221 IF_NZ_SET(out_errstr
, L("malloc"));
226 strcpy(newname
, name
);
229 if (q
->v_num
>= q
->v_len
) {
230 size_t newlen
= q
->v_len
+ 8;
232 if (newlen
>= ((size_t) -1) / newlen
||
233 (newlen
* newlen
) / newlen
!= newlen
) {
234 IF_NZ_SET(out_errstr
, L("too many vertices"));
239 if (!(newmem
= calloc(newlen
* newlen
, sizeof (*q
->e
)))) {
242 IF_NZ_SET(out_errstr
, L("too many vertices"));
247 for (size_t j
= 0; j
< q
->v_num
; ++j
) {
248 for (size_t k
= 0; k
< q
->v_num
; ++k
) {
249 ((struct rational
*) newmem
)[j
* newlen
+ k
] =
250 q
->e
[j
* q
->v_len
+ k
];
253 for (size_t k
= q
->v_num
; k
< newlen
; ++k
) {
254 ((struct rational
*) newmem
)[j
* newlen
+ k
] =
255 (struct rational
) { .p
= 0, .q
= 1 };
259 for (size_t j
= q
->v_num
; j
< newlen
; ++j
) {
260 for (size_t k
= 0; k
< newlen
; ++k
) {
261 ((struct rational
*) newmem
)[j
* newlen
+ k
] =
262 (struct rational
) { .p
= 0, .q
= 1 };
272 if (!(newmem
= (realloc(q
->v
, newlen
* sizeof (*q
->v
))))) {
275 IF_NZ_SET(out_errstr
, L("too many vertices"));
284 for (size_t k
= 0; k
<= q
->v_num
; ++k
) {
285 q
->e
[k
* q
->v_len
+ q
->v_num
] = (struct rational
) { .p
= 0, .q
=
287 q
->e
[q
->v_num
* q
->v_len
+ k
] = (struct rational
) { .p
= 0, .q
=
291 q
->v
[q
->v_num
] = (struct vertex
) { /* */
292 .name
= newname
, .fatness
= fatness
, .x
= x
, .y
= y
, .r
=
293 (argb
& 0xff0000) >> 16, .g
= (argb
& 0xff00) >> 8, .b
=
307 /* Increase or decrease the fatness of a vertex, readjusting edges as well */
309 quiver_adjust_fatness(struct quiver
*q
, size_t i
, int_fast8_t fatness_adj
, const
315 IF_NZ_SET(out_errstr
, L("invalid vertex"));
320 int new_fatness
= q
->v
[i
].fatness
+ fatness_adj
;
322 if (new_fatness
<= 0 ||
323 new_fatness
> UINT_FAST8_MAX
) {
324 IF_NZ_SET(out_errstr
, L("invalid new fatness"));
329 for (size_t j
= 0; j
< q
->v_num
; ++j
) {
334 uint_fast8_t d_orig
= gcd(q
->v
[i
].fatness
, q
->v
[j
].fatness
);
335 uint_fast8_t d_new
= gcd(new_fatness
, q
->v
[j
].fatness
);
336 size_t k
= i
* q
->v_len
+ j
;
338 if ((ret
= reduce_fraction(q
->e
[k
].p
* d_orig
, q
->e
[k
].q
*
339 d_new
, &(q
->e
[k
]), out_errstr
))) {
343 k
= j
* q
->v_len
+ i
;
345 if ((ret
= reduce_fraction(q
->e
[k
].p
* d_orig
* q
->v
[i
].fatness
,
346 q
->e
[k
].q
* d_new
* new_fatness
,
347 &(q
->e
[k
]), out_errstr
))) {
352 q
->v
[i
].fatness
= new_fatness
;
357 /* Free all memory used by this quiver, resetting it to { 0 } */
359 quiver_clean(struct quiver
*q
)
365 for (size_t j
= 0; j
< q
->v_num
; ++j
) {
371 *q
= (struct quiver
) { 0 };
376 /* Delete a vertex (and all associated edges) */
378 quiver_delete_vertex(struct quiver
*q
, size_t i
, const char **out_errstr
)
381 IF_NZ_SET(out_errstr
, L("invalid vertex"));
394 /* First, shift everything over to the left */
395 for (size_t j
= 0; j
< q
->v_num
; ++j
) {
396 for (size_t k
= i
+ 1; k
< q
->v_num
; ++k
) {
397 q
->e
[q
->v_len
* j
+ k
- 1] = q
->e
[q
->v_len
* j
+ k
];
401 /* Second, shift everything below row i up */
402 for (size_t j
= i
+ 1; j
< q
->v_num
; ++j
) {
403 for (size_t k
= 0; k
< q
->v_num
; ++k
) {
404 q
->e
[q
->v_len
* (j
- 1) + k
] = q
->e
[q
->v_len
* j
+ k
];
408 /* Now shift the vertex names/lengths over */
412 for (size_t j
= i
+ 1; j
< q
->v_num
; ++j
) {
413 q
->v
[j
- 1] = q
->v
[j
];
417 q
->v
[q
->v_num
] = (struct vertex
) { 0 };
424 quiver_get_sigma(const struct quiver
* const q
, size_t i
, size_t j
, struct
430 /* No errstr here. This is in a tight loop upstairs */
438 idx
= i
* q
->v_len
+ j
;
439 const struct rational
* const e
= &(q
->e
[idx
]);
443 out_s
->p
*= gcd(q
->v
[i
].fatness
, q
->v
[j
].fatness
);
444 out_s
->q
*= q
->v
[j
].fatness
;
445 d
= gcd(out_s
->p
, out_s
->q
);
452 /* Read quiver from file */
454 quiver_load_from_file(struct quiver
*q
, FILE *f
, const char **out_errstr
)
457 IF_NZ_SET(out_errstr
, L("invalid quiver"));
464 for (size_t k
= 0; k
< q
->v_num
; ++k
) {
465 struct vertex
*v
= &(q
->v
[k
]);
475 if (fscanf(f
, "%zu Vertices\n", &new_vnum
) != 1) {
476 IF_NZ_SET(out_errstr
, L("Invalid file (|V| unspecified)"));
481 for (size_t k
= 0; k
< new_vnum
; ++k
) {
486 unsigned int argb
= 0;
489 if (fscanf(f
, " v[%zu] f:%zu x:%d y:%d c:0x%x l:%zu n:", &j
,
490 &fatness
, &x
, &y
, &argb
, &len
) != 6) {
491 IF_NZ_SET(out_errstr
, L("Invalid file"));
495 IF_NZ_SET(out_errstr
, L(
496 "Invalid file (vertices out of order)"));
499 } else if (fatness
>= UINT_FAST16_MAX
) {
500 IF_NZ_SET(out_errstr
, L(
501 "Invalid file (fatness unsupported)"));
512 IF_NZ_SET(out_errstr
, L(
513 "Invalid file (nameless vertex)"));
519 IF_NZ_SET(out_errstr
, L(
520 "Invalid file (vertex with name too large)"));
525 if (!(p
= malloc(len
+ 1))) {
528 IF_NZ_SET(out_errstr
, L("malloc error"));
533 if (fread(p
, sizeof *p
, len
, f
) != (sizeof *p
) * len
) {
534 IF_NZ_SET(out_errstr
, L(
535 "Invalid file (short vertex name)"));
542 if ((ret
= quiver_add_vertex(q
, 0, p
, (uint_fast16_t) fatness
,
543 x
, y
, (uint32_t) argb
,
555 if (fscanf(f
, " e[%zu][%zu] %ld/%zu ", &i
, &j
, &n
, &d
) != 4) {
556 IF_NZ_SET(out_errstr
, L(
557 "Invalid file (bad edge line)"));
564 IF_NZ_SET(out_errstr
, L(
565 "Invalid file (edge between nonexistent vertices)"));
570 struct rational
*e
= &(q
->e
[i
* q
->v_len
+ j
]);
572 if (n
<= INT_FAST32_MIN
||
573 n
>= INT_FAST32_MAX
||
574 d
> UINT_FAST32_MAX
) {
575 IF_NZ_SET(out_errstr
, L(
576 "Invalid file (edge weight out of range)"));
581 if ((ret
= reduce_fraction((int_fast32_t) n
, (uint_fast32_t) d
,
594 /* Rename a vertex */
596 quiver_rename_vertex(struct quiver
*q
, size_t i
, const char *new_name
, const
603 IF_NZ_SET(out_errstr
, L("invalid quiver"));
609 IF_NZ_SET(out_errstr
, L("vertex out of range"));
614 if (!(newmem
= strdup(new_name
))) {
616 IF_NZ_SET(out_errstr
, L("cannot strdup()"));
621 q
->v
[i
].name
= newmem
;
627 /* Set color of a vertex */
629 quiver_color_vertex(struct quiver
*q
, size_t i
, uint32_t argb
, const
633 IF_NZ_SET(out_errstr
, L("invalid quiver"));
639 IF_NZ_SET(out_errstr
, L("vertex out of range"));
644 q
->v
[i
].r
= (argb
& 0x00ff0000) >> 16;
645 q
->v
[i
].g
= (argb
& 0x0000ff00) >> 8;
646 q
->v
[i
].b
= (argb
& 0x000000ff) >> 0;
651 /* Serialize the quiver */
653 quiver_save_to_file(struct quiver
*q
, FILE *f
, const char **out_errstr
)
656 IF_NZ_SET(out_errstr
, L("invalid quiver"));
663 if (fprintf(f
, "%zu Vertices\n", q
->v_num
) < 0) {
665 IF_NZ_SET(out_errstr
, L("write failed"));
669 for (size_t k
= 0; k
< q
->v_num
; ++k
) {
670 struct vertex
*v
= &(q
->v
[k
]);
671 unsigned int argb
= (v
->r
<< 16) | (v
->g
<< 8) | v
->b
;
673 if (fprintf(f
, "v[%zu] f:%zu x:%d y:%d c:0x%06x l:%zu n:%s\n",
674 k
, (size_t) v
->fatness
, v
->x
, v
->y
, argb
, strlen(
678 IF_NZ_SET(out_errstr
, L("write failed"));
683 for (size_t i
= 0; i
< q
->v_num
; ++i
) {
684 for (size_t j
= 0; j
< q
->v_num
; ++j
) {
685 struct rational
*e
= &(q
->e
[i
* q
->v_len
+ j
]);
691 if (fprintf(f
, "e[%zu][%zu] %d/%u\n", i
, j
, (int) e
->p
,
692 (unsigned) e
->q
) < 0) {
694 IF_NZ_SET(out_errstr
, L("write failed"));
702 IF_NZ_SET(out_errstr
, L("flush failed"));
710 /* Set threshold: if quiver_mutate() creates an edge above, warn_out is set */
712 quiver_set_warning_threshold(struct quiver
*q
, unsigned int warn_threshold
,
713 const char **out_errstr
)
716 IF_NZ_SET(out_errstr
, L("invalid quiver"));
721 q
->edge_weight_threshold
= warn_threshold
;
726 /* Mutate the quiver at vertex k */
728 quiver_mutate(struct quiver
*q
, size_t k
, int *out_warn
, const
734 IF_NZ_SET(out_errstr
, L("invalid quiver"));
739 /* Step one: complete all triangles */
740 for (size_t i
= 0; i
< q
->v_num
; ++i
) {
745 struct rational
*eik
= &(q
->e
[i
* q
->v_len
+ k
]);
747 for (size_t j
= 0; j
< q
->v_num
; ++j
) {
753 struct rational
*eij
= &(q
->e
[i
* q
->v_len
+ j
]);
754 struct rational
*ekj
= &(q
->e
[k
* q
->v_len
+ j
]);
756 if (eik
->p
* ekj
->p
<= 0) {
760 if ((ret
= add_to_fraction(abs(eik
->p
) * ekj
->p
,
761 eik
->q
* ekj
->q
, eij
,
767 q
->edge_weight_threshold
&&
768 eij
->q
* q
->edge_weight_threshold
<= (unsigned
776 /* Step two: invert all edges that touch k */
777 for (size_t i
= 0; i
< q
->v_num
; ++i
) {
782 struct rational
*eik
= &(q
->e
[i
* q
->v_len
+ k
]);
783 struct rational
*eki
= &(q
->e
[k
* q
->v_len
+ i
]);