Code for converting back and forth between pulse vectors and the correspondig
[opus.git] / libcelt / cwrs.c
1 /* (C) 2007 Timothy Terriberry
2 */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7    
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10    
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14    
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18    
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32
33 /*#include <stdio.h>*/
34 #include <stdlib.h>
35
36 #include "cwrs.h"
37
38 /*Returns the numer of ways of choosing _m elements from a set of size _n with
39    replacement when a sign bit is needed for each unique element.*/
40 #if 0
41 static unsigned ncwrs(int _n,int _m){
42   static unsigned c[32][32];
43   if(_n<0||_m<0)return 0;
44   if(!c[_n][_m]){
45     if(_m<=0)c[_n][_m]=1;
46     else if(_n>0)c[_n][_m]=ncwrs(_n-1,_m)+ncwrs(_n,_m-1)+ncwrs(_n-1,_m-1);
47   }
48   return c[_n][_m];
49 }
50
51 #else
52
53 /*Returns the greatest common divisor of _a and _b.*/
54 static unsigned gcd(unsigned _a,unsigned _b){
55   unsigned r;
56   while(_b){
57     r=_a%_b;
58     _a=_b;
59     _b=r;
60   }
61   return _a;
62 }
63
64 /*Returns _a*b/_d, under the assumption that the result is an integer, avoiding
65    overflow.
66   It is assumed, but not required, that _b is smaller than _a.*/
67 static unsigned umuldiv(unsigned _a,unsigned _b,unsigned _d){
68   unsigned d;
69   d=gcd(_b,_d);
70   return (_a/(_d/d))*(_b/d);
71 }
72
73 unsigned ncwrs(int _n,int _m){
74   unsigned ret;
75   unsigned f;
76   unsigned d;
77   int      i;
78   if(_n<0||_m<0)return 0;
79   if(_m==0)return 1;
80   if(_n==0)return 0;
81   ret=0;
82   f=_n;
83   d=1;
84   for(i=1;i<=_m;i++){
85     ret+=f*d<<i;
86 #if 0
87     f=umuldiv(f,_n-i,i+1);
88     d=umuldiv(d,_m-i,i);
89 #else
90     f=(f*(_n-i))/(i+1);
91     d=(d*(_m-i))/i;
92 #endif
93   }
94   return ret;
95 }
96 #endif
97
98 /*Returns the _i'th combination of _m elements chosen from a set of size _n
99    with associated sign bits.
100   _x:      Returns the combination with elements sorted in ascending order.
101   _s:      Returns the associated sign bits.*/
102 void cwrsi(int _n,int _m,unsigned _i,int *_x,int *_s){
103   unsigned pn;
104   int      j;
105   int      k;
106   pn=ncwrs(_n-1,_m);
107   for(k=j=0;k<_m;k++){
108     unsigned pp;
109     unsigned p;
110     unsigned t;
111     pp=0;
112     p=ncwrs(_n-j,_m-k)-pn;
113     if(k>0){
114       t=p>>1;
115       if(t<=_i||_s[k-1])_i+=t;
116     }
117     pn=ncwrs(_n-j-1,_m-k-1);
118     while(p<=_i){
119       pp=p;
120       j++;
121       p+=pn;
122       pn=ncwrs(_n-j-1,_m-k-1);
123       p+=pn;
124     }
125     t=p-pp>>1;
126     _s[k]=_i-pp>=t;
127     _x[k]=j;
128     _i-=pp;
129     if(_s[k])_i-=t;
130   }
131 }
132
133 /*Returns the index of the given combination of _m elements chosen from a set
134    of size _n with associated sign bits.
135   _x:      The combination with elements sorted in ascending order.
136   _s:      The associated sign bits.*/
137 unsigned icwrs(int _n,int _m,const int *_x,const int *_s){
138   unsigned pn;
139   unsigned i;
140   int      j;
141   int      k;
142   i=0;
143   pn=ncwrs(_n-1,_m);
144   for(k=j=0;k<_m;k++){
145     unsigned pp;
146     unsigned p;
147     pp=0;
148     p=ncwrs(_n-j,_m-k)-pn;
149     if(k>0)p>>=1;
150     pn=ncwrs(_n-j-1,_m-k-1);
151     while(j<_x[k]){
152       pp=p;
153       j++;
154       p+=pn;
155       pn=ncwrs(_n-j-1,_m-k-1);
156       p+=pn;
157     }
158     i+=pp;
159     if((k==0||_x[k]!=_x[k-1])&&_s[k])i+=p-pp>>1;
160   }
161   return i;
162 }
163
164 /*Converts a combination _x of _m unit pulses with associated sign bits _s into
165    a pulse vector _y of length _n.
166   _y: Returns the vector of pulses.
167   _x: The combination with elements sorted in ascending order.
168   _s: The associated sign bits.*/
169 void comb2pulse(int _n,int _m,int *_y,const int *_x,const int *_s){
170   int j;
171   int k;
172   int n;
173   for(k=j=0;k<_m;k+=n){
174     for(n=1;k+n<_m&&_x[k+n]==_x[k];n++);
175     while(j<_x[k])_y[j++]=0;
176     _y[j++]=_s[k]?-n:n;
177   }
178   while(j<_n)_y[j++]=0;
179 }
180
181 /*Converts a pulse vector vector _y of length _n into a combination of _m unit
182    pulses with associated sign bits _s.
183   _x: Returns the combination with elements sorted in ascending order.
184   _s: Returns the associated sign bits.
185   _y: The vector of pulses, whose sum of absolute values must be _m.*/
186 void pulse2comb(int _n,int _m,int *_x,int *_s,const int *_y){
187   int j;
188   int k;
189   for(k=j=0;j<_n;j++){
190     if(_y[j]){
191       int n;
192       int s;
193       n=abs(_y[j]);
194       s=_y[j]<0;
195       for(;n-->0;k++){
196         _x[k]=j;
197         _s[k]=s;
198       }
199     }
200   }
201 }
202
203 /*
204 #define NMAX (10)
205 #define MMAX (9)
206
207 int main(int _argc,char **_argv){
208   int n;
209   for(n=0;n<=NMAX;n++){
210     int m;
211     for(m=0;m<=MMAX;m++){
212       unsigned nc;
213       unsigned i;
214       nc=ncwrs(n,m);
215       for(i=0;i<nc;i++){
216         int x[MMAX];
217         int s[MMAX];
218         int x2[MMAX];
219         int s2[MMAX];
220         int y[NMAX];
221         int j;
222         int k;
223         cwrsi(n,m,i,x,s);
224         printf("%6u of %u:",i,nc);
225         for(k=0;k<m;k++){
226           printf(" %c%i",k>0&&x[k]==x[k-1]?' ':s[k]?'-':'+',x[k]);
227         }
228         printf(" ->");
229         if(icwrs(n,m,x,s)!=i){
230           fprintf(stderr,"Combination-index mismatch.\n");
231         }
232         comb2pulse(n,m,y,x,s);
233         for(j=0;j<n;j++)printf(" %c%i",y[j]?y[j]<0?'-':'+':' ',abs(y[j]));
234         printf("\n");
235         pulse2comb(n,m,x2,s2,y);
236         for(k=0;k<m;k++)if(x[k]!=x2[k]||s[k]!=s2[k]){
237           fprintf(stderr,"Pulse-combination mismatch.\n");
238           break;
239         }
240       }
241       printf("\n");
242     }
243   }
244   return 0;
245 }
246 */