|
0
|
1 /*
|
|
|
2 * Copyright (c) 2014 Genome Research Ltd.
|
|
|
3 * Author(s): James Bonfield
|
|
|
4 *
|
|
|
5 * Redistribution and use in source and binary forms, with or without
|
|
|
6 * modification, are permitted provided that the following conditions are met:
|
|
|
7 *
|
|
|
8 * 1. Redistributions of source code must retain the above copyright notice,
|
|
|
9 * this list of conditions and the following disclaimer.
|
|
|
10 *
|
|
|
11 * 2. Redistributions in binary form must reproduce the above
|
|
|
12 * copyright notice, this list of conditions and the following
|
|
|
13 * disclaimer in the documentation and/or other materials provided
|
|
|
14 * with the distribution.
|
|
|
15 *
|
|
|
16 * 3. Neither the names Genome Research Ltd and Wellcome Trust Sanger
|
|
|
17 * Institute nor the names of its contributors may be used to endorse
|
|
|
18 * or promote products derived from this software without specific
|
|
|
19 * prior written permission.
|
|
|
20 *
|
|
|
21 * THIS SOFTWARE IS PROVIDED BY GENOME RESEARCH LTD AND CONTRIBUTORS "AS
|
|
|
22 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
|
|
|
23 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
|
|
|
24 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL GENOME RESEARCH
|
|
|
25 * LTD OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
|
|
|
26 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
|
|
|
27 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
|
|
28 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
|
|
29 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
|
|
30 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
|
|
|
31 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
|
32 */
|
|
|
33
|
|
|
34 /*
|
|
|
35 * Author: James Bonfield, Wellcome Trust Sanger Institute. 2014
|
|
|
36 */
|
|
|
37
|
|
|
38 #include <stdint.h>
|
|
|
39 #include <stdlib.h>
|
|
|
40 #include <stdio.h>
|
|
|
41 #include <unistd.h>
|
|
|
42 #include <assert.h>
|
|
|
43 #include <string.h>
|
|
|
44 #include <sys/time.h>
|
|
|
45
|
|
|
46 #include "cram/rANS_static.h"
|
|
|
47 #include "cram/rANS_byte.h"
|
|
|
48
|
|
|
49 #define TF_SHIFT 12
|
|
|
50 #define TOTFREQ (1<<TF_SHIFT)
|
|
|
51
|
|
|
52 #define ABS(a) ((a)>0?(a):-(a))
|
|
|
53 #ifndef BLK_SIZE
|
|
|
54 # define BLK_SIZE 1024*1024
|
|
|
55 #endif
|
|
|
56
|
|
|
57 // Room to allow for expanded BLK_SIZE on worst case compression.
|
|
|
58 #define BLK_SIZE2 ((int)(1.05*BLK_SIZE))
|
|
|
59
|
|
|
60 /*-----------------------------------------------------------------------------
|
|
|
61 * Memory to memory compression functions.
|
|
|
62 *
|
|
|
63 * These are original versions without any manual loop unrolling. They
|
|
|
64 * are easier to understand, but can be up to 2x slower.
|
|
|
65 */
|
|
|
66
|
|
|
67 unsigned char *rans_compress_O0(unsigned char *in, unsigned int in_size,
|
|
|
68 unsigned int *out_size) {
|
|
|
69 unsigned char *out_buf = malloc(1.05*in_size + 257*257*3 + 9);
|
|
|
70 unsigned char *cp, *out_end;
|
|
|
71 RansEncSymbol syms[256];
|
|
|
72 RansState rans0, rans1, rans2, rans3;
|
|
|
73 uint8_t* ptr;
|
|
|
74 int F[256] = {0}, i, j, tab_size, rle, x, fsum = 0;
|
|
|
75 int m = 0, M = 0;
|
|
|
76 uint64_t tr;
|
|
|
77
|
|
|
78 if (!out_buf)
|
|
|
79 return NULL;
|
|
|
80
|
|
|
81 ptr = out_end = out_buf + (int)(1.05*in_size) + 257*257*3 + 9;
|
|
|
82
|
|
|
83 // Compute statistics
|
|
|
84 for (i = 0; i < in_size; i++) {
|
|
|
85 F[in[i]]++;
|
|
|
86 }
|
|
|
87 tr = ((uint64_t)TOTFREQ<<31)/in_size + (1<<30)/in_size;
|
|
|
88
|
|
|
89 // Normalise so T[i] == TOTFREQ
|
|
|
90 for (m = M = j = 0; j < 256; j++) {
|
|
|
91 if (!F[j])
|
|
|
92 continue;
|
|
|
93
|
|
|
94 if (m < F[j])
|
|
|
95 m = F[j], M = j;
|
|
|
96
|
|
|
97 if ((F[j] = (F[j]*tr)>>31) == 0)
|
|
|
98 F[j] = 1;
|
|
|
99 fsum += F[j];
|
|
|
100 }
|
|
|
101
|
|
|
102 fsum++;
|
|
|
103 if (fsum < TOTFREQ)
|
|
|
104 F[M] += TOTFREQ-fsum;
|
|
|
105 else
|
|
|
106 F[M] -= fsum-TOTFREQ;
|
|
|
107
|
|
|
108 //printf("F[%d]=%d\n", M, F[M]);
|
|
|
109 assert(F[M]>0);
|
|
|
110
|
|
|
111 // Encode statistics.
|
|
|
112 cp = out_buf+9;
|
|
|
113
|
|
|
114 for (x = rle = j = 0; j < 256; j++) {
|
|
|
115 if (F[j]) {
|
|
|
116 // j
|
|
|
117 if (rle) {
|
|
|
118 rle--;
|
|
|
119 } else {
|
|
|
120 *cp++ = j;
|
|
|
121 if (!rle && j && F[j-1]) {
|
|
|
122 for(rle=j+1; rle<256 && F[rle]; rle++)
|
|
|
123 ;
|
|
|
124 rle -= j+1;
|
|
|
125 *cp++ = rle;
|
|
|
126 }
|
|
|
127 //fprintf(stderr, "%d: %d %d\n", j, rle, N[j]);
|
|
|
128 }
|
|
|
129
|
|
|
130 // F[j]
|
|
|
131 if (F[j]<128) {
|
|
|
132 *cp++ = F[j];
|
|
|
133 } else {
|
|
|
134 *cp++ = 128 | (F[j]>>8);
|
|
|
135 *cp++ = F[j]&0xff;
|
|
|
136 }
|
|
|
137 RansEncSymbolInit(&syms[j], x, F[j], TF_SHIFT);
|
|
|
138 x += F[j];
|
|
|
139 }
|
|
|
140 }
|
|
|
141 *cp++ = 0;
|
|
|
142
|
|
|
143 //write(1, out_buf+4, cp-(out_buf+4));
|
|
|
144 tab_size = cp-out_buf;
|
|
|
145
|
|
|
146 RansEncInit(&rans0);
|
|
|
147 RansEncInit(&rans1);
|
|
|
148 RansEncInit(&rans2);
|
|
|
149 RansEncInit(&rans3);
|
|
|
150
|
|
|
151 switch (i=(in_size&3)) {
|
|
|
152 case 3: RansEncPutSymbol(&rans2, &ptr, &syms[in[in_size-(i-2)]]);
|
|
|
153 case 2: RansEncPutSymbol(&rans1, &ptr, &syms[in[in_size-(i-1)]]);
|
|
|
154 case 1: RansEncPutSymbol(&rans0, &ptr, &syms[in[in_size-(i-0)]]);
|
|
|
155 case 0:
|
|
|
156 break;
|
|
|
157 }
|
|
|
158 for (i=(in_size &~3); i>0; i-=4) {
|
|
|
159 RansEncSymbol *s3 = &syms[in[i-1]];
|
|
|
160 RansEncSymbol *s2 = &syms[in[i-2]];
|
|
|
161 RansEncSymbol *s1 = &syms[in[i-3]];
|
|
|
162 RansEncSymbol *s0 = &syms[in[i-4]];
|
|
|
163
|
|
|
164 RansEncPutSymbol(&rans3, &ptr, s3);
|
|
|
165 RansEncPutSymbol(&rans2, &ptr, s2);
|
|
|
166 RansEncPutSymbol(&rans1, &ptr, s1);
|
|
|
167 RansEncPutSymbol(&rans0, &ptr, s0);
|
|
|
168 }
|
|
|
169
|
|
|
170 RansEncFlush(&rans3, &ptr);
|
|
|
171 RansEncFlush(&rans2, &ptr);
|
|
|
172 RansEncFlush(&rans1, &ptr);
|
|
|
173 RansEncFlush(&rans0, &ptr);
|
|
|
174
|
|
|
175 // Finalise block size and return it
|
|
|
176 *out_size = (out_end - ptr) + tab_size;
|
|
|
177
|
|
|
178 cp = out_buf;
|
|
|
179
|
|
|
180 *cp++ = 0; // order
|
|
|
181 *cp++ = ((*out_size-9)>> 0) & 0xff;
|
|
|
182 *cp++ = ((*out_size-9)>> 8) & 0xff;
|
|
|
183 *cp++ = ((*out_size-9)>>16) & 0xff;
|
|
|
184 *cp++ = ((*out_size-9)>>24) & 0xff;
|
|
|
185
|
|
|
186 *cp++ = (in_size>> 0) & 0xff;
|
|
|
187 *cp++ = (in_size>> 8) & 0xff;
|
|
|
188 *cp++ = (in_size>>16) & 0xff;
|
|
|
189 *cp++ = (in_size>>24) & 0xff;
|
|
|
190
|
|
|
191 memmove(out_buf + tab_size, ptr, out_end-ptr);
|
|
|
192
|
|
|
193 return out_buf;
|
|
|
194 }
|
|
|
195
|
|
|
196 typedef struct {
|
|
|
197 struct {
|
|
|
198 int F;
|
|
|
199 int C;
|
|
|
200 } fc[256];
|
|
|
201 unsigned char *R;
|
|
|
202 } ari_decoder;
|
|
|
203
|
|
|
204 unsigned char *rans_uncompress_O0(unsigned char *in, unsigned int in_size,
|
|
|
205 unsigned int *out_size) {
|
|
|
206 /* Load in the static tables */
|
|
|
207 unsigned char *cp = in + 9;
|
|
|
208 int i, j, x, out_sz, in_sz, rle;
|
|
|
209 char *out_buf;
|
|
|
210 ari_decoder D;
|
|
|
211 RansDecSymbol syms[256];
|
|
|
212
|
|
|
213 memset(&D, 0, sizeof(D));
|
|
|
214
|
|
|
215 if (*in++ != 0) // Order-0 check
|
|
|
216 return NULL;
|
|
|
217
|
|
|
218 in_sz = ((in[0])<<0) | ((in[1])<<8) | ((in[2])<<16) | ((in[3])<<24);
|
|
|
219 out_sz = ((in[4])<<0) | ((in[5])<<8) | ((in[6])<<16) | ((in[7])<<24);
|
|
|
220 if (in_sz != in_size-9)
|
|
|
221 return NULL;
|
|
|
222
|
|
|
223 out_buf = malloc(out_sz);
|
|
|
224 if (!out_buf)
|
|
|
225 return NULL;
|
|
|
226
|
|
|
227 //fprintf(stderr, "out_sz=%d\n", out_sz);
|
|
|
228
|
|
|
229 // Precompute reverse lookup of frequency.
|
|
|
230 rle = x = 0;
|
|
|
231 j = *cp++;
|
|
|
232 do {
|
|
|
233 if ((D.fc[j].F = *cp++) >= 128) {
|
|
|
234 D.fc[j].F &= ~128;
|
|
|
235 D.fc[j].F = ((D.fc[j].F & 127) << 8) | *cp++;
|
|
|
236 }
|
|
|
237 D.fc[j].C = x;
|
|
|
238
|
|
|
239 RansDecSymbolInit(&syms[j], D.fc[j].C, D.fc[j].F);
|
|
|
240
|
|
|
241 /* Build reverse lookup table */
|
|
|
242 if (!D.R) D.R = (unsigned char *)malloc(TOTFREQ);
|
|
|
243 memset(&D.R[x], j, D.fc[j].F);
|
|
|
244
|
|
|
245 x += D.fc[j].F;
|
|
|
246
|
|
|
247 if (!rle && j+1 == *cp) {
|
|
|
248 j = *cp++;
|
|
|
249 rle = *cp++;
|
|
|
250 } else if (rle) {
|
|
|
251 rle--;
|
|
|
252 j++;
|
|
|
253 } else {
|
|
|
254 j = *cp++;
|
|
|
255 }
|
|
|
256 } while(j);
|
|
|
257
|
|
|
258 assert(x < TOTFREQ);
|
|
|
259
|
|
|
260 RansState rans0, rans1, rans2, rans3;
|
|
|
261 uint8_t *ptr = cp;
|
|
|
262 RansDecInit(&rans0, &ptr);
|
|
|
263 RansDecInit(&rans1, &ptr);
|
|
|
264 RansDecInit(&rans2, &ptr);
|
|
|
265 RansDecInit(&rans3, &ptr);
|
|
|
266
|
|
|
267 int out_end = (out_sz&~3);
|
|
|
268
|
|
|
269 RansState R[4];
|
|
|
270 R[0] = rans0;
|
|
|
271 R[1] = rans1;
|
|
|
272 R[2] = rans2;
|
|
|
273 R[3] = rans3;
|
|
|
274 uint32_t mask = (1u << TF_SHIFT)-1;
|
|
|
275
|
|
|
276 for (i=0; i < out_end; i+=4) {
|
|
|
277 uint32_t m[4] = {R[0] & mask,
|
|
|
278 R[1] & mask,
|
|
|
279 R[2] & mask,
|
|
|
280 R[3] & mask};
|
|
|
281 uint8_t c[4] = {D.R[m[0]],
|
|
|
282 D.R[m[1]],
|
|
|
283 D.R[m[2]],
|
|
|
284 D.R[m[3]]};
|
|
|
285 out_buf[i+0] = c[0];
|
|
|
286 out_buf[i+1] = c[1];
|
|
|
287 out_buf[i+2] = c[2];
|
|
|
288 out_buf[i+3] = c[3];
|
|
|
289
|
|
|
290 // RansDecAdvanceSymbolStep(&R[0], &syms[c[0]], TF_SHIFT);
|
|
|
291 // RansDecAdvanceSymbolStep(&R[1], &syms[c[1]], TF_SHIFT);
|
|
|
292 // RansDecAdvanceSymbolStep(&R[2], &syms[c[2]], TF_SHIFT);
|
|
|
293 // RansDecAdvanceSymbolStep(&R[3], &syms[c[3]], TF_SHIFT);
|
|
|
294 R[0] = syms[c[0]].freq * (R[0]>>TF_SHIFT);
|
|
|
295 R[1] = syms[c[1]].freq * (R[1]>>TF_SHIFT);
|
|
|
296 R[2] = syms[c[2]].freq * (R[2]>>TF_SHIFT);
|
|
|
297 R[3] = syms[c[3]].freq * (R[3]>>TF_SHIFT);
|
|
|
298
|
|
|
299 R[0] += m[0] - syms[c[0]].start;
|
|
|
300 R[1] += m[1] - syms[c[1]].start;
|
|
|
301 R[2] += m[2] - syms[c[2]].start;
|
|
|
302 R[3] += m[3] - syms[c[3]].start;
|
|
|
303
|
|
|
304 RansDecRenorm(&R[0], &ptr);
|
|
|
305 RansDecRenorm(&R[1], &ptr);
|
|
|
306 RansDecRenorm(&R[2], &ptr);
|
|
|
307 RansDecRenorm(&R[3], &ptr);
|
|
|
308 }
|
|
|
309
|
|
|
310 rans0 = R[0];
|
|
|
311 rans1 = R[1];
|
|
|
312 rans2 = R[2];
|
|
|
313 rans3 = R[3];
|
|
|
314
|
|
|
315 switch(out_sz&3) {
|
|
|
316 unsigned char c;
|
|
|
317 case 0:
|
|
|
318 break;
|
|
|
319 case 1:
|
|
|
320 c = D.R[RansDecGet(&rans0, TF_SHIFT)];
|
|
|
321 RansDecAdvanceSymbol(&rans0, &ptr, &syms[c], TF_SHIFT);
|
|
|
322 out_buf[out_end] = c;
|
|
|
323 break;
|
|
|
324
|
|
|
325 case 2:
|
|
|
326 c = D.R[RansDecGet(&rans0, TF_SHIFT)];
|
|
|
327 RansDecAdvanceSymbol(&rans0, &ptr, &syms[c], TF_SHIFT);
|
|
|
328 out_buf[out_end] = c;
|
|
|
329
|
|
|
330 c = D.R[RansDecGet(&rans1, TF_SHIFT)];
|
|
|
331 RansDecAdvanceSymbol(&rans1, &ptr, &syms[c], TF_SHIFT);
|
|
|
332 out_buf[out_end+1] = c;
|
|
|
333 break;
|
|
|
334
|
|
|
335 case 3:
|
|
|
336 c = D.R[RansDecGet(&rans0, TF_SHIFT)];
|
|
|
337 RansDecAdvanceSymbol(&rans0, &ptr, &syms[c], TF_SHIFT);
|
|
|
338 out_buf[out_end] = c;
|
|
|
339
|
|
|
340 c = D.R[RansDecGet(&rans1, TF_SHIFT)];
|
|
|
341 RansDecAdvanceSymbol(&rans1, &ptr, &syms[c], TF_SHIFT);
|
|
|
342 out_buf[out_end+1] = c;
|
|
|
343
|
|
|
344 c = D.R[RansDecGet(&rans2, TF_SHIFT)];
|
|
|
345 RansDecAdvanceSymbol(&rans2, &ptr, &syms[c], TF_SHIFT);
|
|
|
346 out_buf[out_end+2] = c;
|
|
|
347 break;
|
|
|
348 }
|
|
|
349
|
|
|
350 *out_size = out_sz;
|
|
|
351
|
|
|
352 if (D.R) free(D.R);
|
|
|
353
|
|
|
354 return (unsigned char *)out_buf;
|
|
|
355 }
|
|
|
356
|
|
|
357 unsigned char *rans_compress_O1(unsigned char *in, unsigned int in_size,
|
|
|
358 unsigned int *out_size) {
|
|
|
359 unsigned char *out_buf, *out_end, *cp;
|
|
|
360 unsigned int last_i, tab_size, rle_i, rle_j;
|
|
|
361 RansEncSymbol syms[256][256];
|
|
|
362
|
|
|
363 if (in_size < 4)
|
|
|
364 return rans_compress_O0(in, in_size, out_size);
|
|
|
365
|
|
|
366 out_buf = malloc(1.05*in_size + 257*257*3 + 9);
|
|
|
367 if (!out_buf)
|
|
|
368 return NULL;
|
|
|
369
|
|
|
370 out_end = out_buf + (int)(1.05*in_size) + 257*257*3 + 9;
|
|
|
371 cp = out_buf+9;
|
|
|
372
|
|
|
373 int F[256][256], T[256], i, j;
|
|
|
374 unsigned char c;
|
|
|
375
|
|
|
376 memset(F, 0, 256*256*sizeof(int));
|
|
|
377 memset(T, 0, 256*sizeof(int));
|
|
|
378 //for (last = 0, i=in_size-1; i>=0; i--) {
|
|
|
379 // F[last][c = in[i]]++;
|
|
|
380 // T[last]++;
|
|
|
381 // last = c;
|
|
|
382 //}
|
|
|
383
|
|
|
384 for (last_i=i=0; i<in_size; i++) {
|
|
|
385 F[last_i][c = in[i]]++;
|
|
|
386 T[last_i]++;
|
|
|
387 last_i = c;
|
|
|
388 }
|
|
|
389 F[0][in[1*(in_size>>2)]]++;
|
|
|
390 F[0][in[2*(in_size>>2)]]++;
|
|
|
391 F[0][in[3*(in_size>>2)]]++;
|
|
|
392 T[0]+=3;
|
|
|
393
|
|
|
394 // Normalise so T[i] == TOTFREQ
|
|
|
395 for (rle_i = i = 0; i < 256; i++) {
|
|
|
396 int t2, m, M;
|
|
|
397 unsigned int x;
|
|
|
398
|
|
|
399 if (T[i] == 0)
|
|
|
400 continue;
|
|
|
401
|
|
|
402 //uint64_t p = (TOTFREQ * TOTFREQ) / t;
|
|
|
403 double p = ((double)TOTFREQ)/T[i];
|
|
|
404 for (t2 = m = M = j = 0; j < 256; j++) {
|
|
|
405 if (!F[i][j])
|
|
|
406 continue;
|
|
|
407
|
|
|
408 if (m < F[i][j])
|
|
|
409 m = F[i][j], M = j;
|
|
|
410
|
|
|
411 //if ((F[i][j] = (F[i][j] * p) / TOTFREQ) == 0)
|
|
|
412 if ((F[i][j] *= p) == 0)
|
|
|
413 F[i][j] = 1;
|
|
|
414 t2 += F[i][j];
|
|
|
415 }
|
|
|
416
|
|
|
417 t2++;
|
|
|
418 if (t2 < TOTFREQ)
|
|
|
419 F[i][M] += TOTFREQ-t2;
|
|
|
420 else
|
|
|
421 F[i][M] -= t2-TOTFREQ;
|
|
|
422
|
|
|
423 // Store frequency table
|
|
|
424 // i
|
|
|
425 if (rle_i) {
|
|
|
426 rle_i--;
|
|
|
427 } else {
|
|
|
428 *cp++ = i;
|
|
|
429 // FIXME: could use order-0 statistics to observe which alphabet
|
|
|
430 // symbols are present and base RLE on that ordering instead.
|
|
|
431 if (i && T[i-1]) {
|
|
|
432 for(rle_i=i+1; rle_i<256 && T[rle_i]; rle_i++)
|
|
|
433 ;
|
|
|
434 rle_i -= i+1;
|
|
|
435 *cp++ = rle_i;
|
|
|
436 }
|
|
|
437 }
|
|
|
438
|
|
|
439 int *F_i_ = F[i];
|
|
|
440 x = 0;
|
|
|
441 rle_j = 0;
|
|
|
442 for (j = 0; j < 256; j++) {
|
|
|
443 if (F_i_[j]) {
|
|
|
444 //fprintf(stderr, "F[%d][%d]=%d, x=%d\n", i, j, F_i_[j], x);
|
|
|
445
|
|
|
446 // j
|
|
|
447 if (rle_j) {
|
|
|
448 rle_j--;
|
|
|
449 } else {
|
|
|
450 *cp++ = j;
|
|
|
451 if (!rle_j && j && F_i_[j-1]) {
|
|
|
452 for(rle_j=j+1; rle_j<256 && F_i_[rle_j]; rle_j++)
|
|
|
453 ;
|
|
|
454 rle_j -= j+1;
|
|
|
455 *cp++ = rle_j;
|
|
|
456 }
|
|
|
457 }
|
|
|
458
|
|
|
459 // F_i_[j]
|
|
|
460 if (F_i_[j]<128) {
|
|
|
461 *cp++ = F_i_[j];
|
|
|
462 } else {
|
|
|
463 *cp++ = 128 | (F_i_[j]>>8);
|
|
|
464 *cp++ = F_i_[j]&0xff;
|
|
|
465 }
|
|
|
466
|
|
|
467 RansEncSymbolInit(&syms[i][j], x, F_i_[j], TF_SHIFT);
|
|
|
468 x += F_i_[j];
|
|
|
469 }
|
|
|
470 }
|
|
|
471 *cp++ = 0;
|
|
|
472 }
|
|
|
473 *cp++ = 0;
|
|
|
474
|
|
|
475 //write(1, out_buf+4, cp-(out_buf+4));
|
|
|
476 tab_size = cp - out_buf;
|
|
|
477 assert(tab_size < 257*257*3);
|
|
|
478
|
|
|
479 RansState rans0, rans1, rans2, rans3;
|
|
|
480 RansEncInit(&rans0);
|
|
|
481 RansEncInit(&rans1);
|
|
|
482 RansEncInit(&rans2);
|
|
|
483 RansEncInit(&rans3);
|
|
|
484
|
|
|
485 uint8_t* ptr = out_end;
|
|
|
486
|
|
|
487 int isz4 = in_size>>2;
|
|
|
488 int i0 = 1*isz4-2;
|
|
|
489 int i1 = 2*isz4-2;
|
|
|
490 int i2 = 3*isz4-2;
|
|
|
491 int i3 = 4*isz4-2;
|
|
|
492
|
|
|
493 unsigned char l0 = in[i0+1];
|
|
|
494 unsigned char l1 = in[i1+1];
|
|
|
495 unsigned char l2 = in[i2+1];
|
|
|
496 unsigned char l3 = in[i3+1];
|
|
|
497
|
|
|
498 // Deal with the remainder
|
|
|
499 l3 = in[in_size-1];
|
|
|
500 for (i3 = in_size-2; i3 > 4*isz4-2; i3--) {
|
|
|
501 unsigned char c3 = in[i3];
|
|
|
502 RansEncPutSymbol(&rans3, &ptr, &syms[c3][l3]);
|
|
|
503 l3 = c3;
|
|
|
504 }
|
|
|
505
|
|
|
506 for (; i0 >= 0; i0--, i1--, i2--, i3--) {
|
|
|
507 unsigned char c0, c1, c2, c3;
|
|
|
508 RansEncSymbol *s3 = &syms[c3 = in[i3]][l3];
|
|
|
509 RansEncSymbol *s2 = &syms[c2 = in[i2]][l2];
|
|
|
510 RansEncSymbol *s1 = &syms[c1 = in[i1]][l1];
|
|
|
511 RansEncSymbol *s0 = &syms[c0 = in[i0]][l0];
|
|
|
512
|
|
|
513 RansEncPutSymbol(&rans3, &ptr, s3);
|
|
|
514 RansEncPutSymbol(&rans2, &ptr, s2);
|
|
|
515 RansEncPutSymbol(&rans1, &ptr, s1);
|
|
|
516 RansEncPutSymbol(&rans0, &ptr, s0);
|
|
|
517
|
|
|
518 l0 = c0;
|
|
|
519 l1 = c1;
|
|
|
520 l2 = c2;
|
|
|
521 l3 = c3;
|
|
|
522 }
|
|
|
523
|
|
|
524 RansEncPutSymbol(&rans3, &ptr, &syms[0][l3]);
|
|
|
525 RansEncPutSymbol(&rans2, &ptr, &syms[0][l2]);
|
|
|
526 RansEncPutSymbol(&rans1, &ptr, &syms[0][l1]);
|
|
|
527 RansEncPutSymbol(&rans0, &ptr, &syms[0][l0]);
|
|
|
528
|
|
|
529 RansEncFlush(&rans3, &ptr);
|
|
|
530 RansEncFlush(&rans2, &ptr);
|
|
|
531 RansEncFlush(&rans1, &ptr);
|
|
|
532 RansEncFlush(&rans0, &ptr);
|
|
|
533
|
|
|
534 *out_size = (out_end - ptr) + tab_size;
|
|
|
535
|
|
|
536 cp = out_buf;
|
|
|
537 *cp++ = 1; // order
|
|
|
538
|
|
|
539 *cp++ = ((*out_size-9)>> 0) & 0xff;
|
|
|
540 *cp++ = ((*out_size-9)>> 8) & 0xff;
|
|
|
541 *cp++ = ((*out_size-9)>>16) & 0xff;
|
|
|
542 *cp++ = ((*out_size-9)>>24) & 0xff;
|
|
|
543
|
|
|
544 *cp++ = (in_size>> 0) & 0xff;
|
|
|
545 *cp++ = (in_size>> 8) & 0xff;
|
|
|
546 *cp++ = (in_size>>16) & 0xff;
|
|
|
547 *cp++ = (in_size>>24) & 0xff;
|
|
|
548
|
|
|
549 memmove(out_buf + tab_size, ptr, out_end-ptr);
|
|
|
550
|
|
|
551 return out_buf;
|
|
|
552 }
|
|
|
553
|
|
|
554 unsigned char *rans_uncompress_O1(unsigned char *in, unsigned int in_size,
|
|
|
555 unsigned int *out_size) {
|
|
|
556 /* Load in the static tables */
|
|
|
557 unsigned char *cp = in + 9;
|
|
|
558 int i, j = -999, x, out_sz, in_sz, rle_i, rle_j;
|
|
|
559 char *out_buf;
|
|
|
560 ari_decoder D[256];
|
|
|
561 RansDecSymbol syms[256][256];
|
|
|
562
|
|
|
563 memset(D, 0, 256*sizeof(*D));
|
|
|
564
|
|
|
565 if (*in++ != 1) // Order-1 check
|
|
|
566 return NULL;
|
|
|
567
|
|
|
568 in_sz = ((in[0])<<0) | ((in[1])<<8) | ((in[2])<<16) | ((in[3])<<24);
|
|
|
569 out_sz = ((in[4])<<0) | ((in[5])<<8) | ((in[6])<<16) | ((in[7])<<24);
|
|
|
570 if (in_sz != in_size-9)
|
|
|
571 return NULL;
|
|
|
572
|
|
|
573 out_buf = malloc(out_sz);
|
|
|
574 if (!out_buf)
|
|
|
575 return NULL;
|
|
|
576
|
|
|
577 //fprintf(stderr, "out_sz=%d\n", out_sz);
|
|
|
578
|
|
|
579 //i = *cp++;
|
|
|
580 rle_i = 0;
|
|
|
581 i = *cp++;
|
|
|
582 do {
|
|
|
583 rle_j = x = 0;
|
|
|
584 j = *cp++;
|
|
|
585 do {
|
|
|
586 if ((D[i].fc[j].F = *cp++) >= 128) {
|
|
|
587 D[i].fc[j].F &= ~128;
|
|
|
588 D[i].fc[j].F = ((D[i].fc[j].F & 127) << 8) | *cp++;
|
|
|
589 }
|
|
|
590 D[i].fc[j].C = x;
|
|
|
591
|
|
|
592 //fprintf(stderr, "i=%d j=%d F=%d C=%d\n", i, j, D[i].fc[j].F, D[i].fc[j].C);
|
|
|
593
|
|
|
594 if (!D[i].fc[j].F)
|
|
|
595 D[i].fc[j].F = TOTFREQ;
|
|
|
596
|
|
|
597 RansDecSymbolInit(&syms[i][j], D[i].fc[j].C, D[i].fc[j].F);
|
|
|
598
|
|
|
599 /* Build reverse lookup table */
|
|
|
600 if (!D[i].R) D[i].R = (unsigned char *)malloc(TOTFREQ);
|
|
|
601 memset(&D[i].R[x], j, D[i].fc[j].F);
|
|
|
602
|
|
|
603 x += D[i].fc[j].F;
|
|
|
604 assert(x <= TOTFREQ);
|
|
|
605
|
|
|
606 if (!rle_j && j+1 == *cp) {
|
|
|
607 j = *cp++;
|
|
|
608 rle_j = *cp++;
|
|
|
609 } else if (rle_j) {
|
|
|
610 rle_j--;
|
|
|
611 j++;
|
|
|
612 } else {
|
|
|
613 j = *cp++;
|
|
|
614 }
|
|
|
615 } while(j);
|
|
|
616
|
|
|
617 if (!rle_i && i+1 == *cp) {
|
|
|
618 i = *cp++;
|
|
|
619 rle_i = *cp++;
|
|
|
620 } else if (rle_i) {
|
|
|
621 rle_i--;
|
|
|
622 i++;
|
|
|
623 } else {
|
|
|
624 i = *cp++;
|
|
|
625 }
|
|
|
626 } while (i);
|
|
|
627
|
|
|
628 // Precompute reverse lookup of frequency.
|
|
|
629
|
|
|
630 RansState rans0, rans1, rans2, rans3;
|
|
|
631 uint8_t *ptr = cp;
|
|
|
632 RansDecInit(&rans0, &ptr);
|
|
|
633 RansDecInit(&rans1, &ptr);
|
|
|
634 RansDecInit(&rans2, &ptr);
|
|
|
635 RansDecInit(&rans3, &ptr);
|
|
|
636
|
|
|
637 int isz4 = out_sz>>2;
|
|
|
638 int l0 = 0;
|
|
|
639 int l1 = 0;
|
|
|
640 int l2 = 0;
|
|
|
641 int l3 = 0;
|
|
|
642 int i4[] = {0*isz4, 1*isz4, 2*isz4, 3*isz4};
|
|
|
643
|
|
|
644 RansState R[4];
|
|
|
645 R[0] = rans0;
|
|
|
646 R[1] = rans1;
|
|
|
647 R[2] = rans2;
|
|
|
648 R[3] = rans3;
|
|
|
649
|
|
|
650 for (; i4[0] < isz4; i4[0]++, i4[1]++, i4[2]++, i4[3]++) {
|
|
|
651 uint32_t m[4] = {R[0] & ((1u << TF_SHIFT)-1),
|
|
|
652 R[1] & ((1u << TF_SHIFT)-1),
|
|
|
653 R[2] & ((1u << TF_SHIFT)-1),
|
|
|
654 R[3] & ((1u << TF_SHIFT)-1)};
|
|
|
655
|
|
|
656 uint8_t c[4] = {D[l0].R[m[0]],
|
|
|
657 D[l1].R[m[1]],
|
|
|
658 D[l2].R[m[2]],
|
|
|
659 D[l3].R[m[3]]};
|
|
|
660
|
|
|
661 out_buf[i4[0]] = c[0];
|
|
|
662 out_buf[i4[1]] = c[1];
|
|
|
663 out_buf[i4[2]] = c[2];
|
|
|
664 out_buf[i4[3]] = c[3];
|
|
|
665
|
|
|
666 //RansDecAdvanceSymbolStep(&R[0], &syms[l0][c[0]], TF_SHIFT);
|
|
|
667 //RansDecAdvanceSymbolStep(&R[1], &syms[l1][c[1]], TF_SHIFT);
|
|
|
668 //RansDecAdvanceSymbolStep(&R[2], &syms[l2][c[2]], TF_SHIFT);
|
|
|
669 //RansDecAdvanceSymbolStep(&R[3], &syms[l3][c[3]], TF_SHIFT);
|
|
|
670
|
|
|
671 R[0] = syms[l0][c[0]].freq * (R[0]>>TF_SHIFT);
|
|
|
672 R[1] = syms[l1][c[1]].freq * (R[1]>>TF_SHIFT);
|
|
|
673 R[2] = syms[l2][c[2]].freq * (R[2]>>TF_SHIFT);
|
|
|
674 R[3] = syms[l3][c[3]].freq * (R[3]>>TF_SHIFT);
|
|
|
675
|
|
|
676 R[0] += m[0] - syms[l0][c[0]].start;
|
|
|
677 R[1] += m[1] - syms[l1][c[1]].start;
|
|
|
678 R[2] += m[2] - syms[l2][c[2]].start;
|
|
|
679 R[3] += m[3] - syms[l3][c[3]].start;
|
|
|
680
|
|
|
681 RansDecRenorm(&R[0], &ptr);
|
|
|
682 RansDecRenorm(&R[1], &ptr);
|
|
|
683 RansDecRenorm(&R[2], &ptr);
|
|
|
684 RansDecRenorm(&R[3], &ptr);
|
|
|
685
|
|
|
686 l0 = c[0];
|
|
|
687 l1 = c[1];
|
|
|
688 l2 = c[2];
|
|
|
689 l3 = c[3];
|
|
|
690 }
|
|
|
691
|
|
|
692 rans0 = R[0];
|
|
|
693 rans1 = R[1];
|
|
|
694 rans2 = R[2];
|
|
|
695 rans3 = R[3];
|
|
|
696
|
|
|
697 // Remainder
|
|
|
698 for (; i4[3] < out_sz; i4[3]++) {
|
|
|
699 unsigned char c3 = D[l3].R[RansDecGet(&rans3, TF_SHIFT)];
|
|
|
700 out_buf[i4[3]] = c3;
|
|
|
701 RansDecAdvanceSymbol(&rans3, &ptr, &syms[l3][c3], TF_SHIFT);
|
|
|
702 l3 = c3;
|
|
|
703 }
|
|
|
704
|
|
|
705 *out_size = out_sz;
|
|
|
706
|
|
|
707 for (i = 0; i < 256; i++)
|
|
|
708 if (D[i].R) free(D[i].R);
|
|
|
709
|
|
|
710 return (unsigned char *)out_buf;
|
|
|
711 }
|
|
|
712
|
|
|
713 /*-----------------------------------------------------------------------------
|
|
|
714 * Simple interface to the order-0 vs order-1 encoders and decoders.
|
|
|
715 */
|
|
|
716 unsigned char *rans_compress(unsigned char *in, unsigned int in_size,
|
|
|
717 unsigned int *out_size, int order) {
|
|
|
718 return order
|
|
|
719 ? rans_compress_O1(in, in_size, out_size)
|
|
|
720 : rans_compress_O0(in, in_size, out_size);
|
|
|
721 }
|
|
|
722
|
|
|
723 unsigned char *rans_uncompress(unsigned char *in, unsigned int in_size,
|
|
|
724 unsigned int *out_size) {
|
|
|
725 return in[0]
|
|
|
726 ? rans_uncompress_O1(in, in_size, out_size)
|
|
|
727 : rans_uncompress_O0(in, in_size, out_size);
|
|
|
728 }
|
|
|
729
|
|
|
730
|
|
|
731 #ifdef TEST_MAIN
|
|
|
732 /*-----------------------------------------------------------------------------
|
|
|
733 * Main.
|
|
|
734 *
|
|
|
735 * This is a simple command line tool for testing order-0 and order-1
|
|
|
736 * compression using the rANS codec. Simply compile with
|
|
|
737 *
|
|
|
738 * gcc -DTEST_MAIN -O3 -I. cram/rANS_static.c -o cram/rANS_static
|
|
|
739 *
|
|
|
740 * Usage: cram/rANS_static -o0 < file > file.o0
|
|
|
741 * cram/rANS_static -d < file.o0 > file2
|
|
|
742 *
|
|
|
743 * cram/rANS_static -o1 < file > file.o1
|
|
|
744 * cram/rANS_static -d < file.o1 > file2
|
|
|
745 */
|
|
|
746 int main(int argc, char **argv) {
|
|
|
747 int opt, order = 0;
|
|
|
748 unsigned char in_buf[BLK_SIZE2+257*257*3];
|
|
|
749 int decode = 0;
|
|
|
750 FILE *infp = stdin, *outfp = stdout;
|
|
|
751 struct timeval tv1, tv2;
|
|
|
752 size_t bytes = 0;
|
|
|
753
|
|
|
754 extern char *optarg;
|
|
|
755 extern int optind;
|
|
|
756
|
|
|
757 while ((opt = getopt(argc, argv, "o:d")) != -1) {
|
|
|
758 switch (opt) {
|
|
|
759 case 'o':
|
|
|
760 order = atoi(optarg);
|
|
|
761 break;
|
|
|
762
|
|
|
763 case 'd':
|
|
|
764 decode = 1;
|
|
|
765 break;
|
|
|
766 }
|
|
|
767 }
|
|
|
768
|
|
|
769 order = order ? 1 : 0; // Only support O(0) and O(1)
|
|
|
770
|
|
|
771 if (optind < argc) {
|
|
|
772 if (!(infp = fopen(argv[optind], "rb"))) {
|
|
|
773 perror(argv[optind]);
|
|
|
774 return 1;
|
|
|
775 }
|
|
|
776 optind++;
|
|
|
777 }
|
|
|
778
|
|
|
779 if (optind < argc) {
|
|
|
780 if (!(outfp = fopen(argv[optind], "wb"))) {
|
|
|
781 perror(argv[optind]);
|
|
|
782 return 1;
|
|
|
783 }
|
|
|
784 optind++;
|
|
|
785 }
|
|
|
786
|
|
|
787 gettimeofday(&tv1, NULL);
|
|
|
788
|
|
|
789 if (decode) {
|
|
|
790 // Only used in some test implementations of RC_GetFreq()
|
|
|
791 //RC_init();
|
|
|
792 //RC_init2();
|
|
|
793
|
|
|
794 for (;;) {
|
|
|
795 uint32_t in_size, out_size;
|
|
|
796 unsigned char *out;
|
|
|
797
|
|
|
798 if (4 != fread(&in_size, 1, 4, infp))
|
|
|
799 break;
|
|
|
800 if (in_size != fread(in_buf, 1, in_size, infp)) {
|
|
|
801 fprintf(stderr, "Truncated input\n");
|
|
|
802 exit(1);
|
|
|
803 }
|
|
|
804 out = rans_uncompress(in_buf, in_size, &out_size);
|
|
|
805 if (!out)
|
|
|
806 abort();
|
|
|
807
|
|
|
808 fwrite(out, 1, out_size, outfp);
|
|
|
809 free(out);
|
|
|
810
|
|
|
811 bytes += out_size;
|
|
|
812 }
|
|
|
813 } else {
|
|
|
814 for (;;) {
|
|
|
815 uint32_t in_size, out_size;
|
|
|
816 unsigned char *out;
|
|
|
817
|
|
|
818 in_size = fread(in_buf, 1, BLK_SIZE, infp);
|
|
|
819 if (in_size <= 0)
|
|
|
820 break;
|
|
|
821
|
|
|
822 out = rans_compress(in_buf, in_size, &out_size, order);
|
|
|
823
|
|
|
824 fwrite(&out_size, 1, 4, outfp);
|
|
|
825 fwrite(out, 1, out_size, outfp);
|
|
|
826 free(out);
|
|
|
827
|
|
|
828 bytes += in_size;
|
|
|
829 }
|
|
|
830 }
|
|
|
831
|
|
|
832 gettimeofday(&tv2, NULL);
|
|
|
833
|
|
|
834 fprintf(stderr, "Took %ld microseconds, %5.1f MB/s\n",
|
|
|
835 (long)(tv2.tv_sec - tv1.tv_sec)*1000000 +
|
|
|
836 tv2.tv_usec - tv1.tv_usec,
|
|
|
837 (double)bytes / ((long)(tv2.tv_sec - tv1.tv_sec)*1000000 +
|
|
|
838 tv2.tv_usec - tv1.tv_usec));
|
|
|
839 return 0;
|
|
|
840 }
|
|
|
841 #endif
|
