Symbolic calculations with operator overload

[ambrus]

Here's some code I wrote last summer (and changerd a bit now).

#!/usr/bin/perl -w
#
# We will do operations with multi-variable real polinomials and 3d ve
+ctors 
# of these.  These are slow algorithms, but might be good examples for
# object orientation ond operator overloading in Perl.  It is also not
# commented well.  Many functions are defined but not used, altough I 
+have 
# tested them.  It is not difficult to add some calculations that use
# these functions.
#

# Note that the Vector class knows nothing about the Polinom class, st
+ill
# it can work with vectors of polinomials as well as vectors of number
+s
# transparently.  This is the real advantage of operator overloading, 
+apart
# from easy usage.
#
#

# head -------------------------------------------
use strict;

# Polinom ----------------------------------------

{
package Polinom;


use overload
    '+',    "add",
    '-',    "sub",
    '*',    "mul",
    '**',   "power",
    '""',   "str",
    '0+',    "num",
    'bool',    "bool",
    ;
    

# The internal structure of Polinom objects is like this:
# 6*x^2+3*x*y+8*x-4 ==> {"x*x", 6, "x*y", 3, "x", 8, "", -4}
# The keys of the hash are variables separated with "*", sorted alpha-
# betically, so "y*x" cannot be a key in such a hash.  The values are
# the coeffcients, those are simply a Perl number.

# xtract: this internal function does the magic of converting a number
# to a polinomial when used in a polinomial operation
sub xtract {
    my $a= shift;
    if ( ref($a) ) {
         return %$a
         }
    else {
        return $a ? ("", $a) : ();
        }
    }


# constant (useless, because a number can always be used instd a polin
+om)
sub const {
    return bless {"", $_[1]}, $_[0];
    }


# variable constructor, returns more Polinoms if given more strings
sub var {
    local $_;
    my $t= shift;
    for (@_)
        { m|^[a-zA-Z_][a-zA-Z0-9_\-.\[\]]*$|
            or die "wrong variable name for $t->var"; }
    my @v= ();
    for (@_)
        { push @v, bless {$_, 1}, $t }
    if (wantarray)
        { return @v; }
    else { 
        @_==1 or die 'cannot return list to scalar context';
        return $v[0]; 
        }
    }


# add method
sub add {
    my $t= ref $_[0];
    my %a= %{shift()};
    my %b= xtract(shift());
    for my $k (keys %b) {
        exists $a{$k} or $a{$k}= 0;
        $a{$k}+= $b{$k};
        $a{$k} or delete $a{$k};
        };
    return bless {%a}, $t;
    };


# subtract method
sub sub {
    my $t= ref $_[0];
    my $r= $_[2];
    my %a= xtract($_[$r?1:0]);
    my %b= xtract($_[$r?0:1]);
    for my $k (keys %b) {
        exists $a{$k} or $a{$k}= 0;
        $a{$k}-= $b{$k};
        $a{$k} or delete $a{$k};
        };
    return bless {%a}, $t;
    };


# monom-multiply (internal)
sub mulmonom ($$) { # <---- ELIMINAL @c, STRING HASZNALJ HLYTTE
    my @a= split /\*/, shift;
    my @b= split /\*/, shift;
    my @c= ();
    while ( @a && @b ) 
        { push @c, ( $a[0] le $b[0] ? shift @a : shift @b ); };
    @a ? push (@c, @a) : push (@c, @b);
    return join "*", @c;
    }

# multiply method
# handles mul with number as a separate case
sub mul {
    my $t= ref $_[0];
    my %a= %{shift()};
    my $b= shift;
    if (!ref $b) {
        $b or return bless {}, $t;
        for my $v (@a{keys %a}) 
            { $v*= $b; };
        return bless {%a}, $t;
        }
    else {    
        my %b= %$b;
        my %c= ();
        my $k;
        for my $a (keys %a) {
            for my $b (keys %b) {
                $k= mulmonom($a,$b);
                exists $c{$k} or $c{$k}= 0;
                $c{$k}+= $a{$a}*$b{$b};
                $c{$k} or delete $c{$k};
                }
            }
        return bless {%c}, $t;
        };
    };


