File Coverage

blib/lib/Math/NumSeq/BalancedBinary.pm

Criterion	Covered	Total	%
statement	219	229	95.6
branch	55	68	80.8
condition	20	26	76.9
subroutine	26	26	100.0
pod	10	10	100.0
total	330	359	91.9

# Copyright 2012, 2013, 2014 Kevin Ryde

# This file is part of Math-NumSeq.
#
# Math-NumSeq 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 3, or (at your option) any later
# version.
#
# Math-NumSeq 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 Math-NumSeq.  If not, see .


# http://www.iki.fi/~kartturi/matikka/tab9766.htm
# A009766 Catalan triangle
# A099039 Riordan array

package Math::NumSeq::BalancedBinary;
use 5.004;
use strict;
use List::Util 'max';

use vars '$VERSION', '@ISA';
$VERSION = 72;

use Math::NumSeq;
@ISA = ('Math::NumSeq');
*_is_infinite = \&Math::NumSeq::_is_infinite;
*_to_bigint = \&Math::NumSeq::_to_bigint;

use Math::NumSeq::Repdigits;
*_digit_split_lowtohigh = \&Math::NumSeq::Repdigits::_digit_split_lowtohigh;

use Math::NumSeq::Fibonacci;
*_blog2_estimate = \&Math::NumSeq::Fibonacci::_blog2_estimate;

# uncomment this to run the ### lines
# use Smart::Comments;


# use constant name => Math::NumSeq::__('Balanced Binary');
use constant description => Math::NumSeq::__('Bits 1,0 balanced like parentheses.');
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;
use constant default_i_start => 1;
use constant values_min => 2;

# pred() works, next() doesn't.
# Any merit in bit-reversed form ?
#
# use constant parameter_info_array =>
#   [ { name      => 'direction',
#       share_key => 'pairs_fsb',
#       display   => Math::NumSeq::__('Direction'),
#       type      => 'enum',
#       default   => 'HtoL',
#       choices => ['HtoL','LtoH'],
#       choices_display => [Math::NumSeq::__('HtoL'),
#                           Math::NumSeq::__('LtoH'),
#                          ],
#     },
#   ];

sub new {
  my $self = shift->SUPER::new(@_);
  $self->{'direction'} ||= 'HtoL';
  return $self;
}

#------------------------------------------------------------------------------
# A063171 same in binary
# A080116 predicate 0,1
# A080300 inverse, ranking value->i or 0
# A080237 num trailing zeros
# A085192 first diffs
# A000108 Catalan numbers, count of values in 4^k blocks
# A071152 balanced/2, Lukasiewicz words for the rooted plane binary trees
# A075171 trees by binary runs in integers, coded to Dyck words
# A071153 Lukasiewicz coded digits
#
my %oeis_anum = (HtoL => 'A014486',  # balanced binary
                 # OEIS-Catalogue: A014486

                 LtoH => undef,
                );
sub oeis_anum {
  my ($self) = @_;
  return $oeis_anum{$self->{'direction'}};
}


#------------------------------------------------------------------------------

# When $self->{'wantlen'} becomes equal to _WANTLEN_TO_BIGINT then go to
# Math::BigInt since the bits of a UV would be overflowed.
#
# _WANTLEN_TO_BIGINT is floor(numbits(UV)/2), but calculated by a probe of
# "~0".
#
# On a 32-bit UV balance binary value=0xFFFF0000 is reached at i=48_760_366
# which is fairly big but might be reached by a long running loop.  On a
# 64-bit UV the corresponding limit is i=75*10^15 so would not be reached in
# any reasonable time.
#
use constant 1.02 _WANTLEN_TO_BIGINT => do {
  my $uv = ~0;
  my $limit = 0;
  while ($uv >= 3) {
    $uv >>= 2;
    $limit++;
  }
  ### $limit
  $limit
};

sub rewind {
  my ($self) = @_;
  $self->{'wantlen'} = 0;
  $self->{'bits'} = [];
  $self->{'value'} = 0;
  $self->{'i'} = $self->i_start;
}

sub next {
  my ($self) = @_;
  ### BalancedBinary next() ...
  my $bits = $self->{'bits'};
  my $value;

  if (! @$bits) {
    ### initial 2 ...
    push @$bits, 2;
    $value = 2;

  } else {
    $value = $self->{'value'};
    for (;;) {
      ### at: "bits=".join(',',@{$self->{'bits'}})
      if (scalar(@$bits) < 2) {
        ### extend to wantlen: $self->{'wantlen'}+1
        $value = ($bits->[0] *= 4);
        if (++$self->{'wantlen'} == _WANTLEN_TO_BIGINT) {
          ### promote to BigInt from now on ...
          $bits->[0] = $value = _to_bigint($value);
        }
        last;
      }

      $value -= $bits->[-1];
      my $bit = ($bits->[-1] *= 2);
      if ($bit < $bits->[-2]) {
        ### shifted bit ...
        $value += $bit;
        last;
      }
      ### drop this bit ...
      pop @$bits;
    }

    ### pad for: "value=$value  bits=".join(',',@$bits)
    # trailing bits ...,128, 32, 8, 2
    my $bit = 2 + ($bits->[0]&1); # inherit BigInt
    foreach my $pos (reverse scalar(@$bits) .. $self->{'wantlen'}) {
      ### pad with: $bit
      $bits->[$pos] = $bit;
      $value += $bit;
      $bit *= 4;
    }
  }

  ### return: "value=$value bits=".join(',',@$bits)
  return ($self->{'i'}++,
          $self->{'value'} = $value);
}

