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

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

static CagdCrvStruct *CagdCrvAddSubAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2,
				       CagdBType OperationAdd);

static CagdSrfStruct *CagdSrfAddSubAux(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2,
				       CagdBType OperationAdd);

/******************************************************************************
* Given two curves - add them coordinatewise.				      *
******************************************************************************/
CagdCrvStruct *CagdCrvAdd(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    return CagdCrvAddSubAux(Crv1, Crv2, TRUE);
}

/******************************************************************************
* Given two curves - subtract them coordinatewise.			      *
******************************************************************************/
CagdCrvStruct *CagdCrvSub(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    return CagdCrvAddSubAux(Crv1, Crv2, FALSE);
}

/******************************************************************************
* Given a scalar curve, returns a scalar curve representing reciprocal values.*
******************************************************************************/
CagdCrvStruct *CagdCrvReciprocal(CagdCrvStruct *Crv)
{
    CagdRType *Points;

    if (Crv -> PType == CAGD_PT_E1_TYPE) {
	int i;

	Crv -> Points[0] = Crv -> Points[1];
	Points = Crv -> Points[1] = (CagdRType *)
				IritMalloc(sizeof(CagdRType) * Crv -> Length);
	for (i = 0; i < Crv -> Length; i++)
	    *Points++ = 1.0;
	Crv -> PType = CAGD_PT_P1_TYPE;
	return Crv;
    }
    else if (Crv -> PType == CAGD_PT_P1_TYPE) {
	Crv = CagdCrvCopy(Crv);
	Points = Crv -> Points[0];
	Crv -> Points[0] = Crv -> Points[1];
	Crv -> Points[1] = Points;
	return Crv;
    }
    else {
	FATAL_ERROR(CAGD_ERR_WRONG_PT_TYPE);
	return NULL;
    }
}

/******************************************************************************
* Given two curves - multiply them coordinatewise.			      *
* The two curves are promoted to same point type before the multiplication    *
* can take place.							      *
* Return a curve representing their product.				      *
******************************************************************************/
CagdCrvStruct *CagdCrvMult(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:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
	    break;
    }

    return ProdCrv;
}

/******************************************************************************
* Given two curves - computes their dot product.			      *
* Returned curves is a scalar curves representing the dot product of the      *
* two given curves coordinatewise.					      *
******************************************************************************/
CagdCrvStruct *CagdCrvDotProd(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *PCrvW, *PCrvX, *PCrvY, *PCrvZ, *TCrv1, *TCrv2, *DotProdCrv,
	*ProdCrv = CagdCrvMult(Crv1, Crv2);

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

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

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

    DotProdCrv = CagdCrvMergeScalar(PCrvW, TCrv1, NULL, NULL);
    if (PCrvW != NULL)
	CagdCrvFree(PCrvW);
    CagdCrvFree(TCrv1);

    return DotProdCrv;
}

/******************************************************************************
* Given two curves - computes their cross product.			      *
* Returned curves is a vector field curve representing the cross product      *
* of the two given curves.						      *
******************************************************************************/
CagdCrvStruct *CagdCrvCrossProd(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z,
		  *Crv2W, *Crv2X, *Crv2Y, *Crv2Z,
		  *TCrv1, *TCrv2, *CrossProdCrv,
	*PCrvW = NULL,
	*PCrvX = NULL,
	*PCrvY = NULL,
	*PCrvZ = NULL;

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

    if (Crv1X == NULL || Crv1Y == NULL || Crv2X == NULL || Crv2Y == NULL)
	FATAL_ERROR(CAGD_ERR_NO_CROSS_PROD);

    /* Cross product X axis. */
    TCrv1 = Crv2Z ? CagdCrvMult(Crv1Y, Crv2Z) : NULL;
    TCrv2 = Crv1Z ? CagdCrvMult(Crv2Y, Crv1Z) : NULL;
    if (TCrv1) {
	if (TCrv2) {
	    PCrvX = CagdCrvSub(TCrv1, TCrv2);
	    CagdCrvFree(TCrv2);
	}
	CagdCrvFree(TCrv1);
    }

    /* Cross product Y axis. */
    TCrv1 = Crv1Z ? CagdCrvMult(Crv1Z, Crv2X) : NULL;
    TCrv2 = Crv2Z ? CagdCrvMult(Crv2Z, Crv1X) : NULL;
    if (TCrv1) {
	if (TCrv2) {
	    PCrvY = CagdCrvSub(TCrv1, TCrv2);
	    CagdCrvFree(TCrv2);
	}
	CagdCrvFree(TCrv1);
    }

    /* Cross product Z axis. */
    TCrv1 = CagdCrvMult(Crv1X, Crv2Y);
    TCrv2 = CagdCrvMult(Crv2X, Crv1Y);
    PCrvZ = CagdCrvSub(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 = CagdCrvMult(Crv1W, Crv2W);
    }

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

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

    if (!CagdMakeCrvsCompatible(&PCrvW, &PCrvX, TRUE, TRUE) ||
	!CagdMakeCrvsCompatible(&PCrvW, &PCrvY, TRUE, TRUE) ||
	!CagdMakeCrvsCompatible(&PCrvW, &PCrvZ, TRUE, TRUE))
	FATAL_ERROR(CAGD_ERR_CRV_FAIL_CMPT);

    CrossProdCrv = CagdCrvMergeScalar(PCrvW, PCrvX, PCrvY, PCrvZ);
    if (PCrvX)
	CagdCrvFree(PCrvX);
    if (PCrvY)
	CagdCrvFree(PCrvY);
    if (PCrvZ)
	CagdCrvFree(PCrvZ);
    if (PCrvW)
	CagdCrvFree(PCrvW);

    return CrossProdCrv;
}

