]> git.sesse.net Git - narabu/commitdiff
Encode sign bit directly in rANS, using some symmetry trickery.
authorSteinar H. Gunderson <sgunderson@bigfoot.com>
Sat, 16 Sep 2017 13:10:19 +0000 (15:10 +0200)
committerSteinar H. Gunderson <sgunderson@bigfoot.com>
Sat, 16 Sep 2017 13:58:08 +0000 (15:58 +0200)
decoder.shader
narabu.cpp
qdc.cpp
qdd.cpp

index 6d54e4dc5671173f540a1574666e51870a14b1c2..512d743cb1307b2932f5a198c81f2cb562a4174a 100644 (file)
@@ -48,14 +48,13 @@ layout(std430, binding = 10) buffer layoutName2
 };
 
 struct CoeffStream {
-       uint src_offset, src_len, sign_offset, sign_len, extra_bits;
+       uint src_offset, src_len;
 };
 layout(std430, binding = 0) buffer whatever3
 {
        CoeffStream streams[];
 };
-
-uniform uint src_offset, src_len, sign_offset, sign_len, extra_bits;
+uniform uint sign_bias_per_model[16];
 
 const uint RANS_BYTE_L = (1u << 23);  // lower bound of our normalization interval
 
@@ -73,7 +72,7 @@ uint get_rans_byte(uint offset)
 //     return bitfieldExtract(data_SSBO[offset >> 2], 8 * int(offset & 3u), 8);
 }
 
-void RansDecInit(out uint r, inout uint offset)
+uint RansDecInit(inout uint offset)
 {
        uint x;
 
@@ -83,7 +82,7 @@ void RansDecInit(out uint r, inout uint offset)
        x |= get_rans_byte(offset + 3) << 24;
        offset += 4;
 
-       r = x;
+       return x;
 }
 
 uint RansDecGet(uint r, uint scale_bits)
@@ -206,16 +205,11 @@ void main()
        const uint stream_num = coeff_num * num_blocks + block_row;
        //const uint stream_num = block_row * num_blocks + coeff_num;  // HACK
        const uint model_num = min((coeff_num % 8) + (coeff_num / 8), 7);
+       const uint sign_bias = sign_bias_per_model[model_num];
 
        // Initialize rANS decoder.
        uint offset = streams[stream_num].src_offset;
-       uint rans;
-       RansDecInit(rans, offset);
-
-       // Initialize sign bit decoder. TODO: this ought to be 32-bit-aligned instead!
-       uint soffset = streams[stream_num].sign_offset;
-       uint sign_buf = get_rans_byte(soffset++) >> streams[stream_num].extra_bits;
-       uint sign_bits_left = 8 - streams[stream_num].extra_bits;
+       uint rans = RansDecInit(offset);
 
        float q = (coeff_num == 0) ? 1.0 : (quant_matrix[coeff_num] * quant_scalefac / 128.0 / sqrt(2.0));  // FIXME: fold
        q *= (1.0 / 255.0);
