/src/c-blosc2/plugins/codecs/ndlz/ndlz4x4.c

Source
/*********************************************************************
  Blosc - Blocked Shuffling and Compression Library

  Copyright (c) 2021  Blosc Development Team <blosc@blosc.org>
  https://blosc.org
  License: BSD 3-Clause (see LICENSE.txt)

  See LICENSE.txt for details about copyright and rights to use.
**********************************************************************/

/*********************************************************************
  This codec is meant to leverage multidimensionality for getting
  better compression ratios.  The idea is to look for similarities
  in places that are closer in a euclidean metric, not the typical
  linear one.
**********************************************************************/

#include "ndlz4x4.h"
#include "xxhash.h"
#include "b2nd.h"

#include <stdlib.h>
#include <string.h>

/*
 * Give hints to the compiler for branch prediction optimization.
 */
#if defined(__GNUC__) && (__GNUC__ > 2)
#define NDLZ_EXPECT_CONDITIONAL(c)    (__builtin_expect((c), 1))
#define NDLZ_UNEXPECT_CONDITIONAL(c)  (__builtin_expect((c), 0))
#else
#define NDLZ_EXPECT_CONDITIONAL(c)    (c)
#define NDLZ_UNEXPECT_CONDITIONAL(c)  (c)
#endif

/*
 * Use inlined functions for supported systems.
 */
#if defined(_MSC_VER) && !defined(__cplusplus)   /* Visual Studio */
#define inline __inline  /* Visual C is not C99, but supports some kind of inline */
#endif

#define MAX_COPY 32U
#define MAX_DISTANCE 65535


#ifdef BLOSC_STRICT_ALIGN
#define NDLZ_READU16(p) ((p)[0] | (p)[1]<<8)
#define NDLZ_READU32(p) ((p)[0] | (p)[1]<<8 | (p)[2]<<16 | (p)[3]<<24)
#else
#define NDLZ_READU16(p) *((const uint16_t*)(p))
#define NDLZ_READU32(p) *((const uint32_t*)(p))
#endif

#define HASH_LOG (12)