/******************************************************************************
* Given two curves - multiply them using the quotient product rule:	      *
*  X = X1 W2 +/- X2 W1							      *
* All provided curves are assumed to be non rational scalar curves.	      *
* Returned in a non rational scalar curve (CAGD_PT_E1_TYPE).		      *
******************************************************************************/
CagdCrvStruct *CagdCrvRtnlMult(CagdCrvStruct *Crv1X, CagdCrvStruct *Crv1W,
			       CagdCrvStruct *Crv2X, CagdCrvStruct *Crv2W,
			       CagdBType OperationAdd)
{
    CagdCrvStruct *CTmp1, *CTmp2, *CTmp3;

    /* 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))
	FATAL_ERROR(CAGD_ERR_CRV_FAIL_CMPT);

    CTmp1 = CagdCrvMult(Crv1X, Crv2W);
    CTmp2 = CagdCrvMult(Crv2X, Crv1W);
    CTmp3 = CagdCrvAddSubAux(CTmp1, CTmp2, OperationAdd);
    CagdCrvFree(CTmp1);
    CagdCrvFree(CTmp2);

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

    return CTmp3;
}

/******************************************************************************
* Given two curve - add/sub them coordinatewise.			      *
* The two curves are promoted to same type, point type, and order before the  *
* addition can take place.						      *
* Return a curve representing their sum or difference.			      *
******************************************************************************/
static CagdCrvStruct *CagdCrvAddSubAux(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))
	FATAL_ERROR(CAGD_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. */
	CagdMeshAddSub(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 CagdCrvMult here):			    */
	    /*								    */
	    /*  x1     x2   x1 w2 +/- x2 w1				    */
	    /*  -- +/- -- = ---------------				    */
	    /*  w1     w2        w1 w2					    */
	    /*								    */
	    CagdCrvSplitScalar(Crv1, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z);
	    CagdCrvSplitScalar(Crv2, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z);

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

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

    CagdCrvFree(Crv1);
    CagdCrvFree(Crv2);

    return SumCrv;
}

/******************************************************************************
* Given two surfaces - add them coordinatewise.				      *
* The two surfaces are promoted to same point type and order before the       *
* addition can take place.						      *
* Return a surface representing their sum or difference.		      *
******************************************************************************/
CagdSrfStruct *CagdSrfAdd(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    return CagdSrfAddSubAux(Srf1, Srf2, TRUE);
}

/******************************************************************************
* Given two surfaces - subtract them coordinatewise.			      *
* The two surfaces are promoted to same point type and order before the       *
* subtraction can take place.						      *
* Return a surface representing their sum or difference.		      *
******************************************************************************/
CagdSrfStruct *CagdSrfSub(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    return CagdSrfAddSubAux(Srf1, Srf2, FALSE);
}

/******************************************************************************
* Given a scalar srf, returns a scalar srf representing reciprocal values.    *
******************************************************************************/
CagdSrfStruct *CagdSrfReciprocal(CagdSrfStruct *Srf)
{
    CagdRType *Points;

    if (Srf -> PType == CAGD_PT_E1_TYPE) {
	int i;

	Srf -> Points[0] = Srf -> Points[1];
	Points = Srf -> Points[1] = (CagdRType *)
				IritMalloc(sizeof(CagdRType) *
					   Srf -> ULength * Srf -> VLength);
	for (i = 0; i < Srf -> ULength * Srf -> VLength; i++)
	    *Points++ = 1.0;
	Srf -> PType = CAGD_PT_P1_TYPE;
	return Srf;
    }
    else if (Srf -> PType == CAGD_PT_P1_TYPE) {
	Srf = CagdSrfCopy(Srf);
	Points = Srf -> Points[0];
	Srf -> Points[0] = Srf -> Points[1];
	Srf -> Points[1] = Points;
	return Srf;
    }
    else {
	FATAL_ERROR(CAGD_ERR_WRONG_PT_TYPE);
	return NULL;
    }
}

