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 /* Attempt to store a fraction in out, reducing if possible */
  24 static int reduce_fraction(int_fast32_t n, uint_fast32_t d, struct
  25                            rational *out, const char **out_errstr)
  26 {
  27         size_t j = 0;
  28         uint32_t p = 0;
  29
  30         for (j = 0; j < ((sizeof primes) / (sizeof primes[0])); ++j) {
  31                 p = primes[j];
  32
  33                 if (j &&
  34                     p * 2 > d) {
  35                         break;
  36                 }
  37
  38                 while (n % p == 0 &&
  39                        d % p == 0) {
  40                         n = n / p;
  41                         d = d / p;
  42                 }
  43         }
  44
  45         if (n > INT_FAST8_MAX ||
  46             n < INT_FAST8_MIN ||
  47             d > UINT_FAST8_MAX) {
  48                 IF_NZ_SET(out_errstr, L("unrepresentable fraction"));
  49
  50                 return EDOM;
  51         }
  52
  53         *out = (struct rational) { .p = n, .q = d };
  54
  55         return 0;
  56 }
  57
  58 /* out += a/b */
  59 static int add_to_fraction(int_fast8_t a, uint_fast8_t b, struct rational *out,
  60                            const char **out_errstr)
  61 {
  62         int_fast32_t n = a * out->q + out->p * b;
  63         uint_fast32_t d = b * out->q;
  64
  65         return reduce_fraction(n, d, out, out_errstr);
  66 }
  67
  68 /* Add an arrow of weight a/b from i -> j, affecting e(i,j) and e(j,i) */
  69 int quiver_add_to_edge(struct quiver *q, size_t i, size_t j, int_fast8_t a,
  70                        uint_fast8_t b, const char **out_errstr)
  71 {
  72         int ret = 0;
  73         struct rational *eij = &(q->e[i * q->v_len + j]);
  74         struct rational *eji = &(q->e[j * q->v_len + i]);
  75
  76         if (!q) {
  77                 IF_NZ_SET(out_errstr, L("nonexistant quiver"));
  78
  79                 return EINVAL;
  80         }
  81
  82         if (i >= q->v_num ||
  83             j >= q->v_num) {
  84                 IF_NZ_SET(out_errstr, L("edge includes nonexistant vertex"));
  85
  86                 return EINVAL;
  87         }
  88
  89         if ((ret = add_to_fraction(a * q->v[j].fatness, b, eij, out_errstr))) {
  90                 goto done;
  91         }
  92
  93         if ((ret = add_to_fraction(-1 * a * q->v[i].fatness, b, eji,
  94                                    out_errstr))) {
  95                 goto done;
  96         }
  97
  98 done:
  99
 100         return ret;
 101 }
 102
 103 /* Add a vertex with a name and weight */
 104 int quiver_add_vertex(struct quiver *q, size_t *out_i, const char *name,
 105                       uint_fast16_t fatness, int x, int y, const
 106                       char **out_errstr)
 107 {
 108         void *newmem = 0;
 109         int sv_err = 0;
 110         size_t l = strlen(name);
 111         char *newname;
 112         size_t j;
 113         size_t k;
 114         size_t newlen;
 115
 116         if (!q) {
 117                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 118
 119                 return EINVAL;
 120         }
 121
 122         if (!(newname = malloc(l + 1))) {
 123                 sv_err = errno;
 124                 IF_NZ_SET(out_errstr, L("malloc"));
 125
 126                 return sv_err;
 127         }
 128
 129         strcpy(newname, name);
 130
 131         if (q->v_num >= q->v_len) {
 132                 newlen = q->v_len + 8;
 133
 134                 /* XXX: check for overflow here */
 135                 if (!(newmem = malloc(newlen * newlen * sizeof (*q->e)))) {
 136                         sv_err = errno;
 137                         IF_NZ_SET(out_errstr, L("too many vertices"));
 138
 139                         return sv_err;
 140                 }
 141
 142                 for (j = 0; j < q->v_num; ++j) {
 143                         for (k = 0; k < q->v_num; ++k) {
 144                                 ((struct rational *) newmem)[j * newlen + k] =
 145                                         q->e[j * q->v_len + k];
 146                         }
 147
 148                         for (k = q->v_num; k < newlen; ++k) {
 149                                 ((struct rational *) newmem)[j * newlen + k] =
 150                                         (struct rational) { .p = 0, .q = 1 };
 151                         }
 152                 }
 153
 154                 for (j = q->v_num; j < newlen; ++j) {
 155                         for (k = 0; k < newlen; ++k) {
 156                                 ((struct rational *) newmem)[j * newlen + k] =
 157                                         (struct rational) { .p = 0, .q = 1 };
 158                         }
 159                 }
 160
 161                 q->e = newmem;
 162
 163                 if (!(newmem = (realloc(q->v, newlen * sizeof (*q->v))))) {
 164                         sv_err = errno;
 165                         IF_NZ_SET(out_errstr, L("too many vertices"));
 166
 167                         return sv_err;
 168                 }
 169
 170                 q->v = newmem;
 171                 q->v_len = newlen;
 172         }
 173
 174         for (k = 0; k <= q->v_num; ++k) {
 175                 q->e[k * q->v_len + q->v_num] = (struct rational) { .p = 0, .q =
 176                                                                             1 };
 177                 q->e[q->v_num * q->v_len + k] = (struct rational) { .p = 0, .q =
 178                                                                             1 };
 179         }
 180
 181         q->v[q->v_num] = (struct vertex) { .name = newname, .fatness = fatness,
 182                                            .x = x, .y = y };
 183         *out_i = q->v_num;
 184         q->v_num++;
 185
 186         return 0;
 187 }
 188
 189 /* Mutate the quiver at vertex k */
 190 int quiver_mutate(struct quiver *q, size_t k, const char **out_errstr)
 191 {
 192         size_t i = 0;
 193         size_t j = 0;
 194         struct rational *eik;
 195         struct rational *eij;
 196         struct rational *ekj;
 197         int ret = 0;
 198
 199         if (!q) {
 200                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 201
 202                 return EINVAL;
 203         }
 204
 205         /* Step one: complete all triangles */
 206         for (i = 0; i < q->v_num; ++i) {
 207                 if (i == k) {
 208                         continue;
 209                 }
 210
 211                 eik = &(q->e[i * q->v_len + k]);
 212
 213                 for (j = 0; j < q->v_num; ++j) {
 214                         if (j == k ||
 215                             j == i) {
 216                                 continue;
 217                         }
 218
 219                         eij = &(q->e[i * q->v_len + j]);
 220                         ekj = &(q->e[k * q->v_len + j]);
 221
 222                         if (eik->p * ekj->p <= 0) {
 223                                 continue;
 224                         }
 225
 226                         if ((ret = add_to_fraction(abs(eik->p) * ekj->p,
 227                                                    eik->q * ekj->q, eij,
 228                                                    out_errstr))) {
 229                                 goto done;
 230                         }
 231                 }
 232         }
 233
 234         /* Step two: invert all edges that touch k */
 235         for (i = 0; i < q->v_num; ++i) {
 236                 if (i == k) {
 237                         continue;
 238                 }
 239
 240                 eik = &(q->e[i * q->v_len + k]);
 241                 eik->p = -eik->p;
 242         }
 243
 244 done:
 245
 246         return ret;
 247 }