# multiply method with debug
sub mul_g {
    my $t= ref $_[0];
    my %a= %{shift()};
    my $b= shift;
    if (!ref $b) {
        $b or return bless {}, $t;
        for my $v (@a{keys %a}) 
            { $v*= $b; };
        return bless {%a}, $t;
        }
    else {    
        my %b= %$b;
        my %c= ();
        my $k;
        my $ga= keys(%a);
        print "[poli-mul ".$ga." ".keys(%b).": ";
        for my $a (keys %a) {
            print $ga--." ";
            for my $b (keys %b) {
                $k= mulmonom($a,$b);
                exists $c{$k} or $c{$k}= 0;
                $c{$k}+= $a{$a}*$b{$b};
                $c{$k} or delete $c{$k};
                }
            }
        print "]\n";
        return bless {%c}, $t;
        };
    };


# integer power method
sub power {
    my $t= ref $_[0];
    my ($a, $n, $rev)= @_;
    $rev and die "Wrong op: scalar**Polinom";
    $n==int($n) or die "Wrong op: Polinom**fraction";
    my $x= 1;
    my $c= 1;
    while ($n) {
        $x&$n and $c*= $a, $n-= $x;
        $a*= $a;
        $x+= $x;
        }
    $c;
    };
    

# with: method for applying polinomials
# Usage: $polinomial->with (variable=>value, variable=>value, ...);
# Where variable can be string or polinomial consisting of just a vari
+able;
# values can be numbers or polinomials (but numbers are faster of cour
+se).
# Substitutions are not done simultanously, so do not try to do
# $poli->with ($x->$y, $y->x);

# _with_{n,z,p}: internal functions for setting a variable to {number,
# zero,polinomial} resp.

sub _with_n ($$$) { # <--- TEST
    local $_;
    my %a= %{shift()};
    my $v= shift;
    my $w= shift;
    my %r= ();
    for my $k (keys %a) {
        my @c= split /\*/, $k;
        my $n= "";
        my $y= $a{$k};
        for (@c) {
            if ($v eq $_) { $y*= $w }
            else { $n.= "*$_" }
            }
        $n=~ s/^\*//;
        exists $r{$n} or $r{$n}= 0;
        $r{$n}+= $y;
        }
    return %r;
    }
    
sub _with_z ($$) {
    my %a= %{shift()};
    my $v= shift;
    my $vre= qr!(^|\*)$v(\*|$)!;
    for my $k (keys %a) 
        { $k=~$vre and delete $a{$k}; }
    return %a;
    }
    
sub _with_p ($$$$) {
    local $_;
    my %a=%{shift()};
    my $v= shift;
    my $w= shift;
    my $t= shift;
    my $r= 0;
    my $x;
    for my $k (keys %a) {
        my @c= split /\*/, $k;
        $x= $a{$k};
        for (@c) {
            if ($v eq $_) { $x*= $w }
            else { $x*= bless({$_,1},$t) }
            }
        $r+= $x;
        }
    return %$r;    
    }

# now the definition of with:
sub with {
    my $o= shift;
    my %r= %$o;
    while (@_>=2) {
        my $v= shift;
        my $w= shift;
        if (ref $v) {
            my @c= %$v; 
            @c==2 || !$c[0] || $c[0]=~/\*/ || $c[1]!=1
                or die ref($o).'->with called with wrong polinom as va
+r';
            $v= $c[0];
            };
        if (ref $w) {%r= _with_p \%r, $v, $w, ref $o}
        elsif ($w) {%r= _with_n \%r, $v, $w}
        else {%r= _with_z \%r, $v}
        }
    @_ and die 'odd-sized hash to '.ref $o.'->with';
    return bless {%r}, ref $o;
    }



# overloaded stringify method -- does not print powers of variables.
# A parital solution would be to filter the output through
# perl -wpe's!([a-z])((\*\1)+)!"$1^".(length($2)/2+1)!ge'
# which transforms 3*x*x*x-4*x*x to 3*x^3-4*x^2.

sub str {
    my %a= xtract shift;
    my $r= "";
    my $v;
    for my $k (sort keys %a) {
        $v= $a{$k};
        $r.= $v<0?"-":$r?"+":"";
        if ($k) { $r.= abs($v)."*" unless abs($v)==1; }
            else {$r.= abs($v); }
        $r.= $k;
        };
    $r or $r= "0";
    return $r;
    };


# convert to number method, returns undef if polinom is not constant
sub num {
    my @a= %{shift()};
    if (@a>2)
        { return undef; }
    if (@a==0)
        { return 0; }
    if ($a[0])
        { return undef; }
    return $a[1];
    }


