source/kernel/vector.c

   1 /*
   2     This file is part of PolyLib.
   3
   4     PolyLib is free software: you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation, either version 3 of the License, or
   7     (at your option) any later version.
   8
   9     PolyLib is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with PolyLib.  If not, see <http://www.gnu.org/licenses/>.
  16 */
  17
  18 /* vector.c
  19           COPYRIGHT
  20           Both this software and its documentation are
  21
  22               Copyrighted 1993 by IRISA /Universite de Rennes I - France,
  23               Copyright 1995,1996 by BYU, Provo, Utah
  24                          all rights reserved.
  25
  26 */
  27
  28 #include <stdio.h>
  29 #include <stdlib.h>
  30 #include <polylib/polylib.h>
  31
  32
  33 /* defined by configure script */
  34 /* #define THREAD_SAFE_POLYLIB */
  35
  36 #ifdef MAC_OS
  37   #define abs __abs
  38 #endif
  39
  40 /*
  41  * Compute n!
  42  */
  43 void Factorial(int n, Value *fact) {
  44
  45   int i;
  46   Value tmp;
  47
  48   value_init(tmp);
  49
  50   value_set_si(*fact,1);
  51   for (i=1;i<=n;i++) {
  52     value_set_si(tmp,i);
  53     value_multiply(*fact,*fact,tmp);
  54   }
  55   value_clear(tmp);
  56 } /* Factorial */
  57
  58 /*
  59  * Compute n!/(p!(n-p)!)
  60  */
  61 void Binomial(int n, int p, Value *result) {
  62
  63   int i;
  64   Value tmp;
  65   Value f;
  66
  67   value_init(tmp);value_init(f);
  68
  69   if (n-p<p)
  70     p=n-p;
  71   if (p!=0) {
  72     value_set_si(*result,(n-p+1));
  73     for (i=n-p+2;i<=n;i++) {
  74       value_set_si(tmp,i);
  75       value_multiply(*result,*result,tmp);
  76     }
  77     Factorial(p,&f);
  78     value_division(*result,*result,f);
  79   }
  80   else
  81     value_set_si(*result,1);
  82   value_clear(f);value_clear(tmp);
  83 } /* Binomial */
  84
  85 /*
  86  * Return the number of ways to choose 'b' items from 'a' items, that is,
  87  * return a!/(b!(a-b)!)
  88  */
  89 void CNP(int a,int b, Value *result) {
  90
  91   int i;
  92   Value tmp;
  93   value_init(tmp);
  94
  95   value_set_si(*result,1);
  96
  97   /* If number of items is less than the number to be choosen, return 1 */
  98   if(a <= b) {
  99     value_clear(tmp);
 100     return;
 101   }
 102   for(i=a;i>b;--i) {
 103     value_set_si(tmp,i);
 104     value_multiply(*result,*result,tmp);
 105   }
 106   for(i=1;i<=(a-b);++i) {
 107     value_set_si(tmp,i);
 108     value_division(*result,*result,tmp);
 109   }
 110   value_clear(tmp);
 111 } /* CNP */
 112
 113 /*
 114  * Compute GCD of 'a' and 'b'
 115  */
 116 void Gcd(Value a,Value b,Value *result) {
 117
 118   Value acopy, bcopy;
 119
 120   value_init(acopy);
 121   value_init(bcopy);
 122   value_assign(acopy,a);
 123   value_assign(bcopy,b);
 124   while(value_notzero_p(acopy)) {
 125     value_modulus(*result,bcopy,acopy);
 126     value_assign(bcopy,acopy);
 127     value_assign(acopy,*result);
 128   }
 129   value_absolute(*result,bcopy);
 130   value_clear(acopy);
 131   value_clear(bcopy);
 132 } /* Gcd */
 133
 134 /*
 135  * Return the smallest component index in 'p' whose value is non-zero
 136  */
 137 int First_Non_Zero(Value *p,unsigned length) {
 138
 139   Value *cp;
 140   int i;
 141
 142   cp = p;
 143   for (i=0;i<length;i++) {
 144     if (value_notzero_p(*cp))
 145       break;
 146     cp++;
 147   }
 148   return((i==length) ? -1 : i );
 149 } /* First_Non_Zero */
 150
 151 /*
 152  * Allocate memory space for Vector
 153  */
 154 Vector *Vector_Alloc(unsigned length) {
 155
 156   int i;
 157   Vector *vector;
 158
 159   vector = (Vector *)malloc(sizeof(Vector));
 160   if (!vector) {
 161     errormsg1("Vector_Alloc", "outofmem", "out of memory space");
 162     return 0;
 163   }
 164   vector->Size=length;
 165   vector->p=(Value *)malloc(length * sizeof(Value));
 166   if (!vector->p) {
 167     errormsg1("Vector_Alloc", "outofmem", "out of memory space");
 168     free(vector);
 169     return 0;
 170   }
 171   for(i=0;i<length;i++)
 172     value_init(vector->p[i]);
 173   return vector;
 174 } /* Vector_Alloc */
 175
 176 /*
 177  * Free the memory space occupied by Vector
 178  */
 179 void Vector_Free(Vector *vector) {
 180
 181   int i;
 182
 183   if (!vector) return;
 184   for(i=0;i<vector->Size;i++)
 185     value_clear(vector->p[i]);
 186   free(vector->p);
 187   free(vector);
 188 } /* Vector_Free */
 189
 190 /*
 191  * Print the contents of a Vector
 192  */
 193 void Vector_Print(FILE *Dst, const char *Format, Vector *vector)
 194 {
 195   int i;
 196   Value *p;
 197   unsigned length;
 198
 199   fprintf(Dst, "%d\n", length=vector->Size);
 200   p = vector->p;
 201   for (i=0;i<length;i++) {
 202     if (Format) {
 203       value_print(Dst,Format,*p++);
 204     }
 205     else {
 206       value_print(Dst,P_VALUE_FMT,*p++);
 207     }
 208   }
 209   fprintf(Dst, "\n");
 210 } /* Vector_Print */
 211
 212 /*
 213  * Read the contents of a Vector
 214  */
 215 Vector *Vector_Read() {
 216
 217   Vector *vector;
 218   unsigned length;
 219   int i;
 220   char *s;
 221   Value *p;
 222
 223   scanf("%d", &length);
 224   vector = Vector_Alloc(length);
 225   if (!vector) {
 226     errormsg1("Vector_Read", "outofmem", "out of memory space");
 227     return 0;
 228   }
 229   p = vector->p;
 230   for (i=0;i<length;i++)
 231         {
 232                 int r=scanf("%ms",&s);
 233                 if(r!=1)
 234                 {
 235                         errormsg1("Vector_Read", "hi, Mathematica", "scanf at line 231 failed");
 236                         return 0;
 237                 }
 238                 value_read(*(p++),s);
 239                 free(s);
 240         }
 241   return vector;
 242 } /* Vector_Read */
 243
 244 /*
 245  * Assign 'n' to each component of Vector 'p'
 246  */
 247 void Vector_Set(Value *p,int n,unsigned length) {
 248
 249   Value *cp;
 250   int i;
 251
 252   cp = p;
 253   for (i=0;i<length;i++) {
 254     value_set_si(*cp,n);
 255     cp++;
 256   }
 257   return;
 258 } /* Vector_Set */
 259
 260 /*
 261  * Exchange the components of the vectors 'p1' and 'p2' of length 'length'
 262  */
 263 void Vector_Exchange(Value *p1, Value *p2, unsigned length) {
 264
 265   int i;
 266
 267   for(i=0;i<length;i++) {
 268     value_swap(p1[i],p2[i]);
 269   }
 270   return;
 271 }
 272
 273 /*
 274  * Copy Vector 'p1' to Vector 'p2'
 275  */
 276 void Vector_Copy(Value *p1,Value *p2,unsigned length) {
 277
 278   int i;
 279   Value *cp1, *cp2;
 280
 281   cp1 = p1;
 282   cp2 = p2;
 283
 284   for(i=0;i<length;i++)
 285     value_assign(*cp2++,*cp1++);
 286
 287   return;
 288 }
 289
 290 /*
 291  * Add two vectors 'p1' and 'p2' and store the result in 'p3'
 292  */
 293 void Vector_Add(Value *p1,Value *p2,Value *p3,unsigned length) {
 294
 295   Value *cp1, *cp2, *cp3;
 296   int i;
 297
 298   cp1=p1;
 299   cp2=p2;
 300   cp3=p3;
 301   for (i=0;i<length;i++) {
 302
 303     /* *cp3++ = *cp1++ + *cp2++ */
 304     value_addto(*cp3,*cp1,*cp2);
 305     cp1++; cp2++; cp3++;
 306   }
 307 } /* Vector_Add */
 308
 309 /*
 310  * Subtract two vectors 'p1' and 'p2' and store the result in 'p3'
 311  */
 312 void Vector_Sub(Value *p1,Value *p2,Value *p3,unsigned length) {
 313
 314   Value *cp1, *cp2, *cp3;
 315   int i;
 316
 317   cp1=p1;
 318   cp2=p2;
 319   cp3=p3;
 320   for (i=0;i<length;i++) {
 321
 322     /* *cp3++= *cp1++ - *cp2++ */
 323     value_subtract(*cp3,*cp1,*cp2);
 324     cp1++; cp2++; cp3++;
 325   }
 326 } /* Vector_Sub */
 327
 328 /*
 329  * Compute Bitwise OR of Vectors 'p1' and 'p2' and store in 'p3'
 330  */
 331 void Vector_Or(Value *p1,Value *p2,Value *p3,unsigned length) {
 332
 333   Value *cp1, *cp2, *cp3;
 334   int i;
 335
 336   cp1=p1;
 337   cp2=p2;
 338   cp3=p3;
 339   for (i=0;i<length;i++) {
 340
 341     /* *cp3++=*cp1++ | *cp2++ */
 342     value_orto(*cp3,*cp1,*cp2);
 343     cp1++; cp2++; cp3++;
 344   }
 345 } /* Vector_Or */
 346
 347 /*
 348  * Scale Vector 'p1' lambda times and store in 'p2'
 349  */
 350 void Vector_Scale(Value *p1,Value *p2,Value lambda,unsigned length) {
 351
 352   Value *cp1, *cp2;
 353   int i;
 354
 355   cp1=p1;
 356   cp2=p2;
 357   for (i=0;i<length;i++) {
 358
 359     /* *cp2++=*cp1++ * lambda */
 360     value_multiply(*cp2,*cp1,lambda);
 361     cp1++; cp2++;
 362   }
 363 } /* Vector_Scale */
 364
 365 /*
 366  * Antiscale Vector 'p1' by lambda and store in 'p2'
 367  * Assumes all elements of 'p1' are divisble by lambda.
 368  */
 369 void Vector_AntiScale(Value *p1, Value *p2, Value lambda, unsigned length)
 370 {
 371   int i;
 372
 373   for (i = 0; i < length; i++)
 374     value_divexact(p2[i], p1[i], lambda);
 375 } /* Vector_AntiScale */
 376
 377 /*
 378  * Puts negative of 'p1' in 'p2'
 379  */
 380 void Vector_Oppose(Value *p1, Value *p2, unsigned len)
 381 {
 382   unsigned i;
 383
 384   for (i = 0; i < len; ++i)
 385     value_oppose(p2[i], p1[i]);
 386 }
 387
 388 /*
 389  * Return the inner product of the two Vectors 'p1' and 'p2'
 390  */
 391 void Inner_Product(Value *p1, Value *p2, unsigned length, Value *ip)
 392 {
 393   int i;
 394
 395   if (length != 0)
 396     value_multiply(*ip, p1[0], p2[0]);
 397   else
 398     value_set_si(*ip, 0);
 399   for(i = 1; i < length; i++)
 400     value_addmul(*ip, p1[i], p2[i]);
 401 } /* Inner_Product */
 402
 403 /*
 404  * Return the maximum of the components of 'p'
 405  */
 406 void Vector_Max(Value *p,unsigned length, Value *max) {
 407
 408   Value *cp;
 409   int i;
 410
 411   cp=p;
 412   value_assign(*max,*cp);
 413   cp++;
 414   for (i=1;i<length;i++) {
 415     value_maximum(*max,*max,*cp);
 416     cp++;
 417   }
 418 } /* Vector_Max */
 419
 420 /*
 421  * Return the minimum of the components of Vector 'p'
 422  */
 423 void Vector_Min(Value *p,unsigned length,Value *min) {
 424
 425   Value *cp;
 426   int i;
 427
 428   cp=p;
 429   value_assign(*min,*cp);
 430   cp++;
 431   for (i=1;i<length;i++) {
 432     value_minimum(*min,*min,*cp);
 433     cp++;
 434   }
 435   return;
 436 } /* Vector_Min */
 437
 438 /*
 439  * Given Vectors 'p1' and 'p2', return Vector 'p3' = lambda * p1 + mu * p2.
 440  */
 441 void Vector_Combine(Value *p1, Value *p2, Value *p3, Value lambda, Value mu,
 442                     unsigned length)
 443 {
 444   Value tmp;
 445   int i;
 446
 447   value_init(tmp);
 448   for (i = 0; i < length; i++) {
 449     value_multiply(tmp, lambda, p1[i]);
 450     value_addmul(tmp, mu, p2[i]);
 451     value_assign(p3[i], tmp);
 452   }
 453   value_clear(tmp);
 454   return;
 455 } /* Vector_Combine */
 456
 457 /*
 458  * Return 1 if 'Vec1' equals 'Vec2', otherwise return 0
 459  */
 460 int Vector_Equal(Value *Vec1, Value *Vec2, unsigned n)
 461 {
 462   int i;
 463
 464   for (i = 0; i < n; ++i)
 465     if (value_ne(Vec1[i], Vec2[i]))
 466       return 0;
 467
 468   return 1;
 469 } /* Vector_Equal */
 470
 471 /*
 472  * Return the component of 'p' with minimum non-zero absolute value. 'index'
 473  * points to the component index that has the minimum value. If no such value
 474  * and index is found, Value 1 is returned.
 475  */
 476 void Vector_Min_Not_Zero(Value *p,unsigned length,int *index,Value *min)
 477 {
 478   Value aux;
 479   int i;
 480
 481
 482   i = First_Non_Zero(p, length);
 483   if (i == -1) {
 484     value_set_si(*min,1);
 485     return;
 486   }
 487   *index = i;
 488   value_absolute(*min, p[i]);
 489   value_init(aux);
 490   for (i = i+1; i < length; i++) {
 491     if (value_zero_p(p[i]))
 492       continue;
 493     value_absolute(aux, p[i]);
 494     if (value_lt(aux,*min)) {
 495       value_assign(*min,aux);
 496       *index = i;
 497     }
 498   }
 499   value_clear(aux);
 500 } /* Vector_Min_Not_Zero */
 501
 502 /*
 503  * Return the GCD of components of Vector 'p'
 504  */
 505 void Vector_Gcd(Value *p,unsigned length,Value *min) {
 506
 507   Value *q,*cq, *cp;
 508   int i, Not_Zero, Index_Min=0;
 509
 510   q  = (Value *)malloc(length*sizeof(Value));
 511
 512   /* Initialize all the 'Value' variables */
 513   for(i=0;i<length;i++)
 514     value_init(q[i]);
 515
 516   /* 'cp' points to vector 'p' and cq points to vector 'q' that holds the */
 517   /* absolute value of elements of vector 'p'.                            */
 518   cp=p;
 519   for (cq = q,i=0;i<length;i++) {
 520     value_absolute(*cq,*cp);
 521     cq++;
 522     cp++;
 523   }
 524   do {
 525     Vector_Min_Not_Zero(q,length,&Index_Min,min);
 526
 527     /* if (*min != 1) */
 528     if (value_notone_p(*min)) {
 529
 530       cq=q;
 531       Not_Zero=0;
 532       for (i=0;i<length;i++,cq++)
 533         if (i!=Index_Min) {
 534
 535           /* Not_Zero |= (*cq %= *min) */
 536           value_modulus(*cq,*cq,*min);
 537           Not_Zero |= value_notzero_p(*cq);
 538         }
 539     }
 540     else
 541       break;
 542   } while (Not_Zero);
 543
 544   /* Clear all the 'Value' variables */
 545   for(i=0;i<length;i++)
 546     value_clear(q[i]);
 547   free(q);
 548 } /* Vector_Gcd */
 549
 550 /*
 551  * Given vectors 'p1', 'p2', and a pointer to a function returning 'Value type,
 552  * compute p3[i] = f(p1[i],p2[i]).
 553  */
 554 void Vector_Map(Value *p1,Value *p2,Value *p3,unsigned length,
 555                 Value *(*f)(Value,Value))
 556 {
 557   Value *cp1, *cp2, *cp3;
 558   int i;
 559
 560   cp1=p1;
 561   cp2=p2;
 562   cp3=p3;
 563   for(i=0;i<length;i++) {
 564     value_assign(*cp3,*(*f)(*cp1, *cp2));
 565     cp1++; cp2++; cp3++;
 566   }
 567   return;
 568 } /* Vector_Map */
 569
 570 /*
 571  * Reduce a vector by dividing it by GCD. There is no restriction on
 572  * components of Vector 'p'. Making the last element positive is *not* OK
 573  * for equalities.
 574  */
 575 void Vector_Normalize(Value *p,unsigned length) {
 576
 577   Value gcd;
 578   int i;
 579
 580   value_init(gcd);
 581
 582   Vector_Gcd(p,length,&gcd);
 583
 584   if (value_notone_p(gcd))
 585     Vector_AntiScale(p, p, gcd, length);
 586
 587   value_clear(gcd);
 588 } /* Vector_Normalize */
 589
 590 /*
 591  * Reduce a vector by dividing it by GCD and making sure its pos-th
 592  * element is positive.
 593  */
 594 void Vector_Normalize_Positive(Value *p,int length,int pos) {
 595
 596   Value gcd;
 597   int i;
 598
 599   value_init(gcd);
 600   Vector_Gcd(p,length,&gcd);
 601   if (value_neg_p(p[pos]))
 602     value_oppose(gcd,gcd);
 603   if(value_notone_p(gcd))
 604     Vector_AntiScale(p, p, gcd, length);
 605   value_clear(gcd);
 606 } /* Vector_Normalize_Positive */
 607
 608 /*
 609  * Reduce 'p' by operating binary function on its components successively
 610  */
 611 void Vector_Reduce(Value *p,unsigned length,void(*f)(Value,Value *),Value *r) {
 612
 613   Value *cp;
 614   int i;
 615
 616   cp = p;
 617   value_assign(*r,*cp);
 618   for(i=1;i<length;i++) {
 619     cp++;
 620     (*f)(*cp,r);
 621   }
 622 } /* Vector_Reduce */
 623
 624 /*
 625  * Sort the components of a Vector 'vector' using Heap Sort.
 626  */
 627 void Vector_Sort(Value *vector,unsigned n) {
 628
 629   int i, j;
 630   Value temp;
 631   Value *current_node=(Value *)0;
 632   Value *left_son,*right_son;
 633
 634   value_init(temp);
 635
 636   for (i=(n-1)/2;i>=0;i--) {
 637
 638     /* Phase 1 : build the heap */
 639     j=i;
 640     value_assign(temp,*(vector+i));
 641
 642     /* While not a leaf */
 643     while (j<=(n-1)/2) {
 644       current_node = vector+j;
 645       left_son = vector+(j<<1)+1;
 646
 647       /* If only one son */
 648       if ((j<<1)+2>=n) {
 649         if (value_lt(temp,*left_son)) {
 650           value_assign(*current_node,*left_son);
 651           j=(j<<1)+1;
 652         }
 653         else
 654           break;
 655       }
 656       else {
 657
 658         /* If two sons */
 659         right_son=left_son+1;
 660         if (value_lt(*right_son,*left_son)) {
 661           if (value_lt(temp,*left_son)) {
 662             value_assign(*current_node,*left_son);
 663             j=(j<<1)+1;
 664           }
 665           else
 666             break;
 667         }
 668         else {
 669           if (value_lt(temp,*right_son)) {
 670             value_assign(*current_node,*right_son );
 671             j=(j<<1)+2;
 672           }
 673           else
 674             break;
 675         }
 676       }
 677     }
 678     value_assign(*current_node,temp);
 679   }
 680   for(i=n-1;i>0;i--) {
 681
 682     /* Phase 2 : sort the heap */
 683     value_assign(temp, *(vector+i));
 684     value_assign(*(vector+i),*vector);
 685     j=0;
 686
 687     /* While not a leaf */
 688     while (j<i/2) {
 689       current_node=vector+j;
 690       left_son=vector+(j<<1)+1;
 691
 692       /* If only one son */
 693       if ((j<<1)+2>=i) {
 694         if (value_lt(temp,*left_son)) {
 695           value_assign(*current_node,*left_son);
 696           j=(j<<1)+1;
 697         }
 698         else
 699           break;
 700       }
 701       else {
 702
 703         /* If two sons */
 704         right_son=left_son+1;
 705         if (value_lt(*right_son,*left_son)) {
 706           if (value_lt(temp,*left_son)) {
 707             value_assign(*current_node,*left_son);
 708             j=(j<<1)+1;
 709           }
 710           else
 711             break;
 712         }
 713         else {
 714           if (value_lt(temp,*right_son)) {
 715             value_assign(*current_node,*right_son );
 716             j=(j<<1)+2;
 717           }
 718           else
 719             break;
 720         }
 721       }
 722     }
 723     value_assign(*current_node,temp);
 724   }
 725   value_clear(temp);
 726   return;
 727 } /* Vector_Sort */
 728
 729 /*
 730  * Replaces constraint a x >= c by x >= ceil(c/a)
 731  * where "a" is a common factor in the coefficients
 732  * old is the constraint; v points to an initialized
 733  * value that this procedure can use.
 734  * Return non-zero if something changed.
 735  * Result is placed in newp.
 736  */
 737 int ConstraintSimplify(Value *old, Value *newp, int len, Value* v)
 738 {
 739     /* first remove common factor of all coefficients (including "c") */
 740     Vector_Gcd(old+1, len - 1, v);
 741     if (value_notone_p(*v))
 742         Vector_AntiScale(old+1, newp+1, *v, len-1);
 743
 744     Vector_Gcd(old+1, len - 2, v);
 745
 746     if (value_one_p(*v))
 747         return 0;
 748
 749     Vector_AntiScale(old+1, newp+1, *v, len-2);
 750     value_pdivision(newp[len-1], old[len-1], *v);
 751     return 1;
 752 }
 753
 754 int Vector_IsZero(Value * v, unsigned length) {
 755   unsigned i;
 756   if (value_notzero_p(v[0])) return 0;
 757   else {
 758     value_set_si(v[0], 1);
 759     for (i=length-1; value_zero_p(v[i]); i--);
 760     value_set_si(v[0], 0);
 761     return (i==0);
 762   }
 763 }
 764
 765 typedef struct {
 766   Value *p;
 767   int   size;
 768 } cache_holder;
 769 #define MAX_CACHE_SIZE 20
 770
 771
 772 /************************************************/
 773 /** Vincent's patch for thread safe value cache */
 774 /** each thread has it's own cache and size.    */
 775 /** 02/2012                                     */
 776 /************************************************/
 777
 778 #ifdef THREAD_SAFE_POLYLIB
 779 #include <pthread.h>
 780 #include <assert.h>
 781
 782
 783 static pthread_once_t once_cache = PTHREAD_ONCE_INIT;
 784 static pthread_key_t cache_key;
 785 /* cache_size is stored in the last+1 cache position */
 786 #define cache_size (cache[MAX_CACHE_SIZE].size)
 787
 788 static cache_holder *allocate_local_cache(void)
 789 {
 790         cache_holder *cache;
 791         cache = malloc( sizeof(cache_holder)*(MAX_CACHE_SIZE+1) );
 792         assert( cache!=NULL );
 793         cache_size = 0;
 794         assert( pthread_setspecific( cache_key, cache ) == 0 );
 795         return( cache );
 796 }
 797 static void free_local_cache(void *c)
 798 {
 799         free(c);
 800         assert( pthread_setspecific(cache_key, NULL) == 0 );
 801 }
 802 static void init_value_caches(void)
 803 {
 804         pthread_key_create(&cache_key, free_local_cache);
 805 }
 806
 807 #else
 808 cache_holder cache[MAX_CACHE_SIZE];
 809 static int cache_size = 0;
 810 #endif // THREAD_SAFE_POLYLIB
 811
 812
 813 Value* value_alloc(int want, int *got)
 814 {
 815     int i;
 816     Value *p;
 817
 818 #ifdef THREAD_SAFE_POLYLIB
 819         assert(pthread_once(&once_cache, init_value_caches) == 0);
 820         cache_holder *cache;
 821         if( (cache = pthread_getspecific( cache_key )) == NULL )
 822                 cache = allocate_local_cache();
 823 #endif // THREAD_SAFE_POLYLIB
 824
 825     if (cache_size) {
 826       int best=0;
 827       for (i = 0; i < cache_size; ++i) {
 828         if (cache[i].size >= want) {
 829           Value *p = cache[i].p;
 830           *got = cache[i].size;
 831           if (--cache_size != i)
 832             cache[i] = cache[cache_size];
 833           Vector_Set(p, 0, want);
 834           return p;
 835         }
 836         if (cache[i].size > cache[best].size)
 837           best = i;
 838       }
 839
 840       p = (Value *)realloc(cache[best].p, want * sizeof(Value));
 841       *got = cache[best].size;
 842       if (--cache_size != best)
 843         cache[best] = cache[cache_size];
 844       Vector_Set(p, 0, *got);
 845     }
 846     else {
 847       p = (Value *)malloc(want * sizeof(Value));
 848       *got = 0;
 849     }
 850
 851     if (!p)
 852       return p;
 853
 854     for (i = *got; i < want; ++i)
 855       value_init(p[i]);
 856     *got = want;
 857
 858     return p;
 859 }
 860
 861 void value_free(Value *p, int size)
 862 {
 863     int i;
 864
 865 #ifdef THREAD_SAFE_POLYLIB
 866 /* suppose alloc before free :) */
 867 //      assert(pthread_once(&once_cache, init_value_caches) == 0);
 868         cache_holder *cache;
 869 //      if( (cache = pthread_getspecific( cache_key )) == NULL )
 870 //              cache = allocate_local_cache();
 871         assert( (cache = pthread_getspecific( cache_key )) != NULL );
 872 #endif // THREAD_SAFE_POLYLIB
 873
 874     if (cache_size < MAX_CACHE_SIZE) {
 875       cache[cache_size].p = p;
 876       cache[cache_size].size = size;
 877       ++cache_size;
 878       return;
 879     }
 880
 881     for (i=0; i < size; i++)
 882       value_clear(p[i]);
 883     free(p);
 884 }
 885