Math::Orthonormalize - Gram-Schmidt Orthonormalization of vectors

  use Math::Orthonormalize qw(:all);

  my @base_of_r_2 = (
    [2, 1],
    [1, 3]
  );
  my $vector = [1, 2, 3];

  my @orthonormalized = orthonormalize(@base_of_r_2);
  my @orthogonalized  = orthogonalize(@base_of_r_2);

  my $normalized      = normalize($vector);
  my $scaled          = scale(2, $vector);
  my $scalar          = scalar_product($vector1, $vector2);

package Math::Orthonormalize;

use 5.006;
use strict;
use warnings;

use Math::Symbolic qw/parse_from_string/;
use Carp;

our $VERSION = '1.00';

require Exporter;
our @ISA = qw(Exporter);
our %EXPORT_TAGS = ( 'all' => [ qw(
		orthonormalize
		orthogonalize
		scalar_product
		normalize
		scale
) ] );
our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );
our @EXPORT = ();

sub orthonormalize {
	my @vectors = @_;
	croak "No arguments to orthogonalize"
	  if not @vectors;
	croak "Arguments to orthogonalize must be array refs (vectors)"
	  if grep {ref($_) ne 'ARRAY'} @vectors;
	my $dim = @{$vectors[0]};
	croak "Vectors must have same dimension for orthogonalization"
	  if grep {@$_ != $dim} @vectors;

	my @base = orthogonalize(@vectors);
	@base = map {normalize($_)} @base;
	return @base;
}

sub orthogonalize {
	my @vectors = @_;
	croak "No arguments to orthogonalize"
	  if not @vectors;
	croak "Arguments to orthogonalize must be array refs (vectors)"
	  if grep {ref($_) ne 'ARRAY'} @vectors;
	my $dim = @{$vectors[0]};
	croak "Vectors must have same dimension for orthogonalization"
	  if grep {@$_ != $dim} @vectors;
	
	my @newbase;
	push @newbase, $vectors[0];
	my @squares;
	push @squares, scalar_product($vectors[0], $vectors[0]) if @vectors > 1;
	foreach my $i (1..$#vectors) {
		my $new = $vectors[$i];
		foreach my $j (0..$#newbase) {
			my $vec = scale(
				scalar_product($vectors[$i], $newbase[$j])
				/ $squares[$j],
				$newbase[$j]
			);
			@$new = map {$_ -= shift @$vec} @$new;
		}
		push @newbase, $new;
		push @squares, scalar_product($new, $new) if $i < $#vectors;
	}
	return @newbase;
}

sub normalize {
	my $v = shift;
	croak "Argument to normalize() must be an array ref (vector)"
	  if not ref($v) eq 'ARRAY';
	croak "Cannot normalize 0-dimensional vectors"
	  if @$v == 0;
	my $sum = $v->[0]**2;
	$sum += $v->[$_]**2 for 1..$#$v;
	$sum = sqrt($sum);
	croak "Cannot normalize 0-vector"
	  if $sum == 0;
	return [map {$_ / $sum} @$v];
}

sub scale {
	croak "scale() takes a scalar and a vector (ary ref) as arguments"
	  if not @_ == 2 or not ref($_[1]) eq 'ARRAY';
	croak "Cannot handle 0-dimensional vectors"
	  if @{$_[1]} == 0;
	my $new = [];
	foreach (@{$_[1]}) {
		push @$new, $_[0] * $_;
	}
	return $new;
}

sub scalar_product {
	croak "Invalid number of arguments to scalar_product()"
	  if not @_ == 2;
	my ($v1, $v2) = @_;
	croak "Cannot compute the scalar product of vectors of different ",
		"dimensions"
	  if @$v1 != @$v2;
	croak "Cannot deal with 0-dimensional vectors"
	  if @$v1 == 0;
	my $sum = $v1->[0] * $v2->[0];
	return $sum if @$v1 == 1;
	
	$sum += $v1->[$_] * $v2->[$_] for 1..$#$v1;
	return $sum;
}


1;
__END__