bernoulli.c

   1 #include <assert.h>
   2 #include <barvinok/options.h>
   3 #include <barvinok/util.h>
   4 #include "bernoulli.h"
   5
   6 #define ALLOC(type) (type*)malloc(sizeof(type))
   7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   8 #define REALLOCN(ptr,type,n) (type*)realloc(ptr, (n) * sizeof(type))
   9
  10 static struct bernoulli_coef bernoulli_coef;
  11 static struct poly_list bernoulli;
  12 static struct poly_list faulhaber;
  13
  14 struct bernoulli_coef *bernoulli_coef_compute(int n)
  15 {
  16     int i, j;
  17     Value factor, tmp;
  18
  19     if (n < bernoulli_coef.n)
  20         return &bernoulli_coef;
  21
  22     if (n >= bernoulli_coef.size) {
  23         int size = 3*(n + 5)/2;
  24         Vector *b;
  25
  26         b = Vector_Alloc(size);
  27         if (bernoulli_coef.n) {
  28             Vector_Copy(bernoulli_coef.num->p, b->p, bernoulli_coef.n);
  29             Vector_Free(bernoulli_coef.num);
  30         }
  31         bernoulli_coef.num = b;
  32         b = Vector_Alloc(size);
  33         if (bernoulli_coef.n) {
  34             Vector_Copy(bernoulli_coef.den->p, b->p, bernoulli_coef.n);
  35             Vector_Free(bernoulli_coef.den);
  36         }
  37         bernoulli_coef.den = b;
  38         b = Vector_Alloc(size);
  39         if (bernoulli_coef.n) {
  40             Vector_Copy(bernoulli_coef.lcm->p, b->p, bernoulli_coef.n);
  41             Vector_Free(bernoulli_coef.lcm);
  42         }
  43         bernoulli_coef.lcm = b;
  44
  45         bernoulli_coef.size = size;
  46     }
  47     value_init(factor);
  48     value_init(tmp);
  49     for (i = bernoulli_coef.n; i <= n; ++i) {
  50         if (i == 0) {
  51             value_set_si(bernoulli_coef.num->p[0], 1);
  52             value_set_si(bernoulli_coef.den->p[0], 1);
  53             value_set_si(bernoulli_coef.lcm->p[0], 1);
  54             continue;
  55         }
  56         value_set_si(bernoulli_coef.num->p[i], 0);
  57         value_set_si(factor, -(i+1));
  58         for (j = i-1; j >= 0; --j) {
  59             mpz_mul_ui(factor, factor, j+1);
  60             mpz_divexact_ui(factor, factor, i+1-j);
  61             value_division(tmp, bernoulli_coef.lcm->p[i-1],
  62                            bernoulli_coef.den->p[j]);
  63             value_multiply(tmp, tmp, bernoulli_coef.num->p[j]);
  64             value_multiply(tmp, tmp, factor);
  65             value_addto(bernoulli_coef.num->p[i], bernoulli_coef.num->p[i], tmp);
  66         }
  67         mpz_mul_ui(bernoulli_coef.den->p[i], bernoulli_coef.lcm->p[i-1], i+1);
  68         value_gcd(tmp, bernoulli_coef.num->p[i], bernoulli_coef.den->p[i]);
  69         if (value_notone_p(tmp)) {
  70             value_division(bernoulli_coef.num->p[i],
  71                             bernoulli_coef.num->p[i], tmp);
  72             value_division(bernoulli_coef.den->p[i],
  73                             bernoulli_coef.den->p[i], tmp);
  74         }
  75         value_lcm(bernoulli_coef.lcm->p[i],
  76                   bernoulli_coef.lcm->p[i-1], bernoulli_coef.den->p[i]);
  77     }
  78     bernoulli_coef.n = n+1;
  79     value_clear(factor);
  80     value_clear(tmp);
  81
  82     return &bernoulli_coef;
  83 }
  84
  85 /*
  86  * Compute either Bernoulli B_n or Faulhaber F_n polynomials.
  87  *
  88  * B_n =         sum_{k=0}^n {  n  \choose k } b_k x^{n-k}
  89  * F_n = 1/(n+1) sum_{k=0}^n { n+1 \choose k } b_k x^{n+1-k}
  90  */
  91 static struct poly_list *bernoulli_faulhaber_compute(int n, struct poly_list *pl,
  92                                                      int faulhaber)
  93 {
  94     int i, j;
  95     Value factor;
  96     struct bernoulli_coef *bc;
  97
  98     if (n < pl->n)
  99         return pl;
 100
 101     if (n >= pl->size) {
 102         int size = 3*(n + 5)/2;
 103         Vector **poly;
 104
 105         poly = ALLOCN(Vector *, size);
 106         for (i = 0; i < pl->n; ++i)
 107             poly[i] = pl->poly[i];
 108         free(pl->poly);
 109         pl->poly = poly;
 110
 111         pl->size = size;
 112     }
 113
 114     bc = bernoulli_coef_compute(n);
 115
 116     value_init(factor);
 117     for (i = pl->n; i <= n; ++i) {
 118         pl->poly[i] = Vector_Alloc(i+faulhaber+2);
 119         value_assign(pl->poly[i]->p[i+faulhaber], bc->lcm->p[i]);
 120         if (faulhaber)
 121             mpz_mul_ui(pl->poly[i]->p[i+2], bc->lcm->p[i], i+1);
 122         else
 123             value_assign(pl->poly[i]->p[i+1], bc->lcm->p[i]);
 124         value_set_si(factor, 1);
 125         for (j = 1; j <= i; ++j) {
 126             mpz_mul_ui(factor, factor, i+faulhaber+1-j);
 127             mpz_divexact_ui(factor, factor, j);
 128             value_division(pl->poly[i]->p[i+faulhaber-j],
 129                            bc->lcm->p[i], bc->den->p[j]);
 130             value_multiply(pl->poly[i]->p[i+faulhaber-j],
 131                            pl->poly[i]->p[i+faulhaber-j], bc->num->p[j]);
 132             value_multiply(pl->poly[i]->p[i+faulhaber-j],
 133                            pl->poly[i]->p[i+faulhaber-j], factor);
 134         }
 135         Vector_Normalize(pl->poly[i]->p, i+faulhaber+2);
 136     }
 137     value_clear(factor);
 138     pl->n = n+1;
 139
 140     return pl;
 141 }
 142
 143 struct poly_list *bernoulli_compute(int n)
 144 {
 145     return bernoulli_faulhaber_compute(n, &bernoulli, 0);
 146 }
 147
 148 struct poly_list *faulhaber_compute(int n)
 149 {
 150     return bernoulli_faulhaber_compute(n, &faulhaber, 1);
 151 }
 152
 153 /* shift variables in polynomial one down */
 154 static void shift(evalue *e)
 155 {
 156     int i;
 157     if (value_notzero_p(e->d))
 158         return;
 159     assert(e->x.p->type == polynomial);
 160     assert(e->x.p->pos > 1);
 161     e->x.p->pos--;
 162     for (i = 0; i < e->x.p->size; ++i)
 163         shift(&e->x.p->arr[i]);
 164 }
 165
 166 static evalue *shifted_copy(evalue *src)
 167 {
 168     evalue *e = ALLOC(evalue);
 169     value_init(e->d);
 170     evalue_copy(e, src);
 171     shift(e);
 172     return e;
 173 }
 174
 175 static evalue *power_sums(struct poly_list *faulhaber, evalue *poly,
 176                           Vector *arg, Value denom, int neg, int alt_neg)
 177 {
 178     int i;
 179     evalue *base = affine2evalue(arg->p, denom, arg->Size-1);
 180     evalue *sum = evalue_zero();
 181
 182     for (i = 1; i < poly->x.p->size; ++i) {
 183         evalue *term = evalue_polynomial(faulhaber->poly[i], base);
 184         evalue *factor = shifted_copy(&poly->x.p->arr[i]);
 185         emul(factor, term);
 186         if (alt_neg && (i % 2))
 187             evalue_negate(term);
 188         eadd(term, sum);
 189         evalue_free(factor);
 190         evalue_free(term);
 191     }
 192     if (neg)
 193         evalue_negate(sum);
 194     evalue_free(base);
 195
 196     return sum;
 197 }
 198
 199 /* Given a constraint (cst_affine) a x + b y + c >= 0, compate a constraint (c)
 200  * +- (b y + c) +- a >=,> 0
 201  * ^            ^      ^
 202  * |            |      strict
 203  * sign_affine  sign_cst
 204  */
 205 static void bound_constraint(Value *c, unsigned dim,
 206                              Value *cst_affine,
 207                              int sign_affine, int sign_cst, int strict)
 208 {
 209     if (sign_affine == -1)
 210         Vector_Oppose(cst_affine+1, c, dim+1);
 211     else
 212         Vector_Copy(cst_affine+1, c, dim+1);
 213
 214     if (sign_cst == -1)
 215         value_subtract(c[dim], c[dim], cst_affine[0]);
 216     else if (sign_cst == 1)
 217         value_addto(c[dim], c[dim], cst_affine[0]);
 218
 219     if (strict)
 220         value_decrement(c[dim], c[dim]);
 221 }
 222
 223 struct Bernoulli_data {
 224     unsigned MaxRays;
 225     struct evalue_section *s;
 226     int size;
 227     int ns;
 228     evalue *e;
 229 };
 230
 231 static evalue *compute_poly_u(evalue *poly_u, Value *upper, Vector *row,
 232                               unsigned dim, Value tmp,
 233                               struct poly_list *faulhaber,
 234                               struct Bernoulli_data *data)
 235 {
 236     if (poly_u)
 237         return poly_u;
 238     Vector_Copy(upper+2, row->p, dim+1);
 239     value_oppose(tmp, upper[1]);
 240     value_addto(row->p[dim], row->p[dim], tmp);
 241     return power_sums(faulhaber, data->e, row, tmp, 0, 0);
 242 }
 243
 244 static evalue *compute_poly_l(evalue *poly_l, Value *lower, Vector *row,
 245                               unsigned dim,
 246                               struct poly_list *faulhaber,
 247                               struct Bernoulli_data *data)
 248 {
 249     if (poly_l)
 250         return poly_l;
 251     Vector_Copy(lower+2, row->p, dim+1);
 252     value_addto(row->p[dim], row->p[dim], lower[1]);
 253     return power_sums(faulhaber, data->e, row, lower[1], 0, 1);
 254 }
 255
 256 static void Bernoulli_init(unsigned n, void *cb_data)
 257 {
 258     struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
 259     int cases = 5;
 260
 261     if (cases * n <= data->size)
 262         return;
 263
 264     data->size = cases * (n + 16);
 265     data->s = REALLOCN(data->s, struct evalue_section, data->size);
 266 }
 267
 268 static void Bernoulli_cb(Matrix *M, Value *lower, Value *upper, void *cb_data)
 269 {
 270     struct Bernoulli_data *data = (struct Bernoulli_data *)cb_data;
 271     Matrix *M2;
 272     Polyhedron *T;
 273     evalue *factor = NULL;
 274     evalue *linear = NULL;
 275     int constant = 0;
 276     Value tmp;
 277     unsigned dim = M->NbColumns-2;
 278     Vector *row;
 279     int cases = 5;
 280
 281     assert(lower);
 282     assert(upper);
 283     assert(data->ns + cases <= data->size);
 284
 285     M2 = Matrix_Copy(M);
 286     T = Constraints2Polyhedron(M2, data->MaxRays);
 287     Matrix_Free(M2);
 288
 289     POL_ENSURE_VERTICES(T);
 290     if (emptyQ(T)) {
 291         Polyhedron_Free(T);
 292         return;
 293     }
 294
 295     assert(lower != upper);
 296
 297     row = Vector_Alloc(dim+1);
 298     value_init(tmp);
 299     if (value_notzero_p(data->e->d)) {
 300         factor = data->e;
 301         constant = 1;
 302     } else {
 303         assert(data->e->x.p->type == polynomial);
 304         if (data->e->x.p->pos > 1) {
 305             factor = shifted_copy(data->e);
 306             constant = 1;
 307         } else
 308             factor = shifted_copy(&data->e->x.p->arr[0]);
 309     }
 310     if (!EVALUE_IS_ZERO(*factor)) {
 311         value_absolute(tmp, upper[1]);
 312         /* upper - lower */
 313         Vector_Combine(lower+2, upper+2, row->p, tmp, lower[1], dim+1);
 314         value_multiply(tmp, tmp, lower[1]);
 315         /* upper - lower + 1 */
 316         value_addto(row->p[dim], row->p[dim], tmp);
 317
 318         linear = affine2evalue(row->p, tmp, dim);
 319         emul(factor, linear);
 320     } else
 321         linear = evalue_zero();
 322
 323     if (constant) {
 324         data->s[data->ns].E = linear;
 325         data->s[data->ns].D = T;
 326         ++data->ns;
 327     } else {
 328         evalue *poly_u = NULL, *poly_l = NULL;
 329         Polyhedron *D;
 330         struct poly_list *faulhaber;
 331         assert(data->e->x.p->type == polynomial);
 332         assert(data->e->x.p->pos == 1);
 333         faulhaber = faulhaber_compute(data->e->x.p->size-1);
 334         /* lower is the constraint
 335          *                       b i - l' >= 0          i >= l'/b = l
 336          * upper is the constraint
 337          *                      -a i + u' >= 0          i <= u'/a = u
 338          */
 339         M2 = Matrix_Alloc(3, 2+T->Dimension);
 340         value_set_si(M2->p[0][0], 1);
 341         value_set_si(M2->p[1][0], 1);
 342         value_set_si(M2->p[2][0], 1);
 343         /* Case 1:
 344          *              1 <= l          l' - b >= 0
 345          */
 346         bound_constraint(M2->p[0]+1, T->Dimension, lower+1, -1, -1, 0);
 347         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 348         if (emptyQ2(D))
 349             Polyhedron_Free(D);
 350         else {
 351             evalue *extra;
 352             poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
 353                                         faulhaber, data);
 354             Vector_Oppose(lower+2, row->p, dim+1);
 355             extra = power_sums(faulhaber, data->e, row, lower[1], 1, 0);
 356             eadd(poly_u, extra);
 357             eadd(linear, extra);
 358
 359             data->s[data->ns].E = extra;
 360             data->s[data->ns].D = D;
 361             ++data->ns;
 362         }
 363
 364         /* Case 2:
 365          *              1 <= -u         -u' - a >= 0
 366          */
 367         bound_constraint(M2->p[0]+1, T->Dimension, upper+1, -1, 1, 0);
 368         D = AddConstraints(M2->p_Init, 1, T, data->MaxRays);
 369         if (emptyQ2(D))
 370             Polyhedron_Free(D);
 371         else {
 372             evalue *extra;
 373             poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
 374             Vector_Oppose(upper+2, row->p, dim+1);
 375             value_oppose(tmp, upper[1]);
 376             extra = power_sums(faulhaber, data->e, row, tmp, 1, 1);
 377             eadd(poly_l, extra);
 378             eadd(linear, extra);
 379
 380             data->s[data->ns].E = extra;
 381             data->s[data->ns].D = D;
 382             ++data->ns;
 383         }
 384
 385         /* Case 3:
 386          *              u >= 0          u' >= 0
 387          *              -l >= 0         -l' >= 0
 388          */
 389         bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, 0, 0);
 390         bound_constraint(M2->p[1]+1, T->Dimension, lower+1, 1, 0, 0);
 391         D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
 392         if (emptyQ2(D))
 393             Polyhedron_Free(D);
 394         else {
 395             poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
 396             poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
 397                                         faulhaber, data);
 398
 399             data->s[data->ns].E = ALLOC(evalue);
 400             value_init(data->s[data->ns].E->d);
 401             evalue_copy(data->s[data->ns].E, poly_u);
 402             eadd(poly_l, data->s[data->ns].E);
 403             eadd(linear, data->s[data->ns].E);
 404             data->s[data->ns].D = D;
 405             ++data->ns;
 406         }
 407
 408         /* Case 4:
 409          *              l < 1           -l' + b - 1 >= 0
 410          *              0 < l           l' - 1 >= 0
 411          */
 412         bound_constraint(M2->p[0]+1, T->Dimension, lower+1, 1, 1, 1);
 413         bound_constraint(M2->p[1]+1, T->Dimension, lower+1, -1, 0, 1);
 414         D = AddConstraints(M2->p_Init, 2, T, data->MaxRays);
 415         if (emptyQ2(D))
 416             Polyhedron_Free(D);
 417         else {
 418             poly_u = compute_poly_u(poly_u, upper, row, dim, tmp,
 419                                         faulhaber, data);
 420
 421             eadd(linear, poly_u);
 422             data->s[data->ns].E = poly_u;
 423             poly_u = NULL;
 424             data->s[data->ns].D = D;
 425             ++data->ns;
 426         }
 427
 428         /* Case 5:
 429          *              -u < 1          u' + a - 1 >= 0
 430          *              0 < -u          -u' - 1 >= 0
 431          *              l <= 0          -l' >= 0
 432          */
 433         bound_constraint(M2->p[0]+1, T->Dimension, upper+1, 1, -1, 1);
 434         bound_constraint(M2->p[1]+1, T->Dimension, upper+1, -1, 0, 1);
 435         bound_constraint(M2->p[2]+1, T->Dimension, lower+1, 1, 0, 0);
 436         D = AddConstraints(M2->p_Init, 3, T, data->MaxRays);
 437         if (emptyQ2(D))
 438             Polyhedron_Free(D);
 439         else {
 440             poly_l = compute_poly_l(poly_l, lower, row, dim, faulhaber, data);
 441
 442             eadd(linear, poly_l);
 443             data->s[data->ns].E = poly_l;
 444             poly_l = NULL;
 445             data->s[data->ns].D = D;
 446             ++data->ns;
 447         }
 448
 449         Matrix_Free(M2);
 450         Polyhedron_Free(T);
 451         if (poly_l)
 452             evalue_free(poly_l);
 453         if (poly_u)
 454             evalue_free(poly_u);
 455         evalue_free(linear);
 456     }
 457     if (factor != data->e)
 458         evalue_free(factor);
 459     value_clear(tmp);
 460     Vector_Free(row);
 461 }
 462
 463 /* Looks for variable with integer bounds, i.e., with coefficients 0, 1 or -1.
 464  * Returns 1 if such a variable is found and puts it in the first position,
 465  * possibly changing *P_p and *E_p.
 466  */
 467 static int find_integer_bounds(Polyhedron **P_p, evalue **E_p, unsigned nvar)
 468 {
 469     Polyhedron *P = *P_p;
 470     evalue *E = *E_p;
 471     unsigned dim = P->Dimension;
 472     int i, j;
 473
 474     for (i = 0; i < nvar; ++i) {
 475         for (j = 0; j < P->NbConstraints; ++j) {
 476             if (value_zero_p(P->Constraint[j][1+i]))
 477                 continue;
 478             if (value_one_p(P->Constraint[j][1+i]))
 479                 continue;
 480             if (value_mone_p(P->Constraint[j][1+i]))
 481                 continue;
 482             break;
 483         }
 484         if (j == P->NbConstraints)
 485             break;
 486     }
 487     if (i == nvar)
 488         return 0;
 489     if (i == 0)
 490         return 1;
 491     P = Polyhedron_Copy(P);
 492     Polyhedron_ExchangeColumns(P, 1, 1+i);
 493     *P_p = P;
 494
 495     if (value_zero_p(E->d)) {
 496         evalue **subs;
 497         subs = ALLOCN(evalue *, dim);
 498         for (j = 0; j < dim; ++j)
 499             subs[j] = evalue_var(j);
 500         E = subs[0];
 501         subs[0] = subs[i];
 502         subs[i] = E;
 503         E = evalue_dup(*E_p);
 504         evalue_substitute(E, subs);
 505         for (j = 0; j < dim; ++j)
 506             evalue_free(subs[j]);
 507         free(subs);
 508         *E_p = E;
 509     }
 510
 511     return 1;
 512 }
 513
 514 static evalue *sum_over_polytope(Polyhedron *P, evalue *E, unsigned nvar,
 515                                  struct Bernoulli_data *data,
 516                                  struct barvinok_options *options)
 517 {
 518     unsigned dim = P->Dimension - 1;
 519     evalue *res;
 520
 521     if (value_zero_p(P->Constraint[0][0]) &&
 522             value_notzero_p(P->Constraint[0][1])) {
 523         res = ALLOC(evalue);
 524         value_init(res->d);
 525         value_set_si(res->d, 0);
 526         res->x.p = new_enode(partition, 2, dim);
 527         EVALUE_SET_DOMAIN(res->x.p->arr[0], Polyhedron_Project(P, dim));
 528         evalue_copy(&res->x.p->arr[1], E);
 529         reduce_evalue_in_domain(&res->x.p->arr[1], P);
 530         shift(&res->x.p->arr[1]);
 531     } else {
 532         data->ns = 0;
 533         data->e = E;
 534
 535         for_each_lower_upper_bound(P, Bernoulli_init, Bernoulli_cb, data);
 536
 537         res = evalue_from_section_array(data->s, data->ns);
 538     }
 539
 540     if (nvar > 1) {
 541         evalue *tmp = Bernoulli_sum_evalue(res, nvar-1, options);
 542         evalue_free(res);
 543         res = tmp;
 544     }
 545
 546     return res;
 547 }
 548
 549 evalue *Bernoulli_sum_evalue(evalue *e, unsigned nvar,
 550                              struct barvinok_options *options)
 551 {
 552     struct Bernoulli_data data;
 553     int i, j;
 554     evalue *sum = evalue_zero();
 555
 556     if (EVALUE_IS_ZERO(*e))
 557         return sum;
 558
 559     if (nvar == 0) {
 560         eadd(e, sum);
 561         return sum;
 562     }
 563
 564     assert(value_zero_p(e->d));
 565     assert(e->x.p->type == partition);
 566
 567     data.size = 16;
 568     data.s = ALLOCN(struct evalue_section, data.size);
 569     data.MaxRays = options->MaxRays;
 570
 571     for (i = 0; i < e->x.p->size/2; ++i) {
 572         Polyhedron *D;
 573         for (D = EVALUE_DOMAIN(e->x.p->arr[2*i]); D; D = D->next) {
 574             evalue *E = &e->x.p->arr[2*i+1];
 575             Polyhedron *P = D;
 576             Polyhedron *next = D->next;
 577             evalue *tmp;
 578             int integer_bounds;
 579
 580             P->next = NULL;
 581
 582             integer_bounds = find_integer_bounds(&P, &E, nvar);
 583             if (options->approximation_method == BV_APPROX_NONE &&
 584                 !integer_bounds) {
 585                 evalue_free(sum);
 586                 sum = NULL;
 587             } else {
 588                 evalue *tmp = sum_over_polytope(P, E, nvar, &data, options);
 589                 eadd(tmp, sum);
 590                 evalue_free(tmp);
 591             }
 592
 593             if (P != D)
 594                 Polyhedron_Free(P);
 595             if (E != &e->x.p->arr[2*i+1])
 596                 evalue_free(E);
 597
 598             D->next = next;;
 599
 600             if (!sum)
 601                 break;
 602         }
 603
 604         if (!sum)
 605             break;
 606     }
 607
 608     free(data.s);
 609
 610     if (sum)
 611         reduce_evalue(sum);
 612     return sum;
 613 }
 614
 615 evalue *Bernoulli_sum(Polyhedron *P, Polyhedron *C,
 616                         struct barvinok_options *options)
 617 {
 618     Polyhedron *CA, *D;
 619     evalue e;
 620     evalue *sum;
 621
 622     if (emptyQ(P) || emptyQ(C))
 623         return evalue_zero();
 624
 625     CA = align_context(C, P->Dimension, options->MaxRays);
 626     D = DomainIntersection(P, CA, options->MaxRays);
 627     Domain_Free(CA);
 628
 629     if (emptyQ(D)) {
 630         Domain_Free(D);
 631         return evalue_zero();
 632     }
 633
 634     value_init(e.d);
 635     e.x.p = new_enode(partition, 2, P->Dimension);
 636     EVALUE_SET_DOMAIN(e.x.p->arr[0], D);
 637     evalue_set_si(&e.x.p->arr[1], 1, 1);
 638     sum = Bernoulli_sum_evalue(&e, P->Dimension - C->Dimension, options);
 639     free_evalue_refs(&e);
 640     return sum;
 641 }