File Coverage

blib/lib/Math/Polynomial/Multivariate.pm

Criterion	Covered	Total	%
statement	311	311	100.0
branch	70	70	100.0
condition	15	15	100.0
subroutine	49	49	100.0
pod	20	20	100.0
total	465	465	100.0

# Copyright (c) 2011-2017 Martin Becker.  All rights reserved.
# This package is free software; you can redistribute it and/or modify it
# under the same terms as Perl itself.

package Math::Polynomial::Multivariate;

use 5.008;
use strict;
use warnings;
use Carp qw(croak);

use overload (
    q[neg]      => \&_neg,
    q[+]        => \&_add,
    q[-]        => \&_sub,
    q[*]        => \&_mul,
    q[**]       => \&_pow,
    q[==]       => \&_eq,
    q[!=]       => \&_ne,
    q[!]        => 'is_null',
    q[bool]     => 'is_not_null',
    q[""]       => 'as_string',
    q[fallback] => undef,
);

# ----- object definition -----

# Math::Polynomial::Multivariate=ARRAY(...)

# .......... index ..........   # .......... value ..........
use constant F_TERMS  => 0;     # monomial terms hashref
use constant F_ZERO   => 1;     # zero element of coefficient space
use constant F_ONE    => 2;     # unit element of coefficient space
use constant NFIELDS  => 3;

# a term value is an anonymous arrayref:
# .......... index ..........   # .......... value ..........
use constant T_COEFF  => 0;     # coefficient
use constant T_VARS   => 1;     # variables hashref, mapping names to exponents

use constant K_SEP    => chr 0;
use constant K_SEP_RE => qr/\0/;

use constant _MINUS_INFINITY => - (~0) ** (~0);

# ----- class data -----

our $VERSION = '0.006';

# ----- private subroutines -----

# create a "clean" copy of a variables hashref, without zero exponents
sub _vclean {
    my %vars = %{$_[0]};
    delete @vars{grep {!$vars{$_}} keys %vars};
    die "hash of integer exponents expected"
        if grep {!eval { $_ == int abs $_ }} values %vars;
    return \%vars;
}

# map a clean variables hashref to a key
sub _key {
    my ($vars) = @_;
    return
        join K_SEP,
        map { ($_, chr $vars->{$_}) }
        reverse sort keys %{$vars};
}

# map an arbitrary variables hashref to a key
sub _keyvc {
    my ($vars) = @_;
    return
        join K_SEP,
        map { ($_, chr $vars->{$_}) }
        reverse sort
        grep { $vars->{$_} }
        keys %{$vars};
}

# check if an arbitrary operand is a multivariate polynomial

sub _isa_pol ($) { ref($_[0]) && eval { $_[0]->isa(__PACKAGE__) } }

# ----- private methods -----

# extract coefficient space information from an invocant and some coefficient
sub _space {
    my ($this, $const) = @_;
    my $class = ref $this;
    my ($zero, $one);
    if ($class) {
        ($zero, $one) = @{$this}[F_ZERO, F_ONE];
    }
    else {
        $class = $this;
        $zero = $const - $const;
        $one  = $const ** 0;
    }
    return ($class, $zero, $one);
}

# create an independent duplicate of an object
sub _clone {
    my ($this) = @_;
    my $terms = $this->[F_TERMS];
    return
        bless [
            { map {($_ => [@{$terms->{$_}}])} keys %{$terms} },
            @{$this}[F_ZERO, F_ONE]
        ], ref $this;
}

# multiply a polynomial by variables with exponents
sub _mul_vars {
    my ($this, $vars) = @_;
    return
        bless [
            {
                map {
                    my %tvars = %{$_->[T_VARS]};
                    while (my ($var, $exp) = each %{$vars}) {
                        $tvars{$var} += $exp;
                    }
                    _key(\%tvars) => [$_->[T_COEFF], \%tvars]
                } values %{$this->[F_TERMS]}
            },
            @{$this}[F_ZERO, F_ONE]
        ], ref $this;
}

# overloaded operators

sub _neg {
    my ($this) = @_;
    my $result = $this->_clone;
    foreach my $term (values %{$result->[F_TERMS]}) {
        $term->[T_COEFF] = -$term->[T_COEFF];
    }
    return $result;
}

sub _add {
    my ($this, $that) = @_;
    $that = $this->const($that) if !_isa_pol($that);
    my $result = $this->_clone;
    my $rterms = $result->[F_TERMS];
    my $zero   = $this->[F_ZERO];
    while (my ($key, $term) = each %{$that->[F_TERMS]}) {
        if (exists $rterms->{$key}) {
            if ($zero == ($rterms->{$key}->[T_COEFF] += $term->[T_COEFF])) {
                delete $rterms->{$key};
            }
        }
        else {
            $rterms->{$key} = [@{$term}];
        }
    }
    return $result;
}

sub _sub {
    my ($this, $that, $swap) = @_;
    $that = $this->const($that) if !_isa_pol($that);
    ($this, $that) = ($that, $this) if $swap;
    my $result = $this->_clone;
    my $rterms = $result->[F_TERMS];
    my $zero   = $this->[F_ZERO];
    while (my ($key, $term) = each %{$that->[F_TERMS]}) {
        if (exists $rterms->{$key}) {
            if ($zero == ($rterms->{$key}->[T_COEFF] -= $term->[T_COEFF])) {
                delete $rterms->{$key};
            }
        }
        else {
            my $rterm = $rterms->{$key} = [@{$term}];
            $rterm->[T_COEFF] = -$rterm->[T_COEFF];
        }
    }
    return $result;
}

sub _mul {
    my ($this, $that) = @_;
    $that = $this->const($that) if !_isa_pol($that);
    my $result = $this->null;
    my $bterms = $that->[F_TERMS];
    my $rterms = $result->[F_TERMS];
    my $zero   = $this->[F_ZERO];
    foreach my $aterm (values %{$this->[F_TERMS]}) {
        my @avars = %{$aterm->[T_VARS]};
        foreach my $bterm (values %{$bterms}) {
            my %vars = @avars;
            my $const = $aterm->[T_COEFF] * $bterm->[T_COEFF];
            my $bvars = $bterm->[T_VARS];
            foreach my $bv (keys %{$bvars}) {
                $vars{$bv} += $bvars->{$bv};
            }
            my $key = _key(\%vars);
            if (exists $rterms->{$key}) {
                $rterms->{$key}->[T_COEFF] += $const;
            }
            else {
                $rterms->{$key} = [$const, \%vars];
            }
        }
    }
    delete @{$rterms}{
        grep {$zero == $rterms->{$_}->[T_COEFF]} keys %{$rterms}
    };
    return $result;
}

sub _pow {
    my ($this, $exp, $swap) = @_;
    croak 'illegal exponent' if $swap || $exp != int abs $exp;
    return $this->const($this->[F_ONE]) if !$exp;
    my $result = $this->_clone;
    while (--$exp > 0) {
        $result *= $this;
    }
    return $result;
}

sub _eq_ne {
    my ($this, $that, $eq) = @_;
    $that = $this->const($that) if !_isa_pol($that);
    my $aterms = $this->[F_TERMS];
    my $bterms = $that->[F_TERMS];
    return !$eq if keys(%{$aterms}) != keys(%{$bterms});
    while (my ($key, $term) = each %{$aterms}) {
        return !$eq
            if !exists($bterms->{$key}) ||
                $term->[T_COEFF] != $bterms->{$key}->[T_COEFF];
    }
    return $eq;
}

sub _eq {
    my ($this, $that) = @_;
    return $this->_eq_ne($that, !0);
}

sub _ne {
    my ($this, $that) = @_;
    return $this->_eq_ne($that, !1);
}

# ----- public methods -----

# constructors

sub null {
    my ($this) = @_;
    my ($class, $zero, $one) = $this->_space(1);
    return bless [{}, $zero, $one], $class;
}

sub const {
    my ($this, $value) = @_;
    my ($class, $zero, $one) = $this->_space($value);
    my %terms = $zero == $value? (): (q[] => [$value, {}]);
    return bless [\%terms, $zero, $one], $class;
}

sub var {
    my ($this, $varname) = @_;
    my ($class, $zero, $one) = $this->_space(1);
    my $vars = { $varname => 1 };
    my $key = _key($vars);
    my %terms = ($key => [$one, $vars]);
    return bless [\%terms, $zero, $one], $class;
}

sub monomial {
    my ($this, $const, $vars) = @_;
    my ($class, $zero, $one) = $this->_space($const);
    $vars = _vclean($vars);
    my %terms = $zero == $const? (): (_key($vars) => [$const, $vars]);
    return bless[\%terms, $zero, $one], $class;
}

sub subst {
    my ($this, $varname, $that) = @_;
    my @exps = (-1, $this->exponents_of($varname));
    return $this->null if 1 == @exps;
    my $exp    = pop @exps;
    my $result = $this->factor_of($varname, $exp);
    my $nexp   = pop @exps;
    while ($exp > 0) {
        $result *= $that;
        if (--$exp == $nexp) {
            $result += $this->factor_of($varname, $exp);
            $nexp = pop @exps;
        }
    }
    return $result;
}

