quiver.c

   1 #include <errno.h>
   2 #include <stddef.h>
   3 #include <stdint.h>
   4 #include <stdio.h>
   5 #include <stdlib.h>
   6 #include <string.h>
   7
   8 #include "macros.h"
   9 #include "quiver.h"
  10
  11 #define MAX_SUPPORTED_EDGE_NUM           (((size_t) -1) >> 2)
  12 #define MAX_SUPPORTED_VERTEX_NUM         (((size_t) -1) >> 2)
  13
  14 /* Primes for reducing fractions */
  15 static uint32_t primes[] = {
  16         /* */
  17         2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
  18         71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
  19         149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223,
  20         227, 229, 233, 239, 241, 251
  21 };
  22
  23 /* GCD */
  24 static int gcd(uint_fast8_t x, uint_fast8_t y)
  25 {
  26         uint_fast8_t r = 0;
  27
  28         if (!x &&
  29             !y) {
  30                 return 1;
  31         }
  32
  33         while (y) {
  34                 r = x % y;
  35                 x = y;
  36                 y = r;
  37         }
  38
  39         return x;
  40 }
  41
  42 /* Attempt to store a fraction in out, reducing if possible */
  43 static int reduce_fraction(int_fast32_t n, uint_fast32_t d, struct
  44                            rational *out, const char **out_errstr)
  45 {
  46         if (!n &&
  47             d) {
  48                 *out = (struct rational) { .p = 0, .q = 1 };
  49
  50                 return 0;
  51         }
  52
  53         for (size_t j = 0; j < ((sizeof primes) / (sizeof primes[0])); ++j) {
  54                 uint32_t p = primes[j];
  55
  56                 if (j &&
  57                     p * 2 > d) {
  58                         break;
  59                 }
  60
  61                 while (n % p == 0 &&
  62                        d % p == 0) {
  63                         n = n / p;
  64                         d = d / p;
  65                 }
  66         }
  67
  68         if (n > INT_FAST8_MAX ||
  69             n < INT_FAST8_MIN ||
  70             d == 0 ||
  71             d > UINT_FAST8_MAX) {
  72                 IF_NZ_SET(out_errstr, L("unrepresentable fraction"));
  73
  74                 return EDOM;
  75         }
  76
  77         *out = (struct rational) { .p = n, .q = d };
  78
  79         return 0;
  80 }
  81
  82 /* out += a/b */
  83 static int add_to_fraction(int_fast8_t a, uint_fast8_t b, struct rational *out,
  84                            const char **out_errstr)
  85 {
  86         int_fast32_t n = a * out->q + out->p * b;
  87         uint_fast32_t d = b * out->q;
  88
  89         return reduce_fraction(n, d, out, out_errstr);
  90 }
  91
  92 /* Add an arrow of weight a/b from i -> j, affecting e(i,j) and e(j,i) */
  93 int quiver_add_to_edge(struct quiver *q, size_t i, size_t j, int_fast8_t a,
  94                        uint_fast8_t b, const char **out_errstr)
  95 {
  96         int ret = 0;
  97
  98         if (!q) {
  99                 IF_NZ_SET(out_errstr, L("nonexistant quiver"));
 100
 101                 return EINVAL;
 102         }
 103
 104         if (i >= q->v_num ||
 105             j >= q->v_num) {
 106                 IF_NZ_SET(out_errstr, L("edge includes nonexistant vertex"));
 107
 108                 return EINVAL;
 109         }
 110
 111         if (i == j) {
 112                 IF_NZ_SET(out_errstr, L("loops are not allowed"));
 113
 114                 return EINVAL;
 115         }
 116
 117         struct rational *eij = &(q->e[i * q->v_len + j]);
 118         struct rational *eji = &(q->e[j * q->v_len + i]);
 119         uint_fast8_t d = gcd(q->v[i].fatness, q->v[j].fatness);
 120
 121         if ((ret = add_to_fraction(a * q->v[j].fatness, b * d, eij,
 122                                    out_errstr))) {
 123                 return ret;
 124         }
 125
 126         if ((ret = add_to_fraction(-1 * a * q->v[i].fatness, b * d, eji,
 127                                    out_errstr))) {
 128                 return ret;
 129         }
 130
 131         return 0;
 132 }
 133
 134 /* Add a vertex with a name and weight */
 135 int quiver_add_vertex(struct quiver *q, size_t *out_i, const char *name,
 136                       uint_fast16_t fatness, int x, int y, const
 137                       char **out_errstr)
 138 {
 139         void *newmem = 0;
 140         char *newname;
 141
 142         if (!q) {
 143                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 144
 145                 return EINVAL;
 146         }
 147
 148         if (!name) {
 149                 if (!(newname = malloc(snprintf(0, 0, "%zu", q->v_num) + 1))) {
 150                         int sv_err = errno;
 151
 152                         IF_NZ_SET(out_errstr, L("malloc"));
 153
 154                         return sv_err;
 155                 }
 156
 157                 sprintf(newname, "%zu", q->v_num);
 158         } else {
 159                 if (!(newname = malloc(strlen(name) + 1))) {
 160                         int sv_err = errno;
 161
 162                         IF_NZ_SET(out_errstr, L("malloc"));
 163
 164                         return sv_err;
 165                 }
 166
 167                 strcpy(newname, name);
 168         }
 169
 170         if (q->v_num >= q->v_len) {
 171                 size_t newlen = q->v_len + 8;
 172
 173                 /* XXX: check for overflow here */
 174                 if (!(newmem = malloc(newlen * newlen * sizeof (*q->e)))) {
 175                         int sv_err = errno;
 176
 177                         IF_NZ_SET(out_errstr, L("too many vertices"));
 178
 179                         return sv_err;
 180                 }
 181
 182                 for (size_t j = 0; j < q->v_num; ++j) {
 183                         for (size_t k = 0; k < q->v_num; ++k) {
 184                                 ((struct rational *) newmem)[j * newlen + k] =
 185                                         q->e[j * q->v_len + k];
 186                         }
 187
 188                         for (size_t k = q->v_num; k < newlen; ++k) {
 189                                 ((struct rational *) newmem)[j * newlen + k] =
 190                                         (struct rational) { .p = 0, .q = 1 };
 191                         }
 192                 }
 193
 194                 for (size_t j = q->v_num; j < newlen; ++j) {
 195                         for (size_t k = 0; k < newlen; ++k) {
 196                                 ((struct rational *) newmem)[j * newlen + k] =
 197                                         (struct rational) { .p = 0, .q = 1 };
 198                         }
 199                 }
 200
 201                 q->e = newmem;
 202
 203                 if (!(newmem = (realloc(q->v, newlen * sizeof (*q->v))))) {
 204                         int sv_err = errno;
 205
 206                         IF_NZ_SET(out_errstr, L("too many vertices"));
 207
 208                         return sv_err;
 209                 }
 210
 211                 q->v = newmem;
 212                 q->v_len = newlen;
 213         }
 214
 215         for (size_t k = 0; k <= q->v_num; ++k) {
 216                 q->e[k * q->v_len + q->v_num] = (struct rational) { .p = 0, .q =
 217                                                                             1 };
 218                 q->e[q->v_num * q->v_len + k] = (struct rational) { .p = 0, .q =
 219                                                                             1 };
 220         }
 221
 222         q->v[q->v_num] = (struct vertex) { .name = newname, .fatness = fatness,
 223                                            .x = x, .y = y };
 224
 225         if (out_i) {
 226                 *out_i = q->v_num;
 227         }
 228
 229         q->v_num++;
 230
 231         return 0;
 232 }
 233
 234 /* Increase or decrease the fatness of a vertex, readjusting edges as well */
 235 int quiver_adjust_fatness(struct quiver *q, size_t i, int_fast8_t fatness_adj,
 236                           const char **out_errstr)
 237 {
 238         int ret = 0;
 239
 240         if (i >= q->v_num) {
 241                 IF_NZ_SET(out_errstr, L("invalid vertex"));
 242
 243                 return EINVAL;
 244         }
 245
 246         int new_fatness = q->v[i].fatness + fatness_adj;
 247
 248         if (new_fatness <= 0 ||
 249             new_fatness > UINT_FAST8_MAX) {
 250                 IF_NZ_SET(out_errstr, L("invalid new fatness"));
 251
 252                 return EINVAL;
 253         }
 254
 255         for (size_t j = 0; j < q->v_num; ++j) {
 256                 if (i == j) {
 257                         continue;
 258                 }
 259
 260                 uint_fast8_t d_orig = gcd(q->v[i].fatness, q->v[j].fatness);
 261                 uint_fast8_t d_new = gcd(new_fatness, q->v[j].fatness);
 262                 size_t k = i * q->v_len + j;
 263
 264                 if ((ret = reduce_fraction(q->e[k].p * d_orig, q->e[k].q *
 265                                            d_new, &(q->e[k]), out_errstr))) {
 266                         return ret;
 267                 }
 268
 269                 k = j * q->v_len + i;
 270
 271                 if ((ret = reduce_fraction(q->e[k].p * d_orig * q->v[i].fatness,
 272                                            q->e[k].q * d_new * new_fatness,
 273                                            &(q->e[k]), out_errstr))) {
 274                         return ret;
 275                 }
 276         }
 277
 278         q->v[i].fatness = new_fatness;
 279
 280         return ret;
 281 }
 282
 283 /* Delete a vertex (and all associated edges) */
 284 int quiver_delete_vertex(struct quiver *q, size_t i, const char **out_errstr)
 285 {
 286         if (i >= q->v_num) {
 287                 IF_NZ_SET(out_errstr, L("invalid vertex"));
 288
 289                 return EINVAL;
 290         }
 291
 292         if (q->v_num == 1) {
 293                 q->v_num = 0;
 294                 free(q->v[0].name);
 295
 296                 return 0;
 297         }
 298
 299         /* First, shift everything over to the left */
 300         for (size_t j = 0; j < q->v_num; ++j) {
 301                 for (size_t k = i + 1; k < q->v_num; ++k) {
 302                         q->e[q->v_len * j + k - 1] = q->e[q->v_len * j + k];
 303                 }
 304         }
 305
 306         /* Second, shift everything below row i up */
 307         for (size_t j = i + 1; j < q->v_num; ++j) {
 308                 for (size_t k = 0; k < q->v_num; ++k) {
 309                         q->e[q->v_len * (j - 1) + k] = q->e[q->v_len * j + k];
 310                 }
 311         }
 312
 313         /* Now shift the vertex names/lengths over */
 314         free(q->v[i].name);
 315
 316         for (size_t j = i + 1; j < q->v_num; ++j) {
 317                 q->v[j - 1] = q->v[j];
 318         }
 319
 320         q->v_num--;
 321         q->v[q->v_num] = (struct vertex) { 0 };
 322
 323         return 0;
 324 }
 325
 326 /* Mutate the quiver at vertex k */
 327 int quiver_mutate(struct quiver *q, size_t k, const char **out_errstr)
 328 {
 329         int ret = 0;
 330
 331         if (!q) {
 332                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 333
 334                 return EINVAL;
 335         }
 336
 337         /* Step one: complete all triangles */
 338         for (size_t i = 0; i < q->v_num; ++i) {
 339                 if (i == k) {
 340                         continue;
 341                 }
 342
 343                 struct rational *eik = &(q->e[i * q->v_len + k]);
 344
 345                 for (size_t j = 0; j < q->v_num; ++j) {
 346                         if (j == k ||
 347                             j == i) {
 348                                 continue;
 349                         }
 350
 351                         struct rational *eij = &(q->e[i * q->v_len + j]);
 352                         struct rational *ekj = &(q->e[k * q->v_len + j]);
 353
 354                         if (eik->p * ekj->p <= 0) {
 355                                 continue;
 356                         }
 357
 358                         if ((ret = add_to_fraction(abs(eik->p) * ekj->p,
 359                                                    eik->q * ekj->q, eij,
 360                                                    out_errstr))) {
 361                                 return ret;
 362                         }
 363                 }
 364         }
 365
 366         /* Step two: invert all edges that touch k */
 367         for (size_t i = 0; i < q->v_num; ++i) {
 368                 if (i == k) {
 369                         continue;
 370                 }
 371
 372                 struct rational *eik = &(q->e[i * q->v_len + k]);
 373                 struct rational *eki = &(q->e[k * q->v_len + i]);
 374
 375                 eik->p = -eik->p;
 376                 eki->p = -eki->p;
 377         }
 378
 379         return 0;
 380 }