/*

Copyright (c) 2004 Amir Said (said@ieee.org) & William A. Pearlman (pearlw@ecse.rpi.edu)

All rights reserved.

Redistribution and use in source and binary forms, with or without modification, 

are permitted provided that the following conditions are met:

-   Redistributions of source code must retain the above copyright notice, this list 

    of conditions and the following disclaimer. 

-   Redistributions in binary form must reproduce the above copyright notice, this list of 

    conditions and the following disclaimer in the documentation and/or other materials 

    provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 

EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 

MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL 

THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 

SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 

PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 

OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 

IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING 

IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY 

OF SUCH DAMAGE.

*/

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

//                                                                           -

//                       ****************************                        -

//                        ARITHMETIC CODING EXAMPLES                         -

//                       ****************************                        -

//                                                                           -

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

//                                                                           -

// Fast arithmetic coding implementation                                     -

// -> 32-bit variables, 32-bit product, periodic updates, table decoding     -

//                                                                           -

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

//                                                                           -

// Version 1.00  -  April 25, 2004                                           -

//                                                                           -

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

//                                                                           -

//                                  WARNING                                  -

//                                 =========                                 -

//                                                                           -

// The only purpose of this program is to demonstrate the basic principles   -

// of arithmetic coding. It is provided as is, without any express or        -

// implied warranty, without even the warranty of fitness for any particular -

// purpose, or that the implementations are correct.                         -

//                                                                           -

// Permission to copy and redistribute this code is hereby granted, provided -

// that this warning and copyright notices are not removed or altered.       -

//                                                                           -

// Copyright (c) 2004 by Amir Said (said@ieee.org) &                         -

//                       William A. Pearlman (pearlw@ecse.rpi.edu)           -

//                                                                           -

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

//                                                                           -

// A description of the arithmetic coding method used here is available in   -

//                                                                           -

// Lossless Compression Handbook, ed. K. Sayood                              -

// Chapter 5: Arithmetic Coding (A. Said), pp. 101-152, Academic Press, 2003 -

//                                                                           -

// A. Said, Introduction to Arithetic Coding Theory and Practice             -

// HP Labs report HPL-2004-76  -  http://www.hpl.hp.com/techreports/         -

//                                                                           -

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

// - - Inclusion - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

#include <stdlib.h>

#include "o3dgcArithmeticCodec.h"

namespace o3dgc

{

    // - - Constants - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    const unsigned AC__MinLength = 0x01000000U;   // threshold for renormalization

    const unsigned AC__MaxLength = 0xFFFFFFFFU;      // maximum AC interval length

                                               // Maximum values for binary models

    const unsigned BM__LengthShift = 13;     // length bits discarded before mult.

    const unsigned BM__MaxCount    = 1 << BM__LengthShift;  // for adaptive models

                                              // Maximum values for general models

    const unsigned DM__LengthShift = 15;     // length bits discarded before mult.

    const unsigned DM__MaxCount    = 1 << DM__LengthShift;  // for adaptive models

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    // - - Static functions  - - - - - - - - - - - - - - - - - - - - - - - - - - -

    AI_WONT_RETURN static void AC_Error(const char * msg) AI_WONT_RETURN_SUFFIX;

    static void AC_Error(const char * msg)

