ref: 71c9bbff02edd1b8818b782ebf498bd67e183027
parent: 26c9e1452ae9941e49791dc3acf5dd8d487e876f
author: Timothy B. Terriberry <[email protected]>
date: Tue Jan 1 02:18:06 EST 2008
Updated pulse coding to simpler (slightly faster) code included with http://people.xiph.org/~tterribe/notes/cwrs.html Removed dead code.
--- a/libcelt/cwrs.c
+++ b/libcelt/cwrs.c
@@ -1,21 +1,20 @@
-/* (C) 2007 Timothy Terriberry
-*/
+/* (C) 2007 Timothy B. Terriberry */
/*
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
-
+
- Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
-
+
- Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
-
+
- Neither the name of the Xiph.org Foundation nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
-
+
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
@@ -28,11 +27,7 @@
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
-
-
-/*#include <stdio.h>*/
#include <stdlib.h>
-
#include "cwrs.h"
/*Returns the numer of ways of choosing _m elements from a set of size _n with
@@ -47,29 +42,7 @@
}
return c[_n][_m];
}
-
#else
-
-/*Returns the greatest common divisor of _a and _b.*/
-static unsigned gcd(unsigned _a,unsigned _b){- unsigned r;
- while(_b){- r=_a%_b;
- _a=_b;
- _b=r;
- }
- return _a;
-}
-
-/*Returns _a*b/_d, under the assumption that the result is an integer, avoiding
- overflow.
- It is assumed, but not required, that _b is smaller than _a.*/
-static unsigned umuldiv(unsigned _a,unsigned _b,unsigned _d){- unsigned d;
- d=gcd(_b,_d);
- return (_a/(_d/d))*(_b/d);
-}
-
unsigned ncwrs(int _n,int _m){unsigned ret;
unsigned f;
@@ -83,13 +56,8 @@
d=1;
for(i=1;i<=_m;i++){ret+=f*d<<i;
-#if 0
- f=umuldiv(f,_n-i,i+1);
- d=umuldiv(d,_m-i,i);
-#else
f=(f*(_n-i))/(i+1);
d=(d*(_m-i))/i;
-#endif
}
return ret;
}
@@ -96,27 +64,22 @@
#endif
celt_uint64_t ncwrs64(int _n,int _m){- celt_uint64_t ret;
- celt_uint64_t f;
- celt_uint64_t d;
- int i;
- if(_n<0||_m<0)return 0;
- if(_m==0)return 1;
- if(_n==0)return 0;
- ret=0;
- f=_n;
- d=1;
- for(i=1;i<=_m;i++){- ret+=f*d<<i;
-#if 0
- f=umuldiv(f,_n-i,i+1);
- d=umuldiv(d,_m-i,i);
-#else
- f=(f*(_n-i))/(i+1);
- d=(d*(_m-i))/i;
-#endif
- }
- return ret;
+ celt_uint64_t ret;
+ celt_uint64_t f;
+ celt_uint64_t d;
+ int i;
+ if(_n<0||_m<0)return 0;
+ if(_m==0)return 1;
+ if(_n==0)return 0;
+ ret=0;
+ f=_n;
+ d=1;
+ for(i=1;i<=_m;i++){+ ret+=f*d<<i;
+ f=(f*(_n-i))/(i+1);
+ d=(d*(_m-i))/i;
+ }
+ return ret;
}
/*Returns the _i'th combination of _m elements chosen from a set of size _n
@@ -124,32 +87,29 @@
_x: Returns the combination with elements sorted in ascending order.
_s: Returns the associated sign bits.*/
void cwrsi(int _n,int _m,unsigned _i,int *_x,int *_s){- unsigned pn;
- int j;
- int k;
- pn=ncwrs(_n-1,_m);
+ int j;
+ int k;
for(k=j=0;k<_m;k++){- unsigned pp;
+ unsigned pn;
unsigned p;
unsigned t;
- pp=0;
- p=ncwrs(_n-j,_m-k)-pn;
+ p=ncwrs(_n-j,_m-k-1);
+ pn=ncwrs(_n-j-1,_m-k-1);
+ p+=pn;
if(k>0){t=p>>1;
if(t<=_i||_s[k-1])_i+=t;
}
- pn=ncwrs(_n-j-1,_m-k-1);
while(p<=_i){- pp=p;
+ _i-=p;
j++;
- p+=pn;
+ p=pn;
pn=ncwrs(_n-j-1,_m-k-1);
p+=pn;
}
- t=p-pp>>1;
- _s[k]=_i-pp>=t;
+ t=p>>1;
+ _s[k]=_i>=t;
_x[k]=j;
- _i-=pp;
if(_s[k])_i-=t;
}
}
@@ -159,66 +119,59 @@
_x: The combination with elements sorted in ascending order.
_s: The associated sign bits.*/
unsigned icwrs(int _n,int _m,const int *_x,const int *_s){- unsigned pn;
unsigned i;
int j;
int k;
i=0;
- pn=ncwrs(_n-1,_m);
for(k=j=0;k<_m;k++){- unsigned pp;
+ unsigned pn;
unsigned p;
- pp=0;
- p=ncwrs(_n-j,_m-k)-pn;
- if(k>0)p>>=1;
+ p=ncwrs(_n-j,_m-k-1);
pn=ncwrs(_n-j-1,_m-k-1);
+ p+=pn;
+ if(k>0)p>>=1;
while(j<_x[k]){- pp=p;
+ i+=p;
j++;
- p+=pn;
+ p=pn;
pn=ncwrs(_n-j-1,_m-k-1);
p+=pn;
}
- i+=pp;
- if((k==0||_x[k]!=_x[k-1])&&_s[k])i+=p-pp>>1;
+ if((k==0||_x[k]!=_x[k-1])&&_s[k])i+=p>>1;
}
return i;
}
-
/*Returns the _i'th combination of _m elements chosen from a set of size _n
with associated sign bits.
_x: Returns the combination with elements sorted in ascending order.
_s: Returns the associated sign bits.*/
void cwrsi64(int _n,int _m,celt_uint64_t _i,int *_x,int *_s){- celt_uint64_t pn;
- int j;
- int k;
- pn=ncwrs64(_n-1,_m);
- for(k=j=0;k<_m;k++){- celt_uint64_t pp;
- celt_uint64_t p;
- celt_uint64_t t;
- pp=0;
- p=ncwrs64(_n-j,_m-k)-pn;
- if(k>0){- t=p>>1;
- if(t<=_i||_s[k-1])_i+=t;
- }
+ int j;
+ int k;
+ for(k=j=0;k<_m;k++){+ celt_uint64_t pn;
+ celt_uint64_t p;
+ celt_uint64_t t;
+ p=ncwrs64(_n-j,_m-k-1);
+ pn=ncwrs64(_n-j-1,_m-k-1);
+ p+=pn;
+ if(k>0){+ t=p>>1;
+ if(t<=_i||_s[k-1])_i+=t;
+ }
+ while(p<=_i){+ _i-=p;
+ j++;
+ p=pn;
pn=ncwrs64(_n-j-1,_m-k-1);
- while(p<=_i){- pp=p;
- j++;
- p+=pn;
- pn=ncwrs64(_n-j-1,_m-k-1);
- p+=pn;
- }
- t=p-pp>>1;
- _s[k]=_i-pp>=t;
- _x[k]=j;
- _i-=pp;
- if(_s[k])_i-=t;
- }
+ p+=pn;
+ }
+ t=p>>1;
+ _s[k]=_i>=t;
+ _x[k]=j;
+ if(_s[k])_i-=t;
+ }
}
/*Returns the index of the given combination of _m elements chosen from a set
@@ -226,33 +179,29 @@
_x: The combination with elements sorted in ascending order.
_s: The associated sign bits.*/
celt_uint64_t icwrs64(int _n,int _m,const int *_x,const int *_s){- celt_uint64_t pn;
- celt_uint64_t i;
- int j;
- int k;
- i=0;
- pn=ncwrs64(_n-1,_m);
- for(k=j=0;k<_m;k++){- celt_uint64_t pp;
- celt_uint64_t p;
- pp=0;
- p=ncwrs64(_n-j,_m-k)-pn;
- if(k>0)p>>=1;
+ celt_uint64_t i;
+ int j;
+ int k;
+ i=0;
+ for(k=j=0;k<_m;k++){+ celt_uint64_t pn;
+ celt_uint64_t p;
+ p=ncwrs64(_n-j,_m-k-1);
+ pn=ncwrs64(_n-j-1,_m-k-1);
+ p+=pn;
+ if(k>0)p>>=1;
+ while(j<_x[k]){+ i+=p;
+ j++;
+ p=pn;
pn=ncwrs64(_n-j-1,_m-k-1);
- while(j<_x[k]){- pp=p;
- j++;
- p+=pn;
- pn=ncwrs64(_n-j-1,_m-k-1);
- p+=pn;
- }
- i+=pp;
- if((k==0||_x[k]!=_x[k-1])&&_s[k])i+=p-pp>>1;
- }
- return i;
+ p+=pn;
+ }
+ if((k==0||_x[k]!=_x[k-1])&&_s[k])i+=p>>1;
+ }
+ return i;
}
-
/*Converts a combination _x of _m unit pulses with associated sign bits _s into
a pulse vector _y of length _n.
_y: Returns the vector of pulses.
@@ -292,7 +241,8 @@
}
}
-/*
+#if 0
+#include <stdio.h>
#define NMAX (10)
#define MMAX (9)
@@ -333,53 +283,53 @@
printf("\n");}
}
- return 0;
+ return -1;
}
-*/
+#endif
-/*
+#if 0
#include <stdio.h>
#define NMAX (32)
#define MMAX (16)
int main(int _argc,char **_argv){- int n;
- for(n=0;n<=NMAX;n+=3){- int m;
- for(m=0;m<=MMAX;m++){- celt_uint64_t nc;
- celt_uint64_t i;
- nc=ncwrs64(n,m);
- printf("%d/%d: %llu",n,m, nc);- for(i=0;i<nc;i+=100000){- int x[MMAX];
- int s[MMAX];
- int x2[MMAX];
- int s2[MMAX];
- int y[NMAX];
- int j;
- int k;
- cwrsi64(n,m,i,x,s);
- //printf("%llu of %llu:",i,nc);- for(k=0;k<m;k++){- //printf(" %c%i",k>0&&x[k]==x[k-1]?' ':s[k]?'-':'+',x[k]);- }
- //printf(" ->");- if(icwrs64(n,m,x,s)!=i){- fprintf(stderr,"Combination-index mismatch.\n");
- }
- comb2pulse(n,m,y,x,s);
- //for(j=0;j<n;j++)printf(" %c%i",y[j]?y[j]<0?'-':'+':' ',abs(y[j]));- //printf("\n");- pulse2comb(n,m,x2,s2,y);
- for(k=0;k<m;k++)if(x[k]!=x2[k]||s[k]!=s2[k]){- fprintf(stderr,"Pulse-combination mismatch.\n");
- break;
- }
- }
- printf("\n");+ int n;
+ for(n=0;n<=NMAX;n+=3){+ int m;
+ for(m=0;m<=MMAX;m++){+ celt_uint64_t nc;
+ celt_uint64_t i;
+ nc=ncwrs64(n,m);
+ printf("%d/%d: %llu",n,m, nc);+ for(i=0;i<nc;i+=100000){+ int x[MMAX];
+ int s[MMAX];
+ int x2[MMAX];
+ int s2[MMAX];
+ int y[NMAX];
+ int j;
+ int k;
+ cwrsi64(n,m,i,x,s);
+ /*printf("%llu of %llu:",i,nc);+ for(k=0;k<m;k++){+ printf(" %c%i",k>0&&x[k]==x[k-1]?' ':s[k]?'-':'+',x[k]);+ }
+ printf(" ->");*/+ if(icwrs64(n,m,x,s)!=i){+ fprintf(stderr,"Combination-index mismatch.\n");
+ }
+ comb2pulse(n,m,y,x,s);
+ /*for(j=0;j<n;j++)printf(" %c%i",y[j]?y[j]<0?'-':'+':' ',abs(y[j]));+ printf("\n");*/+ pulse2comb(n,m,x2,s2,y);
+ for(k=0;k<m;k++)if(x[k]!=x2[k]||s[k]!=s2[k]){+ fprintf(stderr,"Pulse-combination mismatch.\n");
+ break;
+ }
}
- }
- return 0;
+ printf("\n");+ }
+ }
+ return 0;
}
-*/
+#endif