File Coverage

blib/lib/Math/NumSeq/NumAronson.pm

Criterion	Covered	Total	%
statement	70	74	94.5
branch	19	22	86.3
condition	6	9	66.6
subroutine	15	15	100.0
pod	4	4	100.0
total	114	124	91.9

# Copyright 2010, 2011, 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 .

package Math::NumSeq::NumAronson;
use 5.004;
use strict;

use vars '$VERSION', '@ISA';
$VERSION = 71;
use Math::NumSeq;
@ISA = ('Math::NumSeq');

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


# use constant name => Math::NumSeq::__('Numerical Aronson');
use constant description => Math::NumSeq::__('Numerical version of Aronson\'s sequence');
use constant values_min => 1;
use constant characteristic_increasing => 1;
use constant characteristic_integer => 1;
use constant i_start => 1;

# cf A080596 - a(1)=1
#    A079253 - even
#    A081023 - lying
#    A014132 - lying opposite parity
use constant oeis_anum => 'A079000';

# a(9*2^k - 3 + j) = 12*2^k - 3 + (3/2)*j + (1/2)*abs(j)
#     where k>=0 and -3*2^k <= j < 3*2^k
# step
#     a(n+1) - 2*a(n) + a(n-1) = 1   if n=9*2^k-3, k>=0
#                              = -1  if n = 2 and 3*2^k-3, k>=1
#                              = 0   otherwise.
#
# lying
# g(3*2^k-1 + j) = 2*2^(k+1)-1 + (3/2)*j + (1/2)*abs(j)
# where -2^k <= j < 2^k  and k>0
#
# then lying is d(n)=g(n+1)-1  n>=1

sub rewind {
  my ($self) = @_;
  $self->{'i'} = $self->i_start;
  $self->{'pow'} = 0;
  $self->{'j'} = -1;
}

sub next {
  my ($self) = @_;
  my $pow = $self->{'pow'};
  my $j = ++ $self->{'j'};

  if ($pow == 0) {
    # low special cases initial 1,4,
    if ($j < 2) {
      return ($self->{'i'}++, $j*3 + 1);
    }
    $pow = $self->{'pow'} = 1;    # 2**k for k=0
    $j   = $self->{'j'}   = -3;   # -3*(2**k) for k=0
  } elsif ($j >= 3 * $pow) {
    $pow = ($self->{'pow'} *= 2);
    $j = $self->{'j'} = -3 * $pow;
  }

  ### assert: -3 * $pow <= $j
  ### assert: $j < 3 * $pow

  return ($self->{'i'}++, 12*$pow - 3 + (3*$j + abs($j))/2);
}

# i = 9*2^k - 3 + j
# base at j=-3*2^k
# is i >= 9*2^k - 3 - 3*2^k
#    i >= 6*2^k - 3
#    i+3 >= 6*2^k
#    (i+3)/6 >= 2^k
#    2^k <= (i+3)/6
# then i = 9*2^k - 3 + j
#      j = i - 9*2^k + 3
#
sub ith {
  my ($self, $i) = @_;
  ### NumAronson ith(): $i

  # special cases ith(1)=1, ith(2)=4
  if ($i <= 2) {
    if ($i < 0) {
      return undef;
    } else {
      return $i*$i;
    }
  }

  my ($pow, $k) = _round_down_pow (int(($i+3)/6), 2);
  my $j = $i - 9*$pow + 3;

  ### round down for k: ($i+3)/6
  ### $k
  ### $pow
  ### $j
  ### assert: $k >= 0
  ### assert: -3 * 2 ** $k <= $j
  ### assert: $j < 3 * 2 ** $k

  return 12*$pow - 3 + (3*$j + abs($j))/2;
}

# value = 12*2^k - 3 + (3/2)*j + (1/2)*abs(j)
# minimum j=-3*2^k
# value = 12*2^k - 3 + (3*j + abs(j))/2
#       = 12*2^k - 3 + (3*-3*2^k + 3*2^k)/2
#       = 12*2^k - 3 + (-9 + 3)*2^k/2
#       = 12*2^k - 3 + -6*2^k/2
#       = 12*2^k - 3 + -3*2^k
#       = 9*2^k - 3
# value >= 9*2^k - 3
# value+3 >= 9*2^k
# 9*2^k <= value+3
# 2^k <= (value+3)/9
#
# from which
#     value = 12*2^k - 3 + (3*j + abs(j))/2
#     (3*j + abs(j))/2 = value - 12*2^k + 3
#     3*j + abs(j) = 2*(value - 12*2^k + 3)
#
# j>=0  4*j = a, j=a/4
# j<0   3*$j-$j = 2*$j = a, j=a/2
#
sub pred {
  my ($self, $value) = @_;
  ### NumAronson pred(): $value

  # special cases pred(1) true, pred(4) true
  if ($value < 6) {
    return ($value == 1 || $value == 4);
  }

  my $k = _round_down_pow (int(($value+3)/9), 2);
  my $pow = 2**$k;
  my $aj = 2*($value - 12*$pow + 3);
  return ($aj % ($aj < 0 ? 2 : 4)) == 0;
}

#------------------------------------------------------------------------------
# generic

# return ($pow, $exp) with $pow = $base**$exp <= $n,
# the next power of $base at or below $n
#
sub _round_down_pow {
  my ($n, $base) = @_;
  ### _round_down_pow(): "$n base $base"

  if ($n < $base) {
    return (1, 0);
  }

  # Math::BigInt and Math::BigRat overloaded log() return NaN, use integer
  # based blog()
  if (ref $n && ($n->isa('Math::BigInt') || $n->isa('Math::BigRat'))) {
    my $exp = $n->copy->blog($base);
    return (Math::BigInt->new(1)->blsft($exp,$base),
            $exp);
  }

  my $exp = int(log($n)/log($base));
  my $pow = $base**$exp;

  ### n:   ref($n)."  $n"
  ### exp: ref($exp)."  $exp"
  ### pow: ref($pow)."  $pow"

  # check how $pow actually falls against $n, not sure should trust float
  # rounding in log()/log($base)
  # Crib: $n as first arg in case $n==BigFloat and $pow==BigInt
  if ($n < $pow) {
    ### hmm, int(log) too big, decrease...
    $exp -= 1;
    $pow = $base**$exp;
  } elsif ($n >= $base*$pow) {
    ### hmm, int(log) too small, increase...
    $exp += 1;
    $pow *= $base;
  }
  return ($pow, $exp);
}

1;
__END__