    {

      fprintf(stderr, "\n\n -> Arithmetic coding error: ");

      fputs(msg, stderr);

      fputs("\n Execution terminated!\n", stderr);

      getchar();

      exit(1);

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    // - - Coding implementations  - - - - - - - - - - - - - - - - - - - - - - - -

    inline void Arithmetic_Codec::propagate_carry(void)

    {

      unsigned char * p;            // carry propagation on compressed data buffer

      for (p = ac_pointer - 1; *p == 0xFFU; p--) *p = 0;

      ++*p;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    inline void Arithmetic_Codec::renorm_enc_interval(void)

    {

      do {                                          // output and discard top byte

        *ac_pointer++ = (unsigned char)(base >> 24);

        base <<= 8;

      } while ((length <<= 8) < AC__MinLength);        // length multiplied by 256

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    inline void Arithmetic_Codec::renorm_dec_interval(void)

    {

      do {                                          // read least-significant byte

        value = (value << 8) | unsigned(*++ac_pointer);

      } while ((length <<= 8) < AC__MinLength);        // length multiplied by 256

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::put_bit(unsigned bit)

    {

    #ifdef _DEBUG

      if (mode != 1) AC_Error("encoder not initialized");

    #endif

      length >>= 1;                                              // halve interval

      if (bit) {

        unsigned init_base = base;

        base += length;                                               // move base

        if (init_base > base) propagate_carry();               // overflow = carry

      }

      if (length < AC__MinLength) renorm_enc_interval();        // renormalization

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    unsigned Arithmetic_Codec::get_bit(void)

    {

    #ifdef _DEBUG

      if (mode != 2) AC_Error("decoder not initialized");

    #endif

      length >>= 1;                                              // halve interval

      unsigned bit = (value >= length);                              // decode bit

      if (bit) value -= length;                                       // move base

      if (length < AC__MinLength) renorm_dec_interval();        // renormalization

      return bit;                                         // return data bit value

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::put_bits(unsigned data, unsigned bits)

    {

    #ifdef _DEBUG

      if (mode != 1) AC_Error("encoder not initialized");

      if ((bits < 1) || (bits > 20)) AC_Error("invalid number of bits");

      if (data >= (1U << bits)) AC_Error("invalid data");

    #endif

      unsigned init_base = base;

      base += data * (length >>= bits);            // new interval base and length

      if (init_base > base) propagate_carry();                 // overflow = carry

      if (length < AC__MinLength) renorm_enc_interval();        // renormalization

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    unsigned Arithmetic_Codec::get_bits(unsigned bits)

    {

    #ifdef _DEBUG

      if (mode != 2) AC_Error("decoder not initialized");

      if ((bits < 1) || (bits > 20)) AC_Error("invalid number of bits");

    #endif

      unsigned s = value / (length >>= bits);      // decode symbol, change length

      value -= length * s;                                      // update interval

      if (length < AC__MinLength) renorm_dec_interval();        // renormalization

      return s;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::encode(unsigned bit,

                                  Static_Bit_Model & M)

    {

    #ifdef _DEBUG

      if (mode != 1) AC_Error("encoder not initialized");

    #endif

      unsigned x = M.bit_0_prob * (length >> BM__LengthShift);   // product l x p0

                                                                // update interval

      if (bit == 0)

        length  = x;

      else {

        unsigned init_base = base;

        base   += x;

        length -= x;

        if (init_base > base) propagate_carry();               // overflow = carry

      }

      if (length < AC__MinLength) renorm_enc_interval();        // renormalization

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    unsigned Arithmetic_Codec::decode(Static_Bit_Model & M)

    {

    #ifdef _DEBUG

      if (mode != 2) AC_Error("decoder not initialized");

    #endif

      unsigned x = M.bit_0_prob * (length >> BM__LengthShift);   // product l x p0

      unsigned bit = (value >= x);                                     // decision

                                                        // update & shift interval

      if (bit == 0)

        length  = x;

      else {

        value  -= x;                                 // shifted interval base = 0

        length -= x;

      }

      if (length < AC__MinLength) renorm_dec_interval();        // renormalization

      return bit;                                         // return data bit value

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::encode(unsigned bit,

                                  Adaptive_Bit_Model & M)

    {

    #ifdef _DEBUG

      if (mode != 1) AC_Error("encoder not initialized");

    #endif

      unsigned x = M.bit_0_prob * (length >> BM__LengthShift);   // product l x p0

                                                                // update interval

      if (bit == 0) {

        length = x;

        ++M.bit_0_count;

      }

      else {

        unsigned init_base = base;

        base   += x;

        length -= x;

        if (init_base > base) propagate_carry();               // overflow = carry

      }

      if (length < AC__MinLength) renorm_enc_interval();        // renormalization

      if (--M.bits_until_update == 0) M.update();         // periodic model update

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    unsigned Arithmetic_Codec::decode(Adaptive_Bit_Model & M)

    {

    #ifdef _DEBUG

      if (mode != 2) AC_Error("decoder not initialized");

    #endif

      unsigned x = M.bit_0_prob * (length >> BM__LengthShift);   // product l x p0

      unsigned bit = (value >= x);                                     // decision

                                                                // update interval

      if (bit == 0) {

        length = x;

        ++M.bit_0_count;

      }

      else {

        value  -= x;

        length -= x;

      }

      if (length < AC__MinLength) renorm_dec_interval();        // renormalization

      if (--M.bits_until_update == 0) M.update();         // periodic model update

      return bit;                                         // return data bit value

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::encode(unsigned data,

                                  Static_Data_Model & M)

    {

    #ifdef _DEBUG

      if (mode != 1) AC_Error("encoder not initialized");

      if (data >= M.data_symbols) AC_Error("invalid data symbol");

    #endif

      unsigned x, init_base = base;

                                                               // compute products

      if (data == M.last_symbol) {

        x = M.distribution[data] * (length >> DM__LengthShift);

        base   += x;                                            // update interval

        length -= x;                                          // no product needed

      }

      else {

        x = M.distribution[data] * (length >>= DM__LengthShift);

        base   += x;                                            // update interval

        length  = M.distribution[data+1] * length - x;

      }

      if (init_base > base) propagate_carry();                 // overflow = carry

      if (length < AC__MinLength) renorm_enc_interval();        // renormalization

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    unsigned Arithmetic_Codec::decode(Static_Data_Model & M)

    {

    #ifdef _DEBUG

      if (mode != 2) AC_Error("decoder not initialized");

    #endif

      unsigned n, s, x, y = length;

      if (M.decoder_table) {              // use table look-up for faster decoding

        unsigned dv = value / (length >>= DM__LengthShift);

        unsigned t = dv >> M.table_shift;

        s = M.decoder_table[t];         // initial decision based on table look-up

        n = M.decoder_table[t+1] + 1;

        while (n > s + 1) {                        // finish with bisection search

          unsigned m = (s + n) >> 1;

          if (M.distribution[m] > dv) n = m; else s = m;

        }

                                                               // compute products

        x = M.distribution[s] * length;

        if (s != M.last_symbol) y = M.distribution[s+1] * length;

      }

      else {                                  // decode using only multiplications

        x = s = 0;

        length >>= DM__LengthShift;

        unsigned m = (n = M.data_symbols) >> 1;

                                                    // decode via bisection search

        do {

          unsigned z = length * M.distribution[m];

          if (z > value) {

            n = m;

            y = z;                                             // value is smaller

          }

          else {

            s = m;

            x = z;                                     // value is larger or equal

          }

        } while ((m = (s + n) >> 1) != s);

      }

      value -= x;                                               // update interval

      length = y - x;

      if (length < AC__MinLength) renorm_dec_interval();        // renormalization

      return s;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::encode(unsigned data,

                                  Adaptive_Data_Model & M)

    {

    #ifdef _DEBUG

      if (mode != 1) AC_Error("encoder not initialized");

      if (data >= M.data_symbols) 

      {

          AC_Error("invalid data symbol");

      }

    #endif

      unsigned x, init_base = base;

                                                               // compute products

      if (data == M.last_symbol) {

        x = M.distribution[data] * (length >> DM__LengthShift);

        base   += x;                                            // update interval

        length -= x;                                          // no product needed

      }

      else {

        x = M.distribution[data] * (length >>= DM__LengthShift);

        base   += x;                                            // update interval

        length  = M.distribution[data+1] * length - x;

      }

      if (init_base > base) propagate_carry();                 // overflow = carry

      if (length < AC__MinLength) renorm_enc_interval();        // renormalization

      ++M.symbol_count[data];

      if (--M.symbols_until_update == 0) M.update(true);  // periodic model update

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    unsigned Arithmetic_Codec::decode(Adaptive_Data_Model & M)

    {

    #ifdef _DEBUG

      if (mode != 2) AC_Error("decoder not initialized");

    #endif

      unsigned n, s, x, y = length;

      if (M.decoder_table) {              // use table look-up for faster decoding

        unsigned dv = value / (length >>= DM__LengthShift);

        unsigned t = dv >> M.table_shift;

        s = M.decoder_table[t];         // initial decision based on table look-up

        n = M.decoder_table[t+1] + 1;

        while (n > s + 1) {                        // finish with bisection search

          unsigned m = (s + n) >> 1;

          if (M.distribution[m] > dv) n = m; else s = m;

        }

                                                               // compute products

        x = M.distribution[s] * length;

        if (s != M.last_symbol) {

            y = M.distribution[s+1] * length;

        }

      }

      else {                                  // decode using only multiplications

        x = s = 0;

        length >>= DM__LengthShift;

        unsigned m = (n = M.data_symbols) >> 1;

                                                    // decode via bisection search

        do {

          unsigned z = length * M.distribution[m];

          if (z > value) {

            n = m;

            y = z;                                             // value is smaller

          }

          else {

            s = m;

            x = z;                                     // value is larger or equal

          }

        } while ((m = (s + n) >> 1) != s);

      }

      value -= x;                                               // update interval

      length = y - x;

      if (length < AC__MinLength) renorm_dec_interval();        // renormalization

      ++M.symbol_count[s];

      if (--M.symbols_until_update == 0) M.update(false);  // periodic model update

      return s;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    // - - Other Arithmetic_Codec implementations  - - - - - - - - - - - - - - - -

    Arithmetic_Codec::Arithmetic_Codec(void)

    {

      mode = buffer_size = 0;

      new_buffer = code_buffer = 0;

    }

    Arithmetic_Codec::Arithmetic_Codec(unsigned max_code_bytes,

                                       unsigned char * user_buffer)

    {

      mode = buffer_size = 0;

      new_buffer = code_buffer = 0;

      set_buffer(max_code_bytes, user_buffer);

    }

    Arithmetic_Codec::~Arithmetic_Codec(void)

    {

      delete [] new_buffer;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::set_buffer(unsigned max_code_bytes,

                                      unsigned char * user_buffer)

    {

                                                      // test for reasonable sizes

      if (!max_code_bytes)// || (max_code_bytes > 0x10000000U)) // updated by K. Mammou

      {

        AC_Error("invalid codec buffer size");

      }

      if (mode != 0) AC_Error("cannot set buffer while encoding or decoding");

      if (user_buffer != 0) {                       // user provides memory buffer

        buffer_size = max_code_bytes;

        code_buffer = user_buffer;               // set buffer for compressed data

        delete [] new_buffer;                 // free anything previously assigned

        new_buffer = 0;

        return;

      }

      if (max_code_bytes <= buffer_size) return;               // enough available

      buffer_size = max_code_bytes;                           // assign new memory

      delete [] new_buffer;                   // free anything previously assigned

      new_buffer = new unsigned char[buffer_size+16];        // 16 extra bytes

      code_buffer = new_buffer;                  // set buffer for compressed data

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::start_encoder(void)

    {

      if (mode != 0) AC_Error("cannot start encoder");

      if (buffer_size == 0) AC_Error("no code buffer set");

      mode   = 1;

      base   = 0;            // initialize encoder variables: interval and pointer

      length = AC__MaxLength;

      ac_pointer = code_buffer;                       // pointer to next data byte

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::start_decoder(void)

    {

      if (mode != 0) AC_Error("cannot start decoder");

      if (buffer_size == 0) AC_Error("no code buffer set");

                      // initialize decoder: interval, pointer, initial code value

      mode   = 2;

      length = AC__MaxLength;

      ac_pointer = code_buffer + 3;

      value = (unsigned(code_buffer[0]) << 24)|(unsigned(code_buffer[1]) << 16) |

              (unsigned(code_buffer[2]) <<  8)| unsigned(code_buffer[3]);

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::read_from_file(FILE * code_file)

    {

      unsigned shift = 0, code_bytes = 0;

      int file_byte;

                          // read variable-length header with number of code bytes

      do {

        if ((file_byte = getc(code_file)) == EOF)

          AC_Error("cannot read code from file");

        code_bytes |= unsigned(file_byte & 0x7F) << shift;

        shift += 7;

      } while (file_byte & 0x80);

                                                           // read compressed data

      if (code_bytes > buffer_size) AC_Error("code buffer overflow");

      if (fread(code_buffer, 1, code_bytes, code_file) != code_bytes)

        AC_Error("cannot read code from file");

      start_decoder();                                       // initialize decoder

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    unsigned Arithmetic_Codec::stop_encoder(void)

    {

      if (mode != 1) AC_Error("invalid to stop encoder");

      mode = 0;

      unsigned init_base = base;            // done encoding: set final data bytes

      if (length > 2 * AC__MinLength) {

        base  += AC__MinLength;                                     // base offset

        length = AC__MinLength >> 1;             // set new length for 1 more byte

      }

      else {

        base  += AC__MinLength >> 1;                                // base offset

        length = AC__MinLength >> 9;            // set new length for 2 more bytes

      }

      if (init_base > base) propagate_carry();                 // overflow = carry

      renorm_enc_interval();                // renormalization = output last bytes

      unsigned code_bytes = unsigned(ac_pointer - code_buffer);

      if (code_bytes > buffer_size) AC_Error("code buffer overflow");

      return code_bytes;                                   // number of bytes used

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    unsigned Arithmetic_Codec::write_to_file(FILE * code_file)

    {

      unsigned header_bytes = 0, code_bytes = stop_encoder(), nb = code_bytes;

                         // write variable-length header with number of code bytes

      do {

        int file_byte = int(nb & 0x7FU);

        if ((nb >>= 7) > 0) file_byte |= 0x80;

        if (putc(file_byte, code_file) == EOF)

          AC_Error("cannot write compressed data to file");

        header_bytes++;

      } while (nb);

                                                          // write compressed data

      if (fwrite(code_buffer, 1, code_bytes, code_file) != code_bytes)

        AC_Error("cannot write compressed data to file");

      return code_bytes + header_bytes;                              // bytes used

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Arithmetic_Codec::stop_decoder(void)

    {

      if (mode != 2) AC_Error("invalid to stop decoder");

      mode = 0;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    // - Static bit model implementation - - - - - - - - - - - - - - - - - - - - -

    Static_Bit_Model::Static_Bit_Model(void)

    {

      bit_0_prob = 1U << (BM__LengthShift - 1);                        // p0 = 0.5

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Static_Bit_Model::set_probability_0(double p0)

    {

      if ((p0 < 0.0001)||(p0 > 0.9999)) AC_Error("invalid bit probability");

      bit_0_prob = unsigned(p0 * (1 << BM__LengthShift));

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    // - Adaptive bit model implementation - - - - - - - - - - - - - - - - - - - -

    Adaptive_Bit_Model::Adaptive_Bit_Model(void)

    {

      reset();

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Adaptive_Bit_Model::reset(void)

    {

                                           // initialization to equiprobable model

      bit_0_count = 1;

      bit_count   = 2;

      bit_0_prob  = 1U << (BM__LengthShift - 1);

      update_cycle = bits_until_update = 4;         // start with frequent updates

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Adaptive_Bit_Model::update(void)

    {

                                       // halve counts when a threshold is reached

      if ((bit_count += update_cycle) > BM__MaxCount) {

        bit_count = (bit_count + 1) >> 1;

        bit_0_count = (bit_0_count + 1) >> 1;

        if (bit_0_count == bit_count) ++bit_count;

      }

                                               // compute scaled bit 0 probability

      unsigned scale = 0x80000000U / bit_count;

      bit_0_prob = (bit_0_count * scale) >> (31 - BM__LengthShift);

                                                 // set frequency of model updates

      update_cycle = (5 * update_cycle) >> 2;

      if (update_cycle > 64) update_cycle = 64;

      bits_until_update = update_cycle;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    // - - Static data model implementation  - - - - - - - - - - - - - - - - - - -

    Static_Data_Model::Static_Data_Model(void)

    {

      data_symbols = 0;

      distribution = 0;

    }

    Static_Data_Model::~Static_Data_Model(void)

    {

      delete [] distribution;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Static_Data_Model::set_distribution(unsigned number_of_symbols,

                                             const double probability[])

    {

      if ((number_of_symbols < 2) || (number_of_symbols > (1 << 11)))

        AC_Error("invalid number of data symbols");

      if (data_symbols != number_of_symbols) {     // assign memory for data model

        data_symbols = number_of_symbols;

        last_symbol = data_symbols - 1;

        delete [] distribution;

                                         // define size of table for fast decoding

        if (data_symbols > 16) {

          unsigned table_bits = 3;

          while (data_symbols > (1U << (table_bits + 2))) ++table_bits;

          table_size  = 1 << table_bits;

          table_shift = DM__LengthShift - table_bits;

          distribution = new unsigned[data_symbols+table_size+2];

          decoder_table = distribution + data_symbols;

        }

        else {                                  // small alphabet: no table needed

          decoder_table = 0;

          table_size = table_shift = 0;

          distribution = new unsigned[data_symbols];

        }

      }

                                 // compute cumulative distribution, decoder table

      unsigned s = 0;

      double sum = 0.0, p = 1.0 / double(data_symbols);

      for (unsigned k = 0; k < data_symbols; k++) {

        if (probability) p = probability[k];

        if ((p < 0.0001) || (p > 0.9999)) AC_Error("invalid symbol probability");

        distribution[k] = unsigned(sum * (1 << DM__LengthShift));

        sum += p;

        if (table_size == 0) continue;

        unsigned w = distribution[k] >> table_shift;

        while (s < w) decoder_table[++s] = k - 1;

      }

      if (table_size != 0) {

        decoder_table[0] = 0;

        while (s <= table_size) decoder_table[++s] = data_symbols - 1;

      }

      if ((sum < 0.9999) || (sum > 1.0001)) AC_Error("invalid probabilities");

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    // - - Adaptive data model implementation  - - - - - - - - - - - - - - - - - -

    Adaptive_Data_Model::Adaptive_Data_Model(void)

    {

      data_symbols = 0;

      distribution = 0;

    }

    Adaptive_Data_Model::Adaptive_Data_Model(unsigned number_of_symbols)

    {

      data_symbols = 0;

      distribution = 0;

      set_alphabet(number_of_symbols);

    }

    Adaptive_Data_Model::~Adaptive_Data_Model(void)

    {

      delete [] distribution;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Adaptive_Data_Model::set_alphabet(unsigned number_of_symbols)

    {

      if ((number_of_symbols < 2) || (number_of_symbols > (1 << 11)))

        AC_Error("invalid number of data symbols");

      if (data_symbols != number_of_symbols) {     // assign memory for data model

        data_symbols = number_of_symbols;

        last_symbol = data_symbols - 1;

        delete [] distribution;

                                         // define size of table for fast decoding

        if (data_symbols > 16) {

          unsigned table_bits = 3;

          while (data_symbols > (1U << (table_bits + 2))) ++table_bits;

          table_size  = 1 << table_bits;

          table_shift = DM__LengthShift - table_bits;

          distribution = new unsigned[2*data_symbols+table_size+2];

          decoder_table = distribution + 2 * data_symbols;

        }

        else {                                  // small alphabet: no table needed

          decoder_table = 0;

          table_size = table_shift = 0;

          distribution = new unsigned[2*data_symbols];

        }

        symbol_count = distribution + data_symbols;

      }

      reset();                                                 // initialize model

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Adaptive_Data_Model::update(bool from_encoder)

    {

                                       // halve counts when a threshold is reached

      if ((total_count += update_cycle) > DM__MaxCount) {

        total_count = 0;

        for (unsigned n = 0; n < data_symbols; n++)

          total_count += (symbol_count[n] = (symbol_count[n] + 1) >> 1);

      }

      assert(total_count > 0);

                                 // compute cumulative distribution, decoder table

      unsigned k, sum = 0, s = 0;

      unsigned scale = 0x80000000U / total_count;

      if (from_encoder || (table_size == 0))

        for (k = 0; k < data_symbols; k++) {

          distribution[k] = (scale * sum) >> (31 - DM__LengthShift);

          sum += symbol_count[k];

        }

      else {

        assert(decoder_table);

        for (k = 0; k < data_symbols; k++) {

          distribution[k] = (scale * sum) >> (31 - DM__LengthShift);

          sum += symbol_count[k];

          unsigned w = distribution[k] >> table_shift;

          while (s < w) decoder_table[++s] = k - 1;

        }

        decoder_table[0] = 0;

        while (s <= table_size) decoder_table[++s] = data_symbols - 1;

      }

                                                 // set frequency of model updates

      update_cycle = (5 * update_cycle) >> 2;

      unsigned max_cycle = (data_symbols + 6) << 3;

      if (update_cycle > max_cycle) update_cycle = max_cycle;

      symbols_until_update = update_cycle;

    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    void Adaptive_Data_Model::reset(void)

    {

      if (data_symbols == 0) return;

                          // restore probability estimates to uniform distribution

      total_count = 0;

      update_cycle = data_symbols;

      for (unsigned k = 0; k < data_symbols; k++) symbol_count[k] = 1;

      update(false);

      symbols_until_update = update_cycle = (data_symbols + 6) >> 1;

    }

}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */