bfzIII.c

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