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 (p * 2 > d) {
  34                         break;
  35                 }
  36
  37                 while (n % p == 0 &&
  38                        d &&
  39                        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         size_t k = 0;
  73         int ret = 0;
  74
  75         if (!q) {
  76                 IF_NZ_SET(out_errstr, L("nonexistant quiver"));
  77
  78                 return EINVAL;
  79         }
  80
  81         if (i >= q->vertices_num ||
  82             j >= q->vertices_num) {
  83                 IF_NZ_SET(out_errstr, L("edge includes nonexistant vertex"));
  84
  85                 return EINVAL;
  86         }
  87
  88         k = (i * q->vertices_len + j);
  89
  90         if ((ret = add_to_fraction(a * q->vertices[i].fatness, b,
  91                                    &(q->edges[k]), out_errstr))) {
  92                 goto done;
  93         }
  94
  95         k = (j * q->vertices_len + i);
  96
  97         if ((ret = add_to_fraction(-1 * a * q->vertices[j].fatness, b,
  98                                    &(q->edges[k]), out_errstr))) {
  99                 goto done;
 100         }
 101
 102 done:
 103
 104         return ret;
 105 }
 106
 107 /* Add a vertex with a name and weight */
 108 int quiver_add_vertex(struct quiver *q, size_t *out_i, const char *name,
 109                       uint_fast16_t fatness, const char **out_errstr)
 110 {
 111         void *newmem = 0;
 112         int sv_err = 0;
 113         size_t l = strlen(name);
 114         char *newname;
 115         size_t k;
 116         size_t newsize;
 117
 118         if (!q) {
 119                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 120
 121                 return EINVAL;
 122         }
 123
 124         if (!(newname = malloc(l + 1))) {
 125                 sv_err = errno;
 126                 IF_NZ_SET(out_errstr, L("malloc"));
 127
 128                 return sv_err;
 129         }
 130
 131         strcat(newname, name);
 132
 133         if (q->vertices_num >= q->vertices_len) {
 134                 newsize = (q->vertices_len + 8) * (q->vertices_len + 8);
 135
 136                 if (!(newmem = realloc(q->edges, newsize *
 137                                        sizeof q->edges[0]))) {
 138                         sv_err = errno;
 139                         IF_NZ_SET(out_errstr, L("too many vertices"));
 140
 141                         return sv_err;
 142                 }
 143
 144                 q->edges = newmem;
 145                 newsize = q->vertices_len + 8;
 146
 147                 if (!(newmem = realloc(q->vertices, newsize *
 148                                        sizeof q->vertices[0]))) {
 149                         sv_err = errno;
 150                         IF_NZ_SET(out_errstr, L("too many vertices"));
 151
 152                         return sv_err;
 153                 }
 154
 155                 q->vertices = newmem;
 156                 q->vertices_len += 8;
 157         }
 158
 159         for (k = 0; k <= q->vertices_num; ++k) {
 160                 q->edges[k * q->vertices_len + q->vertices_num] = (struct
 161                                                                    rational) { .p = 0, .q = 1 };
 162                 q->edges[q->vertices_num * q->vertices_len + k] = (struct
 163                                                                    rational) { .p = 0, .q = 1 };
 164         }
 165
 166         q->vertices[q->vertices_num] = (struct vertex) { .name = newname,
 167                                                          .fatness = fatness };
 168         *out_i = q->vertices_num;
 169         q->vertices_num++;
 170
 171         return 0;
 172 }
 173
 174 /* Mutate the quiver at vertex k */
 175 int quiver_mutate(struct quiver *q, size_t k, const char **out_errstr)
 176 {
 177         size_t i = 0;
 178         size_t j = 0;
 179         size_t ik = 0;
 180         size_t ij = 0;
 181         size_t kj = 0;
 182         struct rational *eik;
 183         struct rational *eij;
 184         struct rational *ekj;
 185         int ret = 0;
 186
 187         if (!q) {
 188                 IF_NZ_SET(out_errstr, L("invalid quiver"));
 189
 190                 return EINVAL;
 191         }
 192
 193         /* Step one: complete all triangles */
 194         for (i = 0; i < q->vertices_num; ++i) {
 195                 if (i == k) {
 196                         continue;
 197                 }
 198
 199                 ik = i * q->vertices_len + k;
 200                 eik = &(q->edges[ik]);
 201
 202                 for (j = 0; j < q->vertices_num; ++j) {
 203                         if (j == k ||
 204                             j == i) {
 205                                 continue;
 206                         }
 207
 208                         ij = i * q->vertices_len + j;
 209                         eij = &(q->edges[ij]);
 210                         kj = k * q->vertices_len + j;
 211                         ekj = &(q->edges[kj]);
 212
 213                         if (eik->p * ekj->p <= 0) {
 214                                 continue;
 215                         }
 216
 217                         if ((ret = add_to_fraction(abs(eik->p) * ekj->p,
 218                                                    eik->q * ekj->q, eij, out_errstr))) {
 219                                 goto done;
 220                         }
 221                 }
 222         }
 223
 224         /* Step two: invert all edges that touch k */
 225         for (i = 0; i < q->vertices_num; ++i) {
 226                 if (i == k) {
 227                         continue;
 228                 }
 229
 230                 ik = i * q->vertices_len + k;
 231                 eik = &(q->edges[ik]);
 232                 eik->p = -eik->p;
 233         }
 234 done:
 235
 236         return ret;
 237 }