/* Lziprecover - Data recovery tool
   Copyright (C) 2009-2026 Antonio Diaz Diaz.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#define _FILE_OFFSET_BITS 64

#include <algorithm>
#include <cerrno>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <string>
#include <vector>
#include <stdint.h>
#include <unistd.h>
#include <sys/mman.h>
#include <sys/stat.h>

#include "lzip.h"
#include "lzip_index.h"


namespace {

const char * format_num( unsigned long long num )
  {
  enum { buffers = 8, bufsize = 32 };
  static char buffer[buffers][bufsize];	// circle of buffers for printf
  static int current = 0;

  unsigned long long den = 1;
  const unsigned factor = 1000;
  char * const buf = buffer[current++]; current %= buffers;
  const char * const prefix = "kMGTPEZYRQ";
  char p[2] = { 0, 0 };

  for( int i = 0; num / den >= factor && den * factor > den && prefix[i]; ++i )
    { den *= factor; *p = prefix[i]; }
  if( num % den == 0 )
    snprintf( buf, bufsize, "%llu %s", num / den, p );
  else
    snprintf( buf, bufsize, "%3.2f %s", (double)num / den, p );
  return buf;
  }

const double accuracy = 3.0e-7;			// relative accuracy

/* Return the incomplete gamma function P(a, x) evaluated by its series
   representation.
*/
double gamma_ser( const double a, const double x )
  {
  if( x <= 0 ) return 0;
  double ap = a;
  double sum = 1 / a;
  double del = sum;
  for( int i = 0; i < 100; ++i )
    {
    ++ap;
    del *= x / ap;
    sum += del;
    if( std::fabs( del ) < std::fabs( sum ) * accuracy )
      return sum * std::exp( -x + a * std::log( x ) - lgamma(a));
    }
  internal_error( "function gamma_ser needs more than 100 iterations to converge." );
  return 0;
  }


/* Return the incomplete gamma function Q(a, x) evaluated by its continued
   fraction representation.
*/
double gamma_cf( const double a, const double x )
  {
  const double epsilon = 1.0e-30;	// Number near the smallest
					// representable floating-point number.
  double b = x + 1 - a;		// Set up for evaluating continued fraction
  double c = 1 / epsilon;	// by modified Lentz's method (§5.2)
  double d = 1 / b;		// with b0 = 0.
  double h = d;
  for( int i = 1; i <= 100; ++i )
    {
    double an = -i * ( i - a );
    b += 2;
    d = an * d + b;
    if( std::fabs( d ) < epsilon ) d = epsilon;
    c = b + an / c;
    if( std::fabs( c ) < epsilon ) c = epsilon;
    d = 1 / d;
    double del = d * c;
    h *= del;
    if( std::fabs( del - 1 ) < accuracy )
      return std::exp( -x + a * std::log( x ) - lgamma( a ) ) * h;
    }
  internal_error( "function gamma_cf needs more than 100 iterations to converge." );
  return 0;
  }


/* Return the incomplete gamma function Q(a, x) = 1 - P(a, x).
   Adapted from Numerical Recipes in C, chapter 6.2.
*/
double gamma_q( const double a, const double x )
  {
  if( x < 0 || a <= 0 ) internal_error( "invalid arguments in routine gamma_q." );
  if( x < a + 1 )	// complement of the series representation
    return 1 - gamma_ser( a, x );
  return gamma_cf( a, x );	// continued fraction representation
  }

} // end namespace


