/*
* 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.
*
*/
/* ------------------------------------------------------------------------- */
/* monomarith.c: monomial arithmetic */
/* ------------------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "mymemory.h"
#include "simple.h"
#include "monomarith.h"
/*****************************************************************************/
/* Global data */
/*****************************************************************************/
/* monxy NULLMONXY = */
/* { */
/* 0.0, */
/* 0, */
/* 0 */
/* }; */
/* monxy ONEMONXY = */
/* { */
/* 1.0, */
/* 0, */
/* 0 */
/* }; */
monxyz NULLMONXYZ =
{
0.0,
0,
0,
0
};
monxyz ONEMONXYZ =
{
1.0,
0,
0,
0
};
/* ------------------------------------------------------------------------- */
/* allocate memory for n monomials and initialize with the nullmonomial */
/* ------------------------------------------------------------------------- */
monxyz *new_monxyz( int n )
{
int i;
monxyz *m;
if( n < 1 )
{
fprintf( stderr,"can't allocate %d monxyz\n",n );
return 0;
}
m = (monxyz*)malloc( n*sizeof(monxyz) );
if( m == (monxyz*)NULL )
{
fprintf( stderr,"can't allocate %d monxyz\n",n );
exit( 1 );
}
for( i = 0; i < n; i++ )
{
m[i] = NULLMONXYZ;
}
return m;
}
/* ------------------------------------------------------------------------- */
/* reallocate memory for n monomials */
/* ------------------------------------------------------------------------- */
monxyz *renew_monxyz( monxyz *m,int n )
{
if( m == (monxyz*)NULL )
{
return new_monxyz( n );
}
else
{
if( n < 1 )
{
fprintf( stderr,"can't reallocate %d monxyz\n",n );
exit( 1 );
}
m = (monxyz*)realloc( (void*)m,n*sizeof(monxyz) );
if( m == (monxyz*)NULL )
{
fprintf( stderr,"can't reallocate %d monxyz\n",n );
exit( 1 );
}
return m;
}
}
/* ------------------------------------------------------------------------- */
/* free memory of the monomials stored at m */
/* ------------------------------------------------------------------------- */
monxyz *delete_monxyz( monxyz *m )
{
if( m != (monxyz*)NULL )
{
free( (void*)m );
}
return (monxyz*)NULL;
}
/* ------------------------------------------------------------------------- */
/* copy monomials */
/* ------------------------------------------------------------------------- */
void copy_monxyz( monxyz *m1,monxyz *m2,int n )
{
memcpy( (void*)m1,(void*)m2,n*sizeof(monxyz) );
}
/* ------------------------------------------------------------------------- */
/* monomial power comparison */
/* ------------------------------------------------------------------------- */
int monxyz_power_equal( monxyz *m1,monxyz *m2 )
{
return m1->kx == m2->kx && m1->ky == m2->ky && m1->kz == m2->kz;
}
/* ------------------------------------------------------------------------- */
/* monomial comparison */
/* ------------------------------------------------------------------------- */
/* int monxyz_equal( monxyz *m1,monxyz *m2 ) */
/* { */
/* return monxyz_power_equal( m1,m2 ) && double_equal( m1->a,m2->a ); */
/* } */
/* ------------------------------------------------------------------------- */
/* monomial output */
/* ------------------------------------------------------------------------- */
void monxyz_print( monxyz *m )
{
fprintf( stderr,"%.4f*x^%d*y^%d*z^%d",m->a,m->kx,m->ky,m->kz );
}
/* ========================================================================= */
/* arithmetics */
/* ========================================================================= */
/* ------------------------------------------------------------------------- */
/* monomial addition: return = m1 + m2 */
/* ------------------------------------------------------------------------- */
/* monxyz monxyz_add( monxyz m1,monxyz *m2 ) */
/* { */
/* m1.a += m2->a; */
/* return m1; */
/* } */
/* ------------------------------------------------------------------------- */
/* monomial addition: m1 += m2 */
/* ------------------------------------------------------------------------- */
void monxyz_add_self( monxyz *m1,monxyz *m2 )
{
m1->a += m2->a;
}
/* ------------------------------------------------------------------------- */
/* monomial subtraction: return = m1 - m2 */
/* ------------------------------------------------------------------------- */
/* monxyz monxyz_sub( monxyz m1,monxyz *m2 ) */
/* { */
/* m1.a -= m2->a; */
/* return m1; */
/* } */
/* ------------------------------------------------------------------------- */
/* monomial subtraction: m1 -= m2 */
/* ------------------------------------------------------------------------- */
/* void monxyz_sub_self( monxyz *m1,monxyz *m2 ) */
/* { */
/* m1->a -= m2->a; */
/* } */
/* ------------------------------------------------------------------------- */
/* monomial multiplication: return = m1*m2 */
/* ------------------------------------------------------------------------- */
monxyz monxyz_mult( monxyz m1,monxyz *m2 )
{
m1.a *= m2->a;
m1.kx += m2->kx;
m1.ky += m2->ky;
m1.kz += m2->kz;
return m1;
}
/* ------------------------------------------------------------------------- */
/* monomial multiplication: m1 *= m2 */
/* ------------------------------------------------------------------------- */
void monxyz_mult_self( monxyz *m1,monxyz *m2 )
{
m1->a *= m2->a;
m1->kx += m2->kx;
m1->ky += m2->ky;
m1->kz += m2->kz;
}
/* ------------------------------------------------------------------------- */
/* monomial scalar multiplication: return = m*mult */
/* ------------------------------------------------------------------------- */
monxyz monxyz_mult_double( monxyz m,double mult )
{
m.a *= mult;
return m;
}
/* ------------------------------------------------------------------------- */
/* monomial scalar multiplication: m *= mult */
/* ------------------------------------------------------------------------- */
void monxyz_mult_double_self( monxyz *m,double mult )
{
m->a *= mult;
}
/* ------------------------------------------------------------------------- */
/* monomial scalar division: return = m/div */
/* ------------------------------------------------------------------------- */
monxyz monxyz_div_double( monxyz m,double div )
{
if( div == 0.0 )
{
fprintf( stderr,"can't divide by zero\n" );
exit( 1 );
}
m.a /= div;
return m;
}
/* ------------------------------------------------------------------------- */
/* monomial scalar division: m /= div */
/* ------------------------------------------------------------------------- */
void monxyz_div_double_self( monxyz *m,double div )
{
if( div == 0.0 )
{
fprintf( stderr,"can't divide by zero\n" );
exit( 1 );
}
m->a /= div;
}
/* ------------------------------------------------------------------------- */
/* monomial power: return = m^n */
/* ------------------------------------------------------------------------- */
monxyz monxyz_power( monxyz m,int n )
{
if( n > 0 )
{
#ifdef SUN
m.a = sunpow( m.a,(double)n );
#else
m.a = pow( m.a,(double)n );
#endif /* SUN */
m.kx *= n;
m.ky *= n;
m.kz *= n;
}
else if( n == 0 )
{
m.kx = m.ky = m.kz = 0;
m.a = 1.0;
}
else
{
fprintf( stderr,"can't take the %d-th power of a monomial\n",n );
exit( 1 );
}
return m;
}
/* ------------------------------------------------------------------------- */
/* monomial power: m ^= n */
/* ------------------------------------------------------------------------- */
/* void monxyz_power_self( monxyz *m,int n ) */
/* { */
/* if( n > 0 ) */
/* { */
/* m->kx *= n; */
/* m->ky *= n; */
/* m->kz *= n; */
/* } */
/* else if( n == 0 ) */
/* { */
/* m->kx = m->ky = m->kz = 0; */
/* m->a = 1.0; */
/* } */
/* else */
/* { */
/* fprintf( stderr,"can't take the %d-th power of a monomial\n",n ); */
/* exit( 1 ); */
/* } */
/* } */
/* ------------------------------------------------------------------------- */
/* monomial negation: return = -m */
/* ------------------------------------------------------------------------- */
monxyz monxyz_neg( monxyz m )
{
m.a = -m.a;
return m;
}
/* ------------------------------------------------------------------------- */
/* monomial negation: m = -m */
/* ------------------------------------------------------------------------- */
void monxyz_neg_self( monxyz *m )
{
m->a = -m->a;
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to x: return = dx( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dx( monxyz m )
{
if( m.kx == 0 )
{
return NULLMONXYZ;
}
else
{
m.a *= m.kx;
m.kx--;
return m;
}
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to x: m = dx( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dx_self( monxyz *m )
{
if( m->kx == 0 )
{
*m = NULLMONXYZ;
}
else
{
m->a *= m->kx;
m->kx--;
}
}
/* ------------------------------------------------------------------------- */
/* monomial n-fold derivative with respect to x: return = d^nx( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dx_n( monxyz m,int n )
{
if( m.kx == 0 || n > m.kx )
{
return NULLMONXYZ;
}
else
{
int i;
for( i=0; i<n; i++ )
{
m.a *= m.kx;
m.kx--;
}
return m;
}
}
/* ------------------------------------------------------------------------- */
/* monomial n-fold derivative with respect to x: m = d^nx( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dx_n_self( monxyz *m,int n )
{
if( m->kx == 0 || n > m->kx )
{
*m = NULLMONXYZ;
}
else
{
int i;
for( i=0; i<n; i++ )
{
m->a *= m->kx;
m->kx--;
}
}
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to y: return = dy( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dy( monxyz m )
{
if( m.ky == 0 )
{
return NULLMONXYZ;
}
else
{
m.a *= m.ky;
m.ky--;
return m;
}
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to y: m = dy( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dy_self( monxyz *m )
{
if( m->ky == 0 )
{
*m = NULLMONXYZ;
}
else
{
m->a *= m->ky;
m->ky--;
}
}
/* ------------------------------------------------------------------------- */
/* monomial n-fold derivative with respect to y: return = d^ny( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dy_n( monxyz m,int n )
{
if( m.ky == 0 || n > m.ky )
{
return NULLMONXYZ;
}
else
{
int i;
for( i=0; i<n; i++ )
{
m.a *= m.ky;
m.ky--;
}
return m;
}
}
/* ------------------------------------------------------------------------- */
/* monomial n-fold derivative with respect to y: m = d^ny( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dy_n_self( monxyz *m,int n )
{
if( m->ky == 0 || n > m->ky )
{
*m = NULLMONXYZ;
}
else
{
int i;
for( i=0; i<n; i++ )
{
m->a *= m->ky;
m->ky--;
}
}
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to z: return = dz( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dz( monxyz m )
{
if( m.kz == 0 )
{
return NULLMONXYZ;
}
else
{
m.a *= m.kz;
m.kz--;
return m;
}
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to x: m = dz( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dz_self( monxyz *m )
{
if( m->kz == 0 )
{
*m = NULLMONXYZ;
}
else
{
m->a *= m->kz;
m->kz--;
}
}
/* ------------------------------------------------------------------------- */
/* monomial n-fold derivative with respect to z: return = d^nz( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dz_n( monxyz m,int n )
{
if( m.kz == 0 || n > m.kz )
{
return NULLMONXYZ;
}
else
{
int i;
for( i=0; i<n; i++ )
{
m.a *= m.kz;
m.kz--;
}
return m;
}
}
/* ------------------------------------------------------------------------- */
/* monomial n-fold derivative with respect to z: m = d^nz( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dz_n_self( monxyz *m,int n )
{
if( m->kz == 0 || n > m->kz )
{
*m = NULLMONXYZ;
}
else
{
int i;
for( i=0; i<n; i++ )
{
m->a *= m->kz;
m->kz--;
}
}
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to w: return = dw( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dw( monxyz m,int deg )
{
int kw = deg - m.kx - m.ky - m.kz;
if( kw <= 0 )
{
return NULLMONXYZ;
}
else
{
m.a *= kw;
return m;
}
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to w: m = dw( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dw_self( monxyz *m,int deg )
{
int kw = deg - m->kx -m->ky - m->kz;
if( kw <= 0 )
{
*m = NULLMONXYZ;
}
else
{
m->a *= kw;
}
}
/* ------------------------------------------------------------------------- */
/* monomial n-fold derivative with respect to w: return = d^nw( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dw_n( monxyz m,int deg,int n )
{
int kw = deg - m.kx - m.ky - m.kz;
if( kw <= 0 || kw < n )
{
return NULLMONXYZ;
}
else
{
int i;
for( i=0; i<n; i++ )
{
m.a *= kw;
kw--;
}
return m;
}
}
/* ------------------------------------------------------------------------- */
/* monomial n-fold derivative with respect to w: m = d^nw( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dw_n_self( monxyz *m,int deg,int n )
{
int kw = deg - m->kx -m->ky - m->kz;
if( kw <= 0 || kw < n )
{
*m = NULLMONXYZ;
}
else
{
int i;
for( i=0; i<n; i++ )
{
m->a *= kw;
kw--;
}
}
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to a monmial x^iy^jz^k: */
/* return = d^ix d^jy d^kz d^lw( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dvector( monxyz m,int deg_m,int i,int j,int k,int l )
{
if( l > 0 ) monxyz_dw_n_self( &m,deg_m,l );
if( i > 0 ) monxyz_dx_n_self( &m,i );
if( j > 0 ) monxyz_dy_n_self( &m,j );
if( k > 0 ) monxyz_dz_n_self( &m,k );
return m;
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to a monmial x^iy^jz^k: */
/* m = d^ix d^jy d^kz d^lw( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dvector_self( monxyz *m,int deg_m,int i,int j,int k,int l )
{
if( l > 0 ) monxyz_dw_n_self( m,deg_m,l );
if( i > 0 ) monxyz_dx_n_self( m,i );
if( j > 0 ) monxyz_dy_n_self( m,j );
if( k > 0 ) monxyz_dz_n_self( m,k );
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to a monmial x^iy^jz^k: */
/* return = d^ix d^jy d^kz d^lw( m ) */
/* ------------------------------------------------------------------------- */
monxyz monxyz_dmonxyz( monxyz m,int deg_m,monxyz *d,int deg_d )
{
int d__kw = deg_d - d->kx -d->ky - d->kz;
if( d__kw > 0 ) monxyz_dw_n_self( &m,deg_m,d__kw );
if( d->kx > 0 ) monxyz_dx_n_self( &m,d->kx );
if( d->ky > 0 ) monxyz_dy_n_self( &m,d->ky );
if( d->kz > 0 ) monxyz_dz_n_self( &m,d->kz );
return m;
}
/* ------------------------------------------------------------------------- */
/* monomial derivative with respect to a monmial x^iy^jz^k: */
/* m = d^ix d^jy d^kz d^lw( m ) */
/* ------------------------------------------------------------------------- */
void monxyz_dmonxyz_self( monxyz *m,int deg_m,monxyz *d,int deg_d )
{
int d__kw = deg_d - d->kx -d->ky - d->kz;
if( d__kw > 0 ) monxyz_dw_n_self( m,deg_m,d__kw );
if( d->kx > 0 ) monxyz_dx_n_self( m,d->kx );
if( d->ky > 0 ) monxyz_dy_n_self( m,d->ky );
if( d->kz > 0 ) monxyz_dz_n_self( m,d->kz );
}
/* ------------------------------------------------------------------------- */
/* monomial conversion: m = (monxyz)i */
/* ------------------------------------------------------------------------- */
monxyz int2monxyz( int i )
{
monxyz m = NULLMONXYZ;
m.a = (double)i;
return m;
}
/* ------------------------------------------------------------------------- */
/* monomial conversion: m = (monxyz)d */
/* ------------------------------------------------------------------------- */
monxyz double2monxyz( double d )
{
monxyz m = NULLMONXYZ;
m.a = d;
return m;
}
/* ------------------------------------------------------------------------- */
/* monomial conversion: m = (monxyz)s */
/* ------------------------------------------------------------------------- */
monxyz atom( char *s )
{
monxyz m = ONEMONXYZ;
switch( s[0] )
{
case 'x' : m.kx++; break;
case 'X' : m.kx++; break;
case 'y' : m.ky++; break;
case 'Y' : m.ky++; break;
case 'z' : m.kz++; break;
case 'Z' : m.kz++; break;
}
return m;
}
/* ------------------------------------------------------------------------- */
/* end of file: monomarith.c */
/* ------------------------------------------------------------------------- */
syntax highlighted by Code2HTML, v. 0.9.1