/******************************************************************************
* Given two surfaces - multiply them coordinatewise.			      *
* The two surfaces are promoted to same point type before the multiplication  *
* can take place.							      *
* Return a surface representing their product.				      *
******************************************************************************/
CagdSrfStruct *CagdSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdSrfStruct
	*ProdSrf = NULL;

    Srf1 = CagdSrfCopy(Srf1);
    Srf2 = CagdSrfCopy(Srf2);
    if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE))
	FATAL_ERROR(CAGD_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:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_SRF);
	    break;
    }

    CagdSrfFree(Srf1);
    CagdSrfFree(Srf2);

    return ProdSrf;
}

/******************************************************************************
* Given two surfaces - computes their dot product.			      *
* Returned surface is a scalar surface representing the dot product of the    *
* two given surfaces coordinatewise.					      *
******************************************************************************/
CagdSrfStruct *CagdSrfDotProd(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdSrfStruct *PSrfW, *PSrfX, *PSrfY, *PSrfZ, *TSrf1, *TSrf2, *DotProdSrf,
	*ProdSrf = CagdSrfMult(Srf1, Srf2);

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

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

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

    DotProdSrf = CagdSrfMergeScalar(PSrfW, TSrf1, NULL, NULL);
    if (PSrfW != NULL)
	CagdSrfFree(PSrfW);
    CagdSrfFree(TSrf1);

    return DotProdSrf;
}

/******************************************************************************
* Given two surfaces - computes their cross product.			      *
* Returned surface is a vector field surface representing the cross product   *
* of the two given surfaces.						      *
******************************************************************************/
CagdSrfStruct *CagdSrfCrossProd(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdSrfStruct *Srf1W, *Srf1X, *Srf1Y, *Srf1Z,
		  *Srf2W, *Srf2X, *Srf2Y, *Srf2Z,
		  *TSrf1, *TSrf2, *CrossProdSrf,
	*PSrfW = NULL,
	*PSrfX = NULL,
	*PSrfY = NULL,
	*PSrfZ = NULL;

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

    if (Srf1X == NULL || Srf1Y == NULL || Srf2X == NULL || Srf2Y == NULL)
	FATAL_ERROR(CAGD_ERR_NO_CROSS_PROD);

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

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

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

    /* Cross product Z axis. */
    TSrf1 = CagdSrfMult(Srf1X, Srf2Y);
    TSrf2 = CagdSrfMult(Srf2X, Srf1Y);
    PSrfZ = CagdSrfSub(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 = CagdSrfMult(Srf1W, Srf2W);
    }

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

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

    if (!CagdMakeSrfsCompatible(&PSrfW, &PSrfX, TRUE, TRUE, TRUE, TRUE) ||
	!CagdMakeSrfsCompatible(&PSrfW, &PSrfY, TRUE, TRUE, TRUE, TRUE) ||
	!CagdMakeSrfsCompatible(&PSrfW, &PSrfZ, TRUE, TRUE, TRUE, TRUE))
	FATAL_ERROR(CAGD_ERR_SRF_FAIL_CMPT);

    CrossProdSrf = CagdSrfMergeScalar(PSrfW, PSrfX, PSrfY, PSrfZ);
    if (PSrfX)
	CagdSrfFree(PSrfX);
    if (PSrfY)
	CagdSrfFree(PSrfY);
    if (PSrfZ)
	CagdSrfFree(PSrfZ);
    if (PSrfW)
	CagdSrfFree(PSrfW);

    return CrossProdSrf;
}

/******************************************************************************
* Given two surfaces - multiply them using the quotient product rule:         *
*  X = X1 W2 +/- X2 W1							      *
* All provided surfaces are assumed to be non rational scalar surfaces.	      *
* Returned in a non rational scalar bezier surface (CAGD_PT_E1_TYPE).	      *
******************************************************************************/
CagdSrfStruct *CagdSrfRtnlMult(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))
	FATAL_ERROR(CAGD_ERR_SRF_FAIL_CMPT);

    CTmp1 = CagdSrfMult(Srf1X, Srf2W);
    CTmp2 = CagdSrfMult(Srf2X, Srf1W);
    CTmp3 = CagdSrfAddSubAux(CTmp1, CTmp2, OperationAdd);
    CagdSrfFree(CTmp1);
    CagdSrfFree(CTmp2);

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

    return CTmp3;
}

