--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <assert.h>
+#include <memory>
+
+#define WIDTH 1280
+#define HEIGHT 720
+#define NUM_SYMS 256
+#define ESCAPE_LIMIT (NUM_SYMS - 1)
+
+#include "ryg_rans/rans_byte.h"
+
+using namespace std;
+
+void fdct_int32(short *const In);
+void idct_int32(short *const In);
+
+unsigned char pix[WIDTH * HEIGHT];
+short coeff[WIDTH * HEIGHT];
+
+struct RansDecoder {
+};
+
+static constexpr uint32_t prob_bits = 12;
+static constexpr uint32_t prob_scale = 1 << prob_bits;
+
+struct RansDecodeTable {
+ int cum2sym[prob_scale];
+ RansDecSymbol dsyms[NUM_SYMS];
+};
+RansDecodeTable decode_tables[16];
+
+static const unsigned char std_luminance_quant_tbl[64] = {
+#if 0
+ 16, 11, 10, 16, 24, 40, 51, 61,
+ 12, 12, 14, 19, 26, 58, 60, 55,
+ 14, 13, 16, 24, 40, 57, 69, 56,
+ 14, 17, 22, 29, 51, 87, 80, 62,
+ 18, 22, 37, 56, 68, 109, 103, 77,
+ 24, 35, 55, 64, 81, 104, 113, 92,
+ 49, 64, 78, 87, 103, 121, 120, 101,
+ 72, 92, 95, 98, 112, 100, 103, 99
+#else
+ // ff_mpeg1_default_intra_matrix
+ 8, 16, 19, 22, 26, 27, 29, 34,
+ 16, 16, 22, 24, 27, 29, 34, 37,
+ 19, 22, 26, 27, 29, 34, 34, 38,
+ 22, 22, 26, 27, 29, 34, 37, 40,
+ 22, 26, 27, 29, 32, 35, 40, 48,
+ 26, 27, 29, 32, 35, 40, 48, 58,
+ 26, 27, 29, 34, 38, 46, 56, 69,
+ 27, 29, 35, 38, 46, 56, 69, 83
+#endif
+};
+
+
+int pick_stats_for(int y, int x)
+{
+ if (x + y >= 7) return 7;
+ return x + y;
+}
+
+uint32_t read_varint(FILE *fp)
+{
+ uint32_t x = 0;
+ int shift = 0;
+ for ( ;; ) {
+ int ch = getc(fp);
+ if (ch == -1) {
+ fprintf(stderr, "Premature EOF\n");
+ exit(1);
+ }
+
+ x |= (ch & 0x7f) << shift;
+ if ((ch & 0x80) == 0) return x;
+ shift += 7;
+ }
+}
+
+static constexpr int dc_scalefac = 8; // Matches the FDCT's gain.
+static constexpr double quant_scalefac = 4.0; // whatever?
+
+static inline int unquantize(int qf, int coeff_idx)
+{
+ if (coeff_idx == 0) {
+ return qf * dc_scalefac;
+ }
+ if (qf == 0) {
+ return 0;
+ }
+
+ const int w = std_luminance_quant_tbl[coeff_idx];
+ const int s = quant_scalefac;
+ return (2 * qf * w * s) / 32;
+}
+
+int main(void)
+{
+ FILE *fp = fopen("coded.dat", "rb");
+ if (fp == nullptr) {
+ perror("coded.dat");
+ exit(1);
+ }
+
+ for (unsigned table = 0; table < 16; ++table) {
+ uint32_t cum_freq = 0;
+ for (unsigned sym = 0; sym < NUM_SYMS; ++sym) {
+ uint32_t freq = read_varint(fp);
+ fprintf(stderr, "sym=%u/%u: freq=%u\n", sym, NUM_SYMS, freq);
+ RansDecSymbolInit(&decode_tables[table].dsyms[sym], cum_freq, freq);
+ for (uint32_t i = 0; i < freq; ++i) {
+ decode_tables[table].cum2sym[cum_freq++] = sym;
+ }
+ }
+ }
+
+ // loop over all coefficients
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 8; ++x) {
+ unsigned tbl = pick_stats_for(x, y);
+
+ RansState rans = 0;
+
+ //unique_ptr<uint8_t[]> rans_bytes(new uint8_t[num_rans_bytes]);
+ //unique_ptr<uint8_t[]> sign_bytes(new uint8_t[num_sign_bytes]);
+ unique_ptr<uint8_t[]> rans_bytes;
+ unique_ptr<uint8_t[]> sign_bytes;
+ uint8_t *rans_ptr = nullptr;
+ uint8_t *sign_ptr = nullptr; // optimize later
+ uint32_t sign_buf = 0, sign_bits_left = 0;
+
+ // loop over all DCT blocks
+ for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
+ if (yb % 16 == 0) {
+ // read a block
+ uint32_t num_rans_bytes = read_varint(fp);
+ rans_bytes.reset(new uint8_t[num_rans_bytes]);
+ fread(rans_bytes.get(), 1, num_rans_bytes, fp);
+
+ uint32_t val = read_varint(fp);
+ uint8_t free_sign_bits = val & 0x7;
+ uint32_t num_sign_bytes = val >> 3;
+ sign_bytes.reset(new uint8_t[num_sign_bytes]);
+ fread(sign_bytes.get(), 1, num_sign_bytes, fp);
+
+ sign_ptr = sign_bytes.get();
+ if (free_sign_bits == 0) {
+ sign_buf = *sign_ptr++;
+ sign_bits_left = 8;
+ } else {
+ sign_buf = *sign_ptr++ >> free_sign_bits;
+ sign_bits_left = 8 - free_sign_bits;
+ }
+
+ printf("%d,%d: read %d rANS bytes, %d sign bytes\n", x, y, num_rans_bytes, num_sign_bytes);
+ //printf("first bytes: %02x %02x %02x %02x %02x %02x %02x %02x\n", rans_bytes[0], rans_bytes[1], rans_bytes[2], rans_bytes[3], rans_bytes[4], rans_bytes[5], rans_bytes[6], rans_bytes[7]);
+
+
+ // init rANS decoder
+ rans_ptr = rans_bytes.get();
+ RansDecInit(&rans, &rans_ptr);
+ }
+ for (unsigned xb = 0; xb < WIDTH; xb += 8) {
+ uint32_t k = decode_tables[tbl].cum2sym[RansDecGet(&rans, prob_bits)];
+ //printf("reading symbol, rans state = %08x\n", rans);
+ RansDecAdvanceSymbol(&rans, &rans_ptr, &decode_tables[tbl].dsyms[k], prob_bits);
+ //printf("done reading symbol, rans state = %08x\n", rans);
+ assert(k <= ESCAPE_LIMIT);
+ if (k == ESCAPE_LIMIT) {
+ k = RansDecGet(&rans, prob_bits);
+ assert(k >= ESCAPE_LIMIT);
+ //printf("reading escape symbol, rans state = %08x\n", rans);
+ RansDecAdvance(&rans, &rans_ptr, k, 1, prob_bits);
+ }
+ if (k != 0) {
+ if (sign_bits_left == 0) {
+ sign_buf = *sign_ptr++;
+ sign_bits_left = 8;
+ }
+ if (sign_buf & 1) k = -k;
+ --sign_bits_left;
+ sign_buf >>= 1;
+ }
+
+ // reverse
+ int reversed_yb = yb ^ 8;
+ int reversed_xb = WIDTH - 8 - xb;
+ coeff[(reversed_yb + y) * WIDTH + (reversed_xb + x)] = k;
+ // printf("coeff %d xb,yb=%d,%d: decoded %d\n", y * 8 + x, reversed_xb, reversed_yb, k);
+ }
+ }
+
+#if 0
+ for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
+ for (unsigned xb = 0; xb < WIDTH; xb += 8) {
+#endif
+ }
+ }
+ fclose(fp);
+
+ // DC coefficient pred from the right to left
+ for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
+ for (int xb = WIDTH - 16; xb >= 0; xb -= 8) {
+ coeff[yb * WIDTH + xb] += coeff[yb * WIDTH + (xb + 8)];
+ }
+ }
+
+ // IDCT
+ for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
+ for (unsigned xb = 0; xb < WIDTH; xb += 8) {
+ // Read one block
+ short in[64];
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 8; ++x) {
+ int k = coeff[(yb + y) * WIDTH + (xb + x)];
+ in[y * 8 + x] = unquantize(k, y * 8 + x);
+ printf("%3d ", in[y * 8 + x]);
+ }
+ printf("\n");
+ }
+ printf("\n");
+
+ idct_int32(in);
+
+ // Clamp and move back
+ for (unsigned y = 0; y < 8; ++y) {
+ for (unsigned x = 0; x < 8; ++x) {
+ int k = in[y * 8 + x];
+ printf("%3d ", k);
+ if (k < 0) k = 0;
+ if (k > 255) k = 255;
+ pix[(yb + y) * WIDTH + (xb + x)] = k;
+ }
+ printf("\n");
+ }
+ printf("\n");
+ }
+ }
+
+ fp = fopen("output.pgm", "wb");
+ fprintf(fp, "P5\n%d %d\n255\n", WIDTH, HEIGHT);
+ fwrite(pix, 1, WIDTH * HEIGHT, fp);
+ fclose(fp);
+}