sub ith {
  my ($self, $i) = @_;
  ### BalancedBinary ith(): $i

  if ($i < 1) {
    return undef;
  }
  if (_is_infinite($i)) {
    return $i;
  }

  $i -= 1;
  ### initial i remainder: $i

  my $zero = ($i*0);

  my @num;
  $num[0][0] = 0;
  $num[1][0] = 1;
  $num[1][1] = 1;
  my $z = 1;

  for ( ; ; $z++) {
    my $prev = $num[$z][0] = 1;  # all zeros, no ones
    if ($z == 16) {
      $prev += $zero;
      if (! ref $zero) {
        ### promote to BigInt ...
        $zero = _to_bigint($zero);
      }
    }
    foreach my $o (1 .. $z) {
      $num[$z][$o]
        = ($prev                        # 1...   $num[$z][$o-1]
           += ($num[$z-1][$o] || 0));   # 0... if $z>=1
    }
    my $catalan = $num[$z][$z];
    if ($i < $catalan) {
      last;
    }
    ### subtract catalan: $catalan
    $i -= $catalan;
  }

  ### i remaining: $i
  my @ret = (1);
  my $o = $z-1;
  while ($o >= 1) {
    ### at: "i=$i  z=$z o=$o  ret=".join('',@ret)
    ### assert: $z >= $o

    if ($z > $o) {
      ### compare: "z=".($z-1).",o=$o  num=".($num[$z-1][$o]||'undef')
      my $znum = $num[$z-1][$o];
      if ($i < $znum) {
        ### 0 ...
        push @ret, 0;
        $z--;
        next;
      }
      $i -= $znum;
    }
    ### 1 ...
    push @ret, 1;
    $o--;
  }

  push @ret, (0) x $z;
  ### final: "ret=".join('',@ret)

  if (! ref $zero && @ret >= 2*_WANTLEN_TO_BIGINT) {
    ### promote to BigInt ...
    $zero = _to_bigint($zero);
  }

  @ret = reverse @ret;
  ### return: _digit_join(\@ret, 2, $zero)
  return _digit_join(\@ret, 2, $zero);
}

# $aref->[0] low digit
sub _digit_join {
  my ($aref, $radix, $zero) = @_;
  my $n = $zero;
  foreach my $digit (reverse @$aref) { # high to low
    $n *= $radix;
    $n += $digit;
  }
  return $n;
}

sub value_to_i {
  my ($self, $value) = @_;
  ### BalancedBinary value_to_i(): $value

  if ($value < 2 || $value != int($value)) {
    return undef;
  }
  if (_is_infinite($value)) {
    return $value;
  }

  my @bits = _digit_split_lowtohigh($value,2)
    or return undef;

  _pred_on_bits($self,\@bits)
    or return undef;

  my @num;
  $num[0][0] = 0;
  $num[1][0] = 1;
  $num[1][1] = 1;
  my $w = scalar(@bits) / 2;
  ### $w

  my $zero = $value*0;
  my $i = 1 + $zero;

  foreach my $z (1 .. $w) {
    my $prev = $num[$z][0] = 1;  # all zeros, no ones
    if ($z > 16) {
      $prev += $zero;
    }
    foreach my $o (1 .. $z) {
      $num[$z][$o]
        = ($prev                        # 1...   $num[$z][$o-1]
           += ($num[$z-1][$o] || 0));   # 0... if $z>=1
    }
    $i += $num[$z-1][$z-1];
  }
  ### base i: $i

  ### bits: join('',@bits)
  shift @bits;  # skip high 1-bit
  my $z = $w;
  my $o = $w-1;
  foreach my $bit (@bits) {   # high to low
    ### at: "z=$z o=$o  bit=$bit"
    ### assert: $o >= 0
    ### assert: $z >= $o
    if ($bit) {
      ### bit 1 add: $num[$z-1][$o]
      $i += $num[$z-1][$o] || 0;
      $o--;
    } else {
      $z--;
    }
  }
  return $i;
}

sub value_to_i_floor {
  my ($self, $value) = @_;
  return _value_to_i_round ($self, $value, 0);
}
sub value_to_i_ceil {
  my ($self, $value) = @_;
  return _value_to_i_round ($self, $value, 1);
}

