/******************************************************************************
* Composit.c - Composition computation (symbolic).			      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, Apr. 93.					      *
******************************************************************************/

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

#define ZERO_SET_EPS 0.0001

static CagdCrvStruct *SymbComposeCrvCrvAux(CagdCrvStruct *Crv1,
					   CagdCrvStruct *Crv2);
static CagdCrvStruct *SymbComposeCrvCrvAux2(CagdCrvStruct *Crv1,
					    CagdCrvStruct *Crv2);
static CagdCrvStruct *SymbComposeSrfCrvAux(CagdSrfStruct *Srf,
					   CagdCrvStruct *Crv);
static CagdCrvStruct **SymbComputeCurvePowers(CagdCrvStruct *Crv, int Order);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two curves, Crv1 and Crv2, computes the composition Crv1(Crv2(t)).   M
*   Crv2 must be a scalar curve completely contained in Crv1's parametric    M
* domain.								     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:   The two curve to compose together.                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   The composed curve.                                   M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbComposeCrvCrv, composition                                           M
*****************************************************************************/
CagdCrvStruct *SymbComposeCrvCrv(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *CmpsCrv,
	*OrigCrv1 = Crv1,
	*OrigCrv2 = Crv2;
    CagdRType TMax, TMin, CTMax, CTMin;
    CagdBType
	DoBezier = Crv1 -> GType == CAGD_CBEZIER_TYPE &&
	           Crv2 -> GType == CAGD_CBEZIER_TYPE;

    switch (Crv1 -> GType) {
	case CAGD_CBEZIER_TYPE:
	    if (!DoBezier)
		Crv1 = CnvrtBezier2BsplineCrv(Crv1);
	    break;
	case CAGD_CBSPLINE_TYPE:
	    break;
	case CAGD_CPOWER_TYPE:
	    SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
	    break;
    }

    switch (Crv2 -> GType) {
	case CAGD_CBEZIER_TYPE:
	    if (!DoBezier)
		Crv2 = CnvrtBezier2BsplineCrv(Crv2);
	    break;
	case CAGD_CBSPLINE_TYPE:
	    break;
	case CAGD_CPOWER_TYPE:
	    SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
	    break;
    }

    CmpsCrv = SymbComposeCrvCrvAux(Crv1, Crv2);

    /* Now map the domain of the curve to be crv2 domain. */
    if (!DoBezier) {
	CagdCrvDomain(Crv2, &TMin, &TMax);
	CagdCrvDomain(CmpsCrv, &CTMin, &CTMax);
	if (CmpsCrv -> GType == CAGD_CBEZIER_TYPE) {
	    CagdCrvStruct
		*CTmp = CnvrtBezier2BsplineCrv(CmpsCrv);

	    CagdCrvFree(CmpsCrv);
	    CmpsCrv = CTmp;
	}
	BspKnotAffineTrans(CmpsCrv -> KnotVector,
			   CmpsCrv -> Length + CmpsCrv -> Order,
			   TMin - CTMin,
			   (TMax - TMin) / (CTMax - CTMin));
    }

    if (Crv1 != OrigCrv1)
	CagdCrvFree(Crv1);
    if (Crv2 != OrigCrv2)
	CagdCrvFree(Crv2);

    return CmpsCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Subdivides Crv2 until it is a Bezier curve, invokes the Bezier composition *
* code on each, and merges them back to a complete curve.		     *
*   At this point, curves can be either Bezier or Bspline only.		     *
*   Curves are both assumed to have open end condition.			     *
*                                                                            *
* PARAMETERS:                                                                *
*   Crv1, Crv2:   The two curve to compose together.                         *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:   The composed curve.                                   *
*****************************************************************************/
static CagdCrvStruct *SymbComposeCrvCrvAux(CagdCrvStruct *Crv1,
					   CagdCrvStruct *Crv2)
{
    CagdCrvStruct *CmpsCrv;

    if (Crv2 -> Length > Crv2 -> Order) {
	int i;
	CagdRType t, *R;
	CagdCrvStruct *Crv1a, *Crv1b, *Crv2a, *Crv2b, *CmpsCrvA, *CmpsCrvB;

	/* Validate that Crv2 is monotone. */
	R = Crv2 -> Points[1];
	for (i = 1; i < Crv2 -> Length; i++)
	    if (R[i-1] > R[i])
		SYMB_FATAL_ERROR(SYMB_ERR_REPARAM_NOT_MONOTONE);

	/* Crv2 is not a Bezier segment. Subdivide, compute for each segment */
	/* and merge back.						     */

	t = Crv2 -> KnotVector[(Crv2 -> Order + Crv2 -> Length) / 2];

	Crv2a = CagdCrvSubdivAtParam(Crv2, t);
	Crv2b = Crv2a -> Pnext;
	Crv2a -> Pnext = NULL;

	R = CagdCrvEval(Crv2, t);
	Crv1a = CagdCrvSubdivAtParam(Crv1, R[1]);
	Crv1b = Crv1a -> Pnext;
	Crv1a -> Pnext = NULL;

	CmpsCrvA = SymbComposeCrvCrvAux(Crv1a, Crv2a);
	CmpsCrvB = SymbComposeCrvCrvAux(Crv1b, Crv2b);
	CagdCrvFree(Crv1a);
	CagdCrvFree(Crv1b);
	CagdCrvFree(Crv2a);
	CagdCrvFree(Crv2b);

	CmpsCrv = CagdMergeCrvCrv(CmpsCrvA, CmpsCrvB, FALSE);
	CagdCrvFree(CmpsCrvA);
	CagdCrvFree(CmpsCrvB);
    }
    else {
    	/* Crv2 is a Bezier segment - compute its composition. */
    	CmpsCrv = SymbComposeCrvCrvAux2(Crv1, Crv2);
    }

    return CmpsCrv;    
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* See SymbComposeCrvCrvAux.                                                  *
*****************************************************************************/
static CagdCrvStruct *SymbComposeCrvCrvAux2(CagdCrvStruct *Crv1,
					    CagdCrvStruct *Crv2)
{
    CagdRType t, TMin, TMax, CTMax, CTMin;
    CagdCrvStruct *CmpsCrv;

    if (Crv1 -> Length > Crv1 -> Order) {
	/* Needs to compose in pieces after subdividing Crv2 at the interior */
	/* knots of Crv1, by finding where Crv2(r) = t, the interior knot.   */
	CagdCrvStruct *Crv1a, *Crv1b, *Crv2a, *Crv2b, *CmpsCrvA, *CmpsCrvB;
	CagdPtStruct *Pts;

	t = Crv1 -> KnotVector[(Crv1 -> Order + Crv1 -> Length) / 2];
	Crv1a = CagdCrvSubdivAtParam(Crv1, t);
	Crv1b = Crv1a -> Pnext;
	Crv1a -> Pnext = NULL;

	Pts = SymbCrvConstSet(Crv2, 1, ZERO_SET_EPS, t);
	if (Pts) {
	    if (Pts -> Pnext != NULL)
		SYMB_FATAL_ERROR(SYMB_ERR_REPARAM_NOT_MONOTONE);

	    t = Pts -> Pt[0];
	    CagdPtFreeList(Pts);
	    Crv2a = CagdCrvSubdivAtParam(Crv2, t);
	    Crv2b = Crv2a -> Pnext;
	    Crv2a -> Pnext = NULL;
	}
	else
	    Crv2a = Crv2b = NULL;

	CagdCrvDomain(Crv1, &TMin, &TMax);
	if (APX_EQ(TMin, t))
	    CmpsCrvA = NULL;
	else
	    CmpsCrvA = SymbComposeCrvCrvAux2(Crv1a, Crv2a ? Crv2a : Crv2);
	if (APX_EQ(TMax, t))
	    CmpsCrvB = NULL;
	else
	    CmpsCrvB = SymbComposeCrvCrvAux2(Crv1b, Crv2b ? Crv2b : Crv2);
	CagdCrvFree(Crv1a);
	CagdCrvFree(Crv1b);
	if (Crv2a)
	    CagdCrvFree(Crv2a);
	if (Crv2b)
	    CagdCrvFree(Crv2b);

	if (CmpsCrvA == NULL)
	    CmpsCrv = CmpsCrvB;
	else if (CmpsCrvB == NULL)
	    CmpsCrv = CmpsCrvA;
	else {
	    CmpsCrv = CagdMergeCrvCrv(CmpsCrvA, CmpsCrvB, FALSE);
	    CagdCrvFree(CmpsCrvA);
	    CagdCrvFree(CmpsCrvB);
	}
    }
    else {
	/* We can compose the curves, but first make sure the domain of */
	/* Crv1 is [0, 1] which is also the range of Crv2.		*/
	CagdCrvDomain(Crv1, &TMin, &TMax);

	if (!APX_EQ(TMin, 0.0) || !APX_EQ(TMax, 1.0)) {
	    CagdRType
		Translate = -TMin;
	    CagdCrvStruct
	        *Crv2Tmp = CagdCrvCopy(Crv2);

	    CagdCrvTransform(Crv2Tmp, &Translate, 1.0 / (TMax - TMin));
	    CmpsCrv = BzrComposeCrvCrv(Crv1, Crv2Tmp);
	    CagdCrvFree(Crv2Tmp);
	}
	else
	    CmpsCrv = BzrComposeCrvCrv(Crv1, Crv2);

	/* Now map the domain of the curve to be domain of Crv2. */
	CagdCrvDomain(Crv2, &TMin, &TMax);

	CagdCrvDomain(CmpsCrv, &CTMin, &CTMax);
	if (CmpsCrv -> GType == CAGD_CBEZIER_TYPE) {
	    CagdCrvStruct
		*CTmp = CnvrtBezier2BsplineCrv(CmpsCrv);

	    CagdCrvFree(CmpsCrv);
	    CmpsCrv = CTmp;
	}
	BspKnotAffineTrans(CmpsCrv -> KnotVector,
			   CmpsCrv -> Length + CmpsCrv -> Order,
			   TMin - CTMin,
			   (TMax - TMin) / (CTMax - CTMin));
    }

    return CmpsCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two Bezier curves, Crv1 and Crv2, computes their composition         M
* Crv1(Crv2(t)).							     M
*   Crv2 must be a scalar curve completely contained in Crv1's parametric    M
* domain.								     M
*   See: "Freeform surfcae analysis using a hybrid of symbolic and numeric   M
* computation" by Gershon Elber, PhD thesis, University of Utah, 1992.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:   The two curve to compose together.                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   The composed curve.                                   M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrComposeCrvCrv, composition                                            M
*****************************************************************************/
CagdCrvStruct *BzrComposeCrvCrv(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdBType
	IsRational = CAGD_IS_RATIONAL_CRV(Crv1);
    int i, j, k, CmpsOrder,
	Order = Crv1 -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv1 -> PType);
    CagdRType
	Translate = 0.0;
    CagdCrvStruct **Crv2Factors,
    	*CmpsCrv = NULL;

    if (CAGD_NUM_OF_PT_COORD(Crv2 -> PType) != 1)
	SYMB_FATAL_ERROR(SYMB_ERR_WRONG_PT_TYPE);

    Crv2Factors = SymbComputeCurvePowers(Crv2, Order);
    CmpsCrv = BzrCrvNew(Crv2Factors[0] -> Length, Crv1 -> PType);
    CmpsOrder = CmpsCrv -> Order;

    for (j = !IsRational; j <= MaxCoord; j++) {
    	CagdRType
	    *CmpsPoints = CmpsCrv -> Points[j],
	    *Crv1Points = Crv1 -> Points[j];

	for (i = 0; i < Order; i++) {
	    CagdCrvStruct
		*CTmp = CagdCrvCopy(Crv2Factors[i]);
	    CagdRType
	        *CTmpPoints = CTmp -> Points[1];

	    CagdCrvTransform(CTmp, &Translate, *Crv1Points++);
	    if (i == 0) {
       		SYMB_GEN_COPY(CmpsPoints, CTmpPoints,
			      CmpsOrder * sizeof(CagdRType));
	    }
	    else {
	    	for (k = 0; k < CmpsOrder; k++)
	    	    CmpsPoints[k] += CTmpPoints[k];
	    }
	}
    }

    for (i = 0; i < Order; i++)
	CagdCrvFree(Crv2Factors[i]);

    if (CAGD_IS_RATIONAL_CRV(Crv2)) {
	CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *CTmp;

	SymbCrvSplitScalar(CmpsCrv, &CrvW, &CrvX, &CrvY, &CrvZ);
	CTmp = SymbCrvMergeScalar(Crv2Factors[Order], CrvX, CrvY, CrvZ);
	CagdCrvFree(CmpsCrv);
	CmpsCrv = CTmp;

	if (CrvX)
	    CagdCrvFree(CrvX);
	if (CrvY)
	    CagdCrvFree(CrvY);
	if (CrvZ)
	    CagdCrvFree(CrvZ);
	
	CagdCrvFree(Crv2Factors[Order]);
    }

    IritFree((VoidPtr) Crv2Factors);

    return CmpsCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a curve Crv and surface Srf, computes the composition Srf(Crv(t)).   M
*   Crv must be a two dimensional curve completely contained in the          M
* parametric domain of Srf.                                                  M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf, Crv:   The curve and surface to compose.                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:    The resulting composition.                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbComposeSrfCrv, composition                                           M
*****************************************************************************/
CagdCrvStruct *SymbComposeSrfCrv(CagdSrfStruct *Srf, CagdCrvStruct *Crv)
{
    CagdCrvStruct *CmpsCrv;

    switch (Srf -> GType) {
	case CAGD_SBEZIER_TYPE:
	    break;
	case CAGD_SBSPLINE_TYPE:
	    SYMB_FATAL_ERROR(SYMB_ERR_BSPLINE_NO_SUPPORT);
	    break;
	case CAGD_SPOWER_TYPE:
	    SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
	    break;
    }

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

    CmpsCrv = SymbComposeSrfCrvAux(Srf, Crv);

    return CmpsCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Aux. function. Subdivides Crv until it is a Bezier curve, invokes the      *
* Bezier composition code on each, and merges them back to a complete curve. *
*   At this point, the curve can be either Bezier or Bspline only.           *
*   Curve is assumed to have open end condition.			     *
*                                                                            *
* PARAMETERS:                                                                *
*   Srf, Crv:   The curve and surface to compose together.                   *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:  The composed curve.                                    *
*****************************************************************************/
static CagdCrvStruct *SymbComposeSrfCrvAux(CagdSrfStruct *Srf,
					   CagdCrvStruct *Crv)
{
    CagdRType t, TMin, TMax;
    CagdCrvStruct *CmpsCrv;

    if (Crv -> Length > Crv -> Order) {
	/* Crv is not a Bezier segment. Subdivide, compute for each segment  */
	/* and merge back.						     */
	CagdCrvStruct *CrvA, *CrvB, *CmpsCrvA, *CmpsCrvB;

	t = Crv -> KnotVector[(Crv -> Order + Crv -> Length) / 2];

	CrvA = CagdCrvSubdivAtParam(Crv, t);
	CrvB = CrvA -> Pnext;

	CmpsCrvA = SymbComposeSrfCrvAux(Srf, CrvA);
	CmpsCrvB = SymbComposeSrfCrvAux(Srf, CrvB);
	CagdCrvFree(CrvA);
	CagdCrvFree(CrvB);

	CmpsCrv = CagdMergeCrvCrv(CmpsCrvA, CmpsCrvB, FALSE);
	CagdCrvFree(CmpsCrvA);
	CagdCrvFree(CmpsCrvB);
    }
    else {
    	/* Crv is a Bezier segment - compute its composition. */
    	CmpsCrv = BzrComposeSrfCrv(Srf, Crv);

	/* Now map the domain of the composed curve to be crv's domain. */
	CagdCrvDomain(Crv, &TMin, &TMax);
	if (CmpsCrv -> GType == CAGD_CBEZIER_TYPE) {
	    CagdCrvStruct
		*CTmp = CnvrtBezier2BsplineCrv(CmpsCrv);

	    CagdCrvFree(CmpsCrv);
	    CmpsCrv = CTmp;
	}
	BspKnotAffineTrans(CmpsCrv -> KnotVector,
			   CmpsCrv -> Length + CmpsCrv -> Order,
			   TMin - CmpsCrv -> KnotVector[0],
			   TMax - TMin);
    }

    return CmpsCrv;    
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a Bezier curve Crv and a Bezier surface Srf, computes their          M
* coomposition Srf(Crv(t)).						     M
*   Crv must be a two dimensional curve completely contained in the          M
* parametric domain of Srf.                                                  M
*   See: "Freeform surfcae analysis using a hybrid of symbolic and numeric   M
* computation" by Gershon Elber, PhD thesis, University of Utah, 1992.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf, Crv:     The curve and surface to compose.                          M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:    The resulting composition.                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrComposeSrfCrv, composition                                            M
*****************************************************************************/
CagdCrvStruct *BzrComposeSrfCrv(CagdSrfStruct *Srf, CagdCrvStruct *Crv)
{
    CagdBType
	IsRational = CAGD_IS_RATIONAL_SRF(Srf);
    int i, j, k, l, CmpsOrder,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdRType
	Translate = 0.0;
    CagdCrvStruct **CrvUFactors, **CrvVFactors, *CrvW, *CrvX, *CrvY, *CrvZ,
	*CrvU, *CrvV,
    	*CmpsCrv = NULL;

    if (CAGD_NUM_OF_PT_COORD(Crv -> PType) != 2)
	SYMB_FATAL_ERROR(SYMB_ERR_WRONG_PT_TYPE);
    SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
    CrvU = CrvW == NULL ? CagdCrvCopy(CrvX)
		        : SymbCrvMergeScalar(CrvW, CrvX, NULL, NULL);
    CrvV = CrvW == NULL ? CagdCrvCopy(CrvY)
			: SymbCrvMergeScalar(CrvW, CrvY, NULL, NULL);
    if (CrvW)
	CagdCrvFree(CrvW);
    CagdCrvFree(CrvX);
    CagdCrvFree(CrvY);

    CrvUFactors = SymbComputeCurvePowers(CrvU, UOrder);
    CrvVFactors = SymbComputeCurvePowers(CrvV, VOrder);

    CmpsCrv = BzrCrvNew(CrvUFactors[0] -> Length + 
			CrvVFactors[0] -> Length - 1, Srf -> PType);
    CmpsOrder = CmpsCrv -> Order;

    for (k = !IsRational; k <= MaxCoord; k++) {
    	CagdRType
	    *CmpsPoints = CmpsCrv -> Points[k],
	    *SPoints = Srf -> Points[k];

	for (j = 0; j < VOrder; j++) {
	    for (i = 0; i < UOrder; i++) {
		CagdCrvStruct
		    *CTmp = SymbCrvMult(CrvUFactors[i], CrvVFactors[j]);
		CagdRType
		    *CTmpPoints = CTmp -> Points[1];

		CagdCrvTransform(CTmp, &Translate, *SPoints++);

		if (i == 0 && j == 0) {
		    SYMB_GEN_COPY(CmpsPoints, CTmpPoints,
				  CmpsOrder * sizeof(CagdRType));
		}
		else {
		    for (l = 0; l < CmpsOrder; l++)
			CmpsPoints[l] += CTmpPoints[l];
		}
	    }
	}
    }

    for (i = 0; i < UOrder; i++)
    	CagdCrvFree(CrvUFactors[i]);
    for (i = 0; i < VOrder; i++)
    	CagdCrvFree(CrvVFactors[i]);

    if (CAGD_IS_RATIONAL_CRV(Crv)) {
	CagdCrvStruct *CTmp,
	    *NewCrvW = SymbCrvMult(CrvUFactors[UOrder], CrvVFactors[VOrder]);

	SymbCrvSplitScalar(CmpsCrv, &CrvW, &CrvX, &CrvY, &CrvZ);
	CTmp = SymbCrvMergeScalar(NewCrvW, CrvX, CrvY, CrvZ);
	CagdCrvFree(NewCrvW);
	CagdCrvFree(CmpsCrv);
	CmpsCrv = CTmp;

	if (CrvX)
	    CagdCrvFree(CrvX);
	if (CrvY)
	    CagdCrvFree(CrvY);
	if (CrvZ)
	    CagdCrvFree(CrvZ);

	CagdCrvFree(CrvUFactors[UOrder]);
	CagdCrvFree(CrvVFactors[VOrder]);
    }

    IritFree((VoidPtr) CrvUFactors);
    IritFree((VoidPtr) CrvVFactors);

    CagdCrvFree(CrvU);
    CagdCrvFree(CrvV);

    return CmpsCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Computes the factors of the Bernstein polynomials where crv is a scalar    *
* curve.								     *
*									     *
*   n            n              n-i         i				     *
*  B (crv(t)) = ( ) (1 - crv(t))    (crv(t))				     *
*   i            i							     *
*									     *
*   The curve crv(t) is a scalar, possibly rational curve.		     *
*   The constant 1 is equal to wcrv(t) if crv(t) is rational.		     *
*   If rational, the returned vector, index Order will contain               *
* wcrv(t)^(Order-1). 							     *
* See: "Freeform surface analysis using a hybrid of symbolic and numeric     *
* computation" by Gershon Elber, PhD thesis, University of Utah, 1992.	     *
*                                                                            *
* PARAMETERS:                                                                *
*   Crv:      Scalar curve to compute factor.                                *
*   Order:    Order is n + 1.                                                *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct **:   A vector of all possibly factors for i equal to      *
*                       zero to i equal to n.                                *
*****************************************************************************/
static CagdCrvStruct **SymbComputeCurvePowers(CagdCrvStruct *Crv, int Order)
{
    CagdBType
	IsRational = CAGD_IS_RATIONAL_CRV(Crv);
    int i;
    CagdCrvStruct *Crv_1,
	**CrvFactors1_Crv = (CagdCrvStruct **)
			    IritMalloc((Order + 1) * sizeof(CagdCrvStruct *)),
	**CrvFactorsCrv = (CagdCrvStruct **)
			    IritMalloc((Order + 1) * sizeof(CagdCrvStruct *)),
	**CrvFactors = (CagdCrvStruct **)
			    IritMalloc((Order + 1) * sizeof(CagdCrvStruct *));
    CagdRType *Points,
	Translate = 0.0;

    if (IsRational) {
	CagdCrvStruct *CrvX, *CrvY, *CrvZ, *CrvW, *CTmp;

	SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
	Crv_1 = SymbCrvSub(CrvW, CrvX);

	/* Set CrvFactors[Order] to CrvW(t)^(Order-1) if curve is rational. */
	CrvFactors[Order] = CagdCrvCopy(CrvW);
	for (i = 1; i < Order - 1; i++) {
	    CTmp = SymbCrvMult(CrvFactors[Order], CrvW);
	    CagdCrvFree(CrvFactors[Order]);
	    CrvFactors[Order] = CTmp;
	}

	CagdCrvFree(CrvW);
	CagdCrvFree(CrvX);
    }
    else {
    	Crv_1 = CagdCrvCopy(Crv);
	Points = Crv_1 -> Points[1];
    	for (i = 0; i < Crv -> Order; i++, Points++)
    	    *Points = 1.0 - *Points;

	CrvFactors[Order] = NULL;
    }

    for (i = 0; i < Order; i++) {
    	if (i == 0) {
    	    CrvFactors1_Crv[0] = NULL;
    	    CrvFactorsCrv[0] = NULL;
	}
    	else if (i == 1) {
    	    CrvFactors1_Crv[1] = Crv_1;
    	    CrvFactorsCrv[1] = CagdCrvCopy(Crv);
    	}
    	else {
    	    CrvFactors1_Crv[i] = SymbCrvMult(CrvFactors1_Crv[i - 1], Crv_1);
    	    CrvFactorsCrv[i] = SymbCrvMult(CrvFactorsCrv[i - 1], Crv);
    	}
    }

    for (i = 0; i < Order; i++) {
	if (i == 0) {
	    CrvFactors[i]= CagdCrvCopy(CrvFactors1_Crv[Order - 1]);
	}
	else if (i == Order - 1) {
	    CrvFactors[i] = CagdCrvCopy(CrvFactorsCrv[Order - 1]);
	}
	else {
	    CrvFactors[i] = SymbCrvMult(CrvFactors1_Crv[Order - 1 - i],
					CrvFactorsCrv[i]);
	}
	CagdCrvTransform(CrvFactors[i],
			 &Translate, CagdIChooseK(i, Order - 1));
    }

    for (i = 1; i < Order; i++) {
    	CagdCrvFree(CrvFactorsCrv[i]);
    	CagdCrvFree(CrvFactors1_Crv[i]);
    }
    IritFree((VoidPtr) CrvFactorsCrv);
    IritFree((VoidPtr) CrvFactors1_Crv);

    return CrvFactors;
}