##############################################################
# VectorDouble.pm
#
#    See POD at the end of the file for documentation.
#
#    This is a demonstration only, though perhaps
#    it'll turn into something like Math::Vector
#    if enough parts get added to make it more useful.
#
##############################################################
package VectorDouble;
require 5.6.0;
our $VERSION = '0.25';

require Exporter;
our @ISA = qw(Exporter);

###############################################
# Implementation Notes
#
#  This package used both overload' and 'tie'.
#
#  The 'overload' capability defines how these objects
#  behave in a variety of perl syntax structures
#    $a+$b
#    print $a
#  and so on.
#
#  It also can defines an array reference to be supplied
#  when the object is treated as an array, via
#    $a->[2] = 3
#  or
#    print $a->[2]
#
#  However, to actually get these subscipts to do anything interesting,
#  the perl array supplied by overload must be 'tied' to a set of
#  routines that implement perl's array reference behavior.
#
##################################################

sub new {
   my $class = shift;
   $class = ref($class) || $class;
   my $self = bless { size => 1, @_ }, $class;
   my @tiedarray;
   tie @tiedarray, VectorDouble, $self;
   my ($setValue, $value, $setSeq, $first, $by) = (0,0,0,0,0);
   if ( exists $self->{value} ) {
       $value = $self->{value};
       $setValue = 1; }
   elsif ( exists $self->{first} ) {
       $first = $self->{first};
       $by = $self->{by} || 1;
       $setSeq = 1;
   }
   $self->{ _tied_array } = \@tiedarray;
   $self->{ _cdata } = vec_allocate( $self->{size},
                                    $setValue, $value,
                                    $setSeq, $first, $by );
   return $self;
};

sub size {
     return shift->{size};
}

# Free the C data structures when the object is garbage collected.
sub DESTROY {
   vec_free( shift->{_cdata} );
}

# --- overload methods -------------------------------------------

use overload
   '+'    => \&add,
   '*'    => \&multiply,
   '""'   => \&stringify,
   '@{}'  => \&as_array,
   '0+'   => \&numify,
   'bool' => \&boolify,
   ;

sub add {
     my ($self, $other, $swap) = @_;
     if ($self->{size} != $other->{size}) {
        # generate warning?
        return;
     }
     my $c = VectorDouble->new( size => $self->{size} );
     if (defined($swap) && $swap) {
        vec_add( $other->{_cdata}, $self->{_cdata}, $c->{_cdata}); }
     else {
        vec_add( $self->{_cdata}, $other->{_cdata}, $c->{_cdata}); }
     return $c;
}

sub multiply {
     my ($self, $other, $swap) = @_;
     if ($self->{size} != $other->{size}) {
        # generate warning?
        return;
     }
     my $c = VectorDouble->new( size => $self->{size} );
     if (defined($swap) && $swap) {
        vec_multiply( $other->{_cdata}, $self->{_cdata}, $c->{_cdata}); }
     else {
        vec_multiply( $self->{_cdata}, $other->{_cdata}, $c->{_cdata}); }
     return $c;
}

# stringify invokes the C routine vec_fetch() to find
# numerical values for the elements of the vector,
# then puts lets perl convert them to strings and concatenate them.
# If we want to do this for long vectors it'd be
# faster (though the string maniupulation would be much messier)
# to do this whole thing in C.
my $SizeToPrint = 10;  # print all elements if size is smaller than this.
sub stringify {
     my $self = shift;
     my $SizeToPrint = 10;
     my $size = $self->{size};
     my $s = '[';
     if ($size < $SizeToPrint) {
        for (my $i=0; $i<$size-1; $i++) {
            $s .= $self->[$i] . ',';
        }
        $s .= $self->[$size-1] . ']';
     }
     else {
        $s .=   $self->[0] . ','
              . $self->[1] . ','
              . $self->[2] . ',...,'
              . $self->[$size-3] . ','
              . $self->[$size-2] . ','
              . $self->[$size-1] . ']';
     }
     return $s;
}

# The array methods are 'tied' to an array reference
# which is returned with this call, invoked via
# the 'overload' mechanism.
sub as_array {
     return shift->{_tied_array};
}

# This object is always true in a boolean context.
sub boolify {
     return 1;
}

# Returned when the object is converted to a single number,
# as in "$a = new VectorDouble size=>10; print $#$a;"
sub numify {
     return shift->{size};
}

# --- TIE-ed methods --------------------------------------------

sub TIEARRAY {
     my ($class, $self) = @_;
     return $self;
}

sub FETCH {
     my ($self, $index) = @_;
     return vec_fetchDouble( $self->{_cdata}, $index );
}

sub STORE {
     my ($self, $index, $value) = @_;
     vec_storeDouble( $self->{_cdata}, $index, $value );
}

sub FETCHSIZE {
     my ($self) = @_;
     return $self->{size};
}

sub EXISTS {
   my ($self,$index) = @_;
   return $index>=0 && $index<$self->{size};
}

# The ther 'tie' array functions which could be
# defined are DELETE, PUSH, POP, SHIFT, UNSHIFT, SPLICE, CLEAR

# --- C methods -----------------------------------------

use Inline C => <<'END_C';

