Updated pulse coding to simpler (slightly faster) code included with

http://people.xiph.org/~tterribe/notes/cwrs.html
Removed dead code.

libcelt/cwrs.c

-/* (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 @@ static unsigned ncwrs(int _n,int _m){
   }
   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,40 +56,30 @@ unsigned ncwrs(int _n,int _m){
   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;
 }
 #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 @@ celt_uint64_t ncwrs64(int _n,int _m){
   _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 @@ void cwrsi(int _n,int _m,unsigned _i,int *_x,int *_s){
   _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 @@ void cwrsi64(int _n,int _m,celt_uint64_t _i,int *_x,int *_s){
   _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 @@ void pulse2comb(int _n,int _m,int *_x,int *_s,const int *_y){
   }
 }
 
-/*
+#if 0
+#include <stdio.h>
 #define NMAX (10)
 #define MMAX (9)
 
@@ -333,53 +283,53 @@ int main(int _argc,char **_argv){
       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