Mercurial > repos > youngkim > ezbamqc

comparison ezBAMQC/src/htslib/cram/rANS_static.c @ 0:dfa3745e5fd8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author youngkim
date Thu, 24 Mar 2016 17:12:52 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:dfa3745e5fd8
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