@@ -234,22 +228,23 @@ void main()
                // rANS decode one coefficient across eight blocks (so 64x8 coefficients).
                for (uint subblock_idx = 8; subblock_idx --> 0; ) {
                        // Read a symbol.
-                       int k = int(cum2sym(RansDecGet(rans, prob_bits), model_num));
+                       uint bottom_bits = RansDecGet(rans, prob_bits + 1);
+                       bool sign = false;
+                       if (bottom_bits >= sign_bias) {
+                               bottom_bits -= sign_bias;
+                               rans -= sign_bias;
+                               sign = true;
+                       }
+                       int k = int(cum2sym(bottom_bits, model_num));  // Can go out-of-bounds; that will return zero.
                        uvec2 sym = get_dsym(k, model_num);
-                       RansDecAdvance(rans, offset, sym.x, sym.y, prob_bits);
+                       RansDecAdvance(rans, offset, sym.x, sym.y, prob_bits + 1);
 
                        if (k == ESCAPE_LIMIT) {
                                k = int(RansDecGet(rans, prob_bits));
                                RansDecAdvance(rans, offset, k, 1, prob_bits);
                        }
-                       if (k != 0) {
-                               if (sign_bits_left == 0) {
-                                       sign_buf = get_rans_byte(soffset++);
-                                       sign_bits_left = 8;
-                               }
-                               if ((sign_buf & 1u) == 1u) k = -k;
-                               --sign_bits_left;
-                               sign_buf >>= 1;
+                       if (sign) {
+                               k = -k;
                        }
 
                        if (coeff_num == 0) {
index 48973664e66a2e3c51f7ed191191186a596e292b..ba1e967e8c95884d36b93e16f6370ad921c1c2f6 100644 (file)
@@ -4,15 +4,18 @@
 #include <SDL2/SDL_error.h>
 #include <SDL2/SDL_video.h>
 #include <epoxy/gl.h>
+#include <movit/util.h>
 #include <string>
 #include <optional>
 #include <algorithm>
 #include <vector>
 #include <memory>
+#include <chrono>
 
 #include "util.h"
 
 using namespace std;
+using namespace std::chrono;
 
 #define WIDTH 1280
 #define HEIGHT 720
@@ -51,7 +54,7 @@ optional<uint32_t> read_varint(const char **ptr, const char *end)
 }
 
 struct CoeffStream {
-        uint src_offset, src_len, sign_offset, sign_len, extra_bits;
+        uint src_offset, src_len;
 };
 CoeffStream streams[45 * 64];  // HACK
 
@@ -83,7 +86,7 @@ int main(int argc, char **argv)
        glGetIntegerv(GL_MAX_COMPUTE_SHARED_MEMORY_SIZE, &size);
        printf("shared_memory_size=%u\n", size);
 
-       string shader_src = read_file("decoder-pre-sign.shader");
+       string shader_src = ::read_file("decoder.shader");
        GLuint shader_num = compile_shader(shader_src, GL_COMPUTE_SHADER);
        GLuint glsl_program_num = glCreateProgram();
        glAttachShader(glsl_program_num, shader_num);
@@ -100,9 +103,10 @@ int main(int argc, char **argv)
 
        glUseProgram(glsl_program_num);
 
-       string coded = read_file(argc >= 2 ? argv[1] : "coded.dat");
+       string coded = ::read_file(argc >= 2 ? argv[1] : "coded.dat");
        const char *ptr = &coded[0];
        const char *end = ptr + coded.size();
+       GLuint sign_bias[NUM_TABLES];
 
 //     printf("first few bytes offs=%zu: %d %d %d %d %d %d %d %d\n", ptr - coded.data(),
 //             (uint8_t)ptr[0], (uint8_t)ptr[1], (uint8_t)ptr[2], (uint8_t)ptr[3],
@@ -118,12 +122,15 @@ int main(int argc, char **argv)
                                exit(1);
                        }
 
-                       decode_tables[table].dsyms[sym].sym_start = cum_freq;
-                       decode_tables[table].dsyms[sym].sym_freq = *freq;
+                       decode_tables[table].dsyms[(sym + 1) & 255].sym_start = cum_freq;
+                       decode_tables[table].dsyms[(sym + 1) & 255].sym_freq = *freq;
                        for (uint32_t i = 0; i < freq; ++i) {
-                               decode_tables[table].cum2sym[cum_freq++] = sym;
+                               if (cum_freq < prob_scale)
+                                       decode_tables[table].cum2sym[cum_freq] = (sym + 1) & 255;
+                               ++cum_freq;
                        }
                }
+               sign_bias[table] = cum_freq;
        }
 
        // Make cum2sym texture.
@@ -141,6 +148,7 @@ int main(int argc, char **argv)
         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
         glTexImage2D(GL_TEXTURE_2D, 0, GL_R8UI, prob_scale, NUM_TABLES, 0, GL_RED_INTEGER, GL_UNSIGNED_BYTE, cum2sym_data.get());
+       check_error();
 
        // Make dsyms texture.
        unique_ptr<pair<uint16_t, uint16_t>[]> dsyms_data(new pair<uint16_t, uint16_t>[NUM_SYMS * NUM_TABLES]);
@@ -158,6 +166,21 @@ int main(int argc, char **argv)
         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
         glTexImage2D(GL_TEXTURE_2D, 0, GL_RG16UI, NUM_SYMS, NUM_TABLES, 0, GL_RG_INTEGER, GL_UNSIGNED_SHORT, dsyms_data.get());
+       check_error();
+
+       GLuint coeff_tex;
+       glGenTextures(1, &coeff_tex);
+        glBindTexture(GL_TEXTURE_2D, coeff_tex);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+       check_error();
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+       check_error();
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
+       check_error();
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
+       check_error();
+        glTexImage2D(GL_TEXTURE_2D, 0, GL_R16I, 1280, 720, 0, GL_RED_INTEGER, GL_SHORT, nullptr);
+       check_error();
 
        GLuint out_tex;
        glGenTextures(1, &out_tex);
@@ -168,22 +191,25 @@ int main(int argc, char **argv)
         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
         glTexImage2D(GL_TEXTURE_2D, 0, GL_R8, 1280, 720, 0, GL_RED, GL_UNSIGNED_BYTE, nullptr);
         //glTexImage2D(GL_TEXTURE_2D, 0, GL_R32F, 1280, 720, 0, GL_RED, GL_FLOAT, nullptr);
