quiver.c

   1 /*
   2  * Copyright (c) 2016, S. Gilles <sgilles@math.umd.edu>
   3  *
   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
   7  * in all copies.
   8  *
   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.
  17  */
  18 #include <errno.h>
  19 #include <stddef.h>
  20 #include <stdint.h>
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24
  25 #include "macros.h"
  26 #include "quiver.h"
  27
  28 #define MAX_SUPPORTED_EDGE_NUM           (((size_t) -1) >> 2)
  29 #define MAX_SUPPORTED_VERTEX_NUM         (((size_t) -1) >> 2)
  30
  31 /* Primes for reducing fractions */
  32 static uint32_t primes[] = {
  33         /* */
  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
  38 };
  39
  40 /* GCD */
  41 static int gcd(uint_fast8_t x, uint_fast8_t y)
  42 {
  43         uint_fast8_t r = 0;
  44
  45         if (!x &&
  46             !y) {
  47                 return 1;
  48         }
  49
  50         while (y) {
  51                 r = x % y;
  52                 x = y;
  53                 y = r;
  54         }
  55
  56         return x;
  57 }
  58
  59 /* Attempt to store a fraction in out, reducing if possible */
  60 static int reduce_fraction(int_fast32_t n, uint_fast32_t d, struct
  61                            rational *out, const char **out_errstr)
  62 {
  63         if (!n &&
  64             d) {
  65                 *out = (struct rational) { .p = 0, .q = 1 };
  66
  67                 return 0;
  68         }
  69
  70         /* We special-case 1-3 since, in those cases, this actually needs to go fast */
  71         if (d == 1) {
  72                 goto calculated;
  73         } else if (d == 2) {
  74                 if (n % 2) {
  75                 } else {
  76                         n /= 2;
  77                         d = 1;
  78                 }
  79
  80                 goto calculated;
  81         } else if (d == 3) {
  82                 if (n % 3) {
  83                 } else {
  84                         n /= 3;
  85                         d = 1;
  86                 }
  87
  88                 goto calculated;
  89         }
  90
  91         for (size_t j = 0; j < ((sizeof primes) / (sizeof primes[0])); ++j) {
  92                 uint32_t p = primes[j];
  93
  94                 if (j &&
  95                     p > d) {
  96                         break;
  97                 }
  98
  99                 while (n % p == 0 &&
 100                        d % p == 0) {
 101                         n = n / p;
 102                         d = d / p;
 103                 }
 104         }
 105
 106 calculated:
 107
 108         if (n > INT_FAST8_MAX ||
 109             n < INT_FAST8_MIN ||
 110             d == 0 ||
 111             d > UINT_FAST8_MAX) {
 112                 IF_NZ_SET(out_errstr, L("unrepresentable fraction"));
 113
 114                 return EDOM;
 115         }
 116
 117         *out = (struct rational) { .p = n, .q = d };
 118
 119         return 0;
 120 }
 121
 122 /* out += a/b */
 123 static int add_to_fraction(int_fast8_t a, uint_fast8_t b, struct rational *out,
 124                            const char **out_errstr)
 125 {
 126         int_fast32_t n = a * out->q + out->p * b;
 127         uint_fast32_t d = b * out->q;
 128
 129         return reduce_fraction(n, d, out, out_errstr);
 130 }
 131
 132 /* Add an arrow of weight a/b from i -> j, affecting e(i,j) and e(j,i) */
 133 int quiver_add_to_edge(struct quiver *q, size_t i, size_t j, int_fast8_t a,
 134                        uint_fast8_t b, const char **out_errstr)
 135 {
 136         int ret = 0;
 137
 138         if (!q) {
 139                 IF_NZ_SET(out_errstr, L("nonexistant quiver"));
 140
 141                 return EINVAL;
 142         }
 143
 144         if (i >= q->v_num ||
 145             j >= q->v_num) {
 146                 IF_NZ_SET(out_errstr, L("edge includes nonexistant vertex"));
 147
 148                 return EINVAL;
 149         }
 150
 151         if (i == j) {
 152                 IF_NZ_SET(out_errstr, L("loops are not allowed"));
 153
 154                 return EINVAL;
 155         }
 156
 157         struct rational *eij = &(q->e[i * q->v_len + j]);
 158         struct rational *eji = &(q->e[j * q->v_len + i]);
 159         uint_fast8_t d = gcd(q->v[i].fatness, q->v[j].fatness);
 160
 161         if ((ret = add_to_fraction(a * q->v[j].fatness, b * d, eij,
 162                                    out_errstr))) {
 163                 return ret;
 164         }
 165
 166         if ((ret = add_to_fraction(-1 * a * q->v[i].fatness, b * d, eji,
 167                                    out_errstr))) {
 168                 return ret;
 169         }
 170
 171         return 0;
 172 }
 173
 174 /* Add a vertex with a name and weight */
 175 int quiver_add_vertex(struct quiver *q, size_t *out_i, const char *name,
 176                       uint_fast16_t fatness, int x, int y, const
 177                       char **out_errstr)
 178 {
 179         void *newmem = 0;
 180         char *newname;
 181
 182         if (!q) {
 183                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 184
 185                 return EINVAL;
 186         }
 187
 188         if (!name) {
 189                 if (!(newname = malloc(snprintf(0, 0, "%zu", q->v_num) + 1))) {
 190                         int sv_err = errno;
 191
 192                         IF_NZ_SET(out_errstr, L("malloc"));
 193
 194                         return sv_err;
 195                 }
 196
 197                 sprintf(newname, "%zu", q->v_num);
 198         } else {
 199                 if (!(newname = malloc(strlen(name) + 1))) {
 200                         int sv_err = errno;
 201
 202                         IF_NZ_SET(out_errstr, L("malloc"));
 203
 204                         return sv_err;
 205                 }
 206
 207                 strcpy(newname, name);
 208         }
 209
 210         if (q->v_num >= q->v_len) {
 211                 size_t newlen = q->v_len + 8;
 212
 213                 /* XXX: check for overflow here */
 214                 if (!(newmem = malloc(newlen * newlen * sizeof (*q->e)))) {
 215                         int sv_err = errno;
 216
 217                         IF_NZ_SET(out_errstr, L("too many vertices"));
 218
 219                         return sv_err;
 220                 }
 221
 222                 for (size_t j = 0; j < q->v_num; ++j) {
 223                         for (size_t k = 0; k < q->v_num; ++k) {
 224                                 ((struct rational *) newmem)[j * newlen + k] =
 225                                         q->e[j * q->v_len + k];
 226                         }
 227
 228                         for (size_t k = q->v_num; k < newlen; ++k) {
 229                                 ((struct rational *) newmem)[j * newlen + k] =
 230                                         (struct rational) { .p = 0, .q = 1 };
 231                         }
 232                 }
 233
 234                 for (size_t j = q->v_num; j < newlen; ++j) {
 235                         for (size_t k = 0; k < newlen; ++k) {
 236                                 ((struct rational *) newmem)[j * newlen + k] =
 237                                         (struct rational) { .p = 0, .q = 1 };
 238                         }
 239                 }
 240
 241                 q->e = newmem;
 242
 243                 if (!(newmem = (realloc(q->v, newlen * sizeof (*q->v))))) {
 244                         int sv_err = errno;
 245
 246                         IF_NZ_SET(out_errstr, L("too many vertices"));
 247
 248                         return sv_err;
 249                 }
 250
 251                 q->v = newmem;
 252                 q->v_len = newlen;
 253         }
 254
 255         for (size_t k = 0; k <= q->v_num; ++k) {
 256                 q->e[k * q->v_len + q->v_num] = (struct rational) { .p = 0, .q =
 257                                                                             1 };
 258                 q->e[q->v_num * q->v_len + k] = (struct rational) { .p = 0, .q =
 259                                                                             1 };
 260         }
 261
 262         q->v[q->v_num] = (struct vertex) { .name = newname, .fatness = fatness,
 263                                            .x = x, .y = y };
 264
 265         if (out_i) {
 266                 *out_i = q->v_num;
 267         }
 268
 269         q->v_num++;
 270
 271         return 0;
 272 }
 273
 274 /* Increase or decrease the fatness of a vertex, readjusting edges as well */
 275 int quiver_adjust_fatness(struct quiver *q, size_t i, int_fast8_t fatness_adj,
 276                           const char **out_errstr)
 277 {
 278         int ret = 0;
 279
 280         if (i >= q->v_num) {
 281                 IF_NZ_SET(out_errstr, L("invalid vertex"));
 282
 283                 return EINVAL;
 284         }
 285
 286         int new_fatness = q->v[i].fatness + fatness_adj;
 287
 288         if (new_fatness <= 0 ||
 289             new_fatness > UINT_FAST8_MAX) {
 290                 IF_NZ_SET(out_errstr, L("invalid new fatness"));
 291
 292                 return EINVAL;
 293         }
 294
 295         for (size_t j = 0; j < q->v_num; ++j) {
 296                 if (i == j) {
 297                         continue;
 298                 }
 299
 300                 uint_fast8_t d_orig = gcd(q->v[i].fatness, q->v[j].fatness);
 301                 uint_fast8_t d_new = gcd(new_fatness, q->v[j].fatness);
 302                 size_t k = i * q->v_len + j;
 303
 304                 if ((ret = reduce_fraction(q->e[k].p * d_orig, q->e[k].q *
 305                                            d_new, &(q->e[k]), out_errstr))) {
 306                         return ret;
 307                 }
 308
 309                 k = j * q->v_len + i;
 310
 311                 if ((ret = reduce_fraction(q->e[k].p * d_orig * q->v[i].fatness,
 312                                            q->e[k].q * d_new * new_fatness,
 313                                            &(q->e[k]), out_errstr))) {
 314                         return ret;
 315                 }
 316         }
 317
 318         q->v[i].fatness = new_fatness;
 319
 320         return ret;
 321 }
 322
 323 /* Delete a vertex (and all associated edges) */
 324 int quiver_delete_vertex(struct quiver *q, size_t i, const char **out_errstr)
 325 {
 326         if (i >= q->v_num) {
 327                 IF_NZ_SET(out_errstr, L("invalid vertex"));
 328
 329                 return EINVAL;
 330         }
 331
 332         if (q->v_num == 1) {
 333                 q->v_num = 0;
 334                 free(q->v[0].name);
 335
 336                 return 0;
 337         }
 338
 339         /* First, shift everything over to the left */
 340         for (size_t j = 0; j < q->v_num; ++j) {
 341                 for (size_t k = i + 1; k < q->v_num; ++k) {
 342                         q->e[q->v_len * j + k - 1] = q->e[q->v_len * j + k];
 343                 }
 344         }
 345
 346         /* Second, shift everything below row i up */
 347         for (size_t j = i + 1; j < q->v_num; ++j) {
 348                 for (size_t k = 0; k < q->v_num; ++k) {
 349                         q->e[q->v_len * (j - 1) + k] = q->e[q->v_len * j + k];
 350                 }
 351         }
 352
 353         /* Now shift the vertex names/lengths over */
 354         free(q->v[i].name);
 355
 356         for (size_t j = i + 1; j < q->v_num; ++j) {
 357                 q->v[j - 1] = q->v[j];
 358         }
 359
 360         q->v_num--;
 361         q->v[q->v_num] = (struct vertex) { 0 };
 362
 363         return 0;
 364 }
 365
 366 /* Read quiver from file */
 367 int quiver_load_from_file(struct quiver *q, FILE *f, const char **out_errstr)
 368 {
 369         if (!q) {
 370                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 371
 372                 return EINVAL;
 373         }
 374
 375         int ret = 0;
 376
 377         for (size_t k = 0; k < q->v_num; ++k) {
 378                 struct vertex *v = &(q->v[k]);
 379
 380                 free(v->name);
 381                 v->name = 0;
 382         }
 383
 384         q->v_num = 0;
 385         size_t new_vnum = 0;
 386         char *p = 0;
 387
 388         if (fscanf(f, "%zu Vertices\n", &new_vnum) != 1) {
 389                 IF_NZ_SET(out_errstr, L("Invalid file (|V| unspecified)"));
 390
 391                 return EINVAL;
 392         }
 393
 394         for (size_t k = 0; k < new_vnum; ++k) {
 395                 size_t j = 0;
 396                 size_t fatness = 0;
 397                 int x = 0;
 398                 int y = 0;
 399                 size_t len = 0;
 400
 401                 if (fscanf(f, " v[%zu] f:%zu x:%d y:%d l:%zu n:", &j, &fatness,
 402                            &x, &y, &len) != 5) {
 403                         IF_NZ_SET(out_errstr, L("Invalid file"));
 404
 405                         return EINVAL;
 406                 } else if (j != k) {
 407                         IF_NZ_SET(out_errstr, L(
 408                                           "Invalid file (vertices out of order)"));
 409
 410                         return EINVAL;
 411                 } else if (fatness >= UINT_FAST16_MAX) {
 412                         IF_NZ_SET(out_errstr, L(
 413                                           "Invalid file (fatness unsupported)"));
 414
 415                         return EINVAL;
 416                 }
 417
 418                 if (p) {
 419                         free(p);
 420                         p = 0;
 421                 }
 422
 423                 if (!(p = malloc(len + 1))) {
 424                         ret = errno;
 425                         perror(L("malloc"));
 426                         IF_NZ_SET(out_errstr, L("malloc error"));
 427
 428                         return ret;
 429                 }
 430
 431                 if (fread(p, sizeof *p, len, f) != (sizeof *p) * len) {
 432                         IF_NZ_SET(out_errstr, L(
 433                                           "Invalid file (short vertex name)"));
 434
 435                         return EINVAL;
 436                 }
 437
 438                 p[len] = '\0';
 439
 440                 if ((ret = quiver_add_vertex(q, 0, p, (uint_fast16_t) fatness,
 441                                              x, y, out_errstr))) {
 442                         goto done;
 443                 }
 444         }
 445
 446         while (!feof(f)) {
 447                 size_t i = 0;
 448                 size_t j = 0;
 449                 long n = 0;
 450                 size_t d = 0;
 451
 452                 if (fscanf(f, " e[%zu][%zu] %ld/%zu ", &i, &j, &n, &d) != 4) {
 453                         IF_NZ_SET(out_errstr, L("Invalid file (bad edge)"));
 454                         ret = EINVAL;
 455                         goto done;
 456                 }
 457
 458                 if (i >= q->v_num ||
 459                     j >= q->v_num) {
 460                         IF_NZ_SET(out_errstr, L(
 461                                           "Invalid file (edge between nonexistent vertices)"));
 462                         ret = EINVAL;
 463                         goto done;
 464                 }
 465
 466                 struct rational *e = &(q->e[i * q->v_len + j]);
 467
 468                 if (n <= INT_FAST32_MIN ||
 469                     n >= INT_FAST32_MAX ||
 470                     d > UINT_FAST32_MAX) {
 471                         IF_NZ_SET(out_errstr, L(
 472                                           "Invalid file (edge weight out of range)"));
 473                         ret = EINVAL;
 474                         goto done;
 475                 }
 476
 477                 if ((ret = reduce_fraction((int_fast32_t) n, (uint_fast32_t) d,
 478                                            e, out_errstr))) {
 479                         goto done;
 480                 }
 481         }
 482
 483         ret = 0;
 484 done:
 485         free(p);
 486
 487         return ret;
 488 }
 489
 490 /* Rename a vertex */
 491 int quiver_rename_vertex(struct quiver *q, size_t i, const char *new_name, const
 492                          char **out_errstr)
 493 {
 494         if (!q) {
 495                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 496
 497                 return EINVAL;
 498         }
 499
 500         if (i >= q->v_num) {
 501                 IF_NZ_SET(out_errstr, L("vertex out of range"));
 502
 503                 return EINVAL;
 504         }
 505
 506         void *newmem = 0;
 507         int ret = 0;
 508
 509         if (!(newmem = strdup(new_name))) {
 510                 ret = errno;
 511                 IF_NZ_SET(out_errstr, L("cannot strdup()"));
 512                 goto done;
 513         }
 514
 515         free(q->v[i].name);
 516         q->v[i].name = newmem;
 517 done:
 518
 519         return ret;
 520 }
 521
 522 /* Serialize the quiver */
 523 int quiver_save_to_file(struct quiver *q, FILE *f, const char **out_errstr)
 524 {
 525         if (!q) {
 526                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 527
 528                 return EINVAL;
 529         }
 530
 531         int ret = 0;
 532
 533         if (fprintf(f, "%zu Vertices\n", q->v_num) < 0) {
 534                 ret = errno;
 535                 IF_NZ_SET(out_errstr, L("write failed"));
 536                 goto done;
 537         }
 538
 539         for (size_t k = 0; k < q->v_num; ++k) {
 540                 struct vertex *v = &(q->v[k]);
 541
 542                 if (fprintf(f, "v[%zu] f:%zu x:%d y:%d l:%zu n:%s\n", k,
 543                             (size_t) v->fatness, v->x, v->y, strlen(v->name),
 544                             v->name) <
 545                     0) {
 546                         ret = errno;
 547                         IF_NZ_SET(out_errstr, L("write failed"));
 548                         goto done;
 549                 }
 550         }
 551
 552         for (size_t i = 0; i < q->v_num; ++i) {
 553                 for (size_t j = 0; j < q->v_num; ++j) {
 554                         struct rational *e = &(q->e[i * q->v_len + j]);
 555
 556                         if (!e->p) {
 557                                 continue;
 558                         }
 559
 560                         if (fprintf(f, "e[%zu][%zu] %d/%u\n", i, j, (int) e->p,
 561                                     (unsigned) e->q) < 0) {
 562                                 ret = errno;
 563                                 IF_NZ_SET(out_errstr, L("write failed"));
 564                                 goto done;
 565                         }
 566                 }
 567         }
 568
 569         if (fflush(f)) {
 570                 ret = errno;
 571                 IF_NZ_SET(out_errstr, L("flush failed"));
 572         }
 573
 574 done:
 575
 576         return ret;
 577 }
 578
 579 /* Set threshold: if quiver_mutate() creates an edge above, warn_out is set */
 580 int quiver_set_warning_threshold(struct quiver *q, unsigned int warn_threshold,
 581                                  char **out_errstr)
 582 {
 583         if (!q) {
 584                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 585
 586                 return EINVAL;
 587         }
 588
 589         q->edge_weight_threshold = warn_threshold;
 590
 591         return 0;
 592 }
 593
 594 /* Mutate the quiver at vertex k */
 595 int quiver_mutate(struct quiver *q, size_t k, int *out_warn, const
 596                   char **out_errstr)
 597 {
 598         int ret = 0;
 599
 600         if (!q) {
 601                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 602
 603                 return EINVAL;
 604         }
 605
 606         /* Step one: complete all triangles */
 607         for (size_t i = 0; i < q->v_num; ++i) {
 608                 if (i == k) {
 609                         continue;
 610                 }
 611
 612                 struct rational *eik = &(q->e[i * q->v_len + k]);
 613
 614                 for (size_t j = 0; j < q->v_num; ++j) {
 615                         if (j == k ||
 616                             j == i) {
 617                                 continue;
 618                         }
 619
 620                         struct rational *eij = &(q->e[i * q->v_len + j]);
 621                         struct rational *ekj = &(q->e[k * q->v_len + j]);
 622
 623                         if (eik->p * ekj->p <= 0) {
 624                                 continue;
 625                         }
 626
 627                         if ((ret = add_to_fraction(abs(eik->p) * ekj->p,
 628                                                    eik->q * ekj->q, eij,
 629                                                    out_errstr))) {
 630                                 return ret;
 631                         }
 632
 633                         if (out_warn &&
 634                             q->edge_weight_threshold &&
 635                             eij->q * q->edge_weight_threshold <= (unsigned
 636                                                                   int) abs(
 637                                     eij->p)) {
 638                                 *out_warn = 1;
 639                         }
 640                 }
 641         }
 642
 643         /* Step two: invert all edges that touch k */
 644         for (size_t i = 0; i < q->v_num; ++i) {
 645                 if (i == k) {
 646                         continue;
 647                 }
 648
 649                 struct rational *eik = &(q->e[i * q->v_len + k]);
 650                 struct rational *eki = &(q->e[k * q->v_len + i]);
 651
 652                 eik->p = -eik->p;
 653                 eki->p = -eki->p;
 654         }
 655
 656         return 0;
 657 }