# CVQM - lightweight Complex numbers, Vectors, Quaternions, and 4x4-Matrices

package CVQM;
our $VERSION = 0.8;

=head1 NAME

CVQM - lightweight Complex numbers, Vectors, Quaternions, and 4x4-Matrices

=head1 SYNOPSIS

    use CVQM;
    my $c = C[1,2];         # 1 + 2i
    my $v = V[1,2,3];       # 1x + 2y + 3z
    my $v = V[1,2,3,4];     # 1x + 2y + 3z + 4w
    my $m = M::identity;    # M[ 1,0,0,0 => 0,1,0,0 => 0,0,1,0 => 0,0,0,1 ]
    my $q = Q(4,[1,2,3]);   # rotate 4 about axis 1,2,3
    my $q = Q[1,2,3,4];     # 4 + 1i + 2j + 3k

=head1 ABSTRACT

This includes four distinct packages: C<C>, C<V>, C<Q>, and C<M>.  These
packages offer the lightest possible pure-Perl implementation of the
typical matrix and vector math used in everyday 3D programming.

C<C> offers basic Complex number capability often used in fractal
geometry.  C<V> offers general Vector support (usually 3-scalar or
4-scalar).  C<Q> offers Quaternions, which are special 4-scalar values
useful in physical modeling.  C<M> offers basic 4x4 Matrix operations.

All three types are implemented as blessed one-dimensional array
references, making it very simple and lightweight to access the scalars
within, or to dump the scalars into argument lists of other routines,
such as C<SDL::OpenGL> routines or other math libraries.  Operator
overloads for many operations make geometrical tasks more reader
friendly.

These packages prefer compactness and usability and pure Perl
implementation, over any raw native power for high intensity computing.

=cut

1;#--------------------------------------------------------------------------

package C;
use strict;
use Exporter;
BEGIN { our @ISA = qw(Exporter); our @EXPORT = qw(C); };

=head1 C - Complex Numbers

    $c = C($original);
    $c = C($x,$yi);         # x + y*i
    $c = C[$x,$yi];         # x + y*i

    ($x,$yi) = @$c;         # extract both scalars
    $c->[0] = $x;
    $yi = $c->[1];          # direct access

    $c = $ca + $cb;         # complex numbers are not ordered

A Vector is a simple array of scalar numerical values.  This C<V> package
does not care how many numbers are included, but most users will use 3-
or 4-dimensional vectors.  Once created, any C<Q> object can be unrolled
back out to its four scalars C<($x, $y, $z, $w)> or the scalars within
can be manipulated directly.

=cut

sub C
{
    return bless [ @{$_[0]} ] if ref $_[0];
    return bless [ $_[0], $_[1] ];
}

use overload
#    '+' => \&add,
#    '-' => \&subtract,
#    '*' => \&multiply,
#    '/' => \&divide,
    '""' => \&stringify;

sub stringify { "C[ " . join(', ', @{$_[0]}) . " ]" }

1;#--------------------------------------------------------------------------

package V;
use strict;
use Exporter;
our @ISA = qw(Exporter);
our @EXPORT = qw(V);
our @EXPORT_OK = qw(magnitude normal normalize urand);
our %EXPORT_TAGS = ( 'all' => [ @EXPORT, @EXPORT_OK ] );

=head1 V - Vectors

    $v = V($original);
    $v = V($x,$y,$z,...);   # 3d or 4d vector most commonly
    $v = V[$x,$y,$z,...];

    ($x,$y,$z,...) = @$v;   # extract all scalars
    $v->[0] = $x;
    $y = $v->[1];           # direct access

    $v = $va + $c;
    $v = $va + $vb;         # real vectors are ordered
    $v = $va * $c;
    $v = $va * $vb;         # real vectors are commutative

A Vector is a simple array of scalar numerical values.  This C<V> package
does not care how many numbers are included, but most users will use 3-
or 4-dimensional vectors.  Once created, any C<Q> object can be unrolled
back out to its four scalars C<($x, $y, $z, $w)> or the scalars within
can be manipulated directly.

=cut

sub V
{
    return bless [ @{$_[0]} ] if ref $_[0];
    return bless [ @_ ];
}

use overload
    '+' => \&add,
    '-' => \&subtract,
    '*' => \&multiply,
    '/' => \&divide,
    'x' => \&cross,
    '.' => \&_dot_overload,
    '""' => \&stringify;

sub stringify { "V[ " . join(', ', @{$_[0]}) . " ]" }