/* Show how well the frequency of sequences of N repeated bytes in LZMA data
   matches the value expected for random data. ( 1 / 2^( 8 * N ) )
   Print cumulative data for all files followed by the name of the first
   file with the longest sequence.
*/
int print_nrep_stats( const std::vector< std::string > & filenames,
                      const Cl_options & cl_opts, const int repeated_byte,
                      const bool ignore_errors )
  {
  unsigned long long hist[256] = {};		// histogram of LZMA data
  std::vector< unsigned long > len_vector;
  unsigned long long udata_size = 0;		// total uncompressed size
  unsigned long long lzma_size = 0;		// total size of LZMA data
  unsigned long best_pos = 0;
  unsigned long members = 0;
  int best_name = -1, retval = 0;
  const bool count_all = repeated_byte < 0 || repeated_byte >= 256;
  bool stdin_used = false;
  for( unsigned i = 0; i < filenames.size(); ++i )
    {
    const bool from_stdin = filenames[i] == "-";
    if( from_stdin ) { if( stdin_used ) continue; else stdin_used = true; }
    const char * const input_filename =
      from_stdin ? "(stdin)" : filenames[i].c_str();
    struct stat in_stats;				// not used
    const int infd = from_stdin ? STDIN_FILENO :
      open_instream( input_filename, &in_stats, false, true );
    if( infd < 0 ) { set_retval( retval, 1 ); continue; }

    const Lzip_index lzip_index( infd, cl_opts, ignore_errors, ignore_errors );
    if( lzip_index.retval() != 0 )
      {
      show_file_error( input_filename, lzip_index.error().c_str() );
      set_retval( retval, lzip_index.retval() );
      close( infd );
      continue;
      }
    const unsigned long long cdata_size = lzip_index.cdata_size();
    if( !fits_in_size_t( cdata_size ) )		// mmap uses size_t
      { show_file_error( input_filename, large_file_msg );
        set_retval( retval, 1 ); close( infd ); continue; }
    const uint8_t * const buffer =
      (const uint8_t *)mmap( 0, cdata_size, PROT_READ, MAP_PRIVATE, infd, 0 );
    close( infd );
    if( buffer == MAP_FAILED )
      { show_file_error( input_filename, mmap_msg, errno );
        set_retval( retval, 1 ); continue; }
    members += lzip_index.members();
    udata_size += lzip_index.udata_size();
    for( long j = 0; j < lzip_index.members(); ++j )
      {
      const Block & mb = lzip_index.mblock( j );
      long pos = mb.pos() + 7;			// skip header (+1 byte) and
      const long end = mb.end() - 24;		// 4 + trailer of each member
      lzma_size += end - pos;
      while( pos < end )
        {
        const uint8_t byte = buffer[pos++];
        ++hist[byte];
        if( buffer[pos] == byte )
          {
          unsigned len = 2;
          ++pos;
          while( pos < end && buffer[pos] == byte ) { ++pos; ++len; }
          hist[byte] += len - 1;
          if( !count_all && repeated_byte != (int)byte ) continue;
          if( len >= len_vector.size() ) { len_vector.resize( len + 1 );
            best_name = i; best_pos = pos - len; }
          ++len_vector[len];
          }
        }
      }
    munmap( (void *)buffer, cdata_size );
    }

  if( lzma_size == 0 ) { show_error( "No LZMA data found." ); return 1; }
  if( verbosity < 0 ) return retval;
  // G-test (improved type of chi-squared test)
  const double expected = 1.0 / 256;
  double gval = 0;			// 2 * sum( oi * ln( oi / ei ) )
  for( unsigned i = 0; i < 256; ++i )
    {
    const double observed = (double)hist[i] / lzma_size;
    if( observed > 0 ) gval += observed * log( observed / expected );
    }
  gval *= 2;

  std::printf( "\nTotal size of LZMA data: %s bytes (%sB) in %s %s\n",
               format_num3( lzma_size ), format_num( lzma_size ),
               format_num3( members ), (members == 1) ? "member" : "members" );
  if( udata_size > 0 && lzma_size > 0 )
    std::printf( "%5.3f:1, %5.2f%% ratio, %5.2f%% saved.\n",
                 (double)udata_size / lzma_size,
                 ( 100.0 * lzma_size ) / udata_size,
                 100.0 - ( ( 100.0 * lzma_size ) / udata_size ) );
  std::printf( "G-test for uniformity: G = %f, p-value = %g\n", gval,
               gamma_q( 0.5, gval / 2 ) );	// chi^2 probability function
  if( len_vector.size() == 0 )
    { std::fputs( "No repeated byte sequences found.\n", stdout );
      return retval; }
  if( count_all )
    std::fputs( "Showing repeated sequences of any byte value:\n", stdout );
  else
    std::printf( "Showing repeated sequences of the byte value 0x%02X:\n",
                 repeated_byte );
  for( unsigned len = 2; len < len_vector.size(); ++len )
    if( len_vector[len] > 0 )
      std::printf( "len %u found %s times, 1 every %s bytes "
                   "(expected 1 every %s B)\n",
                   len, format_num3( len_vector[len] ),
                   format_num3( lzma_size / len_vector[len] ),
                   format_num3( 1ULL << ( 8 * ( len - count_all ) ) ) );
  if( best_name >= 0 )
    std::printf( "Longest sequence found at position %s of '%s'\n",
                 format_num3( best_pos ), filenames[best_name].c_str() );
  return retval;
  }
