http://qs321.pair.com?node_id=132457
Category: Math
Author/Contact Info Mark Stratman (count0) mark@sporkstorms.org
Description: This is a module I was using for basic vector operations. I cleaned it up a bit and documented it for (possible) CPAN submission. I'd be most grateful =) if anyone would be so kind as to provide feedback. Some of the areas I was questionable about were:
  • The namespace.. perhaps Physics::Vector might be more appropriage?
  • The extensive operator overloading.. and more specifically, overloading '.' and 'x' with operations that are completely different than the default '.' and 'x' operators.
  • Anything you'd care to comment on! =)
#
#
# An object oriented module for working with vectors.
# Note that 'scalar' refers to the mathmatical definition, not Perl's 
+;)
#

package Math::Vector;
use strict;
use Carp;
use Math::Trig;

use vars qw($VERSION);
$VERSION = '0.01';

use overload
    "="  => \&assign,
    "+=" => \&add_assign,
    "-=" => \&subtract_assign,
    "+"  => \&add,
    "-"  => \&subtract,
    "x"  => \&cross_product,
    "."  => \&dot_product,
    "*"  => \&multiply_scalar,
    "/"  => \&divide_scalar,
    "*=" => \&multiply_scalar_assign,
    "/=" => \&divide_scalar_assign;


sub new {
    my ($package, %args) = @_;

    return bless {
        _i => $args{i} || 0,
        _j => $args{j} || 0,
        _k => $args{k} || 0,
    }, $package;
}


# $v_1 = $v_2;
sub assign {
    my ($self, $other) = @_;
    carp ('Operation only permitted on Math::Vector objects') && retur
+n undef
        unless (ref $other eq 'Math::Vector');

    $self->set_i($other->i);
    $self->set_j($other->j);
    $self->set_k($other->k);

    return $self;
}


# $v_1 += $v_2;
sub add_assign {
    my ($self, $other) = @_;
    carp ('Operation only permitted on Math::Vector objects') && retur
+n undef
        unless (ref $other eq 'Math::Vector');

    $self->set_i($other->i + $self->i);
    $self->set_j($other->j + $self->j);
    $self->set_k($other->k + $self->k);

    return $self;
}


# $v_1 -= $v_2;
sub subtract_assign {
    my ($self, $other) = @_;
    carp ('Operation only permitted on Math::Vector objects') && retur
+n undef
        unless (ref $other eq 'Math::Vector');

    $self->set_i($other->i - $self->i);
    $self->set_j($other->j - $self->j);
    $self->set_k($other->k - $self->k);

    return $self;
}


# $v_new = $v_1 + $v_2;
# $v_new = $v_1->add($v_2);
sub add {
    my ($self, $other) = @_;
    carp ('Operation only permitted on Math::Vector objects') && retur
+n undef
        unless (ref $other eq 'Math::Vector');

    return new Math::Vector (
        i => $self->i + $other->i,
        j => $self->j + $other->j,
        k => $self->k + $other->k,
    );
}


# $v_new = $v_1 - $v_2;
# $v_new = $v_1->subtract($v_2);
sub subtract {
    my ($self, $other) = @_;
    carp ('Operation only permitted on Math::Vector objects') && retur
+n undef
        unless (ref $other eq 'Math::Vector');

    return new Math::Vector (
        i => $self->i - $other->i,
        j => $self->j - $other->j,
        k => $self->k - $other->k,
    );
}


# $v *= $scalar;
sub multiply_scalar_assign {
    my ($self, $scalar) = @_;

    $self->set_i($self->i * $scalar);
    $self->set_j($self->j * $scalar);
    $self->set_k($self->k * $scalar);

    return $self;
}


# $v /= $scalar;
sub divide_scalar_assign {
    my ($self, $scalar) = @_;
    carp ('Cannot divide by zero') && return undef unless ($scalar);

    $self->set_i($self->i / $scalar);
    $self->set_j($self->j / $scalar);
    $self->set_k($self->k / $scalar);

    return $self;
}