+       check_error();
 
-       //GLint src_offset_pos = glGetUniformLocation(glsl_program_num, "src_offset");
-       //GLint sign_offset_pos = glGetUniformLocation(glsl_program_num, "sign_offset");
-       //GLint extra_bits_pos = glGetUniformLocation(glsl_program_num, "extra_bits");
        GLint cum2sym_tex_pos = glGetUniformLocation(glsl_program_num, "cum2sym_tex");
        GLint dsyms_tex_pos = glGetUniformLocation(glsl_program_num, "dsyms_tex");
        GLint out_tex_pos = glGetUniformLocation(glsl_program_num, "out_tex");
-       printf("%d err=0x%x pos=%d,%d,%d\n", __LINE__, glGetError(), cum2sym_tex_pos, dsyms_tex_pos, out_tex_pos);
+       GLint coeff_tex_pos = glGetUniformLocation(glsl_program_num, "coeff_tex");
+       GLint sign_bias_pos = glGetUniformLocation(glsl_program_num, "sign_bias_per_model");
+       printf("%d err=0x%x pos=%d,%d,%d,%d\n", __LINE__, glGetError(), cum2sym_tex_pos, dsyms_tex_pos, out_tex_pos, sign_bias_pos);
 
        // Bind the textures.
        glUniform1i(cum2sym_tex_pos, 0);
        glUniform1i(dsyms_tex_pos, 1);
        glUniform1i(out_tex_pos, 2);
+       glUniform1i(coeff_tex_pos, 3);
+       glUniform1uiv(sign_bias_pos, 16, sign_bias);
         glBindImageTexture(0, cum2sym_tex, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R8UI);
         glBindImageTexture(1, dsyms_tex, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG16UI);
         glBindImageTexture(2, out_tex, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_R8);
+        glBindImageTexture(3, coeff_tex, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_R16I);
        printf("%d err=0x%x\n", __LINE__, glGetError());
 
        // Decode all luma blocks.
@@ -206,20 +232,6 @@ int main(int argc, char **argv)
                                // TODO: check len
                                ptr += *num_rans_bytes;
 
-                               optional<uint32_t> num_sign_bytes = read_varint(&ptr, end);
-                               if (!num_sign_bytes) {
-                                       fprintf(stderr, "Error parsing varint for block %d rANS bytes\n", yb);
-                                       exit(1);
-                               }
-
-                               stream->sign_offset = ptr - coded.data();
-                               stream->sign_len = *num_sign_bytes >> 3;
-                               stream->extra_bits = *num_sign_bytes & 0x7;
-
-                               // TODO: check len
-                               // TODO: free bits
-                               ptr += *num_sign_bytes >> 3;
-
                                //printf("read %d rANS bytes, %d sign bytes\n", *num_rans_bytes, *num_sign_bytes);
                        }
                }
@@ -244,13 +256,19 @@ int main(int argc, char **argv)
 
        glGenBuffers(1, &ssbo_out);
        glBindBuffer(GL_SHADER_STORAGE_BUFFER, ssbo_out);
-       glBufferData(GL_SHADER_STORAGE_BUFFER, 16384, nullptr, GL_STREAM_DRAW);  // ??
+       glBufferData(GL_SHADER_STORAGE_BUFFER, 65536, nullptr, GL_STREAM_DRAW);  // ??
        glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 10, ssbo_out);
+       check_error();
 
-       for (int i = 0; i < 10000; ++i)
-       glDispatchCompute(1, 45, 1);
+#define PARALLEL_SLICES 1
+       steady_clock::time_point start = steady_clock::now();
+       for (int i = 0; i < 1000; ++i)
+               glDispatchCompute(1, (45+PARALLEL_SLICES-1)/PARALLEL_SLICES, 1);
+       check_error();
+       glFinish();
+       steady_clock::time_point now = steady_clock::now();
 
