OpenXM_contrib2/asir2000/engine/PM.c - view

Return to PM.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / engine

File: [local] / OpenXM_contrib2 / asir2000 / engine / PM.c (download)

Revision 1.1, Fri Dec 3 07:39:08 1999 UTC (24 years, 5 months ago) by noro
Branch: MAIN

Initial revision

/* $OpenXM: OpenXM_contrib2/asir2000/engine/PM.c,v 1.1 1999/12/03 07:39:08 noro Exp $ */
#ifndef MODULAR
#define MODULAR
#endif

#include "b.h"
#include "ca.h"

void D_ADDP(vl,p1,p2,pr)
VL vl;
P p1,p2,*pr;
{
	register DCP dc1,dc2,dcr0,dcr;
	V v1,v2;
	P t;

	if ( !p1 ) 
		*pr = p2;
	else if ( !p2 ) 
		*pr = p1;
	else if ( NUM(p1) )
		if ( NUM(p2) ) 
			ADDNUM(p1,p2,pr);
		else 
			ADDPQ(p2,p1,pr);
	else if ( NUM(p2) ) 
		ADDPQ(p1,p2,pr);
	else if ( ( v1 = VR(p1) ) ==  ( v2 = VR(p2) ) ) {
		for ( dc1 = DC(p1), dc2 = DC(p2), dcr0 = 0; dc1 && dc2; )
			switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
				case 0: 
					ADDP(vl,COEF(dc1),COEF(dc2),&t);
					if ( t )  {
						NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = t;
					}
					dc1 = NEXT(dc1); dc2 = NEXT(dc2); break;
				case 1:
					NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = COEF(dc1);
					dc1 = NEXT(dc1); break;
				case -1:
					NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc2); COEF(dcr) = COEF(dc2);
					dc2 = NEXT(dc2); break;
			}
		if ( !dcr0 ) 
			if ( dc1 )
				dcr0 = dc1;
			else if ( dc2 ) 
				dcr0 = dc2;
			else {
				*pr = 0;
				return;
			}
		else 
			if ( dc1 ) 
				NEXT(dcr) = dc1;
			else if ( dc2 ) 
				NEXT(dcr) = dc2;
			else 
				NEXT(dcr) = 0;
		MKP(v1,dcr0,*pr);
	} else {
		while ( v1 != VR(vl) && v2 != VR(vl) ) 
			vl = NEXT(vl);
		if ( v1 == VR(vl) ) 
			ADDPTOC(vl,p1,p2,pr);
		else 
			ADDPTOC(vl,p2,p1,pr);
	}
}
	
void D_ADDPQ(p,q,pr)
P p,q,*pr;
{
	DCP dc,dcr,dcr0;
	P t;

	if ( NUM(p) ) 
		ADDNUM(p,q,pr);
	else {
		for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {
			NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
		}
		if ( !dc ) {
			NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = q;
		} else {
			ADDPQ(COEF(dc),q,&t);
			if ( t ) {
				NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;
			} 
		}
		NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
	}
}

void D_ADDPTOC(vl,p,c,pr)
VL vl;
P p,c,*pr;
{
	DCP dc,dcr,dcr0;
	P t;

	for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {
		NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
	}
	if ( !dc ) {
		NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = c;
	} else {
		ADDP(vl,COEF(dc),c,&t);
		if ( t ) {
			NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;
		} 
	}
	NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
}

void D_SUBP(vl,p1,p2,pr)
VL vl;
P p1,p2,*pr;
{
	P t;

	if ( !p2 ) 
		*pr = p1;
	else {
		CHSGNP(p2,&t); ADDP(vl,p1,t,pr);
	}
}