int ndlz4_compress(const uint8_t *input, int32_t input_len, uint8_t *output, int32_t output_len,
                   uint8_t meta, blosc2_cparams *cparams) {
  BLOSC_UNUSED_PARAM(meta);
  BLOSC_ERROR_NULL(cparams, BLOSC2_ERROR_NULL_POINTER);
  BLOSC_ERROR_NULL(cparams->schunk, BLOSC2_ERROR_NULL_POINTER);
  uint8_t *smeta;
  int32_t smeta_len;

  if (blosc2_meta_get(cparams->schunk, "b2nd", &smeta, &smeta_len) < 0) {
    BLOSC_TRACE_ERROR("b2nd layer not found!");
    return BLOSC2_ERROR_FAILURE;
  }

  int8_t ndim;
  int64_t *shape = malloc(8 * sizeof(int64_t));
  int32_t *chunkshape = malloc(8 * sizeof(int32_t));
  int32_t *blockshape = malloc(8 * sizeof(int32_t));
  b2nd_deserialize_meta(smeta, smeta_len, &ndim, shape, chunkshape, blockshape, NULL, NULL);
  free(smeta);

  if (ndim != 2) {
    BLOSC_TRACE_ERROR("This codec only works for ndim = 2");
    return BLOSC2_ERROR_FAILURE;
  }

  if (input_len != (blockshape[0] * blockshape[1])) {
    BLOSC_TRACE_ERROR("Length not equal to blocksize");
    return BLOSC2_ERROR_FAILURE;
  }

  if (NDLZ_UNEXPECT_CONDITIONAL(output_len < (int) (1 + ndim * sizeof(int32_t)))) {
    BLOSC_TRACE_ERROR("Output too small");
    return BLOSC2_ERROR_FAILURE;
  }

  uint8_t *ip = (uint8_t *) input;
  uint8_t *op = (uint8_t *) output;
  uint8_t *op_limit;
  uint32_t hval, hash_cell;
  uint32_t hash_triple[2] = {0};
  uint32_t hash_pair[3] = {0};
  uint8_t bufarea[16];
  uint8_t *buf_cell = bufarea;
  uint8_t buf_triple[12];
  uint8_t buf_pair[8];
  uint8_t *buf_aux;
  uint32_t tab_cell[1U << 12U] = {0};
  uint32_t tab_triple[1U << 12U] = {0};
  uint32_t tab_pair[1U << 12U] = {0};
  uint32_t update_triple[2] = {0};
  uint32_t update_pair[3] = {0};

  // Minimum cratios before issuing and _early giveup_
  // Remind that ndlz is not meant for cratios <= 2 (too costly to decompress)

  op_limit = op + output_len;

  // Initialize the hash table to distances of 0
  for (unsigned i = 0; i < (1U << 12U); i++) {
    tab_cell[i] = 0;
  }

  /* input and output buffer cannot be less than 16 and 66 bytes or we can get into trouble */
  int overhead = 17 + (blockshape[0] * blockshape[1] / 16 - 1) * 2;
  if (input_len < 16 || output_len < overhead) {
    BLOSC_TRACE_ERROR("Incorrect length or maxout");
    return 0;
  }

  uint8_t *obase = op;

  /* we start with literal copy */
  *op++ = ndim;
  memcpy(op, &blockshape[0], 4);
  op += 4;
  memcpy(op, &blockshape[1], 4);
  op += 4;

  uint32_t i_stop[2];
  for (int i = 0; i < 2; ++i) {
    i_stop[i] = (blockshape[i] + 3) / 4;
  }

  /* main loop */
  uint32_t padding[2];
  uint32_t ii[2];
  for (ii[0] = 0; ii[0] < i_stop[0]; ++ii[0]) {
    for (ii[1] = 0; ii[1] < i_stop[1]; ++ii[1]) {      // for each cell
      uint8_t token;
      for (int h = 0; h < 2; h++) {         // new cell -> new possible references
        update_triple[h] = 0;
        update_pair[h] = 0;
      }
      update_pair[2] = 0;

      if (NDLZ_UNEXPECT_CONDITIONAL(op + 16 + 1 > op_limit)) {
        free(shape);
        free(chunkshape);
        free(blockshape);
        return 0;
      }

      uint32_t orig = ii[0] * 4 * blockshape[1] + ii[1] * 4;
      if (((blockshape[0] % 4 != 0) && (ii[0] == i_stop[0] - 1)) ||
          ((blockshape[1] % 4 != 0) && (ii[1] == i_stop[1] - 1))) {
        token = 0;                                   // padding -> literal copy
        *op++ = token;
        if (ii[0] == i_stop[0] - 1) {
          padding[0] = (blockshape[0] % 4 == 0) ? 4 : blockshape[0] % 4;
        } else {
          padding[0] = 4;
        }
        if (ii[1] == i_stop[1] - 1) {
          padding[1] = (blockshape[1] % 4 == 0) ? 4 : blockshape[1] % 4;
        } else {
          padding[1] = 4;
        }
        for (uint32_t i = 0; i < padding[0]; i++) {
          memcpy(op, &ip[orig + i * blockshape[1]], padding[1]);
          op += padding[1];
        }
      } else {
        for (uint64_t i = 0; i < 4; i++) {           // fill cell buffer
          uint64_t ind = orig + i * blockshape[1];
          memcpy(buf_cell, &ip[ind], 4);
          buf_cell += 4;
        }
        buf_cell -= 16;

        const uint8_t *ref;
        uint32_t distance;
        uint8_t *anchor = op;    /* comparison starting-point */

        /* find potential match */
        hash_cell = XXH32(buf_cell, 16, 1);        // calculate cell hash
        hash_cell >>= 32U - 12U;
        ref = obase + tab_cell[hash_cell];

        /* calculate distance to the match */
        if (tab_cell[hash_cell] == 0) {
          distance = 0;
        } else {
          bool same = true;
          buf_aux = obase + tab_cell[hash_cell];
          for (int i = 0; i < 16; i++) {
            if (buf_cell[i] != buf_aux[i]) {
              same = false;
              break;
            }
          }
          if (same) {
            distance = (int32_t) (anchor - ref);
          } else {
            distance = 0;
          }
        }

        bool alleq = true;
        for (int i = 1; i < 16; i++) {
          if (buf_cell[i] != buf_cell[0]) {
            alleq = false;
            break;
          }
        }
        if (alleq) {                              // all elements of the cell equal
          token = (uint8_t) (1U << 6U);
          *op++ = token;
          *op++ = buf_cell[0];

        } else if (distance == 0 || (distance >= MAX_DISTANCE)) {   // no cell match
          bool literal = true;

          // 2 rows pairs matches
          for (int j = 1; j < 4; j++) {
            memcpy(buf_pair, buf_cell, 4);
            memcpy(&buf_pair[4], &buf_cell[j * 4], 4);
            hval = XXH32(buf_pair, 8, 1);        // calculate rows pair hash
            hval >>= 32U - 12U;
            ref = obase + tab_pair[hval];
            /* calculate distance to the match */
            bool same = true;
            uint16_t offset;
            if (tab_pair[hval] != 0) {
              buf_aux = obase + tab_pair[hval];
              for (int k = 0; k < 8; k++) {
                if (buf_pair[k] != buf_aux[k]) {
                  same = false;
                  break;
                }
              }
              offset = (uint16_t) (anchor - obase - tab_pair[hval]);
            } else {
              same = false;
            }
            if (same) {
              distance = (int32_t) (anchor - ref);
            } else {
              distance = 0;
            }
            if ((distance != 0) && (distance < MAX_DISTANCE)) {     /* rows pair match */
              int k, m, l = -1;
              for (k = 1; k < 4; k++) {
                if (k != j) {
                  if (l == -1) {
                    l = k;
                  } else {
                    m = k;
                  }
                }
              }
              memcpy(buf_pair, &buf_cell[l * 4], 4);
              memcpy(&buf_pair[4], &buf_cell[m * 4], 4);
              hval = XXH32(buf_pair, 8, 1);        // calculate rows pair hash
              hval >>= 32U - 12U;
              ref = obase + tab_pair[hval];
              same = true;
              if (tab_pair[hval] != 0) {
                buf_aux = obase + tab_pair[hval];
                for (k = 0; k < 8; k++) {
                  if (buf_pair[k] != buf_aux[k]) {
                    same = false;
                    break;
                  }
                }
              } else {
                same = false;
              }
              if (same) {
                distance = (int32_t) (anchor + l * 4 - ref);
              } else {
                distance = 0;
              }
              if ((distance != 0) && (distance < MAX_DISTANCE)) {   /* 2 pair matches */
                literal = false;
                token = (uint8_t) ((1U << 5U) | (j << 3U));
                *op++ = token;
                uint16_t offset_2 = (uint16_t) (anchor - obase - tab_pair[hval]);
                *(uint16_t *) op = offset;
                op += sizeof(offset);
                *(uint16_t *) op = offset_2;
                op += sizeof(offset_2);
                goto match;
              }
            }
          }

          // rows triples
          for (int i = 0; i < 2; i++) {
            memcpy(buf_triple, &buf_cell[i * 4], 4);
            for (int j = i + 1; j < 3; j++) {
              memcpy(&buf_triple[4], &buf_cell[j * 4], 4);
              for (int k = j + 1; k < 4; k++) {
                memcpy(&buf_triple[8], &buf_cell[k * 4], 4);
                hval = XXH32(buf_triple, 12, 1);        // calculate triple hash
                hval >>= 32U - 12U;
                /* calculate distance to the match */
                bool same = true;
                uint16_t offset;
                if (tab_triple[hval] != 0) {
                  buf_aux = obase + tab_triple[hval];
                  for (int l = 0; l < 12; l++) {
                    if (buf_triple[l] != buf_aux[l]) {
                      same = false;
                      break;
                    }
                  }
                  offset = (uint16_t) (anchor - obase - tab_triple[hval]);
                } else {
                  same = false;
                  if ((j - i == 1) && (k - j == 1)) {
                    update_triple[i] = (uint32_t) (anchor + 1 + i * 4 - obase);     /* update hash table */
                    hash_triple[i] = hval;
                  }
                }
                ref = obase + tab_triple[hval];

                if (same) {
                  distance = (int32_t) (anchor + i * 4 - ref);
                } else {
                  distance = 0;
                }
                if ((distance != 0) && (distance < MAX_DISTANCE)) {
                  literal = false;
                  if (i == 1) {
                    token = (uint8_t) (7U << 5U);
                  } else {
                    token = (uint8_t) ((7U << 5U) | ((j + k - 2) << 3U));
                  }
                  *op++ = token;
                  memcpy(op, &offset, 2);
                  op += 2;
                  for (int l = 0; l < 4; l++) {
                    if ((l != i) && (l != j) && (l != k)) {
                      memcpy(op, &buf_cell[4 * l], 4);
                      op += 4;
                      goto match;
                    }
                  }
                }
              }
            }
          }

          // rows pairs
          for (int i = 0; i < 3; i++) {
            memcpy(buf_pair, &buf_cell[i * 4], 4);
            for (int j = i + 1; j < 4; j++) {
              memcpy(&buf_pair[4], &buf_cell[j * 4], 4);
              hval = XXH32(buf_pair, 8, 1);        // calculate rows pair hash
              hval >>= 32U - 12U;
              ref = obase + tab_pair[hval];
              /* calculate distance to the match */
              bool same = true;
              uint16_t offset;
              if (tab_pair[hval] != 0) {
                buf_aux = obase + tab_pair[hval];
                for (int k = 0; k < 8; k++) {
                  if (buf_pair[k] != buf_aux[k]) {
                    same = false;
                    break;
                  }
                }
                offset = (uint16_t) (anchor - obase - tab_pair[hval]);
              } else {
                same = false;
                if (j - i == 1) {
                  update_pair[i] = (uint32_t) (anchor + 1 + i * 4 - obase);     /* update hash table */
                  hash_pair[i] = hval;
                }
              }
              if (same) {
                distance = (int32_t) (anchor + i * 4 - ref);
              } else {
                distance = 0;
              }
              if ((distance != 0) && (distance < MAX_DISTANCE)) {     /* rows pair match */
                literal = false;
                if (i == 2) {
                  token = (uint8_t) (1U << 7U);
                } else {
                  token = (uint8_t) ((1U << 7U) | (i << 5U) | (j << 3U));
                }
                *op++ = token;
                memcpy(op, &offset, 2);
                op += 2;
                for (int k = 0; k < 4; k++) {
                  if ((k != i) && (k != j)) {
                    memcpy(op, &buf_cell[4 * k], 4);
                    op += 4;
                  }
                }
                goto match;
              }
            }
          }

          match:
          if (literal) {
            tab_cell[hash_cell] = (uint32_t) (anchor + 1 - obase);     /* update hash tables */
            if (update_triple[0] != 0) {
              for (int h = 0; h < 2; h++) {
                tab_triple[hash_triple[h]] = update_triple[h];
              }
            }
            if (update_pair[0] != 0) {
              for (int h = 0; h < 3; h++) {
                tab_pair[hash_pair[h]] = update_pair[h];
              }
            }
            token = 0;
            *op++ = token;
            memcpy(op, buf_cell, 16);
            op += 16;
          }

        } else {   // cell match
          token = (uint8_t) ((1U << 7U) | (1U << 6U));
          *op++ = token;
          uint16_t offset = (uint16_t) (anchor - obase - tab_cell[hash_cell]);
          memcpy(op, &offset, 2);
          op += 2;
        }

      }
      if ((op - obase) > input_len) {
        BLOSC_TRACE_ERROR("Compressed data is bigger than input!");
        return 0;
      }
    }
  }

  free(shape);
  free(chunkshape);
  free(blockshape);

  return (int) (op - obase);
}