-       unsigned *timing = (unsigned *)glMapBufferRange(GL_SHADER_STORAGE_BUFFER, 0, 16384, GL_MAP_READ_BIT);
+       unsigned *timing = (unsigned *)glMapBufferRange(GL_SHADER_STORAGE_BUFFER, 0, 65536, GL_MAP_READ_BIT);
        //setlocale(LC_ALL, "nb_NO.UTF-8");
 
        string phases[] = {
@@ -269,7 +287,7 @@ int main(int argc, char **argv)
        for (int i = 0; i < 10; ++i) {
                //printf("%d: %'18.0f  [%s]\n", i, double((uint64_t(timing[i * 2 + 1]) << 32) | timing[i * 2]), phases[i].c_str());
                printf("%d,%s", i, phases[i].c_str());
-               for (int j = 0; j < 64; ++j) {
+               for (int j = 0; j < 512; ++j) {
                        int idx = (j * 10 + i) * 2;
                        uint64_t val = (uint64_t(timing[idx + 1]) << 32) | timing[idx];
                //      printf(" %'18.0f", double(val));
@@ -283,6 +301,7 @@ int main(int argc, char **argv)
 
        unsigned char *data = new unsigned char[1280 * 720];
        glGetTexImage(GL_TEXTURE_2D, 0, GL_RED, GL_UNSIGNED_BYTE, data);
+       check_error();
        printf("%d err=0x%x bufsize=%zu\n", __LINE__, glGetError(), coded.size());
 
 #if 0
@@ -321,8 +340,27 @@ int main(int argc, char **argv)
                }
        }
        fclose(fp);
+
+       int16_t *coeff_data = new int16_t[1280 * 720];
+        glBindTexture(GL_TEXTURE_2D, coeff_tex);
+       check_error();
+       glGetTexImage(GL_TEXTURE_2D, 0, GL_RED_INTEGER, GL_SHORT, coeff_data);
+       check_error();
+       for (int k = 0; k < 4; ++k) {
+               for (int y = 0; y < 8; ++y) {
+                       for (int x = 0; x < 8; ++x) {
+                               printf("%3d ", coeff_data[y * 1280 + x + k*8]);
+                       }
+                       printf("\n");
+               }
+               printf("\n");
+       }
+       printf("\n");
+       
        
+       check_error();
        glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0); // unbind
        
        printf("foo = 0x%x\n", glGetError());
+       printf("Each iteration took %.3f ms.\n", 1e3 * duration<double>(now - start).count() / 1000);
 }
diff --git a/qdc.cpp b/qdc.cpp
index 43fcba0de2dc343fbdee87d762992fa00c9ff95b..a0fe2e7b1a71a36f8f9fe0682cf710218669ca85 100644 (file)
--- a/qdc.cpp
+++ b/qdc.cpp
@@ -1,6 +1,7 @@
 #include <stdio.h>
 #include <stdint.h>
 #include <stdlib.h>
+#include <string.h>
 #include <assert.h>
 #include <math.h>
 
@@ -9,6 +10,7 @@
 #include "ryg_rans/renormalize.h"
 
 #include <memory>
+#include <vector>
 
 #define WIDTH 1280
 #define HEIGHT 720
@@ -67,7 +69,6 @@ struct SymbolStats
     uint32_t cum_freqs[NUM_SYMS + 1];
 
     void clear();
-    void count_freqs(uint8_t const* in, size_t nbytes);
     void calc_cum_freqs();
     void normalize_freqs(uint32_t target_total);
 };
@@ -78,14 +79,6 @@ void SymbolStats::clear()
         freqs[i] = 0;
 }
 
-void SymbolStats::count_freqs(uint8_t const* in, size_t nbytes)
-{
-    clear();
-
-    for (size_t i=0; i < nbytes; i++)
-        freqs[in[i]]++;
-}
-
 void SymbolStats::calc_cum_freqs()
 {
     cum_freqs[0] = 0;
@@ -93,119 +86,6 @@ void SymbolStats::calc_cum_freqs()
         cum_freqs[i+1] = cum_freqs[i] + freqs[i];
 }
 