# $v_new = $v_1 * $scalar;
# $v_new = $v_1->multiply_scalar($scalar);
sub multiply_scalar {
    my ($self, $scalar) = @_;

    return new Math::Vector (
        i => $self->i * $scalar,
        j => $self->j * $scalar,
        k => $self->k * $scalar,
    );
}


# $v_new = $v_1 / $scalar;
# $v_new = $v_1->divide_scalar($scalar);
sub divide_scalar {
    my ($self, $scalar) = @_;
    carp ('Cannot divide by zero') && return undef unless ($scalar);

    return new Math::Vector (
        i => $self->i / $scalar,
        j => $self->j / $scalar,
        k => $self->k / $scalar,
    );
}


# $v_new = $v_1 x $v_2;
# $v_new = $v_1->cross_product($v_2);
sub cross_product {
    my ($self, $other) = @_;
    carp ('Operation only permitted on Math::Vector objects') && retur
+n undef
        unless (ref $other eq 'Math::Vector');

    return new Math::Vector (
        i => $self->j * $other->k - $self->k * $other->j,
        j => $self->k * $other->i - $self->i * $other->k,
        z => $self->i * $other->j - $self->j * $other->i,
    );
}


# $dot_product = $v_1 . $v_2;
# $v_new = $v_1->dot_product($v_2);
sub dot_product {
    my ($self, $other) = @_;
    carp ('Operation only permitted on Math::Vector objects') && retur
+n undef
        unless (ref $other eq 'Math::Vector');

    return $self->i * $other->i + $self->j * $other->j + $self->k * $o
+ther->k;
}


# $magnitude = $v->magnitude();
sub magnitude {
    my $self = shift;

    return sqrt($self->i * $self->i + $self->j * $self->j + $self->k *
+ $self->k);
}


# Converts the vector to its normal
# $v->normalize();
sub normalize {
    my $self = shift;

    my $m = $self->magnitude;
    carp('Cannot normalize a zero vector') && return undef unless ($m)
+;

    $self->set_i( $self->i / $m );
    $self->set_j( $self->j / $m );
    $self->set_k( $self->k / $m );
}


# $v_normalized = $v->normal();
sub normal {
    my $self = shift;

    my $normal = new Math::Vector (i=>$self->i, j=>$self->j, k=>$self-
+>k);
    $normal->normalize();

    return $normal;
}


# $angle_in_radians = $v_1->angle($v_2);
sub angle {
    my ($self, $other) = @_;
    carp ('Operation only permitted on Math::Vector objects') && retur
+n undef
        unless (ref $other eq 'Math::Vector');

    my $self_m  = $self->magnitude;
    my $other_m = $other->magnitude;
    carp('Vectors must be non-zero') && return undef
        unless ($self_m && $other_m);

    return (acos (
        $self->dot_product($other) / ($self_m * $other_m)
    ));
}

sub orthogonal {
    my ($self, $other) = @_;
    carp ('Operation only permitted on Math::Vector objects') && retur
+n undef
        unless (ref $other eq 'Math::Vector');

    carp('Vectors must be non-zero') && return undef
        unless ($self->magnitude && $other->magnitude);

    # v1 is perpendicular to v2 if (v1 . v2 == 0)
    return ($self->dot_product($other) ? 0 : 1);
}

sub i { return $_[0]->{_i}; }
sub j { return $_[0]->{_j}; }
sub k { return $_[0]->{_k}; }

sub set_i {
    my $self = shift;
    $self->{_i} = shift || 0;
}

sub set_j {
    my $self = shift;
    $self->{_j} = shift || 0;
}

sub set_k {
    my $self = shift;
    $self->{_k} = shift || 0;
}

1;
__END__

=head1 NAME

Math::Vector - Object oriented vectors

=head1 SYNOPSIS

  use Math::Vector;
  my $vector = new Math::Vector (i => $i, j => $j, k => $k);

=head1 DESCRIPTION

C<Math::Vector> provides an object oriented approach to working with v
+ectors.

Vectors created with this module are three-dimentional, although
two-dimentional vectors are effectively achieved by omitting the 'k' c
+omponent
(which simply sets it to zero).