sub partial_derivative {
    my ($this, $varname, $cast) = @_;
    my $result = $this->null;
    foreach my $exp ($this->exponents_of($varname)) {
        next if !$exp;
        my $const = $cast? $cast->($exp): $exp;
        $result +=
            $this->factor_of($varname, $exp) *
            $this->monomial($const, { $varname => $exp-1 });
    }
    return $result;
}

sub flatten {
    my ($this) = @_;
    while (_isa_pol(my $aux = $this->[F_ZERO])) {
        foreach my $term (values %{$this->[F_TERMS]}) {
            $aux += $term->[T_COEFF]->_mul_vars($term->[T_VARS]);
        }
        $this = $aux;
    }
    return $this;
}

sub extract {
    my ($this, $var) = @_;
    my %terms =
        map {
            my $vars  = { $_? ($var => $_): () };
            _key($vars) => [$this->factor_of($var, $_), $vars]
        } $this->exponents_of($var);
    return
        bless [
            \%terms,
            $this->null,
            $this->const($this->[F_ONE])
        ], ref $this;
}

# inspection methods

sub is_null     {  !keys %{$_[0]->[F_TERMS]} }
sub is_not_null { !!keys %{$_[0]->[F_TERMS]} }

sub variables {
    my ($this) = @_;
    my %vars = map {%{$_->[T_VARS]}} values %{$this->[F_TERMS]};
    return sort keys %vars;
}

sub exponents_of {
    my ($this, $varname) = @_;
    my $terms = $this->[F_TERMS];
    my %exps = ();
    foreach my $term (values %{$terms}) {
        $exps{$term->[T_VARS]->{$varname} || 0} = undef;
    }
    return sort { $a <=> $b } keys %exps;
}

sub factor_of {
    my ($this, $varname, $exp) = @_;
    my $terms  = $this->[F_TERMS];
    my $result = $this->null;
    my $rterms = $result->[F_TERMS];
    foreach my $term (values %{$terms}) {
        my $aexp = $term->[T_VARS]->{$varname} || 0;
        if ($aexp == $exp) {
            my %vars = %{$term->[T_VARS]};
            delete $vars{$varname};
            $rterms->{_key(\%vars)} = [$term->[T_COEFF], \%vars];
        }
    }
    return $result;
}

sub coefficient {
    my ($this, $vars) = @_;
    my $key = _keyvc($vars);
    my $terms = $this->[F_TERMS];
    return exists($terms->{$key})? $terms->{$key}->[T_COEFF]: $this->[F_ZERO];
}

sub degree {
    my ($this) = @_;
    my $terms = $this->[F_TERMS];
    return _MINUS_INFINITY if !keys %{$terms};
    my $max_degree = 0;
    foreach my $term (values %{$terms}) {
        my $degree = 0;
        foreach my $e (values %{$term->[T_VARS]}) {
            $degree += $e;
        }
        if ($max_degree < $degree) {
            $max_degree = $degree;
        }
    }
    return $max_degree;
}

sub multidegree {
    my ($this) = @_;
    my $terms = $this->[F_TERMS];
    my %result = ();
    foreach my $term (values %{$terms}) {
        while (my ($name, $exp) = each %{$term->[T_VARS]}) {
            if (($result{$name} || 0) < $exp) {
                $result{$name} = $exp;
            }
        }
    }
    return \%result;
}

sub number_of_terms { scalar keys %{$_[0]->[F_TERMS]} }

sub evaluate {
    my ($this, $values) = @_;
    my @vars = $this->variables;
    if (my @miss = grep { !exists $values->{$_} } @vars) {
        my $s = 1 == @miss? q[]: 's';
        croak "missing variable$s: @miss";
    }
    my $result = $this;
    foreach my $varname (@vars) {
        $result = $result->subst($varname, $this->const($values->{$varname}));
    }
    return $result->coefficient({});
}

sub as_monomials {
    my ($this) = @_;
    my $terms = $this->[F_TERMS];
    return scalar keys %{$terms}   if !wantarray;
    return ([$this->[F_ZERO], {}]) if !keys %{$terms};
    return
        map {
            my ($coeff, $vars) = @{$terms->{$_}};
            [$coeff, {%{$vars}}]
        } sort keys %{$terms};
}

sub as_string {
    my ($this) = @_;
    my $one = $this->[F_ONE];
    return
        sprintf '(%s)',
        join q[ + ],
        map {
            my ($const, $vars) = @{$_};
            join q[*],
            keys(%{$vars}) && $one == $const? (): $const,
            map {
                my $exp = $vars->{$_};
                1 == $exp? $_: "$_^$exp"
            }
            sort keys %{$vars}
        }
        $this->as_monomials;
}

1;
__END__