void D_MULP(vl,p1,p2,pr)
VL vl;
P p1,p2,*pr;
{
	register DCP dc,dct,dcr,dcr0;
	V v1,v2;
	P t,s,u;
	int n1,n2;

	if ( !p1 || !p2 ) *pr = 0;
	else if ( NUM(p1) ) 
		MULPQ(p2,p1,pr);
	else if ( NUM(p2) )
		MULPQ(p1,p2,pr);
	else if ( ( v1 = VR(p1) ) ==  ( v2 = VR(p2) ) ) {
		for ( dc = DC(p1), n1 = 0; dc; dc = NEXT(dc), n1++ );
		for ( dc = DC(p2), n2 = 0; dc; dc = NEXT(dc), n2++ );
		if ( n1 > n2 )
			for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
				for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
					NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
					addq(DEG(dct),DEG(dc),&DEG(dcr));
				}
				NEXT(dcr) = 0; MKP(v1,dcr0,t);
				ADDP(vl,s,t,&u); s = u; t = u = 0;
			}
		else
			for ( dc = DC(p1), s = 0; dc; dc = NEXT(dc) ) {
				for ( dcr0 = 0, dct = DC(p2); dct; dct = NEXT(dct) ) {
					NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
					addq(DEG(dct),DEG(dc),&DEG(dcr));
				}
				NEXT(dcr) = 0; MKP(v1,dcr0,t);
				ADDP(vl,s,t,&u); s = u; t = u = 0;
			}
		*pr = s;
	} else {
		while ( v1 != VR(vl) && v2 != VR(vl) ) 
			vl = NEXT(vl);
		if ( v1 == VR(vl) ) 
			MULPC(vl,p1,p2,pr);
		else 
			MULPC(vl,p2,p1,pr);
	}
}

void D_MULPQ(p,q,pr)
P p,q,*pr;
{
	DCP dc,dcr,dcr0;

	if (!p || !q)
		*pr = 0;
	else if ( Uniq(q) ) 
		*pr = p;
	else if ( NUM(p) ) 
		MULNUM(p,q,pr);
	else {
		for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
			NEXTDC(dcr0,dcr); MULPQ(COEF(dc),q,&COEF(dcr)); DEG(dcr) = DEG(dc);
		}
		NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
	}
}

void D_MULPC(vl,p,c,pr)
VL vl;
P p,c,*pr;
{
	DCP dc,dcr,dcr0;

	if ( NUM(c) ) 
		MULPQ(p,c,pr);
	else {
		for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
			NEXTDC(dcr0,dcr); MULP(vl,COEF(dc),c,&COEF(dcr)); DEG(dcr) = DEG(dc);
		}
		NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
	}
}

void D_PWRP(vl,p,q,pr)
VL vl;
P p,*pr;
Q q;
{
	register DCP dc,dcr;
	int n,i;
	P *x,*y;
	P t,s,u;
	DCP dct;
	P *pt;

	if ( !q ) { 
		*pr = (P)One; 
	} else if ( !p ) 
		*pr = 0;
	else if ( UNIQ(q) ) 
		*pr = p;
	else if ( NUM(p) ) 
		PWRNUM(p,q,pr);
	else {
		dc = DC(p);
		if ( !NEXT(dc) ) {
			NEWDC(dcr);
			PWRP(vl,COEF(dc),q,&COEF(dcr)); mulq(DEG(dc),q,&DEG(dcr));
			NEXT(dcr) = 0; MKP(VR(p),dcr,*pr);
		} else if ( !INT(q) ) {
			error("pwrp: can't calculate fractional power."); *pr = 0;
		} else if ( PL(NM(q)) == 1 ) {
			n = QTOS(q); x = (P *)ALLOCA((n+1)*sizeof(pointer));
			NEWDC(dct); DEG(dct) = DEG(dc); COEF(dct) = COEF(dc);
			NEXT(dct) = 0; MKP(VR(p),dct,t);
			for ( i = 0, u = (P)One; i < n; i++ ) {
				x[i] = u; MULP(vl,u,t,&s); u = s;
			}
			x[n] = u; y = (P *)ALLOCA((n+1)*sizeof(pointer));
			if ( DEG(NEXT(dc)) ) {
				dct = NEXT(dc); MKP(VR(p),dct,t);
			} else 
				t = COEF(NEXT(dc));
			for ( i = 0, u = (P)One; i < n; i++ ) {
				y[i] = u; MULP(vl,u,t,&s); u = s;
			}
			y[n] = u; 
			pt = (P *)ALLOCA((n+1)*sizeof(pointer)); MKBC(n,pt);
			for ( i = 0, u = 0; i <= n; i++ ) {
				MULP(vl,x[i],y[n-i],&t); MULP(vl,t,pt[i],&s); 
				ADDP(vl,u,s,&t); u = t;
			}
			*pr = u;
		} else {
			error("exponent too big");
			*pr = 0;
		}		
	}
}

