/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: 2026-01-10 06:17

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	3.83k	#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	1.07k	uint8_t meta, blosc2_dparams *dparams) {
515	1.07k	BLOSC_UNUSED_PARAM(meta);
516	1.07k	BLOSC_UNUSED_PARAM(dparams);
517	1.07k	BLOSC_ERROR_NULL(input, BLOSC2_ERROR_NULL_POINTER);
518	1.07k	BLOSC_ERROR_NULL(output, BLOSC2_ERROR_NULL_POINTER);
519
520	1.07k	uint8_t ip = (uint8_t ) input;
521	1.07k	uint8_t *ip_limit = ip + input_len;
522	1.07k	uint8_t op = (uint8_t ) output;
523	1.07k	uint8_t ndim;
524	1.07k	int32_t blockshape[2];
525	1.07k	int32_t eshape[2];
526	1.07k	uint8_t *buffercpy;
527	1.07k	uint8_t local_buffer[16];
528	1.07k	uint8_t token;
529	1.07k	if (NDLZ_UNEXPECT_CONDITIONAL(input_len < 8)) {
530	8	return 0;
531	8	}
532
533		/* we start with literal copy */
534	1.06k	ndim = *ip;
535	1.06k	ip++;
536	1.06k	if (ndim != 2) {
537	5	BLOSC_TRACE_ERROR("This codec only works for ndim = 2");
538	5	return BLOSC2_ERROR_FAILURE;
539	5	}
540	1.06k	memcpy(&blockshape[0], ip, 4);
541	1.06k	ip += 4;
542	1.06k	memcpy(&blockshape[1], ip, 4);
543	1.06k	ip += 4;
544
545		// Sanity check. See https://www.cve.org/CVERecord?id=CVE-2024-3204
546	1.06k	if (output_len < 0 \|\| blockshape[0] < 0 \|\| blockshape[1] < 0) {
547	65	BLOSC_TRACE_ERROR("Output length or blockshape is negative");
548	65	return BLOSC2_ERROR_FAILURE;
549	65	}
550
551	998	eshape[0] = ((blockshape[0] + 3) / 4) * 4;
552	998	eshape[1] = ((blockshape[1] + 3) / 4) * 4;
553
554	998	if (NDLZ_UNEXPECT_CONDITIONAL((int64_t)output_len < (int64_t)blockshape[0] * (int64_t)blockshape[1])) {
555	108	BLOSC_TRACE_ERROR("The blockshape is bigger than the output buffer");
556	108	return 0;
557	108	}
558	890	memset(op, 0, blockshape[0] * blockshape[1]);
559
560	890	uint32_t i_stop[2];
561	2.67k	for (int i = 0; i < 2; ++i) {
562	1.78k	i_stop[i] = eshape[i] / 4;
563	1.78k	}
564
565		/* main loop */
566	890	uint32_t ii[2];
567	890	uint32_t padding[2] = {0};
568	890	uint32_t ind = 0;
569	890	uint8_t cell_aux[16];
570	431M	for (ii[0] = 0; ii[0] < i_stop[0]; ++ii[0]) {
571	431M	for (ii[1] = 0; ii[1] < i_stop[1]; ++ii[1]) { // for each cell
572	1.76k	if (NDLZ_UNEXPECT_CONDITIONAL(ip > ip_limit)) {
573	16	BLOSC_TRACE_ERROR("Exceeding input length");
574	16	return BLOSC2_ERROR_FAILURE;
575	16	}
576	1.74k	if (ii[0] == i_stop[0] - 1) {
577	1.37k	padding[0] = (blockshape[0] % 4 == 0) ? 4 : blockshape[0] % 4;
578	1.37k	} else {
579	370	padding[0] = 4;
580	370	}
581	1.74k	if (ii[1] == i_stop[1] - 1) {
582	1.18k	padding[1] = (blockshape[1] % 4 == 0) ? 4 : blockshape[1] % 4;
583	1.18k	} else {
584	561	padding[1] = 4;
585	561	}
586	1.74k	token = *ip++;
587	1.74k	if (token == 0) { // no match
588	227	buffercpy = ip;
589	227	ip += padding[0] * padding[1];
590	1.51k	} else if (token == (uint8_t) ((1U << 7U) \| (1U << 6U))) { // cell match
591	71	uint16_t offset = ((uint16_t ) ip);
592	71	buffercpy = ip - offset - 1;
593	71	ip += 2;
594	1.44k	} else if (token == (uint8_t) (1U << 6U)) { // whole cell of same element
595	23	buffercpy = cell_aux;
596	23	memset(buffercpy, *ip, 16);
597	23	ip++;
598	1.42k	} else if (token >= 224) { // three rows match
599	389	buffercpy = local_buffer;
600	389	uint16_t offset = ((uint16_t ) ip);
601	389	offset += 3;
602	389	ip += 2;
603	389	int i, j, k;
604	389	if ((token >> 3U) == 28) {
605	72	i = 1;
606	72	j = 2;
607	72	k = 3;
608	317	} else {
609	317	i = 0;
610	317	if ((token >> 3U) < 30) {
611	76	j = 1;
612	76	k = 2;
613	241	} else {
614	241	k = 3;
615	241	if ((token >> 3U) == 30) {
616	83	j = 1;
617	158	} else {
618	158	j = 2;
619	158	}
620	241	}
621	317	}
622	389	memcpy(&buffercpy[i * 4], ip - offset, 4);
623	389	memcpy(&buffercpy[j * 4], ip - offset + 4, 4);
624	389	memcpy(&buffercpy[k * 4], ip - offset + 8, 4);
625	941	for (int l = 0; l < 4; l++) {
626	941	if ((l != i) && (l != j) && (l != k)) {
627	389	memcpy(&buffercpy[l * 4], ip, 4);
628	389	ip += 4;
629	389	break;
630	389	}
631	941	}
632
633	1.03k	} else if ((token >= 128) && (token <= 191)) { // rows pair match
634	766	buffercpy = local_buffer;
635	766	uint16_t offset = ((uint16_t ) ip);
636	766	offset += 3;
637	766	ip += 2;
638	766	int i, j;
639	766	if (token == 128) {
640	95	i = 2;
641	95	j = 3;
642	671	} else {
643	671	i = (token - 128) >> 5U;
644	671	j = ((token - 128) >> 3U) - (i << 2U);
645	671	}
646	766	memcpy(&buffercpy[i * 4], ip - offset, 4);
647	766	memcpy(&buffercpy[j * 4], ip - offset + 4, 4);
648	3.83k	for (int k = 0; k < 4; k++) {
649	3.06k	if ((k != i) && (k != j)) {
650	1.57k	memcpy(&buffercpy[k * 4], ip, 4);
651	1.57k	ip += 4;
652	1.57k	}
653	3.06k	}
654	766	} else if ((token >= 40) && (token <= 63)) { // 2 rows pair matches
655	245	buffercpy = local_buffer;
656	245	uint16_t offset_1 = ((uint16_t ) ip);
657	245	offset_1 += 5;
658	245	ip += 2;
659	245	uint16_t offset_2 = ((uint16_t ) ip);
660	245	offset_2 += 5;
661	245	ip += 2;
662	245	int i, j, k, l, m;
663	245	i = 0;
664	245	j = ((token - 32) >> 3U);
665	245	l = -1;
666	980	for (k = 1; k < 4; k++) {
667	735	if ((k != i) && (k != j)) {
668	490	if (l == -1) {
669	245	l = k;
670	245	} else {
671	245	m = k;
672	245	}
673	490	}
674	735	}
675	245	memcpy(&buffercpy[i * 4], ip - offset_1, 4);
676	245	memcpy(&buffercpy[j * 4], ip - offset_1 + 4, 4);
677	245	memcpy(&buffercpy[l * 4], ip - offset_2, 4);
678	245	memcpy(&buffercpy[m * 4], ip - offset_2 + 4, 4);
679
680	245	} else {
681	25	BLOSC_TRACE_ERROR("Invalid token: %u at cell [%d, %d]\n", token, ii[0], ii[1]);
682	25	return BLOSC2_ERROR_FAILURE;
683	25	}
684		// fill op with buffercpy
685	1.72k	uint32_t orig = ii[0] * 4 * blockshape[1] + ii[1] * 4;
686	8.60k	for (uint32_t i = 0; i < 4; i++) {
687	6.88k	if (i < padding[0]) {
688	3.61k	ind = orig + i * blockshape[1];
689	3.61k	memcpy(&op[ind], buffercpy, padding[1]);
690	3.61k	}
691	6.88k	buffercpy += padding[1];
692	6.88k	}
693	1.72k	if (ind > (uint32_t) output_len) {
694	0	BLOSC_TRACE_ERROR("Exceeding output size");
695	0	return BLOSC2_ERROR_FAILURE;
696	0	}
697	1.72k	}
698	431M	}
699	849	ind += padding[1];
700
701	849	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	849	if (ind > (uint32_t) output_len) {
706	0	BLOSC_TRACE_ERROR("Exceeding output size");
707	0	return BLOSC2_ERROR_FAILURE;
708	0	}
709
710	849	return (int) ind;
711	849	}