/******************************************************************************
* Symbolic.c - Generic symbolic computation.				      *
*******************************************************************************
* Written by Gershon Elber, Nov. 92.					      *
******************************************************************************/

#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include "symb_loc.h"

static CagdCrvStruct *SymbCrvAddSubAux(CagdCrvStruct *Crv1,
				       CagdCrvStruct *Crv2,
				       CagdBType OperationAdd);
static CagdSrfStruct *SymbSrfAddSubAux(CagdSrfStruct *Srf1,
				       CagdSrfStruct *Srf2,
				       CagdBType OperationAdd);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two curves - add them coordinatewise.				     M
*   The two curves are promoted to same point type before the multiplication M
* can take place. Furthermore, order and continuity are matched as well.     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:  Two curve to add up coordinatewise.                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   The summation of Crv1 + Crv2 coordinatewise.          M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvAdd, addition, symbolic computation                               M
*****************************************************************************/
CagdCrvStruct *SymbCrvAdd(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    return SymbCrvAddSubAux(Crv1, Crv2, TRUE);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two curves - subtract them coordinatewise.			     M
*   The two curves are promoted to same point type before the multiplication M
* can take place. Furthermore, order and continuity are matched as well.     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:  Two curve to subtract coordinatewise.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   The difference of Crv1 - Crv2 coordinatewise.         M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvSub, subtraction, symbolic computation                            M
*****************************************************************************/
CagdCrvStruct *SymbCrvSub(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    return SymbCrvAddSubAux(Crv1, Crv2, FALSE);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two curves - multiply them coordinatewise.			     M
*   The two curves are promoted to same point type before the multiplication M
* can take place.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:  Two curve to multiply coordinatewise.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   The product of Crv1 * Crv2 coordinatewise.            M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvMult, product, symbolic computation                               M
*****************************************************************************/
CagdCrvStruct *SymbCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct
	*ProdCrv = NULL;

    switch (Crv1 -> GType) {
	case CAGD_CBEZIER_TYPE:
	    ProdCrv = BzrCrvMult(Crv1, Crv2);
	    break;
	case CAGD_CBSPLINE_TYPE:
	    ProdCrv = BspCrvMult(Crv1, Crv2);
	    break;
	case CAGD_CPOWER_TYPE:
	    SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
	    break;
    }

    return ProdCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a scalar curve, returns a scalar curve representing the reciprocal   M
* values, by making it rational (if was not one) and flipping the numerator  M
* and the denominator.					       		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       A scalar curve to compute a reciprocal value for.             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A rational scalar curve that is equalto the           M
*                      reciprocal value of Crv.                              M
* KEYWORDS:                                                                  M
*   SymbCrvInvert, division, symbolic computation, reciprocal value          M
*****************************************************************************/
CagdCrvStruct *SymbCrvInvert(CagdCrvStruct *Crv)
{
    int i;
    CagdRType *R, **Points;

    Crv = CagdCrvCopy(Crv);
    Points = Crv -> Points;

    switch (Crv -> PType) {
	case CAGD_PT_E1_TYPE:
	    Points[0] = Points[1];
	    Points[1] = R = (CagdRType *) IritMalloc(sizeof(CagdRType) *
								Crv -> Length);
	    for (i = 0; i < Crv -> Length; i++)
		*R++ = 1.0;
	    Crv -> PType = CAGD_PT_P1_TYPE;
	    break;
	case CAGD_PT_P1_TYPE:
	    R = Points[0];
	    Points[0] = Points[1];
	    Points[1] = R;
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNSUPPORT_PT);
	    break;
    }

    return Crv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a curve, scale it by Scale.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:      A curve to scale by magnitude Scale.                           M
*   Scale:    Scaling factor.                                                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A curves scaled by Scale compared to Crv.             M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvScalarScale, scaling, symbolic computation                        M
*****************************************************************************/
CagdCrvStruct *SymbCrvScalarScale(CagdCrvStruct *Crv, CagdRType Scale)
{
    int i, j;
    CagdRType *R;

    Crv = CagdCrvCopy(Crv);

    for (j = 1; j <= CAGD_NUM_OF_PT_COORD(Crv -> PType); j++)
	for (i = 0, R = Crv -> Points[j]; i < Crv -> Length; i++)
	    *R++ *= Scale;

    return Crv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two curves - computes their dot product.			     M
*   Returned curve is a scalar curve representing the dot product of the     M
* two given curves.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:  Two curve to multiply and compute a dot product for.        M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A scalar curve representing the dot product of        M
*                      Crv1 . Crv2.                                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvDotProd, product, dot product, symbolic computation               M
*****************************************************************************/
CagdCrvStruct *SymbCrvDotProd(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *PCrvW, *PCrvX, *PCrvY, *PCrvZ, *TCrv1, *TCrv2, *DotProdCrv,
	*ProdCrv = SymbCrvMult(Crv1, Crv2);

    SymbCrvSplitScalar(ProdCrv, &PCrvW, &PCrvX, &PCrvY, &PCrvZ);
    CagdCrvFree(ProdCrv);

    if (PCrvY != NULL) {
	TCrv1 = SymbCrvAdd(PCrvX, PCrvY);
	CagdCrvFree(PCrvX);
	CagdCrvFree(PCrvY);
    }
    else
	TCrv1 = PCrvX;

    if (PCrvZ != NULL) {
	TCrv2 = SymbCrvAdd(TCrv1, PCrvZ);
	CagdCrvFree(TCrv1);
	CagdCrvFree(PCrvZ);
	TCrv1 = TCrv2;
    }

    DotProdCrv = SymbCrvMergeScalar(PCrvW, TCrv1, NULL, NULL);
    CagdCrvFree(PCrvW);
    CagdCrvFree(TCrv1);

    return DotProdCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two curves - a vector curve Crv1 and a scalar curve Crv2, multiply   M
* all Crv1's coordinates by the scalar curve Crv2.			     M
*   Returned curve is a curve representing the product of the two given	     M
* curves.								     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:  Two curve to multiply.					     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A curve representing the product of Crv1 and Crv2.    M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvMultScalar, product, symbolic computation  		             M
*****************************************************************************/
CagdCrvStruct *SymbCrvMultScalar(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *PCrv1W, *PCrv1X, *PCrv1Y, *PCrv1Z,
		  *PCrv2W, *PCrv2X, *PCrv2Y, *PCrv2Z,
		  *TCrv, *ProdCrv;

    SymbCrvSplitScalar(Crv1, &PCrv1W, &PCrv1X, &PCrv1Y, &PCrv1Z);
    SymbCrvSplitScalar(Crv2, &PCrv2W, &PCrv2X, &PCrv2Y, &PCrv2Z);

    TCrv = SymbCrvMult(PCrv1X, PCrv2X);
    CagdCrvFree(PCrv1X);
    PCrv1X = TCrv;

    if (PCrv1Y != NULL) {
	TCrv = SymbCrvMult(PCrv1Y, PCrv2X);
	CagdCrvFree(PCrv1Y);
	PCrv1Y = TCrv;
    }

    if (PCrv1Z != NULL) {
	TCrv = SymbCrvMult(PCrv1Z, PCrv2X);
	CagdCrvFree(PCrv1Z);
	PCrv1Z = TCrv;
    }

    if (PCrv1W != NULL && PCrv2W != NULL) {
	TCrv = SymbCrvMult(PCrv1W, PCrv2W);
	CagdCrvFree(PCrv1W);
	PCrv1W = TCrv;
    }
    else if (PCrv2W != NULL) {
	PCrv1W = PCrv2W;
	PCrv2W = NULL;
    }

    ProdCrv = SymbCrvMergeScalar(PCrv1W, PCrv1X, PCrv1Y, PCrv1Z);

    CagdCrvFree(PCrv1W);
    CagdCrvFree(PCrv1X);
    CagdCrvFree(PCrv1Y);
    CagdCrvFree(PCrv1Z);
    CagdCrvFree(PCrv2W);
    CagdCrvFree(PCrv2X);
    CagdCrvFree(PCrv2Y);
    CagdCrvFree(PCrv2Z);

    return ProdCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two curves - computes their cross product.			     M
*   Returned curve is a scalar curve representing the cross product of the   M
* two given curves.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:  Two curve to multiply and compute a cross product for.      M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A scalar curve representing the cross product of      M
*                      Crv1 x Crv2.                                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvCrossProd, product, cross product, symbolic computation           M
*****************************************************************************/
CagdCrvStruct *SymbCrvCrossProd(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z,
		  *Crv2W, *Crv2X, *Crv2Y, *Crv2Z,
		  *TCrv1, *TCrv2, *CrossProdCrv,
		  *PCrvW, *PCrvX, *PCrvY, *PCrvZ;
    CagdBType
	Crv1New = FALSE,
	Crv2New = FALSE;

    if (Crv1 -> GType != CAGD_PT_E3_TYPE && Crv1 -> GType != CAGD_PT_P3_TYPE) {
	Crv1 = CagdCoerceCrvTo(Crv1,
			       CAGD_IS_RATIONAL_CRV(Crv1) ? CAGD_PT_P3_TYPE
							  : CAGD_PT_E3_TYPE);
	Crv1New = TRUE;
    }
    if (Crv2 -> GType != CAGD_PT_E3_TYPE && Crv2 -> GType != CAGD_PT_P3_TYPE) {
	Crv2 = CagdCoerceCrvTo(Crv2,
			       CAGD_IS_RATIONAL_CRV(Crv2) ? CAGD_PT_P3_TYPE
							  : CAGD_PT_E3_TYPE);
	Crv2New = TRUE;
    }

    SymbCrvSplitScalar(Crv1, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z);
    SymbCrvSplitScalar(Crv2, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z);

    if (Crv1New)
	CagdCrvFree(Crv1);
    if (Crv2New)
	CagdCrvFree(Crv2);

    /* Cross product X axis. */
    TCrv1 = SymbCrvMult(Crv1Y, Crv2Z);
    TCrv2 = SymbCrvMult(Crv2Y, Crv1Z);
    PCrvX = SymbCrvSub(TCrv1, TCrv2);
    CagdCrvFree(TCrv2);
    CagdCrvFree(TCrv1);

    /* Cross product Y axis. */
    TCrv1 = SymbCrvMult(Crv1Z, Crv2X);
    TCrv2 = SymbCrvMult(Crv2Z, Crv1X);
    PCrvY = SymbCrvSub(TCrv1, TCrv2);
    CagdCrvFree(TCrv2);
    CagdCrvFree(TCrv1);

    /* Cross product Z axis. */
    TCrv1 = SymbCrvMult(Crv1X, Crv2Y);
    TCrv2 = SymbCrvMult(Crv2X, Crv1Y);
    PCrvZ = SymbCrvSub(TCrv1, TCrv2);
    CagdCrvFree(TCrv1);
    CagdCrvFree(TCrv2);

    /* Cross product W axis. */
    if (Crv1W || Crv2W) {
	if (Crv1W == NULL)
	    PCrvW = CagdCrvCopy(Crv2W);
	else if (Crv2W == NULL)
	    PCrvW = CagdCrvCopy(Crv1W);
	else
	    PCrvW = SymbCrvMult(Crv1W, Crv2W);
    }
    else
	PCrvW = NULL;

    CagdCrvFree(Crv1X);
    CagdCrvFree(Crv1Y);
    CagdCrvFree(Crv1Z);
    CagdCrvFree(Crv1W);

    CagdCrvFree(Crv2X);
    CagdCrvFree(Crv2Y);
    CagdCrvFree(Crv2Z);
    CagdCrvFree(Crv2W);

    if (!CagdMakeCrvsCompatible(&PCrvX, &PCrvY, TRUE, TRUE) ||
	!CagdMakeCrvsCompatible(&PCrvX, &PCrvZ, TRUE, TRUE) ||
	!CagdMakeCrvsCompatible(&PCrvY, &PCrvZ, TRUE, TRUE) ||
	!CagdMakeCrvsCompatible(&PCrvW, &PCrvX, TRUE, TRUE) ||
	!CagdMakeCrvsCompatible(&PCrvW, &PCrvY, TRUE, TRUE) ||
	!CagdMakeCrvsCompatible(&PCrvW, &PCrvZ, TRUE, TRUE))
	SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT);

    CrossProdCrv = SymbCrvMergeScalar(PCrvW, PCrvX, PCrvY, PCrvZ);
    CagdCrvFree(PCrvX);
    CagdCrvFree(PCrvY);
    CagdCrvFree(PCrvZ);
    CagdCrvFree(PCrvW);

    return CrossProdCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two curves - multiply them using the quotient product rule:	     M
*  X = X1 W2 +/- X2 W1							     V
*   All provided curves are assumed to be non rational scalar curves.	     M
*   Returned is a non rational scalar curve (CAGD_PT_E1_TYPE).		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1X:     Numerator of first curve.                                     M
*   Crv1W:     Denominator of first curve.                                   M
*   Crv2X:     Numerator of second curve.                                    M
*   Crv2W:     Denominator of second curve.                                  M
*   OperationAdd:   TRUE for addition, FALSE for subtraction.                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  The result of  Crv1X Crv2W +/- Crv2X Crv1W.            M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvRtnlMult, product, symbolic computation                           M
*****************************************************************************/
CagdCrvStruct *SymbCrvRtnlMult(CagdCrvStruct *Crv1X,
			       CagdCrvStruct *Crv1W,
			       CagdCrvStruct *Crv2X,
			       CagdCrvStruct *Crv2W,
			       CagdBType OperationAdd)
{
    CagdCrvStruct *CTmp1, *CTmp2, *CTmp3;

    if (Crv1X == NULL || Crv2X == NULL)
	return NULL;

    /* Make the two curves - same order and point type. */
    Crv1X = CagdCrvCopy(Crv1X);
    Crv1W = CagdCrvCopy(Crv1W);
    Crv2X = CagdCrvCopy(Crv2X);
    Crv2W = CagdCrvCopy(Crv2W);
    if (!CagdMakeCrvsCompatible(&Crv1X, &Crv2X, FALSE, FALSE) ||
	!CagdMakeCrvsCompatible(&Crv1W, &Crv2W, FALSE, FALSE) ||
	!CagdMakeCrvsCompatible(&Crv1X, &Crv2W, FALSE, FALSE) ||
	!CagdMakeCrvsCompatible(&Crv1W, &Crv2X, FALSE, FALSE))
	SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT);

    CTmp1 = SymbCrvMult(Crv1X, Crv2W);
    CTmp2 = SymbCrvMult(Crv2X, Crv1W);
    CTmp3 = SymbCrvAddSubAux(CTmp1, CTmp2, OperationAdd);
    CagdCrvFree(CTmp1);
    CagdCrvFree(CTmp2);

    CagdCrvFree(Crv1X);
    CagdCrvFree(Crv1W);
    CagdCrvFree(Crv2X);
    CagdCrvFree(Crv2W);

    return CTmp3;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Given two curve - add or subtract them as prescribed by OperationAdd,	     *
* coordinatewise.			      				     *
*   The two curves are promoted to same type, point type, and order before   *
* addition can take place.						     *
*   Returned is a curve representing their sum or difference.		     *
*                                                                            *
* PARAMETERS:                                                                *
*   Crv1, Crv2:    Two curve to subtract coordinatewise.                     *
*   OperationAdd:  TRUE of addition, FALSE for subtraction.                  *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:   The summation or difference of Crv1 and Crv2.         *
*****************************************************************************/
static CagdCrvStruct *SymbCrvAddSubAux(CagdCrvStruct *Crv1,
				       CagdCrvStruct *Crv2,
				       CagdBType OperationAdd)
{
    CagdBType IsNotRational;
    int i, k;
    CagdCrvStruct *SumCrv;
    CagdRType **Points1, **Points2;

    /* Make the two curves - same order and point type. */
    Crv1 = CagdCrvCopy(Crv1);
    Crv2 = CagdCrvCopy(Crv2);
    if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, TRUE, TRUE))
	SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT);

    SumCrv = CagdCrvNew(Crv1 -> GType, Crv1 -> PType, k = Crv1 -> Length);
    SumCrv -> Order = Crv1 -> Order;
    if (CAGD_IS_BSPLINE_CRV(SumCrv)) {
	SumCrv->KnotVector = BspKnotCopy(Crv1 -> KnotVector,
					 Crv1 -> Length + Crv1 -> Order);
    }

    IsNotRational = !CAGD_IS_RATIONAL_CRV(SumCrv);
    Points1 = Crv1 -> Points;
    Points2 = Crv2 -> Points;

    if (IsNotRational) {
	/* Simply add the respective control polygons. */
	SymbMeshAddSub(SumCrv -> Points, Crv1 -> Points, Crv2 -> Points,
		       SumCrv -> PType, SumCrv -> Length, OperationAdd);
    }
    else {
	/* Maybe the weights are identical, in which we can still add the   */
	/* the respective control polygons.				    */
	for (i = 0; i < k; i++)
	    if (!APX_EQ(Points1[0][i], Points2[0][i])) break;

	if (i < k) {
	    CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z,
			  *Crv2W, *Crv2X, *Crv2Y, *Crv2Z,
			  *SumCrvW, *SumCrvX, *SumCrvY, *SumCrvZ;

	    /* Weights are different. Must use the addition of rationals    */
	    /* rule ( we invoke SymbCrvMult here):			    */
	    /*								    */
	    /*  x1     x2   x1 w2 +/- x2 w1				    */
	    /*  -- +/- -- = ---------------				    */
	    /*  w1     w2        w1 w2					    */
	    /*								    */
	    SymbCrvSplitScalar(Crv1, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z);
	    SymbCrvSplitScalar(Crv2, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z);

	    SumCrvW = SymbCrvMult(Crv1W, Crv2W);
	    SumCrvX = SymbCrvRtnlMult(Crv1X, Crv1W, Crv2X, Crv2W, OperationAdd);
	    SumCrvY = SymbCrvRtnlMult(Crv1Y, Crv1W, Crv2Y, Crv2W, OperationAdd);
	    SumCrvZ = SymbCrvRtnlMult(Crv1Z, Crv1W, Crv2Z, Crv2W, OperationAdd);
	    CagdCrvFree(Crv1W);
	    CagdCrvFree(Crv1X);
	    CagdCrvFree(Crv1Y);
	    CagdCrvFree(Crv1Z);
	    CagdCrvFree(Crv2W);
	    CagdCrvFree(Crv2X);
	    CagdCrvFree(Crv2Y);
	    CagdCrvFree(Crv2Z);

	    CagdCrvFree(SumCrv);
	    SumCrv = SymbCrvMergeScalar(SumCrvW, SumCrvX, SumCrvY, SumCrvZ);
	    CagdCrvFree(SumCrvW);
	    CagdCrvFree(SumCrvX);
	    CagdCrvFree(SumCrvY);
	    CagdCrvFree(SumCrvZ);
	}
	else
	{
	    /* Simply add respective control polygons (w is just copied). */
	    SymbMeshAddSub(SumCrv -> Points, Crv1 -> Points, Crv2 -> Points,
			   SumCrv -> PType, SumCrv -> Length, OperationAdd);
	}
    }

    CagdCrvFree(Crv1);
    CagdCrvFree(Crv2);

    return SumCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a planar curve, compute its enclosed area field curve.		     M
*   This has little meaning unless Crv is closed, in which by evaluation     M
* the resulting area filed curve at the end points, the area enclosed by     M
* Crv can be computed.                                                       M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       A curve to compute area filed curve for.			     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   The area field curve.                                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvEnclosedArea, area, symbolic computation                          M
*****************************************************************************/
CagdCrvStruct *SymbCrvEnclosedArea(CagdCrvStruct *Crv)
{
    CagdVType Trans;
    CagdCrvStruct *CrvX, *CrvY, *CrvZ, *CrvW, *DCrvX, *DCrvY,
	*CTmp1, *CTmp2, *CTmp3;

    switch (Crv -> GType) {
	case CAGD_CBEZIER_TYPE:
	case CAGD_CBSPLINE_TYPE:
	    break;
	case CAGD_CPOWER_TYPE:
	    SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
	    return NULL;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
	    return NULL;
    }

    /* Using Green's theorem, if C(t1) = C(t2) the area enclosed is equal to */
    /*									     */
    /*	   t2								     */
    /*	    /								     */
    /*	1  |								     */
    /*  -  | -x'(t) y(t) + x(t) y'(t) dt				     */
    /*	2  |								     */
    /*    /								     */
    /*    t1								     */
    /*									     */
    SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
    if (CrvZ)
        CagdCrvFree(CrvZ);
    if (CrvW) {
	SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_NO_SUPPORT);
	CagdCrvFree(CrvW);
    }

    DCrvX = CagdCrvDerive(CrvX);
    DCrvY = CagdCrvDerive(CrvY);
    CTmp1 = SymbCrvMult(CrvX, DCrvY);
    CTmp2 = SymbCrvMult(DCrvX, CrvY);
    CagdCrvFree(CrvX);
    CagdCrvFree(CrvY);
    CagdCrvFree(DCrvX);
    CagdCrvFree(DCrvY);
    CTmp3 = SymbCrvSub(CTmp1, CTmp2);
    CagdCrvFree(CTmp1);
    CagdCrvFree(CTmp2);
    CTmp1 = CagdCrvIntegrate(CTmp3);
    CagdCrvFree(CTmp3);

    /* Scale by a half. */
    Trans[0] = Trans[1] = Trans[2] = 0.0;
    CagdCrvTransform(CTmp1, Trans, 0.5);

    return CTmp1;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two surfaces - add them coordinatewise.				     M
*   The two surfaces are promoted to same point type before the		     M
* multiplication can take place. Furthermore, order and continuity are	     M
* matched as well.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf1, Srf2:  Two surface to add up coordinatewise.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   The summation of Srf1 + Srf2 coordinatewise.          M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfAdd, addition, symbolic computation                               M
*****************************************************************************/
CagdSrfStruct *SymbSrfAdd(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    return SymbSrfAddSubAux(Srf1, Srf2, TRUE);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two surfaces - subtract them coordinatewise.			     M
*   The two surfaces are promoted to same point type before the		     M
* multiplication can take place. Furthermore, order and continuity are	     M
* matched as well.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf1, Srf2:  Two surface to subtract coordinatewise.                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   The difference of Srf1 - Srf2 coordinatewise.         M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfSub, subtraction, symbolic computation                            M
*****************************************************************************/
CagdSrfStruct *SymbSrfSub(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    return SymbSrfAddSubAux(Srf1, Srf2, FALSE);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two surfaces - multiply them coordinatewise.			     M
*   The two surfaces are promoted to same point type before the 	     M
* multiplication can take place.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf1, Srf2:  Two surface to multiply coordinatewise.                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   The product of Srf1 * Srf2 coordinatewise.            M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfMult, product, symbolic computation                               M
*****************************************************************************/
CagdSrfStruct *SymbSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdSrfStruct
	*ProdSrf = NULL;

    Srf1 = CagdSrfCopy(Srf1);
    Srf2 = CagdSrfCopy(Srf2);
    if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE))
	SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT);

    switch (Srf1 -> GType) {
	case CAGD_SBEZIER_TYPE:
	    ProdSrf = BzrSrfMult(Srf1, Srf2);
	    break;
	case CAGD_SBSPLINE_TYPE:
	    ProdSrf = BspSrfMult(Srf1, Srf2);
	    break;
	case CAGD_SPOWER_TYPE:
	    SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_SRF);
	    break;
    }

    CagdSrfFree(Srf1);
    CagdSrfFree(Srf2);

    return ProdSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a scalar surface, returns a scalar surface representing the	     M
* reciprocal values, by making it rational (if was not one) and flipping     M
* the numerator and the denominator.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       A scalar surface to compute a reciprocal value for.           M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   A rational scalar surface that is equalto the         M
*                      reciprocal value of Srf.                              M
* KEYWORDS:                                                                  M
*   SymbSrfInvert, division, symbolic computation, reciprocal value          M
*****************************************************************************/
CagdSrfStruct *SymbSrfInvert(CagdSrfStruct *Srf)
{
    int i;
    CagdRType *R, **Points;

    Srf = CagdSrfCopy(Srf);
    Points = Srf -> Points;

    switch (Srf -> PType) {
	case CAGD_PT_E1_TYPE:
	    Points[0] = Points[1];
	    Points[1] = R = (CagdRType *) IritMalloc(sizeof(CagdRType) *
					     Srf -> ULength * Srf -> VLength);
	    for (i = 0; i < Srf -> ULength * Srf -> VLength; i++)
		*R++ = 1.0;
	    Srf -> PType = CAGD_PT_P1_TYPE;
	    break;
	case CAGD_PT_P1_TYPE:
	    R = Points[0];
	    Points[0] = Points[1];
	    Points[1] = R;
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNSUPPORT_PT);
	    break;
    }

    return Srf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a surface, scale it by Scale.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:      A surface to scale by magnitude Scale.                         M
*   Scale:    Scaling factor.                                                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   A surfaces scaled by Scale compared to Srf.           M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfScalarScale, scaling, symbolic computation                        M
*****************************************************************************/
CagdSrfStruct *SymbSrfScalarScale(CagdSrfStruct *Srf, CagdRType Scale)
{
    int i;
    CagdRType *R;

    Srf = CagdSrfCopy(Srf);

    switch (Srf -> PType) {
	case CAGD_PT_P1_TYPE:
	case CAGD_PT_E1_TYPE:
	    for (i = 0, R = Srf -> Points[1];
		 i < Srf -> ULength * Srf -> VLength;
		 i++)
		*R++ *= Scale;
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNSUPPORT_PT);
	    break;
    }

    return Srf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two surfaces - computes their dot product.			     M
*   Returned surface is a scalar surface representing the dot product of the M
* two given surfaces.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf1, Srf2:  Two surface to multiply and compute a dot product for.      M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   A scalar surface representing the dot product of      M
*                      Srf1 . Srf2.                                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfDotProd, product, dot product, symbolic computation               M
*****************************************************************************/
CagdSrfStruct *SymbSrfDotProd(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdSrfStruct *PSrfW, *PSrfX, *PSrfY, *PSrfZ, *TSrf1, *TSrf2, *DotProdSrf,
	*ProdSrf = SymbSrfMult(Srf1, Srf2);

    SymbSrfSplitScalar(ProdSrf, &PSrfW, &PSrfX, &PSrfY, &PSrfZ);
    CagdSrfFree(ProdSrf);

    if (PSrfY != NULL) {
	TSrf1 = SymbSrfAdd(PSrfX, PSrfY);
	CagdSrfFree(PSrfX);
	CagdSrfFree(PSrfY);
    }
    else
	TSrf1 = PSrfX;

    if (PSrfZ != NULL) {
	TSrf2 = SymbSrfAdd(TSrf1, PSrfZ);
	CagdSrfFree(TSrf1);
	CagdSrfFree(PSrfZ);
	TSrf1 = TSrf2;
    }

    DotProdSrf = SymbSrfMergeScalar(PSrfW, TSrf1, NULL, NULL);
    CagdSrfFree(PSrfW);
    CagdSrfFree(TSrf1);

    return DotProdSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two surfaces - computes their cross product.			     M
*   Returned surface is a scalar surface representing the cross product of   M
* the two given surfaces.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf1, Srf2:  Two surface to multiply and compute a cross product for.    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   A scalar surface representing the cross product of    M
*                      Srf1 x Srf2.                                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfCrossProd, product, cross product, symbolic computation           M
*****************************************************************************/
CagdSrfStruct *SymbSrfCrossProd(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdSrfStruct *Srf1W, *Srf1X, *Srf1Y, *Srf1Z,
		  *Srf2W, *Srf2X, *Srf2Y, *Srf2Z,
		  *TSrf1, *TSrf2, *CrossProdSrf,
	*PSrfW = NULL,
	*PSrfX = NULL,
	*PSrfY = NULL,
	*PSrfZ = NULL;

    SymbSrfSplitScalar(Srf1, &Srf1W, &Srf1X, &Srf1Y, &Srf1Z);
    SymbSrfSplitScalar(Srf2, &Srf2W, &Srf2X, &Srf2Y, &Srf2Z);

    if (Srf1X == NULL || Srf1Y == NULL || Srf2X == NULL || Srf2Y == NULL)
	SYMB_FATAL_ERROR(SYMB_ERR_NO_CROSS_PROD);

    /* Cross product X axis. */
    TSrf1 = Srf2Z ? SymbSrfMult(Srf1Y, Srf2Z) : NULL;
    TSrf2 = Srf1Z ? SymbSrfMult(Srf2Y, Srf1Z) : NULL;

    if (TSrf1) {
	if (TSrf2) {
	    PSrfX = SymbSrfSub(TSrf1, TSrf2);
	    CagdSrfFree(TSrf2);
	}
	CagdSrfFree(TSrf1);
    }

    /* Cross product Y axis. */
    TSrf1 = Srf1Z ? SymbSrfMult(Srf1Z, Srf2X) : NULL;
    TSrf2 = Srf2Z ? SymbSrfMult(Srf2Z, Srf1X) : NULL;
    if (TSrf1) {
	if (TSrf2) {
	    PSrfY = SymbSrfSub(TSrf1, TSrf2);
	    CagdSrfFree(TSrf2);
	}
	CagdSrfFree(TSrf1);
    }

    /* Cross product Z axis. */
    TSrf1 = SymbSrfMult(Srf1X, Srf2Y);
    TSrf2 = SymbSrfMult(Srf2X, Srf1Y);
    PSrfZ = SymbSrfSub(TSrf1, TSrf2);
    CagdSrfFree(TSrf1);
    CagdSrfFree(TSrf2);

    /* Cross product W axis. */
    if (Srf1W || Srf2W) {
	if (Srf1W == NULL)
	    PSrfW = CagdSrfCopy(Srf2W);
	else if (Srf2W == NULL)
	    PSrfW = CagdSrfCopy(Srf1W);
	else
	    PSrfW = SymbSrfMult(Srf1W, Srf2W);
    }

    CagdSrfFree(Srf1X);
    CagdSrfFree(Srf1Y);
    CagdSrfFree(Srf1Z);
    CagdSrfFree(Srf1W);

    CagdSrfFree(Srf2X);
    CagdSrfFree(Srf2Y);
    CagdSrfFree(Srf2Z);
    CagdSrfFree(Srf2W);

    if (!CagdMakeSrfsCompatible(&PSrfW, &PSrfX, TRUE, TRUE, TRUE, TRUE) ||
	!CagdMakeSrfsCompatible(&PSrfW, &PSrfY, TRUE, TRUE, TRUE, TRUE) ||
	!CagdMakeSrfsCompatible(&PSrfW, &PSrfZ, TRUE, TRUE, TRUE, TRUE))
	SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT);

    CrossProdSrf = SymbSrfMergeScalar(PSrfW, PSrfX, PSrfY, PSrfZ);
    CagdSrfFree(PSrfX);
    CagdSrfFree(PSrfY);
    CagdSrfFree(PSrfZ);
    CagdSrfFree(PSrfW);

    return CrossProdSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two surfaces - multiply them using the quotient product rule:	     M
*  X = X1 W2 +/- X2 W1							     V
*   All provided surfaces are assumed to be non rational scalar surfaces.    M
*   Returned is a non rational scalar surface (CAGD_PT_E1_TYPE).	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf1X:     Numerator of first surface.                                   M
*   Srf1W:     Denominator of first surface.                                 M
*   Srf2X:     Numerator of second surface.                                  M
*   Srf2W:     Denominator of second surface.                                M
*   OperationAdd:   TRUE for addition, FALSE for subtraction.                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  The result of  Srf1X Srf2W +/- Srf2X Srf1W.            M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfRtnlMult, product, symbolic computation                           M
*****************************************************************************/
CagdSrfStruct *SymbSrfRtnlMult(CagdSrfStruct *Srf1X,
			       CagdSrfStruct *Srf1W,
			       CagdSrfStruct *Srf2X,
			       CagdSrfStruct *Srf2W,
			       CagdBType OperationAdd)
{
    CagdSrfStruct *CTmp1, *CTmp2, *CTmp3;

    /* Make the two surfaces - same order and point type. */
    Srf1X = CagdSrfCopy(Srf1X);
    Srf1W = CagdSrfCopy(Srf1W);
    Srf2X = CagdSrfCopy(Srf2X);
    Srf2W = CagdSrfCopy(Srf2W);
    if (!CagdMakeSrfsCompatible(&Srf1X, &Srf2X, FALSE, FALSE, FALSE, FALSE) ||
	!CagdMakeSrfsCompatible(&Srf1W, &Srf2W, FALSE, FALSE, FALSE, FALSE) ||
	!CagdMakeSrfsCompatible(&Srf1X, &Srf2W, FALSE, FALSE, FALSE, FALSE) ||
	!CagdMakeSrfsCompatible(&Srf1W, &Srf2X, FALSE, FALSE, FALSE, FALSE))
	SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT);

    CTmp1 = SymbSrfMult(Srf1X, Srf2W);
    CTmp2 = SymbSrfMult(Srf2X, Srf1W);
    CTmp3 = SymbSrfAddSubAux(CTmp1, CTmp2, OperationAdd);
    CagdSrfFree(CTmp1);
    CagdSrfFree(CTmp2);

    CagdSrfFree(Srf1X);
    CagdSrfFree(Srf1W);
    CagdSrfFree(Srf2X);
    CagdSrfFree(Srf2W);

    return CTmp3;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a surface - compute its unnormalized norrmal vectorfield surface, as M
* the cross product if this partial derivatives.                             M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       To compute an unnormalized normal vector field for.           M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   A vector field representing the unnormalized normal   M
*                      vector field of Srf.                                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfNormalSrf, normal, symbolic computation                           M
*****************************************************************************/
CagdSrfStruct *SymbSrfNormalSrf(CagdSrfStruct *Srf)
{
    CagdSrfStruct
	*SrfDU = CagdSrfDerive(Srf, CAGD_CONST_U_DIR),
	*SrfDV = CagdSrfDerive(Srf, CAGD_CONST_V_DIR),
	*NormalSrf = SymbSrfCrossProd(SrfDV, SrfDU);

    CagdSrfFree(SrfDU);
    CagdSrfFree(SrfDV);

    return NormalSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Given two surface - add or subtract them as prescribed by OperationAdd,    *
* coordinatewise.			      				     *
*   The two surfaces are promoted to same type, point type, and order before *
* addition can take place.						     *
*   Returned is a surface representing their sum or difference.		     *
*                                                                            *
* PARAMETERS:                                                                *
*   Srf1, Srf2:    Two surface to subtract coordinatewise.                   *
*   OperationAdd:  TRUE of addition, FALSE for subtraction.                  *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdSrfStruct *:   The summation or difference of Srf1 and Srf2.         *
*****************************************************************************/
static CagdSrfStruct *SymbSrfAddSubAux(CagdSrfStruct *Srf1,
				       CagdSrfStruct *Srf2,
				       CagdBType OperationAdd)
{

    CagdBType IsNotRational;
    int i, Len;
    CagdSrfStruct *SumSrf;
    CagdRType **Points1, **Points2;

    /* Make the two surfaces - same order and point type. */
    Srf1 = CagdSrfCopy(Srf1);
    Srf2 = CagdSrfCopy(Srf2);
    if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, TRUE, TRUE, TRUE, TRUE))
	SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT);

    Len = Srf1 -> ULength * Srf1 ->VLength;
    SumSrf = CagdSrfNew(Srf1 -> GType, Srf1 -> PType,
			Srf1 -> ULength, Srf1 ->VLength);
    SumSrf -> UOrder = Srf1 -> UOrder;
    SumSrf -> VOrder = Srf1 -> VOrder;
    if (CAGD_IS_BSPLINE_SRF(SumSrf)) {
	SumSrf -> UKnotVector = BspKnotCopy(Srf1 -> UKnotVector,
					    Srf1 -> ULength + Srf1 -> UOrder);
	SumSrf -> VKnotVector = BspKnotCopy(Srf1 -> VKnotVector,
					    Srf1 -> VLength + Srf1 -> VOrder);
    }

    IsNotRational = !CAGD_IS_RATIONAL_SRF(SumSrf);
    Points1 = Srf1 -> Points;
    Points2 = Srf2 -> Points;

    if (IsNotRational) {
	/* Simply add the respective control polygons. */
	SymbMeshAddSub(SumSrf -> Points, Srf1 -> Points, Srf2 -> Points,
		       SumSrf -> PType, Len, OperationAdd);
    }
    else {
	/* Maybe the weights are identical, in which we can still add the   */
	/* the respective control polygons.				    */
	for (i = 0; i < Len; i++)
	    if (!APX_EQ(Points1[0][i], Points2[0][i])) break;

	if (i < Len) {
	    CagdSrfStruct *Srf1W, *Srf1X, *Srf1Y, *Srf1Z,
			  *Srf2W, *Srf2X, *Srf2Y, *Srf2Z,
			  *SumSrfW, *SumSrfX, *SumSrfY, *SumSrfZ;

	    /* Weights are different. Must use the addition of rationals    */
	    /* rule ( we invoke SymbSrfMult here):			    */
	    /*								    */
	    /*  x1     x2   x1 w2 +/- x2 w1				    */
	    /*  -- +/- -- = ---------------				    */
	    /*  w1     w2        w1 w2					    */
	    /*								    */
	    SymbSrfSplitScalar(Srf1, &Srf1W, &Srf1X, &Srf1Y, &Srf1Z);
	    SymbSrfSplitScalar(Srf2, &Srf2W, &Srf2X, &Srf2Y, &Srf2Z);

	    SumSrfW = SymbSrfMult(Srf1W, Srf2W);
	    SumSrfX = SymbSrfRtnlMult(Srf1X, Srf1W, Srf2X, Srf2W, OperationAdd);
	    SumSrfY = SymbSrfRtnlMult(Srf1Y, Srf1W, Srf2Y, Srf2W, OperationAdd);
	    SumSrfZ = SymbSrfRtnlMult(Srf1Z, Srf1W, Srf2Z, Srf2W, OperationAdd);
	    CagdSrfFree(Srf1W);
	    CagdSrfFree(Srf1X);
	    CagdSrfFree(Srf1Y);
	    CagdSrfFree(Srf1Z);
	    CagdSrfFree(Srf2W);
	    CagdSrfFree(Srf2X);
	    CagdSrfFree(Srf2Y);
	    CagdSrfFree(Srf2Z);

	    CagdSrfFree(SumSrf);
	    SumSrf = SymbSrfMergeScalar(SumSrfW, SumSrfX, SumSrfY, SumSrfZ);
	    CagdSrfFree(SumSrfW);
	    CagdSrfFree(SumSrfX);
	    CagdSrfFree(SumSrfY);
	    CagdSrfFree(SumSrfZ);
	}
	else
	{
	    /* Simply add respective control polygons (w is just copied). */
	    SymbMeshAddSub(SumSrf -> Points, Srf1 -> Points, Srf2 -> Points,
			   SumSrf -> PType, Len, OperationAdd);
	}
    }

    CagdSrfFree(Srf1);
    CagdSrfFree(Srf2);

    return SumSrf;
}