/******************************************************************************
* Given a surface - compute is normal surface.				      *
******************************************************************************/
CagdSrfStruct *CagdSrfNormalSrf(CagdSrfStruct *Srf)
{
    CagdSrfStruct
	*SrfDU = CagdSrfDerive(Srf, CAGD_CONST_U_DIR),
	*SrfDV = CagdSrfDerive(Srf, CAGD_CONST_V_DIR),
	*NormalSrf = CagdSrfCrossProd(SrfDV, SrfDU);

    CagdSrfFree(SrfDU);
    CagdSrfFree(SrfDV);

    return NormalSrf;
}

/******************************************************************************
* Given two surfaces - add/sub them coordinatewise.			      *
* The two surfaces are promoted to same point type and order before the       *
* addition/subtraction can take place.					      *
* Return a surface representing their sum or difference.		      *
******************************************************************************/
static CagdSrfStruct *CagdSrfAddSubAux(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))
	FATAL_ERROR(CAGD_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. */
	CagdMeshAddSub(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 CagdSrfMult here):			    */
	    /*								    */
	    /*  x1     x2   x1 w2 +/- x2 w1				    */
	    /*  -- +/- -- = ---------------				    */
	    /*  w1     w2        w1 w2					    */
	    /*								    */
	    CagdSrfSplitScalar(Srf1, &Srf1W, &Srf1X, &Srf1Y, &Srf1Z);
	    CagdSrfSplitScalar(Srf2, &Srf2W, &Srf2X, &Srf2Y, &Srf2Z);

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

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

    CagdSrfFree(Srf1);
    CagdSrfFree(Srf2);

    return SumSrf;
}

/******************************************************************************
* Given two control polygons/meshes - add them coordinatewise.		      *
* If mesh is rational, weights should be identical and are just copied.       *
******************************************************************************/
void CagdMeshAddSub(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))
		FATAL_ERROR(CAGD_ERR_W_NOT_SAME);
	    *DPts++ = *Pts1++;
	    Pts2++;
	}
    }
}

/******************************************************************************
* Given a curve splits it to its scalar component curves.		      *
******************************************************************************/
void CagdCrvSplitScalar(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);

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

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

/******************************************************************************
* Given a set of curve coordinates as scalar curves, assemble one curve.      *
* Assumes at least CrvX is not NULL (in which a scalar curve is returned.).   *
* Assumes CrvX/Y/Z/W are either E1 or P1 (in which the weights are assumed    *
* to be identical and can be ignored if CrvW exists or copied otheriwse).     *
******************************************************************************/
CagdCrvStruct *CagdCrvMergeScalar(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)
		    FATAL_ERROR(CAGD_ERR_SCALAR_EXPECTED);
		else if (CrvW == NULL && WeightCopied == FALSE) {
		    CAGD_GEN_COPY(Crv -> Points[0], Crvs[i] -> Points[0],
				  sizeof(CagdRType) * Length);
		    WeightCopied = TRUE;
		}
	    }

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

    return Crv;
}

/******************************************************************************
* Given a surface splits it to its scalar component surfaces.		      *
******************************************************************************/
void CagdSrfSplitScalar(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);

	CAGD_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];
}

/******************************************************************************
* Given a set of surface coordinates as scalar surfaces, assemble one surface.*
* Assumes at least SrfX is not NULL (in which a scalar surface is returned.). *
* Assumes SrfX/Y/Z/W are either E1 or P1 (in which the weights are assumed    *
* to be identical and can be ignored if SrfW exists or copied otheriwse).     *
******************************************************************************/
CagdSrfStruct *CagdSrfMergeScalar(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)
		    FATAL_ERROR(CAGD_ERR_SCALAR_EXPECTED);
		else if (SrfW == NULL && WeightCopied == FALSE) {
		    CAGD_GEN_COPY(Srf -> Points[0], Srfs[i] -> Points[0],
				  sizeof(CagdRType) * ULength * VLength);
		    WeightCopied = TRUE;
		}
	    }

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

    return Srf;
}

/******************************************************************************
* Promote a scalar curve to two dimensions by moving the scalar axis to be    *
* the Y axis and adding monotone X axis.				      *
******************************************************************************/
CagdCrvStruct *CagdPrmtSclrCrvTo2D(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;
}

/******************************************************************************
* Promote a scalar surface to three dimensions by moving the scalar axis to   *
* be the Z axis and adding monotone X and Y axis.			      *
******************************************************************************/
CagdSrfStruct *CagdPrmtSclrSrfTo3D(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;
}