It is important to note that this module will carp() if operations inv
+olving
division are attempted on zero values.  A basic understanding of vecto
+rs is
expected, and the programmer should check the magnitude of a vector be
+fore
performing operations such as normalization or finding the angle.

=head2 Conventions used in this document

=over 10

The term 'scalar' is used in the mathmatical sense, and is a number re
+presenting
the amount by which a vector is scaled.

Vector objects will always begin with a 'v'. (Example: $v, $v_1, $v_ne
+w)

=back

=head1 OPERATORS

C<Math::Vector> objects override many default operators.

=item '='

$v_1 = $v_2;

=item '+'

$v_new = $v_1 + $v_2;

=item '+='

$v_1 += $v_2

=item '-'

$v_new = $v_1 - $v_2;

=item '-='

$v_1 -= $v_2;

=item 'x'

$v_cross_product = $v_1 x $v_2;

=item '.'

$dot_product = $v_1 . $v_2;

=item '*'

$v_new = $v_1 * $scalar_multiple

=item '*='

$v_1 *= $scalar_multiple; 

=item '/'

$v_new = $v_1 / $scalar_multiple

=item '/='

$v_1 /= $scalar_multiple; 

=head1 METHODS

=head2 Creation

=over 4

=item new Math::Vector (i => $i, j => $j, k => $k)

Creates a new vector object with components ($i, $j, $k).

=back

=head2 Access

=over 4

=item i

Returns the value of the 'i' (first) component of the vector.

=item j

Returns the value of the 'j' (second) component of the vector.

=item k

Returns the value of the 'k' (third) component of the vector.

=item ijk

Returns a list of the three component values.  This is included for co
+nvenience, and is the same as:

    ($v->i, $v->j, $v->k);

=back

=head2 Direct Modifications

=over 4

=item set_i (I)

Sets the value for the 'i' component of the vector.

    $v->set_i($new_i_value)

=item set_j (J)

Sets the value for the 'j' component of the vector.

    $v->set_j($new_j_value)

=item set_k (K)

Sets the value for the 'k' component of the vector.

    $v->set_k($new_k_value)

=back

=head2 Calculations

=over 4

=item dot_product (VECTOR)

Returns the dot product of two vectors.

    $dot_product = $v_1->dot_product( $v_2 );
    # Note this is equivalent:
    $dot_product = $v_1 . $v_2;

=item cross_product (VECTOR)

Returns the cross product of two vectors as a vector object.  Note: yo
+u cannot
take the cross product of 2d vectors.

    $v_cross = $v_1->cross_product( $v_2 );
    # Note this is equivalent:
    $v_cross = $v_1 x $v_2;

=item magnitude

Returns the magnitude of the vector.

    $magnitude = $v->magnitude;

=item normalize

Normalizes the vector.

    $v->normalize;

=item normal

Returns the normal of a vector as a new vector object.

    $v_normalized = $v->normal;
    # Note the following is similar, but does not preserve $v :
    $v->normalize; $v_normalized = $v;

=item angle (VECTOR)

Returns the angle in radians between two vectors.

    $angle = $v_1->angle( $v_2 );

=item orthogonal (VECTOR)

Returns 1 if vectors are orthogonal (or perpendicular), 0 if they are
not, and undef if either is a zero vector.

    $orthogonal = $v_1->orthogonal( $v_2 );
    if ($orthogonal) {
       # $v_1 and $v_2 are perpendicular
    } elsif (defined $orthogonal) {
       # they are not
    } else {
       # you tried to test a zero vector
    }

=back

=head1 SEE ALSO

Math::Trig,
Carp

=head1 BUGS

Please email the author (E<lt>mark@sporkstorms.orgE<gt>) with any.

=head1 FEEDBACK

Feel free to email the author with any questions, comments, or request
+s.
This module was initially written for personal use, and may be lacking
some desired features.  I (Mark) would be glad to add anything to it
that you might need.

=head1 AUTHOR

Mark Stratman, E<lt>mark@sporkstorms.orgE<gt>

=head1 COPYRIGHT