-static double cache[NUM_SYMS + 1][prob_scale + 1];
-static double log2cache[prob_scale + 1];
-static int64_t cachefill = 0;
-
-double find_optimal_cost(const uint32_t *cum_freqs, int sym_to, int available_slots)
-{
-       assert(sym_to >= 0);
-
-       while (sym_to > 0 && cum_freqs[sym_to] == cum_freqs[sym_to - 1]) { --sym_to; }
-       if (cache[sym_to][available_slots] >= 0.0) {
-               //printf("CACHE: %d,%d\n", sym_to, available_slots);
-               return cache[sym_to][available_slots];
-       }
-       if (sym_to == 0) {
-               return 0.0;
-       }
-       if (sym_to == 1) {
-               return cum_freqs[0] * log2cache[available_slots];
-       }
-       if (available_slots == 1) {
-               return cum_freqs[0] * log2cache[1] + find_optimal_cost(cum_freqs, sym_to - 1, 0);
-       }
-
-//     printf("UNCACHE: %d,%d\n", sym_to, available_slots);
-#if 0
-       // ok, test all possible options for the last symbol (TODO: save the choice)
-       double best_so_far = HUGE_VAL;
-       //for (int i = num_syms - 1; i < available_slots; ++i) {
-       double f = freqs[sym_to - 1];
-       for (int i = available_slots; i --> 0; ) {
-               double cost1 = f * log2cache[available_slots - i];
-               double cost2 = find_optimal_cost(freqs, sym_to - 1, i);
-
-               if (sym_to == 3 && available_slots == 838) {
-                       printf("%d %f\n", i, cost1 + cost2);
-               } else
-               if (cost1 + cost2 > best_so_far) {
-                       break;
-               }
-               best_so_far = cost1 + cost2;
-       }
-#elif 1
-       // Minimize the number of total bits spent as a function of how many slots
-       // we assign to this symbol.
-       //
-       // The cost function is convex (I don't know how to prove it, but it makes
-       // intuitively a lot of sense). Find a reasonable guess and see what way
-       // we should search, then iterate until we either hit the end or we start
-       // increasing again.
-       double f = cum_freqs[sym_to - 1] - cum_freqs[sym_to - 2];
-       double start = lrint(available_slots * f / cum_freqs[sym_to - 1]);
-
-       int x1 = std::max<int>(floor(start), 1);
-       int x2 = x1 + 1;
-
-       double f1 = f * log2cache[x1] + find_optimal_cost(cum_freqs, sym_to - 1, available_slots - x1);
-       double f2 = f * log2cache[x2] + find_optimal_cost(cum_freqs, sym_to - 1, available_slots - x2);
-
-       int x, direction;  // -1 or +1
-       double best_so_far = std::min(f1, f2);
-       if (isinf(f1) && isinf(f2)) {
-               // The cost isn't infinite due to the first term, so we need to go downwards
-               // to give the second term more room to breathe.
-               x = x1;
-               direction = -1;
-       } else if (f1 < f2) {
-               x = x1;
-               direction = -1;
-       } else {
-               x = x2;
-               direction = 1;
-       }
-
-       //printf("[%d,%d] freq=%ld cumfreq=%d From %d and %d, chose %d [%f] and direction=%d\n",
-       //      sym_to, available_slots, freqs[sym_to - 1], cum_freqs[sym_to - 1], x1, x2, x, best_so_far, direction);
-
-       while ((x + direction) > 0 && (x + direction) <= available_slots) {
-               x += direction;
-               double fn = f * log2cache[x] + find_optimal_cost(cum_freqs, sym_to - 1, available_slots - x);
-       //      printf("[%d,%d] %d is %f\n", sym_to, available_slots, x, fn);
-               if (fn > best_so_far) {
-                       break;
-               }
-               best_so_far = fn;
-       }
-#endif
-       if (++cachefill % 131072 == 0) {
-       //      printf("%d,%d = %f (cachefill = %.2f%%)\n", sym_to, available_slots, best_so_far,
-       //              100.0 * (cachefill / double((NUM_SYMS + 1) * (prob_scale + 1))));
-       }
-       assert(best_so_far >= 0.0);
-       assert(cache[sym_to][available_slots] < 0.0);
-       cache[sym_to][available_slots] = best_so_far;
-       return best_so_far;
-}
-
-double find_optimal_cost(const uint32_t *cum_freqs, const uint64_t *freqs)
-{
-       for (int j = 0; j <= NUM_SYMS; ++j) {
-               for (unsigned k = 0; k <= prob_scale; ++k) {
-                       cache[j][k] = -1.0;
-               }
-       }
-       for (unsigned k = 0; k <= prob_scale; ++k) {
-               log2cache[k] = -log2(k * (1.0 / prob_scale));
-               //printf("log2cache[%d] = %f\n", k, log2cache[k]);
-       }
-       cachefill = 0;
-       double ret = find_optimal_cost(cum_freqs, NUM_SYMS, prob_scale);
-       printf("Used %ld function invocations\n", cachefill);
-       return ret;
-}
-
 void SymbolStats::normalize_freqs(uint32_t target_total)
 {
     uint64_t real_freq[NUM_SYMS + 1];  // hack
@@ -227,47 +107,6 @@ void SymbolStats::normalize_freqs(uint32_t target_total)
 
     OptimalRenormalize(cum_freqs, NUM_SYMS, prob_scale);
 
-#if 0
-    double optimal_cost = find_optimal_cost(cum_freqs + 1, real_freq + 1);
-
-    // resample distribution based on cumulative freqs
-    for (int i = 1; i <= NUM_SYMS; i++)
-        //cum_freqs[i] = ((uint64_t)target_total * cum_freqs[i])/cur_total;
-        cum_freqs[i] = lrint(cum_freqs[i] * double(target_total) / cur_total);
-
-    // if we nuked any non-0 frequency symbol to 0, we need to steal
-    // the range to make the frequency nonzero from elsewhere.
-    //
-    // this is not at all optimal, i'm just doing the first thing that comes to mind.
-    for (int i=0; i < NUM_SYMS; i++) {
-        if (freqs[i] && cum_freqs[i+1] == cum_freqs[i]) {
-            // symbol i was set to zero freq
-
-            // find best symbol to steal frequency from (try to steal from low-freq ones)
-            uint32_t best_freq = ~0u;
-            int best_steal = -1;
-            for (int j=0; j < NUM_SYMS; j++) {
-                uint32_t freq = cum_freqs[j+1] - cum_freqs[j];
-                if (freq > 1 && freq < best_freq) {
-                    best_freq = freq;
-                    best_steal = j;
-                }
-            }
-            assert(best_steal != -1);
-
-            // and steal from it!
-            if (best_steal < i) {
-                for (int j = best_steal + 1; j <= i; j++)
-                    cum_freqs[j]--;
-            } else {
-                assert(best_steal > i);
-                for (int j = i + 1; j <= best_steal; j++)
-                    cum_freqs[j]++;
-            }
-        }
-    }
-#endif
-
     // calculate updated freqs and make sure we didn't screw anything up
     assert(cum_freqs[0] == 0 && cum_freqs[NUM_SYMS] == target_total);
     for (int i=0; i < NUM_SYMS; i++) {
@@ -302,6 +141,7 @@ SymbolStats stats[64];
 
 int pick_stats_for(int y, int x)
 {
+       //return 0;
        //return std::min<int>(hypot(x, y), 7);
        return std::min<int>(x + y, 7);
        //if (x + y >= 7) return 7;
@@ -331,26 +171,23 @@ public:
        RansEncoder()
        {
                out_buf.reset(new uint8_t[out_max_size]);
-               sign_buf.reset(new uint8_t[max_num_sign]);
                clear();
        }
 
-       void init_prob(const SymbolStats &s)
+       void init_prob(SymbolStats &s)
        {
                for (int i = 0; i < NUM_SYMS; i++) {
-                       printf("%d: cumfreqs=%d freqs=%d prob_bits=%d\n", i, s.cum_freqs[i], s.freqs[i], prob_bits);
-                       RansEncSymbolInit(&esyms[i], s.cum_freqs[i], s.freqs[i], prob_bits);
+                       printf("%d: cumfreqs=%d freqs=%d prob_bits=%d\n", i, s.cum_freqs[i], s.freqs[i], prob_bits + 1);
+                       RansEncSymbolInit(&esyms[i], s.cum_freqs[i], s.freqs[i], prob_bits + 1);
                }
+               sign_bias = s.cum_freqs[256];
        }
 
        void clear()
        {
                out_end = out_buf.get() + out_max_size;
-               sign_end = sign_buf.get() + max_num_sign;
                ptr = out_end; // *end* of output buffer
-               sign_ptr = sign_end; // *end* of output buffer
                RansEncInit(&rans);
-               free_sign_bits = 0;
        }
 
        uint32_t save_block(FILE *codedfp)  // Returns number of bytes.
@@ -359,33 +196,41 @@ public:
                //printf("post-flush = %08x\n", rans);
 
                uint32_t num_rans_bytes = out_end - ptr;
+#if 0
+               if (num_rans_bytes == 4) {
+                       uint32_t block;
+                       memcpy(&block, ptr, 4);
+
+                       if (block == last_block) {
+                               write_varint(0, codedfp);
+                               clear();
+                               return 1;
+                       }
+
+                       last_block = block;
+               } else {
+                       last_block = 0;
+               }
+#endif
+
                write_varint(num_rans_bytes, codedfp);
                //fwrite(&num_rans_bytes, 1, 4, codedfp);
                fwrite(ptr, 1, num_rans_bytes, codedfp);
 
                //printf("first rANS bytes: %02x %02x %02x %02x %02x %02x %02x %02x\n", ptr[0], ptr[1], ptr[2], ptr[3], ptr[4], ptr[5], ptr[6], ptr[7]);
 
-               if (free_sign_bits > 0) {
-                       *sign_ptr <<= free_sign_bits;
-               }
-
-#if 1
-               uint32_t num_sign_bytes = sign_end - sign_ptr;
-               write_varint((num_sign_bytes << 3) | free_sign_bits, codedfp);
-               fwrite(sign_ptr, 1, num_sign_bytes, codedfp);
-#endif
 
                clear();
 
-               printf("Saving block: %d rANS bytes, %d sign bytes\n", num_rans_bytes, num_sign_bytes);
-               return num_rans_bytes + num_sign_bytes;
+               printf("Saving block: %d rANS bytes\n", num_rans_bytes);
+               return num_rans_bytes;
                //return num_rans_bytes;
        }
 
        void encode_coeff(short signed_k)
        {
-               printf("encoding coeff %d\n", signed_k);
-               short k = abs(signed_k);
+               //printf("encoding coeff %d (sym %d), rans before encoding = %08x\n", signed_k, ((abs(signed_k) - 1) & 255), rans);
+               unsigned short k = abs(signed_k);
                if (k >= ESCAPE_LIMIT) {
                        // Put the coefficient as a 1/(2^12) symbol _before_
                        // the 255 coefficient, since the decoder will read the
@@ -393,33 +238,24 @@ public:
                        RansEncPut(&rans, &ptr, k, 1, prob_bits);
                        k = ESCAPE_LIMIT;
                }
-               if (k != 0) {
-#if 1
-                       if (free_sign_bits == 0) {
-                               --sign_ptr;
-                               *sign_ptr = 0;
-                               free_sign_bits = 8;
-                       }
-                       *sign_ptr <<= 1;
-                       *sign_ptr |= (signed_k < 0);
-                       --free_sign_bits;
-#else
-                       RansEncPut(&rans, &ptr, (k < 0) ? prob_scale / 2 : 0, prob_scale / 2, prob_bits);
-#endif
+               RansEncPutSymbol(&rans, &ptr, &esyms[(k - 1) & 255]);
+               if (signed_k < 0) {
+                       rans += sign_bias;
                }
-               RansEncPutSymbol(&rans, &ptr, &esyms[k]);
        }
 
 private:
        static constexpr size_t out_max_size = 32 << 20; // 32 MB.
        static constexpr size_t max_num_sign = 1048576;  // Way too big. And actually bytes.
 
-       unique_ptr<uint8_t[]> out_buf, sign_buf;
-       uint8_t *out_end, *sign_end;
-       uint8_t *ptr, *sign_ptr;
+       unique_ptr<uint8_t[]> out_buf;
+       uint8_t *out_end;
+       uint8_t *ptr;
        RansState rans;
-       size_t free_sign_bits;
        RansEncSymbol esyms[NUM_SYMS];
+       uint32_t sign_bias;
+
+       uint32_t last_block = 0;  // Not a valid 4-byte rANS block (?)
 };
 
 static constexpr int dc_scalefac = 8;  // Matches the FDCT's gain.
@@ -778,7 +614,7 @@ int main(int argc, char **argv)
        fclose(fp);
 
        // For each coefficient, make some tables.
-       size_t extra_bits = 0, sign_bits = 0;
+       size_t extra_bits = 0;
        for (unsigned i = 0; i < 64; ++i) {
                stats[i].clear();
        }
@@ -791,20 +627,19 @@ int main(int argc, char **argv)
                        // Luma
                        for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
                                for (unsigned xb = 0; xb < WIDTH; xb += 8) {
-                                       short k = abs(coeff_y[(yb + y) * WIDTH + (xb + x)]);
+                                       unsigned short k = abs(coeff_y[(yb + y) * WIDTH + (xb + x)]);
                                        if (k >= ESCAPE_LIMIT) {
                                                k = ESCAPE_LIMIT;
                                                extra_bits += 12;  // escape this one
                                        }
-                                       if (k != 0) ++sign_bits;
-                                       ++s_luma.freqs[k];
+                                       ++s_luma.freqs[(k - 1) & 255];
                                }
                        }
                        // Chroma
                        for (unsigned yb = 0; yb < HEIGHT; yb += 8) {
                                for (unsigned xb = 0; xb < WIDTH/2; xb += 8) {
-                                       short k_cb = abs(coeff_cb[(yb + y) * WIDTH/2 + (xb + x)]);
-                                       short k_cr = abs(coeff_cr[(yb + y) * WIDTH/2 + (xb + x)]);
+                                       unsigned short k_cb = abs(coeff_cb[(yb + y) * WIDTH/2 + (xb + x)]);
+                                       unsigned short k_cr = abs(coeff_cr[(yb + y) * WIDTH/2 + (xb + x)]);
                                        if (k_cb >= ESCAPE_LIMIT) {
                                                k_cb = ESCAPE_LIMIT;
                                                extra_bits += 12;  // escape this one
@@ -813,16 +648,17 @@ int main(int argc, char **argv)
                                                k_cr = ESCAPE_LIMIT;
                                                extra_bits += 12;  // escape this one
                                        }
-                                       if (k_cb != 0) ++sign_bits;
-                                       if (k_cr != 0) ++sign_bits;
-                                       ++s_chroma.freqs[k_cb];
-                                       ++s_chroma.freqs[k_cr];
+                                       ++s_chroma.freqs[(k_cb - 1) & 255];
+                                       ++s_chroma.freqs[(k_cr - 1) & 255];
                                }
                        }
                }
        }
        for (unsigned i = 0; i < 64; ++i) {
+               stats[i].freqs[255] /= 2;  // zero, has no sign bits (yes, this is trickery)
                stats[i].normalize_freqs(prob_scale);
+               stats[i].cum_freqs[256] += stats[i].freqs[255];
+               stats[i].freqs[255] *= 2;
        }
 
        FILE *codedfp = fopen("coded.dat", "wb");
