### Symbolic calculations with operator overload

by ambrus (Abbot)
 on Oct 09, 2003 at 06:21 UTC ( #297829=note: print w/replies, xml ) Need Help??

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
+apart
# from easy usage.
#
#

use strict;

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

{
package Polinom;

'-',    "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];
}
}

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;

'*',    "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__

Create A New User
Node Status?
node history
Node Type: note [id://297829]
help
Chatterbox?
 [marinersk]: I've often marvelled over the past few decades at just how little it takes to mess me up with this spacial dependence thing. I wonder if I'm mildly autistic or something. [Eily]: can't the nodelet be moved though? Maybe you could put one that doesn't change first ("Sections" or "Find Nodes" for example) [Eily]: "Other Users" seems like a poor choice :P

How do I use this? | Other CB clients
Other Users?
Others pondering the Monastery: (8)
As of 2017-05-29 14:05 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
My favorite model of computation is ...

Results (192 votes). Check out past polls.