// See https://habr.com/en/company/yandex/blog/457612/
#ifdef __AVX2__

#if defined(_MSC_VER)
#define ALIGNED_(x) __declspec(align(x))
#else
#if defined(__GNUC__)
#define ALIGNED_(x) __attribute__ ((aligned(x)))
#endif
#endif
#define ALIGNED_TYPE_(t, x) t ALIGNED_(x)

static unsigned char* copy_match_16(unsigned char *op, const unsigned char *match, int32_t len)
{
  size_t offset = op - match;
  while (len >= 16) {

    static const ALIGNED_TYPE_(uint8_t, 16) masks[] =
      {
                0,  1,  2,  1,  4,  1,  4,  2,  8,  7,  6,  5,  4,  3,  2,  1, // offset = 0, not used as mask, but for shift
                0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, // offset = 1
                0,  1,  0,  1,  0,  1,  0,  1,  0,  1,  0,  1,  0,  1,  0,  1,
                0,  1,  2,  0,  1,  2,  0,  1,  2,  0,  1,  2,  0,  1,  2,  0,
                0,  1,  2,  3,  0,  1,  2,  3,  0,  1,  2,  3,  0,  1,  2,  3,
                0,  1,  2,  3,  4,  0,  1,  2,  3,  4,  0,  1,  2,  3,  4,  0,
                0,  1,  2,  3,  4,  5,  0,  1,  2,  3,  4,  5,  0,  1,  2,  3,
                0,  1,  2,  3,  4,  5,  6,  0,  1,  2,  3,  4,  5,  6,  0,  1,
                0,  1,  2,  3,  4,  5,  6,  7,  0,  1,  2,  3,  4,  5,  6,  7,
                0,  1,  2,  3,  4,  5,  6,  7,  8,  0,  1,  2,  3,  4,  5,  6,
                0,  1,  2,  3,  4,  5,  6,  7,  8,  9,  0,  1,  2,  3,  4,  5,
                0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10,  0,  1,  2,  3,  4,
                0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11,  0,  1,  2,  3,
                0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12,  0,  1,  2,
                0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13,  0,  1,
                0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14,  0,
                0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14,  15, // offset = 16
      };

    _mm_storeu_si128((__m128i *)(op),
                     _mm_shuffle_epi8(_mm_loadu_si128((const __m128i *)(match)),
                                      _mm_load_si128((const __m128i *)(masks) + offset)));

    match += masks[offset];

    op += 16;
    len -= 16;
  }
  // Deal with remainders
  for (; len > 0; len--) {
    *op++ = *match++;
  }
  return op;
}
#endif


