File Coverage

blib/lib/LCS/BV.pm

Criterion	Covered	Total	%
statement	161	161	100.0
branch	30	30	100.0
condition	50	51	100.0
subroutine	11	11	100.0
pod	5	5	100.0
total	257	258	100.0

package LCS::BV;

use 5.010001;
use strict;
use warnings;
our $VERSION = '0.13';
#use utf8;

our $width = int 0.999+log(~0)/log(2);

use integer;
no warnings 'portable'; # for 0xffffffffffffffff

sub new {
  my $class = shift;
  # uncoverable condition false
  bless @_ ? @_ > 1 ? {@_} : {%{$_[0]}} : {}, ref $class || $class;
}

# H. Hyyroe. A Note on Bit-Parallel Alignment Computation. In
# M. Simanek and J. Holub, editors, Stringology, pages 79-87. Department
# of Computer Science and Engineering, Faculty of Electrical
# Engineering, Czech Technical University, 2004.

sub LLCS {
  my ($self,$a,$b) = @_;

  #use integer;
  #no warnings 'portable'; # for 0xffffffffffffffff

  # TODO: maybe faster, if we have fewer expensive iterations
  #if (@$a < @$b) {
  #  my $temp = $a;
  #  $a = $b;
  #  $b = $temp;
  #}

  my ($amin, $amax, $bmin, $bmax) = (0, $#$a, 0, $#$b);

  while ($amin <= $amax and $bmin <= $bmax and $a->[$amin] eq $b->[$bmin]) {
    $amin++;
    $bmin++;
  }
  while ($amin <= $amax and $bmin <= $bmax and $a->[$amax] eq $b->[$bmax]) {
    $amax--;
    $bmax--;
  }

  #my %positions;

  if ($amax < $width ) {
    my %positions;
    #$positions{$a->[$_]}  |= 1 << ($_ % $width) for $amin..$amax;
    for ($amin..$amax) { $positions{$a->[$_]} |= 1 << ($_ % $width); }

    my $v = ~0;
    my ($p,$u);

    for ($bmin..$bmax) {
      $p = $positions{$b->[$_]} // 0;
      $u = $v & $p;
      $v = ($v + $u) | ($v - $u);
    }
    return $amin + _count_bits(~$v) + $#$a - $amax;
  }
  else {
    my %positions;
    #$positions{$a->[$_]}->[$_ / $width] |= 1 << ($_ % $width) for $amin..$amax;
    for ($amin..$amax) { $positions{$a->[$_]}->[$_ / $width] |= 1 << ($_ % $width); }

    my $S;
    my @Vs = ();      # $Vs->[$k] = bits;

    my ($p, $u, $carry);

    my $kmax = ($amax+1) / $width;
    $kmax++ if (($amax+1) % $width);

    for (my $k=0; $k < $kmax; $k++ ) { $Vs[$k] = ~0; }

    for my $j ($bmin..$bmax) {
      $carry = 0;
      for (my $k=0; $k < $kmax; $k++ ) {
        $S = $Vs[$k];
        $p = $positions{$b->[$j]}->[$k] // 0;
        $u = $S & $p;             # [Hyy04]
        $Vs[$k] = ($S + $u + $carry) | ($S - $u);
        $carry = (($S & $u) | (($S | $u) & ~($S + $u + $carry))) >> ($width-1) & 1;
      }
    }

    my $bitcount = 0;

    for my $k ( @Vs ) {
      $bitcount += _count_bits(~$k);
    }
    return $amin + $bitcount + $#$a - $amax;
  }
}

sub LLCS_prepared {
  my ($self,$positions,$b) = @_;

  #if ( ref($a) eq 'HASH' )    # $a contains $positions

  my ($bmin, $bmax) = (0, $#$b);
  my ($amin, $amax) = (0, $self->{'length_a'} - 1);

  if ($amax < $width ) {

    my $v = ~0;
    my ($p,$u);

    for ($bmin..$bmax) {
      $p = $positions->{$b->[$_]} // 0;
      $u = $v & $p;
      $v = ($v + $u) | ($v - $u);
    }
    return _count_bits(~$v);
  }
  else {

    my $S;
    my @Vs = ();      # $Vs->[$k] = bits;

    my ($p, $u, $carry);

    my $kmax = ($amax+1) / $width;
    $kmax++ if (($amax+1) % $width);

    for (my $k=0; $k < $kmax; $k++ ) { $Vs[$k] = ~0; }

    for my $j ($bmin..$bmax) {
      $carry = 0;
      for (my $k=0; $k < $kmax; $k++ ) {
        $S = $Vs[$k];
        $p = $positions->{$b->[$j]}->[$k] // 0;
        $u = $S & $p;             # [Hyy04]
        $Vs[$k] = ($S + $u + $carry) | ($S - $u);
        $carry = (($S & $u) | (($S | $u) & ~($S + $u + $carry))) >> ($width-1) & 1;
      }
    }

    my $bitcount = 0;

    for my $k ( @Vs ) {
      $bitcount += _count_bits(~$k);
    }
    return $bitcount;
  }
}


sub LCS {
  my ($self, $a, $b) = @_;

  #use integer;
  #no warnings 'portable'; # for 0xffffffffffffffff

  my ($amin, $amax, $bmin, $bmax) = (0, $#$a, 0, $#$b);

  while ($amin <= $amax and $bmin <= $bmax and $a->[$amin] eq $b->[$bmin]) {
    $amin++;
    $bmin++;
  }
  while ($amin <= $amax and $bmin <= $bmax and $a->[$amax] eq $b->[$bmax]) {
    $amax--;
    $bmax--;
  }

  my %positions;
  my @lcs;

  if ($amax < $width ) {
    $positions{$a->[$_]} |= 1 << ($_ % $width) for $amin..$amax;

    my $S = ~0;
    my @Vs = (~0);

    my ($y,$u);

    # outer loop
    for my $j ($bmin..$bmax) {
      $y = $positions{$b->[$j]} // 0;
      $u = $S & $y;               # [Hyy04]
      $S = ($S + $u) | ($S - $u); # [Hyy04]
      $Vs[$j] = $S;
    }

    # recover alignment
    my $i = $amax;
    my $j = $bmax;

    while ($i >= $amin && $j >= $bmin) {
      if ($Vs[$j] & (1<<$i)) {
        $i--;
      }
      else {
        unless (
           $j && ~$Vs[$j-1] & (1<<$i)
        ) {
           unshift @lcs, [$i,$j];
           $i--;
        }
        $j--;
      }
    }
  }
  else {
    $positions{$a->[$_]}->[$_ / $width] |= 1 << ($_ % $width) for $amin..$amax;

    my $S;
    my @Vs = ([~0]);
    my ($y,$u,$carry);

    my $kmax = ($amax+1) / $width;
    $kmax++ if (($amax+1) % $width);

    # outer loop
    for my $j ($bmin..$bmax) {
      for (my $k=0; $k < $kmax; $k++ ) { $Vs[$j]->[$k] = ~0; }
      $carry = 0;

      for (my $k=0; $k < $kmax; $k++ ) {
        $S = ($j > $bmin) ? $Vs[$j-1]->[$k] : ~0;
        $y = $positions{$b->[$j]}->[$k] // 0;
        $u = $S & $y;             # [Hyy04]

        $Vs[$j]->[$k] = ($S + $u + $carry) | ($S - $u);

        $carry = (($S & $u) | (($S | $u) & ~($S + $u + $carry))) >> ($width-1) & 1;
      }
    }

    # recover alignment
    my $i = $amax;
    my $j = $bmax;

    while ($i >= $amin && $j >= $bmin) {
      my $k = $i / $width;
      if ($Vs[$j]->[$k] & (1<<($i % $width))) {
        $i--;
      }
      else {
        unless (
           $j && ~$Vs[$j-1]->[$k] & (1<<($i % $width))
        ) {
           unshift @lcs, [$i,$j];
           $i--;
        }
        $j--;
      }
    }
  }

  return [
    map([$_ => $_], 0 .. ($bmin-1)), ## no critic qw(BuiltinFunctions::RequireBlockMap)
    @lcs,
    map([++$amax => $_], ($bmax+1) .. $#$b) ## no critic qw(BuiltinFunctions::RequireBlockMap)
  ];
}

sub prepare {
  my ($self, $a) = @_;

  $self->{'length_a'} = @{$a};

  my $positions;

  if ($#$a < $width ) {
    $positions->{$a->[$_]} |= 1 << ($_ % $width) for 0..$#$a;
  }
  else {
    $positions->{$a->[$_]}->[$_ / $width] |= 1 << ($_ % $width) for 0..$#$a;
  }
  return $positions;
}

sub _count_bits {
  my $v = shift;

  #use integer;
  #no warnings 'portable'; # for 0xffffffffffffffff

  if ($width == 64) {
    $v = $v - (($v >> 1) & 0x5555555555555555);
    $v = ($v & 0x3333333333333333) + (($v >> 2) & 0x3333333333333333);
    # (bytesof($v) -1) * bitsofbyte = (8-1)*8 = 56 ----------------------vv
    $v = (($v + ($v >> 4) & 0x0f0f0f0f0f0f0f0f) * 0x0101010101010101) >> 56;
    return $v;
  }
  else {
    #$v = $v - (($v >> 1) & 0x55555555);
    #$v = ($v & 0x33333333) + (($v >> 2) & 0x33333333);
    ## (bytesof($v) -1) * bitsofbyte = (4-1)*8 = 24 ------vv
    #$v = (($v + ($v >> 4) & 0x0f0f0f0f) * 0x01010101) >> 24

    my $c; # count
    for ($c = 0; $v; $c++) {
      $v &= $v - 1; # clear the least significant bit set
    }
    return $c;
  }
}

1;

__END__