#!/usr/bin/perl

use PDL;
use List::Util;

sub candidate_recurrence {
    my ($d, $lin, @values) = @_;
    my $N = $d+$lin;
    my $samp = 2*$N -$lin; ## number of samples we'll need
    
    ## some base cases
    return if $N == 0;
    return Recurrence->new( 0, 1, $values[0] ) if $d == 0 and $lin == 1;
    return if @values < $samp; # need enough points
    
    if ($d == 1 and $lin == 0) {
      return $values[0] ? Recurrence->new( 1, 0, $values[1] / $values[0] ) : undef;
    }

    my @samples = @values[ 0 .. $samp-1 ];

    my @matrix;
    for (0 .. $N-1) {
        push @matrix, [ @samples[ $_ .. $_+$d-1 ] ];
    }
    if ($lin) {
      push @$_, 1 for @matrix;
    }

    my $M = pdl @matrix;
    my $v = pdl( [ @samples[ $d .. $samp-1 ] ] )->transpose;
    
    return if $M->det == 0;
    return Recurrence->new( $d, $lin, list( $M->inv() x $v ) );
}

sub check_recurrence {
    my $R = shift;
    for my $n ($R->d .. $#_-1) {
        return if $R->next( @_[ 0 .. $n ] ) != $_[$n+1];
    }
    return 1;
}

sub identify {
    for my $d (0 .. 2) {
        for my $lin (0 .. 1) {
            my $R = candidate_recurrence($d, $lin, @_);
            return $R if $R and check_recurrence($R, @_);
        }
    }
}


sub series {
    my $R = identify(@_) or return;
    return $R->next(@_);
}

############
## cute recurrence class to make things simpler

{
 package Recurrence;
 sub new {
    my ($pkg, $d, $lin, @coeff) = @_;
    my $const = $lin ? pop @coeff : 0;
    return bless { d=>$d, lin=>$lin, coeff=>\@coeff, const=>$const }, $pkg;
 }
  
 sub next {
    my $self = shift;
    return if @_ < $self->{d};
    my @window =  @_[ -$self->{d} .. -1 ];
    
    $self->{const}
      + List::Util::sum map { $self->{coeff}[$_] * $window[$_] } 0 .. $#window;    
 }
  
 sub print {
    my $self = shift;
    my $d = $self->d;
    my @terms = grep { !/^0\*/ }
                map { $self->{coeff}[$d-$_] . "*a(n-$_)" } 1 .. $self->{d};
    push @terms, $self->{const} if $self->{const} || @terms == 0;
    my $str = "a(n) = " . join " + ", @terms;
    $str =~ s/\b 1\*//gx;
    return $str; }
  
 sub d { return $_[0]->{d}; }
}