int ndlz4_decompress(const uint8_t *input, int32_t input_len, uint8_t *output, int32_t output_len,
                     uint8_t meta, blosc2_dparams *dparams) {
  BLOSC_UNUSED_PARAM(meta);
  BLOSC_UNUSED_PARAM(dparams);
  BLOSC_ERROR_NULL(input, BLOSC2_ERROR_NULL_POINTER);
  BLOSC_ERROR_NULL(output, BLOSC2_ERROR_NULL_POINTER);

  uint8_t *ip = (uint8_t *) input;
  uint8_t *ip_limit = ip + input_len;
  uint8_t *op = (uint8_t *) output;
  uint8_t ndim;
  int32_t blockshape[2];
  int32_t eshape[2];
  uint8_t *buffercpy;
  uint8_t local_buffer[16];
  uint8_t token;
  if (NDLZ_UNEXPECT_CONDITIONAL(input_len < 8)) {
    return 0;
  }

  /* we start with literal copy */
  ndim = *ip;
  ip++;
  if (ndim != 2) {
    BLOSC_TRACE_ERROR("This codec only works for ndim = 2");
    return BLOSC2_ERROR_FAILURE;
  }
  memcpy(&blockshape[0], ip, 4);
  ip += 4;
  memcpy(&blockshape[1], ip, 4);
  ip += 4;

  // Sanity check.  See https://www.cve.org/CVERecord?id=CVE-2024-3204
  if (output_len < 0 || blockshape[0] < 0 || blockshape[1] < 0) {
    BLOSC_TRACE_ERROR("Output length or blockshape is negative");
    return BLOSC2_ERROR_FAILURE;
  }

  eshape[0] = ((blockshape[0] + 3) / 4) * 4;
  eshape[1] = ((blockshape[1] + 3) / 4) * 4;

  if (NDLZ_UNEXPECT_CONDITIONAL((int64_t)output_len < (int64_t)blockshape[0] * (int64_t)blockshape[1])) {
    BLOSC_TRACE_ERROR("The blockshape is bigger than the output buffer");
    return 0;
  }
  memset(op, 0, blockshape[0] * blockshape[1]);

  uint32_t i_stop[2];
  for (int i = 0; i < 2; ++i) {
    i_stop[i] = eshape[i] / 4;
  }

  /* main loop */
  uint32_t ii[2];
  uint32_t padding[2] = {0};
  uint32_t ind = 0;
  uint8_t cell_aux[16];
  for (ii[0] = 0; ii[0] < i_stop[0]; ++ii[0]) {
    for (ii[1] = 0; ii[1] < i_stop[1]; ++ii[1]) {      // for each cell
      if (NDLZ_UNEXPECT_CONDITIONAL(ip > ip_limit)) {
        BLOSC_TRACE_ERROR("Exceeding input length");
        return BLOSC2_ERROR_FAILURE;
      }
      if (ii[0] == i_stop[0] - 1) {
        padding[0] = (blockshape[0] % 4 == 0) ? 4 : blockshape[0] % 4;
      } else {
        padding[0] = 4;
      }
      if (ii[1] == i_stop[1] - 1) {
        padding[1] = (blockshape[1] % 4 == 0) ? 4 : blockshape[1] % 4;
      } else {
        padding[1] = 4;
      }
      token = *ip++;
      if (token == 0) {    // no match
        buffercpy = ip;
        ip += padding[0] * padding[1];
      } else if (token == (uint8_t) ((1U << 7U) | (1U << 6U))) {  // cell match
        uint16_t offset = *((uint16_t *) ip);
        buffercpy = ip - offset - 1;
        ip += 2;
      } else if (token == (uint8_t) (1U << 6U)) { // whole cell of same element
        buffercpy = cell_aux;
        memset(buffercpy, *ip, 16);
        ip++;
      } else if (token >= 224) { // three rows match
        buffercpy = local_buffer;
        uint16_t offset = *((uint16_t *) ip);
        offset += 3;
        ip += 2;
        int i, j, k;
        if ((token >> 3U) == 28) {
          i = 1;
          j = 2;
          k = 3;
        } else {
          i = 0;
          if ((token >> 3U) < 30) {
            j = 1;
            k = 2;
          } else {
            k = 3;
            if ((token >> 3U) == 30) {
              j = 1;
            } else {
              j = 2;
            }
          }
        }
        memcpy(&buffercpy[i * 4], ip - offset, 4);
        memcpy(&buffercpy[j * 4], ip - offset + 4, 4);
        memcpy(&buffercpy[k * 4], ip - offset + 8, 4);
        for (int l = 0; l < 4; l++) {
          if ((l != i) && (l != j) && (l != k)) {
            memcpy(&buffercpy[l * 4], ip, 4);
            ip += 4;
            break;
          }
        }

      } else if ((token >= 128) && (token <= 191)) { // rows pair match
        buffercpy = local_buffer;
        uint16_t offset = *((uint16_t *) ip);
        offset += 3;
        ip += 2;
        int i, j;
        if (token == 128) {
          i = 2;
          j = 3;
        } else {
          i = (token - 128) >> 5U;
          j = ((token - 128) >> 3U) - (i << 2U);
        }
        memcpy(&buffercpy[i * 4], ip - offset, 4);
        memcpy(&buffercpy[j * 4], ip - offset + 4, 4);
        for (int k = 0; k < 4; k++) {
          if ((k != i) && (k != j)) {
            memcpy(&buffercpy[k * 4], ip, 4);
            ip += 4;
          }
        }
      } else if ((token >= 40) && (token <= 63)) {  // 2 rows pair matches
        buffercpy = local_buffer;
        uint16_t offset_1 = *((uint16_t *) ip);
        offset_1 += 5;
        ip += 2;
        uint16_t offset_2 = *((uint16_t *) ip);
        offset_2 += 5;
        ip += 2;
        int i, j, k, l, m;
        i = 0;
        j = ((token - 32) >> 3U);
        l = -1;
        for (k = 1; k < 4; k++) {
          if ((k != i) && (k != j)) {
            if (l == -1) {
              l = k;
            } else {
              m = k;
            }
          }
        }
        memcpy(&buffercpy[i * 4], ip - offset_1, 4);
        memcpy(&buffercpy[j * 4], ip - offset_1 + 4, 4);
        memcpy(&buffercpy[l * 4], ip - offset_2, 4);
        memcpy(&buffercpy[m * 4], ip - offset_2 + 4, 4);

      } else {
        BLOSC_TRACE_ERROR("Invalid token: %u at cell [%d, %d]\n", token, ii[0], ii[1]);
        return BLOSC2_ERROR_FAILURE;
      }
      // fill op with buffercpy
      uint32_t orig = ii[0] * 4 * blockshape[1] + ii[1] * 4;
      for (uint32_t i = 0; i < 4; i++) {
        if (i < padding[0]) {
          ind = orig + i * blockshape[1];
          memcpy(&op[ind], buffercpy, padding[1]);
        }
        buffercpy += padding[1];
      }
      if (ind > (uint32_t) output_len) {
        BLOSC_TRACE_ERROR("Exceeding output size");
        return BLOSC2_ERROR_FAILURE;
      }
    }
  }
  ind += padding[1];

  if ((int32_t)ind != (blockshape[0] * blockshape[1])) {
    BLOSC_TRACE_ERROR("Output size is not compatible with embedded blockshape");
    return BLOSC2_ERROR_FAILURE;
  }
  if (ind > (uint32_t) output_len) {
    BLOSC_TRACE_ERROR("Exceeding output size");
    return BLOSC2_ERROR_FAILURE;
  }

  return (int) ind;
}

Coverage Report

Created: 2025-12-11 06:39