@@ -853,6 +689,7 @@ int main(int argc, char **argv)
                        rans_encoder.init_prob(s_luma);
 
                        // Luma
+                       std::vector<int> lens;
 
                        // need to reverse later
                        rans_encoder.clear();
@@ -864,7 +701,9 @@ int main(int argc, char **argv)
                                        rans_encoder.encode_coeff(k);
                                }
                                if (yb % 16 == 8) {
-                                       num_bytes += rans_encoder.save_block(codedfp);
+                                       int l = rans_encoder.save_block(codedfp);
+                                       num_bytes += l;
+                                       lens.push_back(l);
                                }
                        }
                        if (HEIGHT % 16 != 0) {
@@ -872,6 +711,19 @@ int main(int argc, char **argv)
                        }
                        tot_bytes += num_bytes;
                        printf("coeff %d Y': %ld bytes\n", y * 8 + x, num_bytes);
+
+                       double sum_l = 0.0;
+                       for (int l : lens) {
+                               sum_l += l;
+                       }
+                       double avg_l = sum_l / lens.size();
+
+                       double sum_sql = 0.0;
+                       for (int l : lens) {
+                               sum_sql += (l - avg_l) * (l - avg_l);
+                       }
+                       double stddev_l = sqrt(sum_sql / (lens.size() - 1));
+                       printf("coeff %d: avg=%.2f bytes, stddev=%.2f bytes\n", y*8+x, avg_l, stddev_l);
                }
        }
 
