bfzIII.c

   1 /*
   2  * Copyright (c) 2016-2021, S. Gilles <sgilles@sgilles.net>
   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 <limits.h>
  20 #include <stdint.h>
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include <unistd.h>
  25
  26 #include "macros.h"
  27 #include "quiver.h"
  28
  29 /*
  30  * This is an implementation of the algorithm of Berenstein, Fomin,
  31  * and Zelevinsky from "Cluster Algebras III", which produces quivers
  32  * for the double Bruhat cell G^{u,v}.
  33  */
  34
  35 #define zabs(x)    (((x) < 0) ? -(x) : (x))
  36 #define zabsdiff(x, y) ((x) < (y) ? ((y) - (x)) : ((x) - (y)))
  37 #define zmin(x, y) (((x) < (y)) ? (x) : (y))
  38 #define zmax(x, y) (((x) > (y)) ? (x) : (y))
  39 #define zsgn(x) ((x) < 0 ? -1 : (x) > 0 ? 1 : 0)
  40
  41 /* Which of the Lie classifications we're looking at */
  42 enum dynkin { INVALID, An, Bn, Cn, Dn, G2, F4, E6, E7, E8 };
  43
  44 /* Cartan matrix entry a_{i,j}, given dynkin type. 1-based. */
  45 static int
  46 aij(enum dynkin typ, size_t n, size_t i, size_t j)
  47 {
  48         i--;
  49         j--;
  50         size_t max_idx = (i > j) ? i : j;
  51         size_t k = 0;
  52
  53         switch (typ) {
  54         case INVALID:
  55
  56                 return INT_MAX;
  57                 break;
  58         case An:
  59
  60                 if (max_idx >= n) {
  61                         return INT_MAX;
  62                 } else if (i == j) {
  63                         return 2;
  64                 } else if (zabsdiff(i, j) == 1) {
  65                         return -1;
  66                 }
  67
  68                 return 0;
  69                 break;
  70         case Bn:
  71         case Cn:
  72
  73                 /* TODO: do you have Bn and Cn switched? You always do that */
  74                 if (max_idx >= n) {
  75                         return INT_MAX;
  76                 } else if (i == j) {
  77                         return 2;
  78                 } else if (zabsdiff(i, j) == 1) {
  79                         if (typ == Cn &&
  80                             i == n - 2 &&
  81                             j == n - 1) {
  82                                 return -2;
  83                         } else if (typ == Bn &&
  84                                    i == n - 1 &&
  85                                    j == n - 2) {
  86                                 return -2;
  87                         }
  88
  89                         return -1;
  90                 }
  91
  92                 return 0;
  93                 break;
  94         case Dn:
  95
  96                 if (max_idx >= n) {
  97                         return INT_MAX;
  98                 } else if (i == j) {
  99                         return 2;
 100                 } else if (zabsdiff(i, j) == 1 &&
 101                            i <= n - 2 &&
 102                            j <= n - 2) {
 103                         return -1;
 104                 } else if ((i == n - 3 &&
 105                             j == n - 1) ||
 106                            (i == n - 1 &&
 107                             j == n - 3)) {
 108                         return -1;
 109                 }
 110
 111                 return 0;
 112                 break;
 113         case G2:
 114
 115                 if (max_idx >= 2) {
 116                         return INT_MAX;
 117                 } else if (i == j) {
 118                         return 2;
 119                 } else if (i == 0 &&
 120                            j == 1) {
 121                         return -3;
 122                 } else if (i == 1 &&
 123                            j == 0) {
 124                         return -1;
 125                 }
 126
 127                 return 0;
 128                 break;
 129         case F4:
 130
 131                 if (max_idx >= 4) {
 132                         return INT_MAX;
 133                 } else if (i == j) {
 134                         return 2;
 135                 } else if (i == 1 &&
 136                            j == 2) {
 137                         return -2;
 138                 } else if (zabsdiff(i, j) == 1) {
 139                         return -1;
 140                 }
 141
 142                 return 0;
 143                 break;
 144         case E6:
 145         case E7:
 146         case E8:
 147                 k = (typ == E6) ? 6 : (typ == E7) ? 7 : 8;
 148
 149                 if (max_idx >= k) {
 150                         return INT_MAX;
 151                 } else if (i == j) {
 152                         return 2;
 153                 } else if (zabsdiff(i, j) == 1 &&
 154                            zmax(i, j) >= 3) {
 155                         return -1;
 156                 } else if (zabsdiff(i, j) == 2 &&
 157                            zmin(i, j) <= 1) {
 158                         return -1;
 159                 }
 160
 161                 return 0;
 162                 break;
 163         }
 164
 165         return INT_MAX;
 166 }
 167
 168 const int x_mult = 25;
 169 const int y_mult = 25;
 170
 171 /* Given vertex idx of k, return vertex idx of k⁺ */
 172 static int
 173 find_plus(int k, size_t *u, size_t lu, size_t *v, size_t lv)
 174 {
 175         size_t vk = (size_t) -1;
 176
 177         if (k < 0) {
 178                 vk = -1 * (k + 1);
 179
 180                 for (int j = 1; j <= (int) (lu + lv); ++j) {
 181                         size_t vj = (size_t) -1;
 182
 183                         vj = (j <= (int) lu) ? u[j - 1] : v[j - 1 - lu];
 184
 185                         if (vk == vj) {
 186                                 return j;
 187                         }
 188                 }
 189         } else if (1 <= k &&
 190                    k <= (int) lu) {
 191                 vk = u[k - 1];
 192
 193                 for (int j = k + 1; j <= (int) (lu + lv); ++j) {
 194                         size_t vj = (size_t) -1;
 195
 196                         vj = (j <= (int) lu) ? u[j - 1] : v[j - 1 - lu];
 197
 198                         if (vk == vj) {
 199                                 return j;
 200                         }
 201                 }
 202         } else if ((int) lu + 1 <= k &&
 203                    k <= (int) (lu + lv)) {
 204                 vk = v[k - 1 - lu];
 205
 206                 for (int j = k + 1; j <= (int) (lu + lv); ++j) {
 207                         size_t vj = v[j - 1 - lu];
 208
 209                         if (vk == vj) {
 210                                 return j;
 211                         }
 212                 }
 213         }
 214
 215         return lu + lv + 1;
 216 }
 217
 218 /* Given k in { -1, ..., -r } U { 1, 2, ..., l(u) + l(v) }, return the vertex idx for it */
 219 static size_t
 220 v_of(int k, int r, size_t lu, size_t lv)
 221 {
 222         if (-1 >= k &&
 223             k >= -r) {
 224                 return (-1 * k) - 1;
 225         } else if (1 <= k &&
 226                    k <= (int) (lu + lv)) {
 227                 return k - 1 + r;
 228         }
 229
 230         return (size_t) -1;
 231 }
 232
 233 /* [BFZIII] */
 234 static int
 235 run_bfzIII(struct quiver *qq, enum dynkin typ, size_t r, size_t *u, size_t lu,
 236            size_t *v, size_t lv, int ymult)
 237 {
 238         int ret = -1;
 239         size_t nl = 1 + snprintf(0, 0, "%d%zu", -1 * (int) r, lu + lv);
 240         char *nbuf = 0;
 241         const char *errstr = 0;
 242         size_t x = 0;
 243         int *i = 0;
 244         int *e = 0;
 245         uint32_t c_froz = 0x4ce7ff;
 246         uint32_t c_nonf = 0xad7fa8;
 247         int need_fatness_run = 1;
 248
 249         if (!(nbuf = malloc(nl))) {
 250                 goto done;
 251         }
 252
 253         if (!(i = calloc(r + lu + lv + 1, sizeof *i))) {
 254                 goto done;
 255         }
 256
 257         if (!(e = calloc(r + lu + lv + 1, sizeof *e))) {
 258                 goto done;
 259         }
 260
 261         /* First, add { -1, -2, ..., -r } */
 262         for (int k = -1; k >= -1 * (int) r; k--) {
 263                 size_t vv = v_of(k, r, lu, lv);
 264
 265                 i[vv] = -1 * k;
 266                 sprintf(nbuf, "%d", k);
 267
 268                 if (quiver_add_vertex(qq, 0, nbuf, 1, 50 * x, ymult * 50 *
 269                                       i[vv], c_froz, &errstr) < 0) {
 270                         fprintf(stderr, "%s\n", errstr);
 271                         fprintf(stderr, L("quiver_add_vertex failed\n"));
 272                         goto done;
 273                 }
 274
 275                 x++;
 276         }
 277
 278         /* Now add { 1, 2, ..., u_seq_len + v_seq_len } */
 279         for (int k = 1; k <= (int) (lu + lv); ++k) {
 280                 size_t vv = v_of(k, r, lu, lv);
 281                 int kplus = find_plus(k, u, lu, v, lv);
 282                 uint32_t c = (kplus <= (int) (lu + lv)) ? c_nonf : c_froz;
 283
 284                 i[vv] = 1 + ((k <= (int) lu) ? u[k - 1] : v[k - lu - 1]);
 285                 sprintf(nbuf, "%d", k);
 286
 287                 if (quiver_add_vertex(qq, 0, nbuf, 1, 50 * x, ymult * 50 *
 288                                       i[vv], c, &errstr) < 0) {
 289                         fprintf(stderr, "%s\n", errstr);
 290                         fprintf(stderr, L("quiver_add_vertex failed\n"));
 291                         goto done;
 292                 }
 293
 294                 x++;
 295         }
 296
 297         /* Now add edges by remark 2.4 */
 298         for (int k = -r; k <= (int) (lu + lv); ++k) {
 299                 int kplus = find_plus(k, u, lu, v, lv);
 300                 size_t vk = v_of(k, r, lu, lv);
 301
 302                 if (k == 0) {
 303                         continue;
 304                 }
 305
 306                 for (int l = -r; l <= (int) (lu + lv); ++l) {
 307                         int lplus = find_plus(l, u, lu, v, lv);
 308                         int p = zmax(k, l);
 309                         int eip = (1 <= p &&
 310                                    p <= (int) lu) ? -1 : 1;
 311                         int q = zmin((int) kplus, (int) lplus);
 312                         int eiq = (1 <= q &&
 313                                    q <= (int) lu) ? -1 : 1;
 314                         size_t vl = v_of(l, r, lu, lv);
 315                         size_t idx = vk * qq->v_len + vl;
 316                         struct rational *e = &qq->e[idx];
 317
 318                         if (l == 0 ||
 319                             (k < 0 &&
 320                              l < 0)) {
 321                                 continue;
 322                         }
 323
 324                         if (p == q) {
 325                                 e->p = zsgn(k - l) * eip;
 326                                 e->q = 1;
 327                         } else if (p < q &&
 328                                    eip * eiq * (k - l) * (kplus - lplus) > 0) {
 329                                 e->p = zsgn(k - l) * eip * aij(typ, r, i[vk],
 330                                                                i[vl]);
 331                                 e->q = 1;
 332                         }
 333                 }
 334         }
 335
 336         /*
 337          * Now modify the fatness of various vertices to normalize |sigma|.
 338          *
 339          * We could just compute this by looking at the y-position
 340          * and the Dynkin type, but this serves as a nice sanity check.
 341          */
 342         while (need_fatness_run) {
 343                 need_fatness_run = 0;
 344
 345                 for (size_t j = 0; j < qq->v_num; ++j) {
 346                         struct vertex *v = &qq->v[j];
 347
 348                         for (size_t k = 0; k < qq->v_num; ++k) {
 349                                 struct vertex *w = &qq->v[k];
 350                                 struct rational *e = &qq->e[j * qq->v_len + k];
 351                                 struct rational *f = &qq->e[k * qq->v_len + j];
 352
 353                                 if (zabs(e->p * v->fatness) > zabs(f->p *
 354                                                                    w->fatness))
 355                                 {
 356                                         w->fatness++;
 357                                         need_fatness_run = 1;
 358                                 } else if (zabs(e->p * v->fatness) < zabs(f->p *
 359                                                                           w->
 360                                                                           fatness))
 361                                 {
 362                                         v->fatness++;
 363                                         need_fatness_run = 1;
 364                                 }
 365                         }
 366                 }
 367         }
 368
 369         quiver_save_to_file(qq, stdout, 0);
 370         ret = 0;
 371 done:
 372         free(nbuf);
 373         free(i);
 374         free(e);
 375
 376         return ret;
 377 }
 378
 379 /* "e" -> { }, "1,2" -> { 1, 2 } */
 380 static int
 381 parse_word(const char *s, size_t n, size_t **out_word, size_t *out_len)
 382 {
 383         size_t *seq = 0;
 384         size_t len = 1;
 385         size_t j = 0;
 386         int ret = -1;
 387         char *err = 0;
 388
 389         if (!(strcmp(s, "e"))) {
 390                 *out_word = 0;
 391                 *out_len = 0;
 392
 393                 return 0;
 394         }
 395
 396         for (const char *p = s; *p; p++) {
 397                 len += (*p == ',');
 398         }
 399
 400         if (!(seq = calloc(len, sizeof *seq))) {
 401                 perror(L("calloc"));
 402                 goto done;
 403         }
 404
 405         errno = 0;
 406         seq[j] = strtoll(s, 0, 0) - 1;
 407
 408         if (seq[j] >= n ||
 409             errno) {
 410                 goto done;
 411         }
 412
 413         j++;
 414
 415         for (const char *p = s; *p; p++) {
 416                 if (*p != ',') {
 417                         continue;
 418                 }
 419
 420                 errno = 0;
 421                 seq[j] = strtoll(p + 1, &err, 0) - 1;
 422
 423                 if (seq[j] >= n ||
 424                     errno) {
 425                         goto done;
 426                 }
 427
 428                 j++;
 429         }
 430
 431         ret = 0;
 432 done:
 433
 434         if (ret < 0) {
 435                 free(seq);
 436         } else {
 437                 *out_word = seq;
 438                 *out_len = len;
 439         }
 440
 441         return ret;
 442 }
 443
 444 static void
 445 usage(char *argv0)
 446 {
 447         printf(
 448                 "Usage: %s -c <Dynkin classification> [ -u <word>]  [-v <word>] [ -U ]\n",
 449                 argv0);
 450         printf(
 451                 "       Dynkin classification: Something like \u2018F4\u2019; permitted are\n");
 452         printf("           A\u2081, A\u2082, A\u2083, \u2026\n");
 453         printf("           B\u2081, B\u2082, B\u2083, \u2026\n");
 454         printf("           C\u2081, C\u2082, C\u2083, \u2026\n");
 455         printf("           D\u2081, D\u2082, D\u2083, \u2026\n");
 456         printf("           G\u2082, F\u2084, E\u2086, E\u2087, E\u2088\n");
 457         printf(
 458                 "       Word: something like \u20181,2,1,2\u2019 (\u2018e\u2019 is permitted)\n");
 459         printf("             if none given, \u2018e\u2019 is assumed\n");
 460         printf("\n");
 461         printf("        With -U, -r is above -1. Without, -1 is above -r\n");
 462 }
 463
 464 /* Run the thing */
 465 int
 466 main(int argc, char **argv)
 467 {
 468         int ret = 0;
 469         int opt = 0;
 470         enum dynkin typ = INVALID;
 471         int n = -1;
 472         char *p = 0;
 473         char *uword = "e";
 474         char *vword = "e";
 475         int ymult = 1;
 476         size_t *u_seq = 0;
 477         size_t *v_seq = 0;
 478         size_t u_seq_len = 0;
 479         size_t v_seq_len = 0;
 480         struct quiver q = { 0 };
 481
 482         while ((opt = getopt(argc, argv, "hUc:u:v:")) != -1) {
 483                 switch (opt) {
 484                 case 'h':
 485                         usage(argv[0]);
 486                         ret = 0;
 487                         goto done;
 488                 case 'u':
 489                         uword = optarg;
 490                         break;
 491                 case 'v':
 492                         vword = optarg;
 493                         break;
 494                 case 'U':
 495                         ymult = -1;
 496                         break;
 497                 case 'c':
 498
 499                         switch (optarg[0]) {
 500                         case 'a':
 501                         case 'A':
 502                                 typ = An;
 503                                 n = strtoll(optarg + 1, &p, 0);
 504                                 break;
 505                         case 'b':
 506                         case 'B':
 507                                 typ = Bn;
 508                                 n = strtoll(optarg + 1, &p, 0);
 509                                 break;
 510                         case 'c':
 511                         case 'C':
 512                                 typ = Cn;
 513                                 n = strtoll(optarg + 1, &p, 0);
 514                                 break;
 515                         case 'd':
 516                         case 'D':
 517                                 typ = Dn;
 518                                 n = strtoll(optarg + 1, &p, 0);
 519                                 break;
 520                         case 'e':
 521                         case 'E':
 522                                 n = strtoll(optarg + 1, &p, 0);
 523
 524                                 switch (n) {
 525                                 case 6:
 526                                         typ = E6;
 527                                         break;
 528                                 case 7:
 529                                         typ = E7;
 530                                         break;
 531                                 case 8:
 532                                         typ = E8;
 533                                         break;
 534                                 default:
 535                                         fprintf(stderr,
 536                                                 "Type E may only be E\u2086, E\u2087, or E\u2088\n");
 537                                         n = -1;
 538                                         break;
 539                                 }
 540
 541                                 break;
 542                         case 'f':
 543                         case 'F':
 544
 545                                 if (optarg[1] == '4' &&
 546                                     optarg[2] == 0) {
 547                                         typ = F4;
 548                                         n = 4;
 549                                 } else {
 550                                         fprintf(stderr,
 551                                                 "Type F may only be F\u2084\n");
 552                                 }
 553
 554                                 break;
 555                         case 'g':
 556                         case 'G':
 557
 558                                 if (optarg[1] == '2' &&
 559                                     optarg[2] == 0) {
 560                                         typ = G2;
 561                                         n = 2;
 562                                 } else {
 563                                         fprintf(stderr,
 564                                                 "Type G may only be G\u2082\n");
 565                                 }
 566
 567                                 break;
 568                         default:
 569                                 fprintf(stderr,
 570                                         "Invalid Dynkin classification: \u2018%s\u2019\n",
 571                                         optarg);
 572                                 goto done;
 573                                 break;
 574                         }
 575
 576                         if (n <= 0 ||
 577                             (p &&
 578                              *p)) {
 579                                 fprintf(stderr, "Invalid n: \u2018%s\u2019\n",
 580                                         optarg + 1);
 581                                 ret = EINVAL;
 582                                 goto done;
 583                         }
 584                 }
 585         }
 586
 587         switch (typ) {
 588         case INVALID:
 589                 fprintf(stderr,
 590                         "Dynkin classification is required (ex: \u2018-c F4\u2019)\n");
 591                 ret = EINVAL;
 592                 goto done;
 593                 break;
 594         default:
 595                 break;
 596         }
 597
 598         if (parse_word(uword, n, &u_seq, &u_seq_len) < 0) {
 599                 fprintf(stderr, "Invalid u-word: \u2018%s\u2019\n", uword);
 600                 ret = EINVAL;
 601                 goto done;
 602         }
 603
 604         if (parse_word(vword, n, &v_seq, &v_seq_len) < 0) {
 605                 fprintf(stderr, "Invalid v-word: \u2018%s\u2019\n", vword);
 606                 ret = EINVAL;
 607                 goto done;
 608         }
 609
 610         if ((ret = run_bfzIII(&q, typ, n, u_seq, u_seq_len, v_seq, v_seq_len,
 611                               ymult)) < 0) {
 612                 fprintf(stderr, "BFZIII algorithm failed\n");
 613                 goto done;
 614         }
 615
 616 done:
 617         quiver_clean(&q);
 618         free(u_seq);
 619         free(v_seq);
 620
 621         return ret;
 622 }