/******************************************************************************
* DESCRIPTION:                                                               M
* Given two control polygons/meshes - add them coordinatewise.		     M
* If mesh is rational, weights are assumed identical and are just copied.    M
*                                                                            *
* PARAMETERS:                                                                M
*   DestPoints:    Where addition or difference result should go to.         M
*   Points1:       First control polygon/mesh.                               M
*   Points2:       Second control polygon/mesh.                              M
*   PType:         Type of points we are dealing with.                       M
*   Size:          Length of each vector in Points1/2.                       M
*   OperationAdd:  TRUE of addition, FALSE for subtraction.                  M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbMeshAddSub, addition, subtraction, symbolic computation              M
*****************************************************************************/
void SymbMeshAddSub(CagdRType **DestPoints,
		    CagdRType **Points1,
		    CagdRType **Points2,
		    CagdPointType PType,
		    int Size,
		    CagdBType OperationAdd)
{
    CagdBType
	IsRational = CAGD_IS_RATIONAL_PT(PType);
    int i, j,
	NumCoords = CAGD_NUM_OF_PT_COORD(PType);

    for (i = 1; i <= NumCoords; i++) {
	CagdRType
	    *DPts = DestPoints[i],
	    *Pts1 = Points1[i],
	    *Pts2 = Points2[i];

	for (j = 0; j < Size; j++)
	    *DPts++ = OperationAdd ? *Pts1++ + *Pts2++ : *Pts1++ - *Pts2++;
    }

    if (IsRational) {     /* Copy the weights (should be identical in both). */
	CagdRType
	    *DPts = DestPoints[0],
	    *Pts1 = Points1[0],
	    *Pts2 = Points2[0];

	for (j = 0; j < Size; j++) {
	    if (!APX_EQ(*Pts1, *Pts2))
		SYMB_FATAL_ERROR(SYMB_ERR_W_NOT_SAME);
	    *DPts++ = *Pts1++;
	    Pts2++;
	}
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a curve splits it to its scalar component curves. Ignores all	     M
* dimensions beyond the third, Z, dimension.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:      Curve to split.                                                M
*   CrvW:     The weight component of Crv, if have any.                      M
*   CrvX:     The X component of Crv.                                        M
*   CrvY:     The Y component of Crv, if have any.                           M
*   CrvZ:     The Z component of Crv, if have any.                           M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvSplitScalar, split, symbolic computation                          M
*****************************************************************************/
void SymbCrvSplitScalar(CagdCrvStruct *Crv,
			CagdCrvStruct **CrvW,
			CagdCrvStruct **CrvX,
			CagdCrvStruct **CrvY,
			CagdCrvStruct **CrvZ)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i,
	Length = Crv -> Length,
	NumCoords = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct
	*Crvs[CAGD_MAX_PT_SIZE];

    for (i = 0; i < CAGD_MAX_PT_SIZE; i++)
	Crvs[i] = NULL;

    for (i = IsNotRational; i <= NumCoords; i++) {
	Crvs[i] = CagdCrvNew(Crv -> GType, CAGD_PT_E1_TYPE, Length);
	Crvs[i] -> Order = Crv -> Order;
	if (Crv -> KnotVector != NULL)
	    Crvs[i] -> KnotVector = BspKnotCopy(Crv -> KnotVector,
						Crv -> Length + Crv -> Order);

	SYMB_GEN_COPY(Crvs[i] -> Points[1], Crv -> Points[i],
		      sizeof(CagdRType) * Length);
    }

    *CrvW = Crvs[0];
    *CrvX = Crvs[1];
    *CrvY = Crvs[2];
    *CrvZ = Crvs[3];
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a set of scalar curves, treat them as coordinates into a new curve.  M
*   Assumes at least CrvX is not NULL in which a scalar curve is returned.   M
*   Assumes CrvX/Y/Z/W are either E1 or P1 in which the weights are assumed  M
* to be identical and can be ignored if CrvW exists or copied otheriwse.     M
*                                                                            *
* PARAMETERS:                                                                M
*   CrvW:     The weight component of new constructed curve, if have any.    M
*   CrvX:     The X component of new constructed curve.                      M
*   CrvY:     The Y component of new constructed curve, if have any.         M
*   CrvZ:     The Z component of new constructed curve, if have any.         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A new curve constructed from given scalar curves.     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvMergeScalar, merge, symbolic computation                          M
*****************************************************************************/
CagdCrvStruct *SymbCrvMergeScalar(CagdCrvStruct *CrvW,
				  CagdCrvStruct *CrvX,
				  CagdCrvStruct *CrvY,
				  CagdCrvStruct *CrvZ)
{
    CagdBType
	WeightCopied = FALSE,
	IsRational = CrvW != NULL;
    int i,
	Length = CrvX -> Length,
	NumCoords = (CrvX != NULL) + (CrvY != NULL) + (CrvZ != NULL);
    CagdPointType
	PType = CAGD_MAKE_PT_TYPE(IsRational, NumCoords);
    CagdCrvStruct *Crvs[CAGD_MAX_PT_SIZE],
	*Crv = CagdCrvNew(CrvX -> GType, PType, Length);

    Crvs[0] = CrvW;
    Crvs[1] = CrvX;
    Crvs[2] = CrvY;
    Crvs[3] = CrvZ;

    Crv -> Order = CrvX -> Order;
    if (CrvX -> KnotVector != NULL)
	Crv -> KnotVector = BspKnotCopy(CrvX -> KnotVector,
					Length + CrvX -> Order);

    for (i = !IsRational; i <= NumCoords; i++) {
	if (Crvs[i] != NULL) {
	    if (Crvs[i] -> PType != CAGD_PT_E1_TYPE) {
		if (Crvs[i] -> PType != CAGD_PT_P1_TYPE)
		    SYMB_FATAL_ERROR(SYMB_ERR_SCALAR_EXPECTED);
		else if (CrvW == NULL && WeightCopied == FALSE) {
		    SYMB_GEN_COPY(Crv -> Points[0], Crvs[i] -> Points[0],
				  sizeof(CagdRType) * Length);
		    WeightCopied = TRUE;
		}
	    }

	    SYMB_GEN_COPY(Crv -> Points[i], Crvs[i] -> Points[1],
			  sizeof(CagdRType) * Length);
	}
    }

    return Crv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a surface splits it to its scalar component surfaces. Ignores all    M
* dimensions beyond the third, Z, dimension.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:      Surface to split.                                              M
*   SrfW:     The weight component of Srf, if have any.                      M
*   SrfX:     The X component of Srf.                                        M
*   SrfY:     The Y component of Srf, if have any.                           M
*   SrfZ:     The Z component of Srf, if have any.                           M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfSplitScalar, split, symbolic computation                          M
*****************************************************************************/
void SymbSrfSplitScalar(CagdSrfStruct *Srf,
			CagdSrfStruct **SrfW,
			CagdSrfStruct **SrfX,
			CagdSrfStruct **SrfY,
			CagdSrfStruct **SrfZ)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int i,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength,
	NumCoords = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdSrfStruct
	*Srfs[CAGD_MAX_PT_SIZE];

    for (i = 0; i < CAGD_MAX_PT_SIZE; i++)
	Srfs[i] = NULL;

    for (i = IsNotRational; i <= NumCoords; i++) {
	Srfs[i] = CagdSrfNew(Srf -> GType, CAGD_PT_E1_TYPE, ULength, VLength);
	Srfs[i] -> UOrder = Srf -> UOrder;
	Srfs[i] -> VOrder = Srf -> VOrder;
	if (Srf -> UKnotVector != NULL)
	    Srfs[i] -> UKnotVector = BspKnotCopy(Srf -> UKnotVector,
						 ULength + Srf -> UOrder);
	if (Srf -> VKnotVector != NULL)
	    Srfs[i] -> VKnotVector = BspKnotCopy(Srf -> VKnotVector,
						 VLength + Srf -> VOrder);

	SYMB_GEN_COPY(Srfs[i] -> Points[1], Srf -> Points[i],
		      sizeof(CagdRType) * Srf -> ULength * Srf -> VLength);
    }

    *SrfW = Srfs[0];
    *SrfX = Srfs[1];
    *SrfY = Srfs[2];
    *SrfZ = Srfs[3];
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a set of scalar surfaces, treat them as coordinates into a new	     M
* surface.								     M
*   Assumes at least SrfX is not NULL in which a scalar surface is returned.  M
*   Assumes SrfX/Y/Z/W are either E1 or P1 in which the weights are assumed  M
* to be identical and can be ignored if SrfW exists or copied otheriwse.     M
*                                                                            *
* PARAMETERS:                                                                M
*   SrfW:     The weight component of new constructed surface, if have any.  M
*   SrfX:     The X component of new constructed surface.                    M
*   SrfY:     The Y component of new constructed surface, if have any.       M
*   SrfZ:     The Z component of new constructed surface, if have any.       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   A new surface constructed from given scalar surfaces. M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfMergeScalar, merge, symbolic computation                          M
*****************************************************************************/
CagdSrfStruct *SymbSrfMergeScalar(CagdSrfStruct *SrfW,
				  CagdSrfStruct *SrfX,
				  CagdSrfStruct *SrfY,
				  CagdSrfStruct *SrfZ)
{
    CagdBType
	WeightCopied = FALSE,
	IsRational = SrfW != NULL;
    int i,
	ULength = SrfX -> ULength,
	VLength = SrfX -> VLength,
	NumCoords = (SrfX != NULL) + (SrfY != NULL) + (SrfZ != NULL);
    CagdPointType
	PType = CAGD_MAKE_PT_TYPE(IsRational, NumCoords);
    CagdSrfStruct *Srfs[CAGD_MAX_PT_SIZE],
	*Srf = CagdSrfNew(SrfX -> GType, PType, ULength, VLength);

    Srfs[0] = SrfW;
    Srfs[1] = SrfX;
    Srfs[2] = SrfY;
    Srfs[3] = SrfZ;

    Srf -> UOrder = SrfX -> UOrder;
    Srf -> VOrder = SrfX -> VOrder;
    if (SrfX -> UKnotVector != NULL)
	Srf -> UKnotVector = BspKnotCopy(SrfX -> UKnotVector,
					 ULength + SrfX -> UOrder);
    if (SrfX -> VKnotVector != NULL)
	Srf -> VKnotVector = BspKnotCopy(SrfX -> VKnotVector,
					 VLength + SrfX -> VOrder);

    for (i = !IsRational; i <= NumCoords; i++) {
	if (Srfs[i] != NULL) {
	    if (Srfs[i] -> PType != CAGD_PT_E1_TYPE) {
		if (Srfs[i] -> PType != CAGD_PT_P1_TYPE)
		    SYMB_FATAL_ERROR(SYMB_ERR_SCALAR_EXPECTED);
		else if (SrfW == NULL && WeightCopied == FALSE) {
		    SYMB_GEN_COPY(Srf -> Points[0], Srfs[i] -> Points[0],
				  sizeof(CagdRType) * ULength * VLength);
		    WeightCopied = TRUE;
		}
	    }

	    SYMB_GEN_COPY(Srf -> Points[i], Srfs[i] -> Points[1],
			  sizeof(CagdRType) * ULength * VLength);
	}
    }

    return Srf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Promote a scalar curve to two dimensions by moving the scalar axis to be   M
* the Y axis and adding monotone X axis.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:      Scalar curve to promose to a two dimensional one.              M
*   Min:      Minimum of new monotone X axis.                                M
*   Max:      Maximum of new monotone X axis.                                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A two dimensional curve.                               M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbPrmtSclrCrvTo2D, promotion, conversion, symbolic computation         M
*****************************************************************************/
CagdCrvStruct *SymbPrmtSclrCrvTo2D(CagdCrvStruct *Crv,
				   CagdRType Min,
				   CagdRType Max)
{
    int i,
	Length = Crv -> Length;
    CagdBType
	IsRational = CAGD_IS_RATIONAL_CRV(Crv);
    CagdRType *R, *Points, *WPoints,
	Step = (Max - Min) / (Length - 1);
    CagdCrvStruct
	*PrmtCrv = CagdCoerceCrvTo(Crv, IsRational ? CAGD_PT_P2_TYPE
						   : CAGD_PT_E2_TYPE);

    R = PrmtCrv -> Points[1];
    PrmtCrv -> Points[1] = PrmtCrv -> Points[2];
    PrmtCrv -> Points[2] = R;

    Points = PrmtCrv -> Points[1];
    WPoints = IsRational ? PrmtCrv -> Points[0] : NULL;
    for (i = 0; i < Length; i++)
	*Points++ = (Min + Step * i) * (IsRational ? *WPoints++ : 1.0);

    return PrmtCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Promote a scalar surface to three dimensions by moving the scalar axis to  M
* be the Z axis and adding monotone X and Y axes.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:                                                                     M
*   UMin:     Minimum of new monotone X axis.                                M
*   UMax:     Maximum of new monotone X axis.                                M
*   VMin:     Minimum of new monotone Y axis.                                M
*   VMax:     Maximum of new monotone X axis.                                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:    A three dimensional surface.                         M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbPrmtSclrSrfTo3D, promotion, conversion, symbolic computation         M
*****************************************************************************/
CagdSrfStruct *SymbPrmtSclrSrfTo3D(CagdSrfStruct *Srf,
				   CagdRType UMin,
				   CagdRType UMax,
				   CagdRType VMin,
				   CagdRType VMax)
{
    int i, j,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength;
    CagdBType
	IsRational = CAGD_IS_RATIONAL_SRF(Srf);
    CagdRType *R, *Points, *WPoints,
	UStep = (UMax - UMin) / (ULength - 1),
	VStep = (VMax - VMin) / (VLength - 1);
    CagdSrfStruct
	*PrmtSrf = CagdCoerceSrfTo(Srf, IsRational ? CAGD_PT_P3_TYPE
						   : CAGD_PT_E3_TYPE);

    R = PrmtSrf -> Points[1];
    PrmtSrf -> Points[1] = PrmtSrf -> Points[3];
    PrmtSrf -> Points[3] = R;

    Points = PrmtSrf -> Points[1];
    WPoints = IsRational ? PrmtSrf -> Points[0] : NULL;
    for (j = 0; j < VLength; j++)
	for (i = 0; i < ULength; i++)
	    *Points++ = (UMin + UStep * i) *
		(IsRational ? *WPoints++ : 1.0);

    Points = PrmtSrf -> Points[2];
    WPoints = IsRational ? PrmtSrf -> Points[0] : NULL;
    for (j = 0; j < VLength; j++)
	for (i = 0; i < ULength; i++)
	    *Points++ = (VMin + VStep * j) *
		(IsRational ? *WPoints++ : 1.0);

    return PrmtSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a controlpolygon/mesh, computes the extremum values of them all.     M
*                                                                            *
* PARAMETERS:                                                                M
*   Points:       To scan for extremum values.				     M
*   Length:       Length of each vector in Points.                           M
*   FindMinimum:  TRUE for minimum, FALSE for maximum.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:  A vector holding PType point with the extremum values of   M
*                 each axis independently.                                   M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbExtremumCntPtVals, extremum values, symbolic computation             M
*****************************************************************************/
CagdRType *SymbExtremumCntPtVals(CagdRType **Points,
				 int Length,
				 CagdBType FindMinimum)
{
    static CagdRType Extremum[CAGD_MAX_PT_SIZE];
    int i, j;

    for (i = 1; Points[i] != NULL && i <= CAGD_MAX_PT_COORD; i++) {
	CagdRType
	    *R = Points[0],
	    *P = Points[i];

	Extremum[i] = FindMinimum ? INFINITY : -INFINITY;

	for (j = 0; j < Length; j++) {
	    CagdRType
		Val = R == NULL ? *P++ : *P++ / *R++;

	    if (FindMinimum) {
		if (Extremum[i] > Val)
		    Extremum[i] = Val;
	    }
	    else {
		if (Extremum[i] < Val)
		    Extremum[i] = Val;
	    }
	}
    }

    return Extremum;
}

