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 /* Set threshold: if quiver_mutate() creates an edge above, warn_out is set */
 367 int quiver_set_warning_threshold(struct quiver *q, unsigned int warn_threshold,
 368                                  char **out_errstr)
 369 {
 370         if (!q) {
 371                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 372
 373                 return EINVAL;
 374         }
 375
 376         q->edge_weight_threshold = warn_threshold;
 377
 378         return 0;
 379 }
 380
 381 /* Mutate the quiver at vertex k */
 382 int quiver_mutate(struct quiver *q, size_t k, int *out_warn, const
 383                   char **out_errstr)
 384 {
 385         int ret = 0;
 386
 387         if (!q) {
 388                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 389
 390                 return EINVAL;
 391         }
 392
 393         /* Step one: complete all triangles */
 394         for (size_t i = 0; i < q->v_num; ++i) {
 395                 if (i == k) {
 396                         continue;
 397                 }
 398
 399                 struct rational *eik = &(q->e[i * q->v_len + k]);
 400
 401                 for (size_t j = 0; j < q->v_num; ++j) {
 402                         if (j == k ||
 403                             j == i) {
 404                                 continue;
 405                         }
 406
 407                         struct rational *eij = &(q->e[i * q->v_len + j]);
 408                         struct rational *ekj = &(q->e[k * q->v_len + j]);
 409
 410                         if (eik->p * ekj->p <= 0) {
 411                                 continue;
 412                         }
 413
 414                         if ((ret = add_to_fraction(abs(eik->p) * ekj->p,
 415                                                    eik->q * ekj->q, eij,
 416                                                    out_errstr))) {
 417                                 return ret;
 418                         }
 419
 420                         if (out_warn &&
 421                             q->edge_weight_threshold &&
 422                             eij->q * q->edge_weight_threshold <= (unsigned int)abs(eij->p)) {
 423                                 *out_warn = 1;
 424                         }
 425                 }
 426         }
 427
 428         /* Step two: invert all edges that touch k */
 429         for (size_t i = 0; i < q->v_num; ++i) {
 430                 if (i == k) {
 431                         continue;
 432                 }
 433
 434                 struct rational *eik = &(q->e[i * q->v_len + k]);
 435                 struct rational *eki = &(q->e[k * q->v_len + i]);
 436
 437                 eik->p = -eik->p;
 438                 eki->p = -eki->p;
 439         }
 440
 441         return 0;
 442 }