/******************************************************************************
* CBzr_Pwr.c - Bezier to pawer basis conversion.			      *
*******************************************************************************
* Written by Gershon Elber, Jun. 90.					      *
******************************************************************************/

#ifdef __MSDOS__
#include <stdlib.h>
#endif /* __MSDOS__ */

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

static CagdRType BinomCoef(int n, int i);

/******************************************************************************
* Convert the given curve from Bezier basis function to Power basis.	      *
* Using:								      *
*									      *
*	  n								      *
*	  __								      *
*  n	  \	j-i n	j   j						      *
* B (t) = /  (-1)  ( ) ( ) t						      *
*  i	  --	    j   i						      *
*	 j=i								      *
*					     n-i			      *
* Which can be derived by expanding the (1-t)    term in bezier basis	      *
* function definition as:						      *
*	  								      *
* 	   n-i  							      *
*           __  							      *
*      n-i  \  n-i      j						      *
* (1-t)   = / (   ) (-t)	using binomial expansion.		      *
*	    --  j							      *
*	   j=0								      *
*                                               			      *
* This routine simply take the weight of each Bezier basis function B(t) and  *
* spread it into the different power basis t^j function scaled by:	      *
*                                               			      *
*	j-i n	j							      *
*   (-1)   ( ) ( )                                             		      *
*           j   i                                 			      *
******************************************************************************/
CagdCrvStruct *CnvrtBezier2PowerCrv(CagdCrvStruct *Crv)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j, l,
	n = Crv -> Length,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType *PwrP, *BzrP;
    CagdCrvStruct
	*NewCrv = CagdCrvNew(CAGD_CPOWER_TYPE, Crv -> PType, n);

    NewCrv -> Order = n;

    for (l = IsNotRational; l <= MaxCoord; l++) {
	PwrP = NewCrv -> Points[l];
	BzrP = Crv -> Points[l];
	ZAP_MEM(PwrP, sizeof(CagdRType) * n);

	for (i = 0; i < n; i++) {
	    for (j = i; j < n; j++) {
		PwrP[j] += BzrP[i] * BinomCoef(n, j) * BinomCoef(j, i) *
						(((j - i) & 0x01) ? -1 : 1);
	    }
	}
    }

    return NewCrv;
}

/******************************************************************************
* Convert the given curve from Power basis function to Bezier basis.	      *
* Using:								      *
*									      *
*      n    j								      *
*      __  ( )								      *
*   i  \    i	 n							      *
*  t = /  ----- B (t)							      *
*      --   n	 j							      *
*      j=i ( )								      *
*	    i								      *
*                                               			      *
* This routine simply take the weight of each Power basis function t^i and    *
* spread it into the different basis basis function B(t) scaled by:	      *
*                                               			      *
*     j   / n								      *
*    ( ) / ( )                                             		      *
*     i /   i	                                 			      *
******************************************************************************/
CagdCrvStruct *CnvrtPower2BezierCrv(CagdCrvStruct *Crv)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j, l,
	n = Crv -> Length,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType *PwrP, *BzrP;
    CagdCrvStruct
	*NewCrv = BzrCrvNew(n, Crv -> PType);

    for (l = IsNotRational; l <= MaxCoord; l++) {
	PwrP = Crv -> Points[l];
	BzrP = NewCrv -> Points[l];
	ZAP_MEM(BzrP, sizeof(CagdRType) * n);

	for (i = 0; i < n; i++) {
	    for (j = i; j < n; j++) {
		BzrP[j] += PwrP[i] * BinomCoef(j, i) / BinomCoef(n, i);
	    }
	}
    }

    return NewCrv;
}

/******************************************************************************
* Evaluate the following:						      *
*			 n         n!					      *
*			( ) = -------------				      *
*			 i    i! * (n - i)!				      *
******************************************************************************/
static CagdRType BinomCoef(int n, int i)
{
    int j;
    CagdRType c = 1.0;

    if ((n >> 1) > i) {				/* i is less than half of n: */
	for (j = n - i + 1; j <= n; j++) c *= j;
	for (j = 2; j <= i; j++) c /= j;
    }
    else {
	for (j = i + 1; j <= n; j++) c *= j;
	for (j = 2; j <= n - i; j++) c /= j;
    }

    return c;
}