/*
* surf - visualizing algebraic curves and algebraic surfaces
* Copyright (C) 1996-1997 Friedrich-Alexander-Universitaet
* Erlangen-Nuernberg
* 1997-2000 Johannes Gutenberg-Universitaet Mainz
* Authors: Stephan Endrass, Hans Huelf, Ruediger Oertel,
* Kai Schneider, Ralf Schmitt, Johannes Beigel
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "def.h"
#include "degree.h"
#include "simple.h"
#include "MultiVariatePolynom.h"
#include "MultiIndex.h" // for power function
#include "DrawFunc.h" // for errnums in check
#include "Monomial.h"
// ----------------------------------------------------------------------------
// ---------------- construktors of base - template ---------------------------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon>::MultiPoly()
: Number(1), deg(0), Monomial(new Mon [1])
{
Monomial[0] = 0.0;
}
template<class Mon>
MultiPoly<Mon>::MultiPoly( const int num )
: Number(num), deg(DEG_UNSPEC), Monomial(new Mon [num])
{
}
template<class Mon>
MultiPoly<Mon>::MultiPoly(const polyxyz alt)
: Number(alt.n), deg(alt.deg), Monomial(new Mon [alt.n])
{
for( int i =0; i < Number; i++)
Monomial[i]=alt.m[i];
}
template<class Mon>
MultiPoly<Mon>::MultiPoly( const MultiPoly &alt )
: Number(alt.Number), deg(alt.deg), Monomial(new Mon [alt.Number])
{
for( int i = 0; i < Number; i++ )
Monomial[i] = alt.Monomial[i];
}
template<class Mon>
MultiPoly<Mon>::MultiPoly( const double value )
: Number(1), deg(0), Monomial(new Mon [1])
{
Monomial[0] = value;
}
template<class Mon>
MultiPoly<Mon>::MultiPoly( const Mon &m )
: Number(1), deg(m.Degree()), Monomial(new Mon [1])
{
Monomial[0] = m;
}
template<class Mon>
MultiPoly<Mon>::~MultiPoly()
{
delete [] Monomial;
}
// ----------------------------------------------------------------------------
// -------------------- operators of template ------------------------------
// -------------------- reflexive ----------------------------------
// ----------------------------------------------------------------------------
// ----------------------- multiplications ------------------------------------
template<class Mon>
MultiPoly<Mon>& MultiPoly<Mon>::operator*=(const double Mult)
{
for( int i =0; i< Number; i++)
Monomial[i] *= Mult;
return *this;
}
template<class Mon>
MultiPoly<Mon>& MultiPoly<Mon>::operator*=(const Mon &mon)
{
for( int i = 0; i < Number; i++ )
Monomial[i] *= mon;
Collect();
return *this;
}
template<class Mon>
MultiPoly<Mon>& MultiPoly<Mon>::operator*=(const MultiPoly& Pol2)
{
Mon* Temp = new Mon [ Number * Pol2.Number ];
int i = 0, j = 0, k = 0, l = 0;
for( i = 0; i < Number; i++ )
for( j = 0; j < Pol2.Number; j++ ) {
// calculate next monomial
Mon Monom = Monomial[i] * Pol2.Monomial[j];
for( l = 0; l < k && !(Temp[l]>=Monom); l++ )
;
// see if monomial with same exponents exists
// if not first monomial or found -- add up
if ( k != l )
Temp[l] += Monom;
else
Temp[k++] = Monom; // otherwise add to chain
}
delete [] Monomial;
Monomial = Temp;
Number = k;
Renew( Number );
Collect();
return *this;
}
// ----------------------------------------------------------------------------
// ------------------- additions ----------------------------------------------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon>& MultiPoly<Mon>::operator+=( const MultiPoly& Pol )
{
Mon* Temp = new Mon [ Number + Pol.Number ];
int i = 0, j = 0;
for( i = 0; i < Number; i++ )
Temp[i] = Monomial[i];
for( i = 0; i < Pol.Number; i++ ) {
for( j = 0; j < Number && !(Pol.Monomial[i]>=Temp[j]); j++ );
// suche ob schon ein monom mit gleichen exponenten existiert
if ( j < Number )
Temp[j] += Pol.Monomial[i]; // falls gefunden aufaddieren
else
Temp[Number++] = Pol.Monomial[i]; // falls nicht anhängen
}
delete [] Monomial;
Monomial = Temp;
Renew( Number );
Collect();
return *this;
}
template<class Mon>
MultiPoly<Mon>& MultiPoly<Mon>::operator+=( const Mon& m )
{
int i = 0;
for( i = 0; i < Number && !(Monomial[i]>=m); i++ );
if ( i == Number ) {
Mon *Temp = new Mon [Number+1];
for( i = 0; i < Number; i++ )
Temp[i] = Monomial[i];
Temp[Number++] = m;
delete [] Monomial;
Monomial = Temp;
if( m.Degree() > deg ) deg = m.Degree();
} else
Monomial[i] += m;
return *this;
}
// ----------------------------------------------------------------------------
// ----------------------- assignments ----------------------------------------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon>& MultiPoly<Mon>::operator=( const double var )
{
delete [] Monomial;
Monomial = new Mon [1];
Monomial[0] = var;
deg = 0;
Number = 1;
return *this;
}
template<class Mon>
MultiPoly<Mon>& MultiPoly<Mon>::operator=(const MultiPoly& Pol)
{
Mon* Temp = new Mon [Pol.Number];
deg = Pol.deg; Number = Pol.Number;
for( int i = 0; i < Number; i++ )
Temp[i] = Pol.Monomial[i];
delete [] Monomial;
Monomial = Temp;
return *this;
}
// ----------------------------------------------------------------------------
// ------------------- normal multiplication ----------------------------------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon> MultiPoly<Mon>::operator*(const Mon &m) const
{
MultiPoly Pol = *this; Pol *= m; return Pol;
}
template<class Mon>
MultiPoly<Mon> MultiPoly<Mon>::operator*(const MultiPoly& Pol2) const
{
MultiPoly Pol = *this; Pol *= Pol2; return Pol;
}
// ----------------------------------------------------------------------------
// ------------------- comparators --------------------------------------------
// ----------------------------------------------------------------------------
template<class Mon>
int MultiPoly<Mon>::operator==( const MultiPoly& Pol ) const
{
if( Number != Pol.Number || deg != Pol.deg )
return FALSE;
int diff = 0;
for( int i = 0; i < Number && !diff; i++ )
diff = !( Monomial[i] == Pol.Monomial[i] );
return diff;
}
template<class Mon>
int MultiPoly<Mon>::operator!=( const MultiPoly& Pol ) const
{
return !( *this == Pol );
}
// ----------------------------------------------------------------------------
// ---------------- element functions of base template ------------------------
// ----------------------------------------------------------------------------
template<class Mon>
int MultiPoly<Mon>::Check(void) const
{
if( Number == 0 )
return CURVE_ZERO;
if( deg < 1 )
return CURVE_CONST;
if( deg >= MAX_DEGREE )
return CURVE_DEG_HIGH;
return FALSE;
}
// ----------------------------------------------------------------------------
// ---------- printout for debugging ------------------------------------------
// ----------------------------------------------------------------------------
template<class Mon>
void MultiPoly<Mon>::Print (ostream &os) const
{
for( int i=0; i< Number; i++) {
if( Monomial[i].Coeff() >= 0.0 )
os << " +";
else
os << " ";
//os << Monomial[i];
Monomial[i].Print(os);
os << endl;
}
}
// ----------------------------------------------------------------------------
// ------------ garbage collect -----------------------------------------------
// ----------------------------------------------------------------------------
template<class Mon>
void MultiPoly<Mon>::Collect(void)
{
int i,j;
for( i = Number - 1; i > 0 ; i--) {
for( j = i-1; j>=0 && !(Monomial[i]>=Monomial[j]); j-- )
;
if( j>=0 ) {
Monomial[j]+=Monomial[i];
if( (--Number) > i)
Monomial[i]=Monomial[Number];
}
}
while( i < Number ) {
if( Monomial[i].Coeff() )
i++;
else if( (--Number) > i )
Monomial[i] = Monomial[Number];
}
Renew( Number );
SetDegree();
}
// ----------------------------------------------------------------------------
// ------------ translate polynom -- part of positioning ----------------------
// ----------------------------------------------------------------------------
template<class Mon>
void MultiPoly<Mon>::Shift(double sx, double sy, double sz)
{
MultiPoly *SubPol = new MultiPoly [3];
SubPol[0] = -sx;
SubPol[1] = -sy;
SubPol[2] = -sz;
for( int i = 0; i < 3; i++ ) {
// Substitutionspolynome ' x-shiftx', 'y-shifty', 'z-shiftz' erzeugen
Mon Monom = 1.0; Monom.Exponent( i, 1 ); // monom ist aktuelle variable
SubPol[i] += Monom;
}
Subst( SubPol ); // und substituieren
delete [] SubPol;
Collect();
}
// ----------------------------------------------------------------------------
// ------------- substitute polynom for variable nr. var ----------------------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon> MultiPoly<Mon>::Subst( const MultiPoly<Mon> &SubstPoly, const int Var ) const
{
MultiPoly Result, Part; int i,j;
for( i = 0; i < Number; i++ ) {
Part = Monomial[i];
j = Part.Monomial[0].ExtractExponent( Var );
Part *= pow( SubstPoly, j );
Result += Part;
}
return Result;
}
// ----------------------------------------------------------------------------
// -------------- substitute polyfield for each variable ----------------------
// ----------------------------------------------------------------------------
template<class Mon>
void MultiPoly<Mon>::Subst( const MultiPoly<Mon> *SubstPoly )
{
MultiPoly<Mon> Result, Part; int i,j;
SetDegree();
// precalculate all powers 0...deg for all substitution polynoms
MultiPoly<Mon> **PowPoly = new MultiPoly<Mon>* [3];
for( i = 0; i < 3; i++ ) {
PowPoly[i] = new MultiPoly<Mon> [deg+1];
for( j = 0; j <= deg; j++ )
PowPoly[i][j] = pow( SubstPoly[i], j );
}
// place the adequately powered substitution polynoms in every monomial
for( i = 0; i < Number ; i++ ) {
Part = Monomial[i].Coeff();
for( j = 0; j < 3; j++ )
Part *= PowPoly[j][Monomial[i].Exponent(j)];
Result += Part;
}
for( j = 0; j < 3; j++ )
delete [] PowPoly[j];
delete [] PowPoly;
*this = Result;
}
// ----------------------------------------------------------------------------
// ------------ swap variables x and y ----------------------------------------
// ----------------------------------------------------------------------------
template<class Mon>
void MultiPoly<Mon>::SwapXY()
{
for( int i = 0; i < Number; i++ )
Monomial[i].SwapXY();
}
// ----------------------------------------------------------------------------
// ------------ solve a polynomial for one variable ---------------------------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon> MultiPoly<Mon>::Solve(const int Var)
{
int Location = -1;
int error = 0;
Collect();
for( int i = 0; i < Number && !error; i++ ) {
if( Monomial[i].Exponent( Var ) > 1 )
error = 1;
if( Monomial[i].Exponent( Var ) == 1 ) {
if( Monomial[i].Degree() != 1 )
error = 1;
Location = i;
}
}
// if anything is strange return zero
if( error == 1
|| Location == -1
|| Number == 1
|| Monomial[Location].Coeff() == 0.0 )
return 0.0;
double CoeffRez = - 1.0 / Monomial[Location].Coeff();
MultiPoly Result = *this;
Result.Monomial[Location] = 0.0;
Result *= CoeffRez;
Result.Collect();
return Result;
}
// ----------------------------------------------------------------------------
// ------------ sort monomials ------------------------------------------------
// ----------------------------------------------------------------------------
template<class Mon>
void MultiPoly<Mon>::Sort( const int Var )
{
for( int j = 1; j < Number; j++ ) {
for( int i = 1; i < Number; i++ ) {
if( Monomial[i-1].Compare( Monomial[i], Var ) > 0 ) {
Mon Temp = Monomial[i-1];
Monomial[i-1] = Monomial[i];
Monomial[i] = Temp;
}
}
}
}
// ----------------------------------------------------------------------------
// ------------------- power function for polynoms -- declared friend ---------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon> MultiPoly<Mon>::computePower (const int s) const
{
if( s < 0 ) {
fprintf( stderr, "Trying to get %d-th power of a MultiPoly\n",s);
exit(1);
}
if( s == 0 )
return 1.0;
if( s == 1 )
return *this;
if( Number == 0 )
return 0.0;
if( Number == 1 )
return (Mon) pow( Monomial[0], s );
MultiPoly<Mon> Pneu = 0.0;
Index MultiIndex( Number, s );
do {
Mon Monom = MultiIndex.Coeff();
for( int i = 0; i < Number; i++ )
Monom *= pow( Monomial[i], MultiIndex.Element(i) );
Pneu += Monom;
} while( ++MultiIndex );
return Pneu;
}
// template<class Mon>
// MultiPoly<Mon> pow( const MultiPoly<Mon>& Poly, const int s )
// {
// return Poly.computePower(s);
// }
// template<class Mon>
// MultiPoly<Mon> pow( const MultiPoly<Mon>& Poly, const int s )
// {
// if( s < 0 )
// {
// fprintf( stderr, "Trying to get %d-th power of a MultiPoly\n",s);
// exit(1);
// }
// if( s == 0 ) return 1.0;
// if( s == 1 ) return Poly;
// if( Poly.Number == 0 ) return 0.0;
// if( Poly.Number == 1 )
// return (Mon) pow( Poly.Monomial[0], s );
// MultiPoly<Mon> Pneu = 0.0;
// Index MultiIndex( Poly.Number, s );
// do{
// Mon Monom = MultiIndex.Coeff();
// for( int i = 0; i < Poly.Number; i++ )
// Monom *= pow( Poly.Monomial[i], MultiIndex.Element(i) );
// Pneu += Monom;
// }while( ++MultiIndex );
// return Pneu;
// }
// ----------------------------------------------------------------------------
// --------------------- special functions for derived classes ----------------
// - transformation from polyxyz to hornerxyz called by hornerxyz::operator= --
// ----------------------------------------------------------------------------
void Polyxyz::Transform( Polyx *PolY, int *pYepX, int *pYeX,
Polyx *PolX, int *pXeZ, Polyx *PolZ )
{
int pxn = 0, pyn = 0;
PolX[0].NewCoeff( Monomial[0].Exponent(0)+1 );
PolY[0].NewCoeff( Monomial[0].Exponent(1)+1 );
PolY[0].WrtCoeff( Monomial[0].Exponent(1), Monomial[0].Coeff() );
pYepX[0] = 0;
pYeX[0] = Monomial[0].Exponent(0);
pXeZ[0] = Monomial[0].Exponent(2);
for( int i = 1; i < Number; i++ ) {
if( !Monomial[i].CompareDeg( Monomial[i-1], 2 ) ) {
PolX[++pxn].NewCoeff( Monomial[i].Exponent(0)+1 );
pXeZ[pxn] = Monomial[i].Exponent(2);
}
if( !Monomial[i].CompareDeg( Monomial[i-1], 2 )
|| !Monomial[i].CompareDeg( Monomial[i-1], 0 ) ) {
PolY[++pyn].NewCoeff( Monomial[i].Exponent(1)+1 );
pYepX[pyn] = pxn;
pYeX[pyn] = Monomial[i].Exponent(0);
}
PolY[pyn].WrtCoeff( Monomial[i].Exponent(1) , Monomial[i].Coeff() );
}
PolZ->NewCoeff( Monomial[0].Exponent(2)+1 );
}
// ----------------------------------------------------------------------------
// --------- perspective trivariate polynom -----------------------------------
// ----------------------------------------------------------------------------
void Polyxyz::Perspective( const double z )
{
double z1 = -1.0/z;
MultiPoly<MonomXYZ> *SubPol = new MultiPoly<MonomXYZ>[3];
for( int i = 0; i < 2; i++ ) {
MonomXYZ Monom1 = z1;
Monom1.Exponent( i,1 );
Monom1.Exponent( 2,1 );
MonomXYZ Monom2 = 1.0;
Monom2.Exponent( i,1 );
SubPol[i] = Monom1; SubPol[i] += Monom2;
}
// substitute x by z1*xz + x , and y by z1*y z + y
MonomXYZ Monom3 = 1.0;
Monom3.Exponent( 2, 1 );
SubPol[2] = Monom3; // substitute z by z
Subst( SubPol );
delete [] SubPol;
Collect();
}
// ----------------------------------------------------------------------------
// - transformation from polyxy to hornerxy called by hornerxy::transform -----
// ------------------ two possible directions for ''hornering'' ---------------
// ----------------------------------------------------------------------------
void Polyxy::Transform( Polyx *Polynom, int *Exponent, Polyx* PolynomX,
int dir ) const
{
dir = dir & 1;
int Counter = 0; Polynom[0].NewCoeff( Monomial[0].Exponent(1-dir)+1 );
Polynom[0].WrtCoeff( Monomial[0].Exponent(1-dir) , Monomial[0].Coeff() );
Exponent[0] = Monomial[0].Exponent(dir);
for( int i = 1; i < Number; i++ ) {
if( !Monomial[i].CompareDeg( Monomial[i-1], dir ) ) {
Polynom[ ++Counter ].NewCoeff( Monomial[i].Exponent(1-dir)+1 );
Exponent[ Counter ] = Monomial[i].Exponent(dir);
}
Polynom[ Counter ].WrtCoeff( Monomial[i].Exponent(1-dir),
Monomial[i].Coeff() );
}
PolynomX->NewCoeff( Monomial[0].Exponent(dir)+1 );
}
// ----------------------------------------------------------------------------
// solve polynom cut for variable z and subst result in this, return as polyxy
// ----------------------------------------------------------------------------
// -------- return denominatorof polynom solved for variable nr. var ----------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon> MultiPoly<Mon>::CounterOfSolve( int Var ) const
{
MultiPoly Counter = 0.0;
for( int i = 0; i < Number; i++ ) {
if( Monomial[i].Exponent( Var ) == 0 ) {
Mon Monom = Monomial[i]; // copy monom
Monom.SetCoeff( -Monom.Coeff( ) ); // invert sign of coefficient
Counter += Monom; // add to denomiantor
}
}
Counter.Collect();
return Counter;
}
// ----------------------------------------------------------------------------
// -------- return nominator of polynom solved for variable nr. var -----------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon> MultiPoly<Mon>::NominatorOfSolve( int Var ) const
{
MultiPoly Nominator = 0.0;
for( int i = 0; i < Number; i++ ) {
if( Monomial[i].Exponent( Var ) == 1 ) {
Mon Monom = Monomial[i]; // copy monom
Monom.Exponent( Var, 0 ); // clear z exponent
Nominator += Monom; // add to nominator
}
}
Nominator.Collect();
return Nominator;
}
// ----------------------------------------------------------------------------
// ------------- substitute rational polynom for variable nr. var -------------
// ----------------------------------------------------------------------------
template<class Mon>
MultiPoly<Mon> MultiPoly<Mon>::Subst( const MultiPoly& Counter,
const MultiPoly& Nominator,
const int Var ) const
{
int i = 0;
int ZmaxDeg = MaxDegreeOfVariable( 2 ); // calc max z degree of surface
MultiPoly *PowPoly = new MultiPoly [ZmaxDeg+1];
for( i = 0; i <= ZmaxDeg; i++ ) {
PowPoly[i] = pow( Counter,i )*pow( Nominator,ZmaxDeg-i );
}
MultiPoly Result = 0.0, Part = 0.0;
for( i = 0; i < Number; i++ ) {
Mon Monom = Monomial[i]; // copy monom
int zDeg = Monom.ExtractExponent( Var );
Part = Monom;
Part *= PowPoly[zDeg];
Result += Part;
}
delete []PowPoly;
return Result;
}
double Polyxyz::evaluate(double x, double y, double z)
{
double result=0;
int i;
for (i=0; i<Number; i++) {
result += Monomial[i].a*pow(x, (double)Monomial[i].k[0])*pow(y, (double)Monomial[i].k[1])*pow(z, (double)Monomial[i].k[2]);
}
return result;
}
// ----------------------------------------------------------------------------
// ------------ declare templates -- at end for gcc 2.6.3 ---------------------
// ----------------------------------------------------------------------------
#ifdef HAVE_EXPLICIT_TEMPLATE_INSTANTIATION
template class MultiPoly<MonomXY>;
template class MultiPoly<MonomXYZ>;
template class Monom<2>;
template class Monom<3>;
#endif
// ----------------------------------------------------------------------------
// ------------- end of file MultiVariatePolynom.cc ---------------------------
// ----------------------------------------------------------------------------
syntax highlighted by Code2HTML, v. 0.9.1