sub add
{
    my @p = @{+shift};
    my $q = shift;
    if (ref($q)) { $p[$_] += $q->[$_] for (0..$#p); }
    else { $_ += $q for @p; }
    return bless [ @p ];
}

sub subtract
{
    my @p = @{+shift};
    my $q = shift;
    if (ref($q)) { $p[$_] -= $q->[$_] for (0..$#p); }
    else { $_ -= $q for @p; }
    return bless [ @p ];
}

sub multiply # vector * scalar
{
    my @p = @{+shift};
    my $q = shift; $_ *= $q for @p;
    return bless [ @p ];
}

sub divide # vector / scalar
{
    my @p = @{+shift};
    my $q = 1 / (+shift); $_ *= $q for @p;
    return bless [ @p ];
}

sub dot
{
    my ($p, $q) = @_;
    my $s = 0; $s += ($p->[$_] * $q->[$_]) for (0..$#{$p});
    return $s;
}

sub cross # specific to 3D
{
    my ($p, $q) = @_;
    return bless [ $p->[1] * $q->[2] - $p->[2] * $q->[1],
                   $p->[2] * $q->[0] - $p->[0] * $q->[2],
                   $p->[0] * $q->[1] - $p->[1] * $q->[0] ];
}

sub magnitude
{
    my $p = shift;
    my $s = 0; $s += ($_ * $_) for @$p;
    return sqrt( $s );
}

sub normalize { multiply($_[0], 1 / magnitude($_[0])) }

sub normal # specific to 3D
{
    return normal([0,0,0], $_[0], $_[1]) if @_ == 2;
    my ($p, $q, $r) = @_;
    my @Q = @$q; $Q[$_] -= $p->[$_] for (0..$#Q);
    my @R = @$r; $R[$_] -= $p->[$_] for (0..$#R);
    return normalize(cross(\@Q, \@R));
}

sub centroid
{
    my $c = @_;
    my $v = V(shift());
    $v += shift() while @_;
    return $v / $c;
}

sub urand { normalize( [ rand()*2 - 1, rand()*2 - 1, rand()*2 - 1 ] ) }

sub _dot_overload
{
    my ($obj, $arg, $flip) = @_;
    return dot($obj, $arg)
        if $arg and UNIVERSAL::isa($arg, __PACKAGE__);
    return $arg . stringify($obj) if $flip;
    return stringify($obj) . $arg;
}

if (0)
{
    my $x = V[ 1, 2, 3 ];
    my $y = V[ 2, 3, 4 ];
    my $z;

    for ( '$z = $x + $y',
          '$z = $y - $x',
          '$z = $x - $y',
          '$z = $x + 8',
          '$z = $x - 8',
          '$z = $x * 4',
          '$z = $x . $y',
          '$z = $x x $y',
          '$z = magnitude($x)',
          '$z = normalize($x)',
          '$z = V[0,0,0]',
          '$z = normal($x, $y, $z)',
          '$z = normal($z, $y, $x)',
          '$z = centroid($x, $y)',
          )
    {
        eval $_;
        print $_, " \t|   ", eval("qq{$_}"), $/;
    }
}

1;

1;#--------------------------------------------------------------------------

package Q;
use strict;
use Exporter;
BEGIN { our @ISA = qw(Exporter); our @EXPORT = qw(Q); };

=head1 Q - Quaternions

    $q = Q($original);
    $q = Q($w, [$x,$y,$z]);
    $q = Q($x,$y,$z,$w);
    $q = Q[$x,$y,$z,$w];

    ($x,$y,$z,$w) = @$q;

    $q = $qa * $qb;

A Quaternion is a special four-dimensional vector which is useful for
describing rotations and angular momentum in 3D space.  Once created, any
C<Q> object can be unrolled back out to its four scalars C<($x, $y, $z,
$w)> or the scalars within can be manipulated directly.  Quaternions
resemble normal 4-scalar vectors but differ in some operations; they are
the 4D progression of Complex numbers.  (The C<CVQM> package does not
implement the rarely used but similar progression: the 8-dimensional
Octonian.)

=cut

sub Q
{
    return bless [ @{$_[1]}, $_[0] ]
        if ref $_[1] and not ref $_[0];
    return bless [ @{$_[0]} ] if ref $_[0];
    return bless [ @_ ];
}

use overload
#    '+' => \&add,
#    '-' => \&subtract,
#    '*' => \&multiply,
#    '/' => \&divide,
    '""' => \&stringify;

sub stringify { "Q[ " . join(', ', @{$_[0]}) . " ]" }

sub conjugate($)
{
    my @q = @{+shift};
    $q[$_] = -$q[$_] for 0..2;
    return bless [ @q ];
}

sub normalize($) { V::normalize(@{+shift}) }

sub magnitude($) { V::magnitude(@{+shift}) }

sub multiply($$)
{
    my ($x1,$y1,$z1,$w1) = @{+shift};
    my ($x2,$y2,$z2,$w2) = @{+shift};
    return bless
        [ $w1*$x2 + $x1*$w2 + $y1*$z2 - $z1*$y2,
          $w1*$y2 + $y1*$w2 + $z1*$x2 - $x1*$z2,
          $w1*$z2 + $z1*$w2 + $x1*$y2 - $y1*$x2,
          $w1*$w2 - $x1*$x2 - $y1*$y2 - $z1*$z2 ];
}

1;
#----------------------------------------------------------------------------

package M;
use strict;
use Exporter;
BEGIN { our @ISA = qw(Exporter); our @EXPORT = qw(M); };

=head1 M - 4x4 Matrices

    $mm = M($original);
    $mm = M( $a, $b, $c, $d,
             $e, $f, $g, $h,
             $i, $j, $k, $l,
             $m, $n, $o, $p );
    $mm = M[ @s16 ];

    ( $a, $b, $c, $d,
      $e, $f, $g, $h,
      $i, $j, $k, $l,
      $m, $n, $o, $p ) = @$mm;

    $mm = $mm * $nn;

A matrix is a rectangular grid of scalar numbers; this package only deals
with one certain kind of matrix C<M> which has four rows and four
columns.  The 4x4 matrix is very useful for describing almost any linear
transformation in 3D space.

To keep the implementation simple, an C<M> object is kept as a single
16-element array.  It can be unrolled back out to its sixteen scalars
C<($x, $y, $z, $w)> or the scalars within can be manipulated directly.
While this C<M> package does offer pure Perl matrix multiplication
features, most C<OpenGL> applications probably want to rely on native
processors to do actual transformation work and use C<M> for a few basic
bookkeeping tasks.

=cut

sub M
{
    return bless [ (0)x($_[0]*$_[1]) ] if @_ == 2 and not ref($_[0]);
    return bless [ (0)x($_[0]*$_[0]) ] if @_ == 1 and not ref($_[0]);
    return bless [ @{$_[0]} ] if ref($_[0]);
    return bless [ @_ ];
}

use overload
#    '+' => \&add,
#    '-' => \&subtract,
#    '*' => \&multiply,
#    '/' => \&divide,
    '""' => \&stringify;

sub stringify
{
    return join('',
                "M[ ",
                join(',', @{$_[0]}[0..3]), " => ",
                join(',', @{$_[0]}[4..7]), " => ",
                join(',', @{$_[0]}[8..11]), " => ",
                join(',', @{$_[0]}[12..15]),
                " ]");
}

sub identity { bless [ 1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1 ] }

sub multiply_mm
{
    my ($a,$b,$c,$d,$e,$f,$g,$h,$i,$j,$k,$l,$m,$n,$o,$p) = @{+shift};
    my ($A,$B,$C,$D,$E,$F,$G,$H,$I,$J,$K,$L,$M,$N,$O,$P) = @{+shift};
    return bless 
        [ $a*$A+$b*$E+$c*$I+$d*$M,
          $a*$B+$b*$F+$c*$J+$d*$N,
          $a*$C+$b*$G+$c*$K+$d*$O,
          $a*$D+$b*$H+$c*$L+$d*$P,
          $e*$A+$f*$E+$g*$I+$h*$M,
          $e*$B+$f*$F+$g*$J+$h*$N,
          $e*$C+$f*$G+$g*$K+$h*$O,
          $e*$D+$f*$H+$g*$L+$h*$P,
          $i*$A+$j*$E+$k*$I+$l*$M,
          $i*$B+$j*$F+$k*$J+$l*$N,
          $i*$C+$j*$G+$k*$K+$l*$O,
          $i*$D+$j*$H+$k*$L+$l*$P,
          $m*$A+$n*$E+$o*$I+$p*$M,
          $m*$B+$n*$F+$o*$J+$p*$N,
          $m*$C+$n*$H+$o*$K+$p*$O,
          $m*$D+$n*$H+$o*$L+$p*$P ];
}

sub multiply_vm
{
    my ($x,$y,$z,$w) = @{+shift};
    my ($A,$B,$C,$D,$E,$F,$G,$H,$I,$J,$K,$L,$M,$N,$O,$P) = @{+shift};
    return 
        V[ $x*$A+$y*$E+$z*$I+$w*$M,
           $x*$B+$y*$F+$z*$J+$w*$N,
           $x*$C+$y*$G+$z*$K+$w*$O,
           $x*$D+$y*$H+$z*$L+$w*$P ];
}

sub multiply_mv
{
    my ($a,$b,$c,$d,$e,$f,$g,$h,$i,$j,$k,$l,$m,$n,$o,$p) = @{+shift};
    my ($X,$Y,$Z,$W) = @{+shift};
    return 
        V[ $a*$X+$b*$Y+$c*$Z+$d*$W,
           $e*$X+$f*$Y+$g*$Z+$h*$W,
           $i*$X+$j*$Y+$k*$Z+$l*$W,
           $m*$X+$n*$Y+$o*$Z+$p*$W ];
}

if (1)
{
    my $mm = M::identity;
    print $mm, $/;
}

1;#--------------------------------------------------------------------------
__END__

=head1 SEE ALSO

More generalized matrix and vector routines are implemented in the CPAN
modules C<Math::MatrixReal> and C<Math::VectorReal>.  These offer support
for grabbing vector slices from matrices, doing various pivots,
decompositions, eigensystems, and many other capabilities, at the expense
of compactness.

Some basic quick-reference resources on these concepts and operations can
be found at many places on the web, but this one served very helpful for
quick implementations.  http://www.j3d.org/matrix_faq/matrfaq_latest.html

=head1 AUTHOR

Ed Halley <C<< ed@halley.cc >>>

=head1 COPYRIGHT AND LICENSE

Copyright 2003-2004 by Ed Halley

This library is free software; you can redistribute it and/or modify it
under the same terms as Perl itself.

=cut