Copyright (c) 2001, Mark Stratman.  All Rights Reserved.
This modules is free software. It may be used, redistributed and/or
modified under the same terms as Perl itself.

=cut
#
Replies are listed 'Best First'.
Re: Math::Vector - request for comments
by Zaxo (Archbishop) on Dec 17, 2001 at 08:37 UTC

    A few points on implementation:

    1. The sub orthogonal is not very stable. Floating point numbers should be compared in terms of a small tolerance. A proper choice for two vectors could be pseudocoded as magnitude(v1-v2) < $EPS * sqrt(magnitude(v1)**2 + magnitude(v2)**2); $EPS is epsilon, the builtin precision of floats.
    2. Sub angle loses precision for angles near integer multiples of pi. An implementation in terms of atan2 would be better, extracting sin of the angle from the cross product.
    3. You may prefer to simply return the zero vector from normalize, instead of carping out.
    4. You'll find some help in the Math::Trig module. For high performance, take a look at Math::GSL and Math::Pari. If you decide to follow Masem's suggestion to generalize, PDL is another high performance library specializing in arrays of values.
    I implemented some of the sexier parts of this stuff in the snippet Cartesian 3-Vectors.

    Update: Added #4. U2: The pseudocode in #1 is for equality ( == operator). Other comparisons are similarly made.

    After Compline,
    Zaxo

      Floating point numbers should be compared in terms of a small tolerance.
      That is a wonderful point! Many thanks. I have been so used to Perl taking care of the nitty gritty numeric stuff, that I didn't even think of that.

      An implementation in terms of atan2 would be better...
      I had originally used atan2(), simply for the lack of an inverse-cosine function.. until I was pointed to Math::Trig and told to use it ;) So i just haphazardly switched, never considering accuracy issues.

      Thanks a lot for the pointers. They're greatly appreciated =)
Re: Math::Vector - request for comments
by Masem (Monsignor) on Dec 17, 2001 at 07:43 UTC
    One thing that I would consider, but that would require a bit of a rewrite, is to keep "Math::Vector" general to n dimensions, and create "Math::Vector::3D" for the specific 3-dimensional extentions that you have (eg x(), y(), z()). You'd have to add some size checks on the operations that take two vectors, of course.

    -----------------------------------------------------
    Dr. Michael K. Neylon - mneylon-pm@masemware.com || "You've left the lens cap of your mind on again, Pinky" - The Brain
    "I can see my house from here!"
    It's not what you know, but knowing how to find it if you don't know that's important

      Thanks for the input. =)

      I considered that when rewriting bits of this. The main reason it isn't implemented is because personally I had only needed 3d vectors (physics fun =).

      I think I will definitely consider this before making it more widely available.
Re: Math::Vector - request for comments
by Starky (Chaplain) on Dec 17, 2001 at 10:32 UTC
    I would concur with the comment regarding generalizing to arbitrary dimensionality. I don't the the generalization would be all that difficult and the payoff would be worth the effort. In terms of the changes that would be required, other than the implications for your algorithms
    • I think that the arguments to new() would most appropriately be an array ref
    • A size() method would be useful both internally and externally
    • You would have to have some internal exception handling for your overloaded operators in case, for example, someone tried to subtract vectors of different dimension

    My only other comment is that I think the namespace should stay Math::Vector rather than your alternative, Physics::Vector. After all, the concept of a vector is defined precisely in mathematics and only used by physicists. (Athough the mathematical definition does not necessarily concern vectors of numbers, I don't think you'll confuse anyone.) So the designation in the Math::* namespace is entirely appropriate.

    Excellent contribution! I look forward to using it.

Re: Math::Vector - request for comments
by chromatic (Archbishop) on Dec 18, 2001 at 01:58 UTC
    Your last six methods could be consolidated:
    no strict 'refs'; foreach my $axis ('i', 'j', 'k') { *{ $axis } = sub { my $self = shift; if (@_) { $self->{"_$axis"} = shift; } return $self->{"_$axis"}; } }
    If you're willing to limit this to Perl 5.6 and newer, you could also look into lvaluable subroutines. Class::Accessor also handles much of the same thing for you nearly automagically.