/* ####################################################################### */ /* PROGOL random sampling */ /* --------------------- */ /* ####################################################################### */ #include #include "progol.h" /* * hibit/1 - finds the highest bit set in an integer by binary chop. * Fails when given 0. */ LONG hibit(n) register unsigned long int n; { register long int hi=INT_SZ-1,lo=0,mid; if (!n) d_error("hibit - given 0 as input"); while (hi != lo) { mid=((hi-lo)>>1)+lo; if (!(n>>(mid+1))) hi=mid; else lo=mid+1; } return(hi); } /* * b_smpair/3 - efficient method for picking a random pair * where both x and y are in the range 0 to sz-1 and * x=(*y=val%sz)) { swap= *y; *y= *x+1; *x= swap; } } /* isqrt(x) - returns the largest integer less than the sqrt of x. * Successive approximation calculated for y^2=x as follows * yi+1 = yi+((x/yi)-yi)/2 = (yi+(x/yi))/2 * where y0 = 2^(log2(x)/2) * Converges in time O(log x) */ LONG isqrt(sq) LONG sq; { LONG oldi,newi; double oldf,newf,sqf=sq; if(sq==0) return(0); else if(sq<0) d_error("isqrt - negative argument"); newi=1<<(LOG2(sq)>>1); newf=newi; do { oldf=newf; newf=(oldf+sqf/oldf)/2; oldi=oldf; newi=newf; } while(oldi!=newi); return(newi); } /* * prencode/3 - given 0<=x to a unique number in the range [1,n*(n-1)/2] * The result k= x((2n-3)-x)/2 + y */ LONG prencode(x,y,n) LONG x,y,n; { if(n>30000) d_error("prencode - set too large"); else if(x<0||y<=x||n<=y) d_error("prencode - bad values"); return(((x*((n<<1)-3-x))>>1)+y); } /* * prdecode/4 - given k=prencode(x,y,n) this function decodes * x and y from k and n. Decoding is done by binary * chop. Upper and lower bounds are computed using * (2n-3-sqrt(abs((2n-3)^2-8k+8(n-1))))/2 <=x<= * (2n-3-sqrt(abs((2n-3)^2-8k+8)))/2. y is * computed using y=(x(x-(2n-3))+2k)/2 */ prdecode(k,n,x,y) LONG k,n,*x,*y; { LONG lo,hi,mid,n2=(n<<1)-3,k2=k<<1,sq=(n2*n2)-(k<<3); if(n>30000) d_error("prdecode - set too large"); else if(k<=0||k>(n*(n-1))>>1) d_error("prdecode - bad values"); lo=(n2-isqrt(sq+((n-1)<<3)))>>1; hi=sq+8>0?(n2-isqrt(sq+8))>>1:n2>>1; do { mid=((hi-lo)>>1)+lo; *y=(mid*(mid-n2)+k2)>>1; if(*y<=mid) hi=mid; else if(*ylo); *x=lo; *y=(lo*(lo-n2)+k2)>>1; } struct otree * ot_create(v,nl) LONG v,nl; { struct otree *result; if (!(result = (struct otree *)PROGOL_CALLOC(1, sizeof(struct otree)))) d_error("ot_create - calloc failure"); result->val=v; result->nleft=nl; result->left=(struct otree *)NULL; result->right=(struct otree *)NULL; return(result); } ot_delete(o) struct otree *o; { if(!o) return; ot_delete(o->left); ot_delete(o->right); PROGOL_CFREE(o,sizeof(struct otree)); } struct otree * ot_ins(v,lo,tree) LONG v,lo; struct otree *tree; { struct otree **ptree= &tree; while(*ptree) if(v<= (*ptree)->nleft) {(*ptree)->nleft--;ptree= &((*ptree)->left);} else {v-=(*ptree)->nleft;lo=(*ptree)->val;ptree= &((*ptree)->right);} *ptree=ot_create(v+lo,v-1); return(tree); } ot_print(indent,tree) LONG indent; struct otree *tree; { LONG cnt; if(tree) { ot_print(indent+1,tree->left); for(cnt=indent;cnt>0;cnt--) printf(" "); printf("<%d,%d>\n",tree->val,tree->nleft); ot_print(indent+1,tree->right); } else { for(cnt=indent;cnt>0;cnt--) printf(" "); printf("_\n"); } } /* * ot_sample(m,n,tree) - randomly choose m different pairs from * a set size n and insert into tree. Return tree. */ struct otree * ot_sample(m,n) LONG m,n; { struct otree *tree=(struct otree *)NULL; if(m>n) m=n; for(;m>0;m--,n--) tree=ot_ins(MYRAND(1,n),0,tree); return(tree); } ot_prprint(tree,n) struct otree *tree; LONG n; { LONG x,y; if(tree) { ot_prprint(tree->left,n); prdecode(tree->val,n,&x,&y); printf("<%d,%d>\n",x,y); ot_prprint(tree->right,n); } } /* * b_bitsum/2 - returns an array giving for each word in the bitset the sum of the * number of bits up to that word. */ BLOCK b_bitsum(b1,total) BLOCK b1; LONG *total; { register LONG cnt,sum=0l; register unsigned long int val; register BLOCK bp1,be1,b2=b_create(B_SIZE(b1)),bp2,be2; BLOCK_LOOP2(bp1,bp2,be1,be2,b1,b2) { *bp2=sum; if (val= *bp1) for (cnt=0l;cnt>=BYTE_SZ) sum += byte_sz[val&0xff]; } *total=sum; return(b2); } /* * b_samprem/4 - carries out a binary split search for the * integer containing the given element and returns * the bit number. */ LONG w_findnth(); LONG b_findnth(elem,b,bitsum) LONG elem; BLOCK b,bitsum; { register LONG lo,hi,mid; hi=B_SIZE(b)-1; b++; bitsum++; for(lo=0l;lo>1); if (elem<=bitsum[mid]) hi=mid-1; else lo=mid; } return((lo<>=BYTE_SZ) if((nset=byte_sz[index=val&0xff])left,b,bsum,result); b_add(b_findnth(tree->val,b,bsum),result); b_sample1(tree->right,b,bsum,result); } } /* normal(x,mu,sg) - gives the tail of the normal distribution * with mean mu and std. dvn. sg. Calculation is done using * the following Taylor series for mu=0, sg=1. * * Integral_0^x p(x)dx = Sum_(i=0)^(inf)((-1^i)x^(2i+1)/((2i+1)i!(2^i))) * sqrt(2 PI) */ #define PI5 3.1416 #define DELTA 0.001 double s_normal(z) double z; { PREDICATE neg; double fy; register double i,sign,y,dy,ti,ifact; y=0.0; dy=RINF; sign=1.0; ifact=1.0; if(neg=(z<0.0)) z= -z; if(z>5.0) y=0.5; else for(i=0.0;dy>=DELTA;) { ti=2.0*i+1.0; dy=(pow(z,ti))/(sqrt(2.0*PI5)*ti*ifact*pow(2.0,i)); y+=sign*dy; sign= -sign; ifact*=(i+=1.0); } fy=(neg?0.5-y:y+0.5); return(fy); } /* s_nmsample - randomly samples according to the density function * of the normal distribution. This is done by transforming * the choice to the [0-1] interval, and carrying out a binary * chop procedure. */ #define GN(z,mu,sg) ((z)*(sg)+(mu)) double s_nmsample(mu,sg) double mu,sg; { double r=RAND,lo= -5.0,hi=5.0,mid; double result; while(lo+DELTA>1l,cnt; sprintf(fileO,"%sO",fileI); if(len<=1l) { sprintf(command,"cp %s %s",fileI,fileO); system(command); return; } if(len<=INMEM) {s_kpermute(fileI,len,fileO,depth+1l); return;} for(cnt=depth;cnt>0l;cnt--) printf(" "); g_message(1l,"Permute %s (%d)",fileI,len); s_randsplit(fileI,file0I,file1I,len0,len); s_permute(file0I,len0,file0O,depth+1l); s_delfile(file0I); s_permute(file1I,len-len0,file1O,depth+1l); s_delfile(file1I); s_appfile(file0O,file1O,fileO); s_delfile(file0O); s_delfile(file1O); } s_randsplit(file,file0,file1,len0,len) STRING file,file0,file1; LONG len0,len; { FILEREC *f,*f0,*f1,*fout; ITEM i; DOUBLE dlen0=len0,dlen=len; char mess[MAXMESS],mess1[MAXMESS]; if (!(f=frecopen(file,"r"))) { printf("[Cannot find %s]\n",file); return; } sprintf(file0,"%s0",file); sprintf(file1,"%s1",file); f0=frecopen(file0,"w"); f1=frecopen(file1,"w"); while((i=ccl_fread(f))!=(ITEM)I_TERM) { if(i==(ITEM)I_ERROR) continue; else { if(drand48()<(dlen0/dlen)) {fout=f0; (dlen0)--;} else fout=f1; dlen--; ccl_fwrite(fout,i,(BIND)NULL); i_fnl(fout); } i_delete(i); } freclose(f); freclose(f0); freclose(f1); } s_delfile(f) STRING f; { char command[MAXMESS]; sprintf(command,"rm %s",f); system(command); } s_appfile(f0,f1,f) STRING f0,f1,f; { char command[MAXMESS]; sprintf(command,"cat %s %s > %s",f0,f1,f); system(command); } LONG s_length(file) STRING file; { LONG res=0l; FILEREC *f; ITEM i; char mess[MAXMESS],mess1[MAXMESS]; if (!(f=frecopen(file,"r"))) { printf("[Cannot find %s]\n",file); return; } g_message(1l,"Counting clauses in <%s>",file); while((i=ccl_fread(f))!=(ITEM)I_TERM) { if(i==(ITEM)I_ERROR) continue; else res++; i_delete(i); } freclose(f); return(res); } /* s_kpermute - Knuth's algorithm for in memory random permutation * of an array. */ s_kpermute(file,len,fileO,depth) STRING file,fileO; LONG len,depth; { ITEM perms=F_EMPTY,i,swap; LONG cnt,r; FILEREC *f,*fout; DOUBLE dlen,dcnt,dr; if (!(f=frecopen(file,"r"))) { printf("[Cannot find %s]\n",file); return; } /* Read clauses into array */ for(cnt=depth;cnt>0l;cnt--) printf(" "); g_message(1l,"IN MEMORY permute %s (%d)",file,len); cnt=0l; while((i=ccl_fread(f))!=(ITEM)I_TERM) { if(i==(ITEM)I_ERROR) continue; else *f_ins(cnt++,perms)=i; } if(cnt!=len) { printf("cnt=%d; len=%d\n",cnt,len); d_error("Miscount of clauses in file"); } freclose(f); fout=frecopen(fileO,"w"); dlen=len; for(cnt=0l;cnt