void D_CHSGNP(p,pr)
P p,*pr;
{
	register DCP dc,dcr,dcr0;
	P t;

	if ( !p ) 
		*pr = NULL;
	else if ( NUM(p) ) {
#if defined(THINK_C) || defined(_PA_RISC1_1) || defined(__alpha) || defined(mips)
#ifdef FBASE
		chsgnnum((Num)p,(Num *)pr);
#else
		MQ mq;

		NEWMQ(mq); CONT(mq)=mod-CONT((MQ)p); *pr = (P)mq;
#endif
#else
		CHSGNNUM(p,t); *pr = t;
#endif
	} else {
		for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
			NEXTDC(dcr0,dcr); CHSGNP(COEF(dc),&COEF(dcr)); DEG(dcr) = DEG(dc);
		}
		NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
	}
}


#ifdef FBASE
void ptozp(p,sgn,c,pr)
P p;
int sgn;
Q *c;
P *pr;
{
	N nm,dn;

	if ( !p ) {
		*c = 0; *pr = 0;
	} else {
		lgp(p,&nm,&dn);
		if ( UNIN(dn) )
			NTOQ(nm,sgn,*c);
		else 
			NDTOQ(nm,dn,sgn,*c);
		divcp(p,*c,pr);
	}
}

void lgp(p,g,l)
P p;
N *g,*l;
{
	N g1,g2,l1,l2,l3,l4,tmp;
	DCP dc;

	if ( NUM(p) ) {
		*g = NM((Q)p);
		if ( INT((Q)p) ) 
			*l = ONEN;
		else 
			*l = DN((Q)p);
	} else {
		dc = DC(p); lgp(COEF(dc),g,l);
		for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
			lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;
			gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;
		}
	} 
}

void divcp(f,q,rp)
P f;
Q q;
P *rp;
{
	DCP dc,dcr,dcr0;

	if ( !f ) 
		*rp = 0;
	else if ( NUM(f) )
		DIVNUM(f,q,rp);
	else {
		for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
			NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
			divcp(COEF(dc),q,&COEF(dcr));
		}
		NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
	}
}

void diffp(vl,p,v,r)
VL vl;
P p;
V v;
P *r;
{
	P t;
	DCP dc,dcr,dcr0;

	if ( !p || NUM(p) ) 
		*r = 0;
	else {
		if ( v == VR(p) ) {
			for ( dc = DC(p), dcr0 = 0; 
				dc && DEG(dc); dc = NEXT(dc) ) {
				NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));
				MULPQ(COEF(dc),(P)DEG(dc),&COEF(dcr));
			}
			NEXT(dcr) = 0; MKP(v,dcr0,*r);
		} else {
			for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
				diffp(vl,COEF(dc),v,&t);
				if ( t ) {
					NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
				}
			}
			if ( !dcr0 )
				*r = 0;
			else {
				NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
			}
		}
	}
}

void coefp(p,d,pr)
P p;
int d;
P *pr;
{
	DCP dc;
	int sgn;
	Q dq;

	if ( NUM(p) )
		if ( d == 0 ) 
			*pr = p;
		else 
			*pr = 0;
	else {
		for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) ) 
			if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
				continue;
			else if ( sgn == 0 ) {
				*pr = COEF(dc);
				return;
			} else {
				*pr = 0;
				break;	
			}
		*pr = 0;
	}
}

int compp(vl,p1,p2)
VL vl;
P p1,p2;
{
	DCP dc1,dc2;
	V v1,v2;

	if ( !p1 )
		return p2 ? -1 : 0;
	else if ( !p2 )
		return 1;
	else if ( NUM(p1) )
		return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
	else if ( NUM(p2) )
		return 1;
	if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
		for ( dc1 = DC(p1), dc2 = DC(p2); 
			dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
			switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
				case 1:
					return 1; break;
				case -1:
					return -1; break;
				default:
					switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
						case 1:
							return 1; break;
						case -1:
							return -1; break;
						default:
							break;
					}
					break;
			}
		return dc1 ? 1 : (dc2 ? -1 : 0);
	} else {
		for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
		return v1 == VR(vl) ? 1 : -1;
	}
}

void csump(vl,p,c)
VL vl;
P p;
Q *c;
{
	DCP dc;
	Q s,s1,s2;

	if ( !p || NUM(p) )
		*c = (Q)p;
	else {
		for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
			csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
		}
		*c = s;
	}
}

void degp(v,p,d)
V v;
P p;
Q *d;
{
	DCP dc;
	Q m,m1;
	
	if ( !p || NUM(p) ) 
		*d = 0;
	else if ( v == VR(p) ) 
		*d = DEG(DC(p));
	else {
		for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
			degp(v,COEF(dc),&m1);
			if ( cmpq(m,m1) < 0 )
				m = m1;
		}
		*d = m;
	}
}
#endif