# convert to boolean method
sub bool {
    return scalar keys(%{shift()});
    }


}



# Vector ------------------------------------------
# Vector objects are vectors of 3 perl scalars.  See note at top.

{
package Vector;


use overload
    '*',    "mul",
    'x',    "cross",
    '""',    "str",
    ;

# the xi+yj+zk constructor
sub new {
    @_==4 or die "a Vector should have 3 coordinates";
    my $o= [@_[1..3]];
    return bless $o, $_[0];
    }

# the 0-vector constructor
sub zero {
    my $o= [0,0,0];
    return bless $o, $_[0];
    }

# multiply method, overloads "*" operator
# muliplies with scalar or calculates dot product as appropriate
sub mul {
    local $_;
    my $a= shift;
    my $b= shift;
    if ( ref($b) && $b->isa("Vector") ) {
        my $x= $$a[0]*$$b[0]+$$a[1]*$$b[1]+$$a[2]*$$b[2];
        return $x;
        }
    else {
        my @x= ($$a[0]*$b, $$a[1]*$b, $$a[2]*$b);
        return bless [@x], ref $a;
        }
    }

# multiply method with debug (to see how much you have to wait)
sub mul_g {
    local $_;
    my $a= shift;
    my $b= shift;
    if ( ref($b) && $b->isa("Vector") ) {
        print "[vector-inmul ";
        my $x= 0;
        $x+= $$a[0]->mul_g($$b[0]);
        $x+= $$a[1]->mul_g($$b[1]);
        $x+= $$a[2]->mul_g($$b[2]);
        print "]\n";
        return $x;
        }
    else {
        my @x= ($$a[0]*$b, $$a[1]*$b, $$a[2]*$b);
        return bless [@x], ref $a;
        }
    }

# cross product method, overloads "x"
sub cross {
    local $_;
    my $a= shift;
    my $b= shift;
    my @x= (
        $$a[1]*$$b[2]-$$a[2]*$$b[1],
        $$a[2]*$$b[0]-$$a[0]*$$b[2],
        $$a[0]*$$b[1]-$$a[1]*$$b[0]);
    return bless [@x], ref $a;
    }

# stringify method
sub str {
    my @o= @{shift()};
    return "(".$o[0].",".$o[1].",".$o[2].")";    
    }


}


# Test --------------------------------------------

{ # empty-subclass test
package Poland;
use vars qw{@ISA};
@ISA= "Polinom";
}

{ # empty-subclass test
package Hector;
use vars qw{@ISA};
@ISA= "Vector";
}


no strict "vars";


# This test will calculate Chebyshev-polinomials defined by the recurs
+ion
# t(0,x)= 1; t(1,x)= x;
# t(n+1,x)= 2*x*t(n,x)-t(n-1,x)
print "Chebyshev-polinomials\n";

# define the variables x
$x= var Polinom qw!x!;

($n, $u, $t)= (1, 1, $x);
while ($n<12) {
print "t($n,$x)= ", $t, ";\n";
print "2*t($n,$x/2)= ", 2*$t->with ($x, 1/2*$x), ";\n";
print "t($n,-1)= ", $t->with ($x, -1), ";\n";
print "t($n,1/2)= ", $t->with ($x, 1/2), ";\n";
($n, $u, $t)= ($n+1, $t, 2*$x*$t-$u);
};

print "\n";

# This test will verify the Pappos-Pascal-theorem. The final result
# should be zero.  This test is quite slow.

# variables
$a= Hector->new(Poland->var(qw(a1 a2 a3)));
$b= Hector->new(Poland->var(qw(b1 b2 b3)));
$c= Hector->new(Poland->var(qw(c1 c2 c3)));
$e= Hector->new(Poland->var(qw(d1 d2 d3)));
$f= Hector->new(Poland->var(qw(e1 e2 e3)));

$|= 1;
print "Verifying the Pappos-Pascal theorem\n";
$v=    (($c x $e)x($a x $f)) x (($a x $e)x($c x $f));
print "v= ";
# if you see the equation sign, the calculation's been done, 
# only stringification remains
print $v, "\n";
$w=    ( (($a x $e)x($b x $f)) x (($b x $e)x($a x $f)) )
      x ( (($b x $e)x($c x $f)) x (($c x $e)x($b x $f)) );
print "w= ";
print $w, "\n";
$r= $v->mul_g($w);
print "v.w=";
print $r, "\n". "The final result above should be zero.\n";


__END__