Line	Count	Source
1		/*********************************************************************
2		Blosc - Blocked Shuffling and Compression Library
3
4		Copyright (c) 2021 Blosc Development Team <blosc@blosc.org>
5		https://blosc.org
6		License: BSD 3-Clause (see LICENSE.txt)
7
8		See LICENSE.txt for details about copyright and rights to use.
9		**********************************************************************/
10
11		/*********************************************************************
12		This codec is meant to leverage multidimensionality for getting
13		better compression ratios. The idea is to look for similarities
14		in places that are closer in a euclidean metric, not the typical
15		linear one.
16		**********************************************************************/
17
18		#include "ndlz4x4.h"
19		#include "xxhash.h"
20		#include "b2nd.h"
21
22		#include <stdlib.h>
23		#include <string.h>
24
25		/*
26		* Give hints to the compiler for branch prediction optimization.
27		*/
28		#if defined(__GNUC__) && (__GNUC__ > 2)
29		#define NDLZ_EXPECT_CONDITIONAL(c) (__builtin_expect((c), 1))
30	0	#define NDLZ_UNEXPECT_CONDITIONAL(c) (__builtin_expect((c), 0))
31		#else
32		#define NDLZ_EXPECT_CONDITIONAL(c) (c)
33		#define NDLZ_UNEXPECT_CONDITIONAL(c) (c)
34		#endif
35
36		/*
37		* Use inlined functions for supported systems.
38		*/
39		#if defined(_MSC_VER) && !defined(__cplusplus) /* Visual Studio */
40		#define inline __inline /* Visual C is not C99, but supports some kind of inline */
41		#endif
42
43		#define MAX_COPY 32U
44	0	#define MAX_DISTANCE 65535
45
46
47		#ifdef BLOSC_STRICT_ALIGN
48		#define NDLZ_READU16(p) ((p)[0] \| (p)[1]<<8)
49		#define NDLZ_READU32(p) ((p)[0] \| (p)[1]<<8 \| (p)[2]<<16 \| (p)[3]<<24)
50		#else
51		#define NDLZ_READU16(p) ((const uint16_t)(p))
52		#define NDLZ_READU32(p) ((const uint32_t)(p))
53		#endif
54
55		#define HASH_LOG (12)
56
57
58		int ndlz4_compress(const uint8_t input, int32_t input_len, uint8_t output, int32_t output_len,
59	0	uint8_t meta, blosc2_cparams *cparams) {
60	0	BLOSC_UNUSED_PARAM(meta);
61	0	BLOSC_ERROR_NULL(cparams, BLOSC2_ERROR_NULL_POINTER);
62	0	BLOSC_ERROR_NULL(cparams->schunk, BLOSC2_ERROR_NULL_POINTER);
63	0	uint8_t *smeta;
64	0	int32_t smeta_len;
65
66	0	if (blosc2_meta_get(cparams->schunk, "b2nd", &smeta, &smeta_len) < 0) {
67	0	BLOSC_TRACE_ERROR("b2nd layer not found!");
68	0	return BLOSC2_ERROR_FAILURE;
69	0	}
70
71	0	int8_t ndim;
72	0	int64_t shape = malloc(8 sizeof(int64_t));
73	0	int32_t chunkshape = malloc(8 sizeof(int32_t));
74	0	int32_t blockshape = malloc(8 sizeof(int32_t));
75	0	b2nd_deserialize_meta(smeta, smeta_len, &ndim, shape, chunkshape, blockshape, NULL, NULL);
76	0	free(smeta);
77
78	0	if (ndim != 2) {
79	0	BLOSC_TRACE_ERROR("This codec only works for ndim = 2");
80	0	return BLOSC2_ERROR_FAILURE;
81	0	}
82
83	0	if (input_len != (blockshape[0] * blockshape[1])) {
84	0	BLOSC_TRACE_ERROR("Length not equal to blocksize");
85	0	return BLOSC2_ERROR_FAILURE;
86	0	}
87
88	0	if (NDLZ_UNEXPECT_CONDITIONAL(output_len < (int) (1 + ndim * sizeof(int32_t)))) {
89	0	BLOSC_TRACE_ERROR("Output too small");
90	0	return BLOSC2_ERROR_FAILURE;
91	0	}
92
93	0	uint8_t ip = (uint8_t ) input;
94	0	uint8_t op = (uint8_t ) output;
95	0	uint8_t *op_limit;
96	0	uint32_t hval, hash_cell;
97	0	uint32_t hash_triple[2] = {0};
98	0	uint32_t hash_pair[3] = {0};
99	0	uint8_t bufarea[16];
100	0	uint8_t *buf_cell = bufarea;
101	0	uint8_t buf_triple[12];
102	0	uint8_t buf_pair[8];
103	0	uint8_t *buf_aux;
104	0	uint32_t tab_cell[1U << 12U] = {0};
105	0	uint32_t tab_triple[1U << 12U] = {0};
106	0	uint32_t tab_pair[1U << 12U] = {0};
107	0	uint32_t update_triple[2] = {0};
108	0	uint32_t update_pair[3] = {0};
109
110		// Minimum cratios before issuing and _early giveup_
111		// Remind that ndlz is not meant for cratios <= 2 (too costly to decompress)
112
113	0	op_limit = op + output_len;
114
115		// Initialize the hash table to distances of 0
116	0	for (unsigned i = 0; i < (1U << 12U); i++) {
117	0	tab_cell[i] = 0;
118	0	}
119
120		/* input and output buffer cannot be less than 16 and 66 bytes or we can get into trouble */
121	0	int overhead = 17 + (blockshape[0] * blockshape[1] / 16 - 1) * 2;
122	0	if (input_len < 16 \|\| output_len < overhead) {
123	0	BLOSC_TRACE_ERROR("Incorrect length or maxout");
124	0	return 0;
125	0	}
126
127	0	uint8_t *obase = op;
128
129		/* we start with literal copy */
130	0	*op++ = ndim;
131	0	memcpy(op, &blockshape[0], 4);
132	0	op += 4;
133	0	memcpy(op, &blockshape[1], 4);
134	0	op += 4;
135
136	0	uint32_t i_stop[2];
137	0	for (int i = 0; i < 2; ++i) {
138	0	i_stop[i] = (blockshape[i] + 3) / 4;
139	0	}
140
141		/* main loop */
142	0	uint32_t padding[2];
143	0	uint32_t ii[2];
144	0	for (ii[0] = 0; ii[0] < i_stop[0]; ++ii[0]) {
145	0	for (ii[1] = 0; ii[1] < i_stop[1]; ++ii[1]) { // for each cell
146	0	uint8_t token;
147	0	for (int h = 0; h < 2; h++) { // new cell -> new possible references
148	0	update_triple[h] = 0;
149	0	update_pair[h] = 0;
150	0	}
151	0	update_pair[2] = 0;
152
153	0	if (NDLZ_UNEXPECT_CONDITIONAL(op + 16 + 1 > op_limit)) {
154	0	free(shape);
155	0	free(chunkshape);
156	0	free(blockshape);
157	0	return 0;
158	0	}
159
160	0	uint32_t orig = ii[0] * 4 * blockshape[1] + ii[1] * 4;
161	0	if (((blockshape[0] % 4 != 0) && (ii[0] == i_stop[0] - 1)) \|\|
162	0	((blockshape[1] % 4 != 0) && (ii[1] == i_stop[1] - 1))) {
163	0	token = 0; // padding -> literal copy
164	0	*op++ = token;
165	0	if (ii[0] == i_stop[0] - 1) {
166	0	padding[0] = (blockshape[0] % 4 == 0) ? 4 : blockshape[0] % 4;
167	0	} else {
168	0	padding[0] = 4;
169	0	}
170	0	if (ii[1] == i_stop[1] - 1) {
171	0	padding[1] = (blockshape[1] % 4 == 0) ? 4 : blockshape[1] % 4;
172	0	} else {
173	0	padding[1] = 4;
174	0	}
175	0	for (uint32_t i = 0; i < padding[0]; i++) {
176	0	memcpy(op, &ip[orig + i * blockshape[1]], padding[1]);
177	0	op += padding[1];
178	0	}
179	0	} else {
180	0	for (uint64_t i = 0; i < 4; i++) { // fill cell buffer
181	0	uint64_t ind = orig + i * blockshape[1];
182	0	memcpy(buf_cell, &ip[ind], 4);
183	0	buf_cell += 4;
184	0	}
185	0	buf_cell -= 16;
186
187	0	const uint8_t *ref;
188	0	uint32_t distance;
189	0	uint8_t anchor = op; / comparison starting-point */
190
191		/* find potential match */
192	0	hash_cell = XXH32(buf_cell, 16, 1); // calculate cell hash
193	0	hash_cell >>= 32U - 12U;
194	0	ref = obase + tab_cell[hash_cell];
195
196		/* calculate distance to the match */
197	0	if (tab_cell[hash_cell] == 0) {
198	0	distance = 0;
199	0	} else {
200	0	bool same = true;
201	0	buf_aux = obase + tab_cell[hash_cell];
202	0	for (int i = 0; i < 16; i++) {
203	0	if (buf_cell[i] != buf_aux[i]) {
204	0	same = false;
205	0	break;
206	0	}
207	0	}
208	0	if (same) {
209	0	distance = (int32_t) (anchor - ref);
210	0	} else {
211	0	distance = 0;
212	0	}
213	0	}
214
215	0	bool alleq = true;
216	0	for (int i = 1; i < 16; i++) {
217	0	if (buf_cell[i] != buf_cell[0]) {
218	0	alleq = false;
219	0	break;
220	0	}
221	0	}
222	0	if (alleq) { // all elements of the cell equal
223	0	token = (uint8_t) (1U << 6U);
224	0	*op++ = token;
225	0	*op++ = buf_cell[0];
226
227	0	} else if (distance == 0 \|\| (distance >= MAX_DISTANCE)) { // no cell match
228	0	bool literal = true;
229
230		// 2 rows pairs matches
231	0	for (int j = 1; j < 4; j++) {
232	0	memcpy(buf_pair, buf_cell, 4);
233	0	memcpy(&buf_pair[4], &buf_cell[j * 4], 4);
234	0	hval = XXH32(buf_pair, 8, 1); // calculate rows pair hash
235	0	hval >>= 32U - 12U;
236	0	ref = obase + tab_pair[hval];
237		/* calculate distance to the match */
238	0	bool same = true;
239	0	uint16_t offset;
240	0	if (tab_pair[hval] != 0) {
241	0	buf_aux = obase + tab_pair[hval];
242	0	for (int k = 0; k < 8; k++) {
243	0	if (buf_pair[k] != buf_aux[k]) {
244	0	same = false;
245	0	break;
246	0	}
247	0	}
248	0	offset = (uint16_t) (anchor - obase - tab_pair[hval]);
249	0	} else {
250	0	same = false;
251	0	}
252	0	if (same) {
253	0	distance = (int32_t) (anchor - ref);
254	0	} else {
255	0	distance = 0;
256	0	}
257	0	if ((distance != 0) && (distance < MAX_DISTANCE)) { /* rows pair match */
258	0	int k, m, l = -1;
259	0	for (k = 1; k < 4; k++) {
260	0	if (k != j) {
261	0	if (l == -1) {
262	0	l = k;
263	0	} else {
264	0	m = k;
265	0	}
266	0	}
267	0	}
268	0	memcpy(buf_pair, &buf_cell[l * 4], 4);
269	0	memcpy(&buf_pair[4], &buf_cell[m * 4], 4);
270	0	hval = XXH32(buf_pair, 8, 1); // calculate rows pair hash
271	0	hval >>= 32U - 12U;
272	0	ref = obase + tab_pair[hval];
273	0	same = true;
274	0	if (tab_pair[hval] != 0) {
275	0	buf_aux = obase + tab_pair[hval];
276	0	for (k = 0; k < 8; k++) {
277	0	if (buf_pair[k] != buf_aux[k]) {
278	0	same = false;
279	0	break;
280	0	}
281	0	}
282	0	} else {
283	0	same = false;
284	0	}
285	0	if (same) {
286	0	distance = (int32_t) (anchor + l * 4 - ref);
287	0	} else {
288	0	distance = 0;
289	0	}
290	0	if ((distance != 0) && (distance < MAX_DISTANCE)) { /* 2 pair matches */
291	0	literal = false;
292	0	token = (uint8_t) ((1U << 5U) \| (j << 3U));
293	0	*op++ = token;
294	0	uint16_t offset_2 = (uint16_t) (anchor - obase - tab_pair[hval]);
295	0	(uint16_t ) op = offset;
296	0	op += sizeof(offset);
297	0	(uint16_t ) op = offset_2;
298	0	op += sizeof(offset_2);
299	0	goto match;
300	0	}
301	0	}
302	0	}
303
304		// rows triples
305	0	for (int i = 0; i < 2; i++) {
306	0	memcpy(buf_triple, &buf_cell[i * 4], 4);
307	0	for (int j = i + 1; j < 3; j++) {
308	0	memcpy(&buf_triple[4], &buf_cell[j * 4], 4);
309	0	for (int k = j + 1; k < 4; k++) {
310	0	memcpy(&buf_triple[8], &buf_cell[k * 4], 4);
311	0	hval = XXH32(buf_triple, 12, 1); // calculate triple hash
312	0	hval >>= 32U - 12U;
313		/* calculate distance to the match */
314	0	bool same = true;
315	0	uint16_t offset;
316	0	if (tab_triple[hval] != 0) {
317	0	buf_aux = obase + tab_triple[hval];
318	0	for (int l = 0; l < 12; l++) {
319	0	if (buf_triple[l] != buf_aux[l]) {
320	0	same = false;
321	0	break;
322	0	}
323	0	}
324	0	offset = (uint16_t) (anchor - obase - tab_triple[hval]);
325	0	} else {
326	0	same = false;
327	0	if ((j - i == 1) && (k - j == 1)) {
328	0	update_triple[i] = (uint32_t) (anchor + 1 + i * 4 - obase); /* update hash table */
329	0	hash_triple[i] = hval;
330	0	}
331	0	}
332	0	ref = obase + tab_triple[hval];
333
334	0	if (same) {
335	0	distance = (int32_t) (anchor + i * 4 - ref);
336	0	} else {
337	0	distance = 0;
338	0	}
339	0	if ((distance != 0) && (distance < MAX_DISTANCE)) {
340	0	literal = false;
341	0	if (i == 1) {
342	0	token = (uint8_t) (7U << 5U);
343	0	} else {
344	0	token = (uint8_t) ((7U << 5U) \| ((j + k - 2) << 3U));
345	0	}
346	0	*op++ = token;
347	0	memcpy(op, &offset, 2);
348	0	op += 2;
349	0	for (int l = 0; l < 4; l++) {
350	0	if ((l != i) && (l != j) && (l != k)) {
351	0	memcpy(op, &buf_cell[4 * l], 4);
352	0	op += 4;
353	0	goto match;
354	0	}
355	0	}
356	0	}
357	0	}
358	0	}
359	0	}
360
361		// rows pairs
362	0	for (int i = 0; i < 3; i++) {
363	0	memcpy(buf_pair, &buf_cell[i * 4], 4);
364	0	for (int j = i + 1; j < 4; j++) {
365	0	memcpy(&buf_pair[4], &buf_cell[j * 4], 4);
366	0	hval = XXH32(buf_pair, 8, 1); // calculate rows pair hash
367	0	hval >>= 32U - 12U;
368	0	ref = obase + tab_pair[hval];
369		/* calculate distance to the match */
370	0	bool same = true;
371	0	uint16_t offset;
372	0	if (tab_pair[hval] != 0) {
373	0	buf_aux = obase + tab_pair[hval];
374	0	for (int k = 0; k < 8; k++) {
375	0	if (buf_pair[k] != buf_aux[k]) {
376	0	same = false;
377	0	break;
378	0	}
379	0	}
380	0	offset = (uint16_t) (anchor - obase - tab_pair[hval]);
381	0	} else {
382	0	same = false;
383	0	if (j - i == 1) {
384	0	update_pair[i] = (uint32_t) (anchor + 1 + i * 4 - obase); /* update hash table */
385	0	hash_pair[i] = hval;
386	0	}
387	0	}
388	0	if (same) {
389	0	distance = (int32_t) (anchor + i * 4 - ref);
390	0	} else {
391	0	distance = 0;
392	0	}
393	0	if ((distance != 0) && (distance < MAX_DISTANCE)) { /* rows pair match */
394	0	literal = false;
395	0	if (i == 2) {
396	0	token = (uint8_t) (1U << 7U);
397	0	} else {
398	0	token = (uint8_t) ((1U << 7U) \| (i << 5U) \| (j << 3U));
399	0	}
400	0	*op++ = token;
401	0	memcpy(op, &offset, 2);
402	0	op += 2;
403	0	for (int k = 0; k < 4; k++) {
404	0	if ((k != i) && (k != j)) {
405	0	memcpy(op, &buf_cell[4 * k], 4);
406	0	op += 4;
407	0	}
408	0	}
409	0	goto match;
410	0	}
411	0	}
412	0	}
413
414	0	match:
415	0	if (literal) {
416	0	tab_cell[hash_cell] = (uint32_t) (anchor + 1 - obase); /* update hash tables */
417	0	if (update_triple[0] != 0) {
418	0	for (int h = 0; h < 2; h++) {
419	0	tab_triple[hash_triple[h]] = update_triple[h];
420	0	}
421	0	}
422	0	if (update_pair[0] != 0) {
423	0	for (int h = 0; h < 3; h++) {
424	0	tab_pair[hash_pair[h]] = update_pair[h];
425	0	}
426	0	}
427	0	token = 0;
428	0	*op++ = token;
429	0	memcpy(op, buf_cell, 16);
430	0	op += 16;
431	0	}
432
433	0	} else { // cell match
434	0	token = (uint8_t) ((1U << 7U) \| (1U << 6U));
435	0	*op++ = token;
436	0	uint16_t offset = (uint16_t) (anchor - obase - tab_cell[hash_cell]);
437	0	memcpy(op, &offset, 2);
438	0	op += 2;
439	0	}
440
441	0	}
442	0	if ((op - obase) > input_len) {
443	0	BLOSC_TRACE_ERROR("Compressed data is bigger than input!");
444	0	return 0;
445	0	}
446	0	}
447	0	}
448
449	0	free(shape);
450	0	free(chunkshape);
451	0	free(blockshape);
452
453	0	return (int) (op - obase);
454	0	}
455
456
457		// See https://habr.com/en/company/yandex/blog/457612/
458		#ifdef __AVX2__
459
460		#if defined(_MSC_VER)
461		#define ALIGNED_(x) __declspec(align(x))
462		#else
463		#if defined(__GNUC__)
464		#define ALIGNED_(x) __attribute__ ((aligned(x)))
465		#endif
466		#endif
467		#define ALIGNED_TYPE_(t, x) t ALIGNED_(x)
468
469		static unsigned char* copy_match_16(unsigned char op, const unsigned char match, int32_t len)
470		{
471		size_t offset = op - match;
472		while (len >= 16) {
473
474		static const ALIGNED_TYPE_(uint8_t, 16) masks[] =
475		{
476		0, 1, 2, 1, 4, 1, 4, 2, 8, 7, 6, 5, 4, 3, 2, 1, // offset = 0, not used as mask, but for shift
477		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // offset = 1
478		0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
479		0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0,
480		0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3,
481		0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0,
482		0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3,
483		0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1,
484		0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
485		0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6,
486		0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5,
487		0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4,
488		0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 1, 2, 3,
489		0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 0, 1, 2,
490		0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 0, 1,
491		0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 0,
492		0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, // offset = 16
493		};
494
495		_mm_storeu_si128((__m128i *)(op),
496		_mm_shuffle_epi8(_mm_loadu_si128((const __m128i *)(match)),
497		_mm_load_si128((const __m128i *)(masks) + offset)));
498
499		match += masks[offset];
500
501		op += 16;
502		len -= 16;
503		}
504		// Deal with remainders
505		for (; len > 0; len--) {
506		op++ = match++;
507		}
508		return op;
509		}
510		#endif
511
512
513		int ndlz4_decompress(const uint8_t input, int32_t input_len, uint8_t output, int32_t output_len,
514	0	uint8_t meta, blosc2_dparams *dparams) {
515	0	BLOSC_UNUSED_PARAM(meta);
516	0	BLOSC_UNUSED_PARAM(dparams);
517	0	BLOSC_ERROR_NULL(input, BLOSC2_ERROR_NULL_POINTER);
518	0	BLOSC_ERROR_NULL(output, BLOSC2_ERROR_NULL_POINTER);
519
520	0	uint8_t ip = (uint8_t ) input;
521	0	uint8_t *ip_limit = ip + input_len;
522	0	uint8_t op = (uint8_t ) output;
523	0	uint8_t ndim;
524	0	int32_t blockshape[2];
525	0	int32_t eshape[2];
526	0	uint8_t *buffercpy;
527	0	uint8_t local_buffer[16];
528	0	uint8_t token;
529	0	if (NDLZ_UNEXPECT_CONDITIONAL(input_len < 8)) {
530	0	return 0;
531	0	}
532
533		/* we start with literal copy */
534	0	ndim = *ip;
535	0	ip++;
536	0	if (ndim != 2) {
537	0	BLOSC_TRACE_ERROR("This codec only works for ndim = 2");
538	0	return BLOSC2_ERROR_FAILURE;
539	0	}
540	0	memcpy(&blockshape[0], ip, 4);
541	0	ip += 4;
542	0	memcpy(&blockshape[1], ip, 4);
543	0	ip += 4;
544
545		// Sanity check. See https://www.cve.org/CVERecord?id=CVE-2024-3204
546	0	if (output_len < 0 \|\| blockshape[0] < 0 \|\| blockshape[1] < 0) {
547	0	BLOSC_TRACE_ERROR("Output length or blockshape is negative");
548	0	return BLOSC2_ERROR_FAILURE;
549	0	}
550
551	0	eshape[0] = ((blockshape[0] + 3) / 4) * 4;
552	0	eshape[1] = ((blockshape[1] + 3) / 4) * 4;
553
554	0	if (NDLZ_UNEXPECT_CONDITIONAL((int64_t)output_len < (int64_t)blockshape[0] * (int64_t)blockshape[1])) {
555	0	BLOSC_TRACE_ERROR("The blockshape is bigger than the output buffer");
556	0	return 0;
557	0	}
558	0	memset(op, 0, blockshape[0] * blockshape[1]);
559
560	0	uint32_t i_stop[2];
561	0	for (int i = 0; i < 2; ++i) {
562	0	i_stop[i] = eshape[i] / 4;
563	0	}
564
565		/* main loop */
566	0	uint32_t ii[2];
567	0	uint32_t padding[2] = {0};
568	0	uint32_t ind = 0;
569	0	uint8_t cell_aux[16];
570	0	for (ii[0] = 0; ii[0] < i_stop[0]; ++ii[0]) {
571	0	for (ii[1] = 0; ii[1] < i_stop[1]; ++ii[1]) { // for each cell
572	0	if (NDLZ_UNEXPECT_CONDITIONAL(ip > ip_limit)) {
573	0	BLOSC_TRACE_ERROR("Exceeding input length");
574	0	return BLOSC2_ERROR_FAILURE;
575	0	}
576	0	if (ii[0] == i_stop[0] - 1) {
577	0	padding[0] = (blockshape[0] % 4 == 0) ? 4 : blockshape[0] % 4;
578	0	} else {
579	0	padding[0] = 4;
580	0	}
581	0	if (ii[1] == i_stop[1] - 1) {
582	0	padding[1] = (blockshape[1] % 4 == 0) ? 4 : blockshape[1] % 4;
583	0	} else {
584	0	padding[1] = 4;
585	0	}
586	0	token = *ip++;
587	0	if (token == 0) { // no match
588	0	buffercpy = ip;
589	0	ip += padding[0] * padding[1];
590	0	} else if (token == (uint8_t) ((1U << 7U) \| (1U << 6U))) { // cell match
591	0	uint16_t offset = ((uint16_t ) ip);
592	0	buffercpy = ip - offset - 1;
593	0	ip += 2;
594	0	} else if (token == (uint8_t) (1U << 6U)) { // whole cell of same element
595	0	buffercpy = cell_aux;
596	0	memset(buffercpy, *ip, 16);
597	0	ip++;
598	0	} else if (token >= 224) { // three rows match
599	0	buffercpy = local_buffer;
600	0	uint16_t offset = ((uint16_t ) ip);
601	0	offset += 3;
602	0	ip += 2;
603	0	int i, j, k;
604	0	if ((token >> 3U) == 28) {
605	0	i = 1;
606	0	j = 2;
607	0	k = 3;
608	0	} else {
609	0	i = 0;
610	0	if ((token >> 3U) < 30) {
611	0	j = 1;
612	0	k = 2;
613	0	} else {
614	0	k = 3;
615	0	if ((token >> 3U) == 30) {
616	0	j = 1;
617	0	} else {
618	0	j = 2;
619	0	}
620	0	}
621	0	}
622	0	memcpy(&buffercpy[i * 4], ip - offset, 4);
623	0	memcpy(&buffercpy[j * 4], ip - offset + 4, 4);
624	0	memcpy(&buffercpy[k * 4], ip - offset + 8, 4);
625	0	for (int l = 0; l < 4; l++) {
626	0	if ((l != i) && (l != j) && (l != k)) {
627	0	memcpy(&buffercpy[l * 4], ip, 4);
628	0	ip += 4;
629	0	break;
630	0	}
631	0	}
632
633	0	} else if ((token >= 128) && (token <= 191)) { // rows pair match
634	0	buffercpy = local_buffer;
635	0	uint16_t offset = ((uint16_t ) ip);
636	0	offset += 3;
637	0	ip += 2;
638	0	int i, j;
639	0	if (token == 128) {
640	0	i = 2;
641	0	j = 3;
642	0	} else {
643	0	i = (token - 128) >> 5U;
644	0	j = ((token - 128) >> 3U) - (i << 2U);
645	0	}
646	0	memcpy(&buffercpy[i * 4], ip - offset, 4);
647	0	memcpy(&buffercpy[j * 4], ip - offset + 4, 4);
648	0	for (int k = 0; k < 4; k++) {
649	0	if ((k != i) && (k != j)) {
650	0	memcpy(&buffercpy[k * 4], ip, 4);
651	0	ip += 4;
652	0	}
653	0	}
654	0	} else if ((token >= 40) && (token <= 63)) { // 2 rows pair matches
655	0	buffercpy = local_buffer;
656	0	uint16_t offset_1 = ((uint16_t ) ip);
657	0	offset_1 += 5;
658	0	ip += 2;
659	0	uint16_t offset_2 = ((uint16_t ) ip);
660	0	offset_2 += 5;
661	0	ip += 2;
662	0	int i, j, k, l, m;
663	0	i = 0;
664	0	j = ((token - 32) >> 3U);
665	0	l = -1;
666	0	for (k = 1; k < 4; k++) {
667	0	if ((k != i) && (k != j)) {
668	0	if (l == -1) {
669	0	l = k;
670	0	} else {
671	0	m = k;
672	0	}
673	0	}
674	0	}
675	0	memcpy(&buffercpy[i * 4], ip - offset_1, 4);
676	0	memcpy(&buffercpy[j * 4], ip - offset_1 + 4, 4);
677	0	memcpy(&buffercpy[l * 4], ip - offset_2, 4);
678	0	memcpy(&buffercpy[m * 4], ip - offset_2 + 4, 4);
679
680	0	} else {
681	0	BLOSC_TRACE_ERROR("Invalid token: %u at cell [%d, %d]\n", token, ii[0], ii[1]);
682	0	return BLOSC2_ERROR_FAILURE;
683	0	}
684		// fill op with buffercpy
685	0	uint32_t orig = ii[0] * 4 * blockshape[1] + ii[1] * 4;
686	0	for (uint32_t i = 0; i < 4; i++) {
687	0	if (i < padding[0]) {
688	0	ind = orig + i * blockshape[1];
689	0	memcpy(&op[ind], buffercpy, padding[1]);
690	0	}
691	0	buffercpy += padding[1];
692	0	}
693	0	if (ind > (uint32_t) output_len) {
694	0	BLOSC_TRACE_ERROR("Exceeding output size");
695	0	return BLOSC2_ERROR_FAILURE;
696	0	}
697	0	}
698	0	}
699	0	ind += padding[1];
700
701	0	if ((int32_t)ind != (blockshape[0] * blockshape[1])) {
702	0	BLOSC_TRACE_ERROR("Output size is not compatible with embedded blockshape");
703	0	return BLOSC2_ERROR_FAILURE;
704	0	}
705	0	if (ind > (uint32_t) output_len) {
706	0	BLOSC_TRACE_ERROR("Exceeding output size");
707	0	return BLOSC2_ERROR_FAILURE;
708	0	}
709
710	0	return (int) ind;
711	0	}