# Return the i corresponding to $value.
# If $value is not balanced binary then round according to $ceil.
# $ceil=1 to round up, or $ceil=0 to round down.
#
sub _value_to_i_round {
  my ($self, $value, $ceil) = @_;
  ### _value_to_i_round(): $value

  if ($value < 2) {
    return $self->i_start();
  }
  if (_is_infinite($value)) {
    return $value;
  }

  {
    my $int = int($value);
    if ($value != $int) {
      $value = $int + $ceil;  # +1 if not integer and want ceil
    }
  }
  my @bits = reverse _digit_split_lowtohigh($value,2);

  if (scalar(@bits) & 1) {
    # ENHANCE-ME: this is Catalan cumulative, or cumulative+1
    ### odd num bits ...
    @bits = ((0) x scalar(@bits), 1);
  }
  ### assert: (scalar(@bits) % 2) == 0

  my @num;
  $num[0][0] = 0;
  $num[1][0] = 1;
  $num[1][1] = 1;
  my $w = scalar(@bits)/2;
  ### $w

  my $zero = $value*0;
  my $i = 1 + $zero;

  foreach my $z (1 .. $w) {
    my $prev = $num[$z][0] = 1;  # all zeros, no ones
    if ($z > 16) {
      $prev += $zero;
    }
    foreach my $o (1 .. $z) {
      $num[$z][$o]
        = ($prev                        # 1...   $num[$z][$o-1]
           += ($num[$z-1][$o] || 0));   # 0... if $z>=1
    }
    $i += $num[$z-1][$z-1];
  }
  ### base i: $i

  ### bits: join('',@bits)
  shift @bits;  # skip high 1-bit
  my $z = $w;
  my $o = $w-1;

  foreach my $bit (@bits) {   # high to low
    ### at: "z=$z o=$o  bit=$bit"
    ### assert: $o >= 0
    ### assert: $z >= $o
    if ($bit) {
      if ($o == 0) {
        ### all 1s used, rest round down to zeros, so done: $i + $ceil
        return $i + $ceil;
      }
      ### bit 1 add: $num[$z-1][$o]
      $i += $num[$z-1][$o] || 0;
      $o--;

    } else {
      if ($z == $o) {
        ### 0 out of place, round up to 101010 ...
        while ($o) {
          $i += $num[$o-1][$o] || 0;
          $o--;
        }
        return $i - (1-$ceil);
      }

      $z--;
    }
  }
  return $i;
}



# 1       10,
# 2     1010,
# 3     1100,
# 4   101010,
# 5   101100,
# 6   110010,
# 7   110100,
# 8   111000,
# 9 10101010,   170
#
# 10xxxx
# 1xxxx0
# num(width) = 2*num(width-1) + extra
# num(1) = 1
# num(2) = 2*1 = 2
# num(3) = 2*2 + 1 = 5
# total(width)
#   = num(1) + num(2) + num(3) + ... + num(width)
#   = 1 + 2*1 + 2*2+1 +

sub pred {
  my ($self, $value) = @_;
  ### BalancedBinary pred(): $value

  if ($value != int($value) || _is_infinite($value)) {
    return 0;
  }

  my @bits = _digit_split_lowtohigh($value,2)
    or return 0;
  return _pred_on_bits($self,\@bits);
}
sub _pred_on_bits {
  my ($self, $bits) = @_;
  ### _pred_on_bits(): $bits

  if (scalar(@$bits) & 1) {
    return 0;
  }

  if ($self->{'direction'} eq 'HtoL') {
    @$bits = reverse @$bits;
    ### reversed bits: @$bits
  }

  my @count = (0,0);
  foreach my $bit (@$bits) {
    $count[$bit]++;
    if ($count[0] > $count[1]) {
      return 0;
    }
  }
  ### @count
  return ($count[0] == $count[1]);
}

# w = log2(value) / 2
# Catalan = 4^w / (sqrt(Pi * w) * (w + 1))
#         = 4^(log2(v)/2) / (sqrt(pi*log2(v)/2) * (log2(v)/2+1))
#         = v / (sqrt(pi*log2(v)/2) * (log2(v)/2+1))
#         = v / (sqrt(pi*log2(v)/2) * (log2(v)+2)/2)
#         = v / (sqrt(pi*log2(v)/2)/2 * (log2(v)+2))
#         = v / (sqrt(pi/8*log2(v)) * (log2(v)+2))
#
# if value has 2*w bits then cumulative catalan
# i = sum w=0 to log2(value)/2
#       of value / (sqrt(pi/8*log2(value)) * (log2(value)+2))
#
# is cumulative close enough to the plain ?
#
# 2*w bits, ipart=Catalan(w) value=4^w
# so w = log4(value)
#    ipart = Catalan(log4(value))
#          =

sub value_to_i_estimate {
  my ($self, $value) = @_;
  ### value_to_i_estimate: $value

  if ($value <= 2) {
    return 1;
  }

  my $log2 = _blog2_estimate($value);
  if (! defined $log2) {
    $log2 = log($value) * (1/log(2));
  }
  $log2 = max($log2,1);

  if (ref $value && $value->isa('Math::BigInt')) {
    # oldish BigInt doesn't like to divide BigInt/NV
    require Math::BigFloat;
    $value = Math::BigFloat->new($value);
  }

  return max(1,
             int($value / (sqrt((3.141592/8)*$log2) * ($log2+1))));
}

1;
__END__