@@ -927,10 +779,8 @@ int main(int argc, char **argv)
                }
        }
 
-       printf("%ld bytes + %ld sign bits (%ld) + %ld escape bits (%ld) = %ld total bytes\n",
-               tot_bytes - sign_bits / 8 - extra_bits / 8,
-               sign_bits,
-               sign_bits / 8,
+       printf("%ld bytes + %ld escape bits (%ld) = %ld total bytes\n",
+               tot_bytes - extra_bits / 8,
                extra_bits,
                extra_bits / 8,
                tot_bytes);
diff --git a/qdd.cpp b/qdd.cpp
index 8c8196dfdcab8bc537901c27d5ff872d08581898..a4fa907a297549f03d0860bb016d41fb7bf83420 100644 (file)
--- a/qdd.cpp
+++ b/qdd.cpp
@@ -103,16 +103,21 @@ int main(void)
                exit(1);
        }
 
+       uint32_t sign_bias[16];
        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);
+                       RansDecSymbolInit(&decode_tables[table].dsyms[(sym + 1) & 255], cum_freq, freq);
                        for (uint32_t i = 0; i < freq; ++i) {
-                               decode_tables[table].cum2sym[cum_freq++] = sym;
+                               if (cum_freq < prob_scale)
+                                       decode_tables[table].cum2sym[cum_freq] = (sym + 1) & 255;
+                               ++cum_freq;
                        }
                }
+               sign_bias[table] = cum_freq;
+               printf("sign_bias=%u (of %d)\n", sign_bias[table], prob_scale * 2);
        }
 
        // loop over all coefficients
@@ -122,13 +127,8 @@ int main(void)
                
                        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) {
@@ -138,22 +138,7 @@ int main(void)
                                        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("%d,%d: read %d rANS bytes\n", x, y, num_rans_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]);
 
 
@@ -162,25 +147,24 @@ int main(void)
                                        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);
+                                       uint32_t bottom_bits = RansDecGet(&rans, prob_bits + 1);
+                                       uint32_t sign = 0;
+                                       if (bottom_bits >= sign_bias[tbl]) {
+                                               bottom_bits -= sign_bias[tbl];
+                                               rans -= sign_bias[tbl];
+                                               sign = 1;
+                                       }
+                                       uint32_t k = decode_tables[tbl].cum2sym[std::min(bottom_bits, prob_scale - 1)];
+                                       RansDecAdvanceSymbol(&rans, &rans_ptr, &decode_tables[tbl].dsyms[k], prob_bits + 1);
                                        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;
+                                       if (sign) {
+                                               assert(k != 0);
+                                               k = -k;
                                        }
 
                                        // reverse