#define MinSize 1

  typedef struct {
     int       size;
     double*  numbers;
   } cVectorDouble;

   /* Return a perl SV who's IV (integer value) is a pointer
      to a cVectorDouble struct.  This technique (from Inline::C-Cookbook)
      embeds an arbitrary C structure with a single perl scalar that
      we can then pass back and forth.  Memory management of
      these internal structures must be coded explicity.  */
   /* Not yet done:
        Add error tests for malloc failure, if size is too big.
        Or bypass mallocs by using perl strings (SV's) for memory ? */
   SV* vec_allocate(int size_in,
                   bool setValue, double value,
                   bool setSeq, double first, double by  ) {
     cVectorDouble*   vector  = malloc(sizeof(cVectorDouble));
     SV*              data    = newSViv(0);
     double*          nums;
     int              i;
     if ( size_in < MinSize)  { size_in = MinSize; }
     nums = malloc( size_in * sizeof(double) );
     vector->size  = size_in;
     vector->numbers = nums;
     sv_setiv( data, (IV)vector );
     SvREADONLY_on(data);
     /* --- initialize the vector values --- */
     if (setValue) {
        for (i=0; i < vector->size; i++){
            nums[i] = value;
        }
     }
     else if (setSeq) {
        for (i=0; i < vector->size; i++){
            nums[i] = first;
            first += by;
        }
     }
     return data;
   }

   /* Returns j'th data  */
   double vec_fetchDouble( SV* data, int j ){
     cVectorDouble* vector = (cVectorDouble*)SvIV(data);
     j = j % vector->size;   /* make sure 0 <= j <= (size-1) */
     return vector->numbers[j];
   }

   /* Sets j'th data = value  */
   double vec_storeDouble( SV* data, int j, double value){
       cVectorDouble* vector = (cVectorDouble*)SvIV(data);
       j = j % vector->size;   /* make sure 0 <= j <= (size-1) */
       vector->numbers[j] = value;
       return value;
   }

   void vec_free( SV* data ){
     cVectorDouble* vector = (cVectorDouble*)SvIV(data);
     free(vector->numbers);
     free(vector);
   }

   /* c = a*b.  a,b,c should have the same size before we get here. */
   void vec_multiply( SV* a, SV* b, SV* c ){
       cVectorDouble* va = (cVectorDouble*)SvIV(a);
       cVectorDouble* vb = (cVectorDouble*)SvIV(b);
       cVectorDouble* vc = (cVectorDouble*)SvIV(c);
       int i;
        for (i=0; i<va->size; i++) {
           vc->numbers[i] = va->numbers[i] * vb->numbers[i];
        }
   }

   /* c = a+b.  a,b,c should have the same size before we get here. */
   void vec_add( SV* a, SV* b, SV* c ){
       cVectorDouble* va = (cVectorDouble*)SvIV(a);
       cVectorDouble* vb = (cVectorDouble*)SvIV(b);
       cVectorDouble* vc = (cVectorDouble*)SvIV(c);
       int i;
        for (i=0; i<va->size; i++) {
           vc->numbers[i] = va->numbers[i] + vb->numbers[i];
        }
   }

END_C

if (not caller) {     # Do this only if we weren't called
                       # from another module.  See the "SYNOPSIS" below.
     print "  --  VectorDouble.pm -- \n";
     print "  Arrays of packed doubles implemented in C \n";
     print "  Demo : \n\n";

     print " \$a = new VectorDouble size=>5000, value=>10.1 ; \n";
     my $a = new VectorDouble size => 5000, value => 10.1 ;
     print " \$a is " . $a . "\n\n";

     print " \$b = new VectorDouble size=>5000, first=>1, by=>3 ; \n";
     my $b = VectorDouble->new( size => 5000, first => 1, by => 3 );
     print " \$b is " . $b . "\n\n";

     print " Setting \$a->[1] to 22 :  \n";
     $a->[1] = 22;
     print " \$a->[1] is " . $a->[1] . "\n\n";

     my $c = $a * $b + $a;
     print " \$a * \$b + \$a is " . $c . "\n\n";
}

1;

##########################################################################

=head1 NAME

VectorDouble - a perl scalar reference to a fixed packed array of
doubles with arithmetic operations on the scalars, implemented in C.

=head1 SYNOPSIS

  use VectorDouble;

  my $a = new VectorDouble size => 5000, value => 10.1 ;
  print " \$a is " . $a . "\n";

  my $b = new VectorDouble size => 5000, first => 1, by => 3 ;
  print " \$b is " . $b . "\n";

  $a->[1] = 22;
  print " \$a->[1] is " . $a->[1] . "\n";

  my $c = $a * $b + $a;
  print " \$a * \$b + \$a is " . $c . "\n";

=head1 DESCRIPTION

A demonstration of Inline for the IPL class in the MSIE program
at the Gradcenter of Marlboro College.

To see a demo, type

     perl VectorDouble.pm

at the shell prompt.

The currently implemented operations are

  * creating a constant array
     $a = new VectorDouble size=>10, value=>2;

  * creating a sequence
     $b = new VectorDouble size=>10, first=>1, by=>2

  * addition and multiplication
     $a + $b ; $a * $b

  * setting and fetching individual values
      $a->[2] = $b->[5]

  * testing for the existance of a value
      if ( exists $a->[101] ) { ... }

  * string representation
      print $a;

The size of a VectorDouble is fixed when it is created,
and can be found from either $a->size or the perlism $#$a

Like perl arrays, the first element is always $a->[0] and the
last is $a->[size-1].  All subscripts outside this range
are converted to ones within this range with j=j % size.

A full implementation would also include things like

  * more binary operations (-, /, **, etc)

  * arithmetic with perl scalars,
    i.e. "$a=new VectorDouble; $b = $a + 1;"

  * mathematical functions like sin, cos, ...

  * some kind iteration,
    i.e. functions (perl or C) applied to each element.

  * better error testing and return values

=head1 AUTHOR

Jim Mahoney, Marlboro College (mahoney@marlboro.edu)

=head1 LICENSE

Same terms as Perl.

=head1 SEE ALSO

perl(1), PDL

=cut