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