src/mesa/math/m_eval.c - platform/external/mesa3d - Gitiles

 /* $Id: m_eval.c,v 1.2 2001/03/07 05:06:12 brianp Exp $ */

 /*
  * Mesa 3-D graphics library
  * Version:  3.5
  *
  * Copyright (C) 1999-2000  Brian Paul   All Rights Reserved.
  *
  * Permission is hereby granted, free of charge, to any person obtaining a
  * copy of this software and associated documentation files (the "Software"),
  * to deal in the Software without restriction, including without limitation
  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  * and/or sell copies of the Software, and to permit persons to whom the
  * Software is furnished to do so, subject to the following conditions:
  *
  * The above copyright notice and this permission notice shall be included
  * in all copies or substantial portions of the Software.
  *
  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  * BRIAN PAUL BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN
  * AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  */


 /*
  * eval.c was written by
  * Bernd Barsuhn (bdbarsuh@cip.informatik.uni-erlangen.de) and
  * Volker Weiss (vrweiss@cip.informatik.uni-erlangen.de).
  *
  * My original implementation of evaluators was simplistic and didn't
  * compute surface normal vectors properly.  Bernd and Volker applied
  * used more sophisticated methods to get better results.
  *
  * Thanks guys!
  */


 #include "glheader.h"
 #include "config.h"
 #include "m_eval.h"

 static GLfloat inv_tab[MAX_EVAL_ORDER];


 /*
  * Horner scheme for Bezier curves
  *
  * Bezier curves can be computed via a Horner scheme.
  * Horner is numerically less stable than the de Casteljau
  * algorithm, but it is faster. For curves of degree n
  * the complexity of Horner is O(n) and de Casteljau is O(n^2).
  * Since stability is not important for displaying curve
  * points I decided to use the Horner scheme.
  *
  * A cubic Bezier curve with control points b0, b1, b2, b3 can be
  * written as
  *
  *        (([3]        [3]     )     [3]       )     [3]
  * c(t) = (([0]*s*b0 + [1]*t*b1)*s + [2]*t^2*b2)*s + [3]*t^2*b3
  *
  *                                           [n]
  * where s=1-t and the binomial coefficients [i]. These can
  * be computed iteratively using the identity:
  *
  * [n]               [n  ]             [n]
  * [i] = (n-i+1)/i * [i-1]     and     [0] = 1
  */


 void
 _math_horner_bezier_curve(const GLfloat *cp, GLfloat *out, GLfloat t,
 			  GLuint dim, GLuint order)
 {
    GLfloat s, powert;
    GLuint i, k, bincoeff;

    if(order >= 2)
    {
       bincoeff = order-1;
       s = 1.0-t;

       for(k=0; k<dim; k++)
 	 out[k] = s*cp[k] + bincoeff*t*cp[dim+k];

       for(i=2, cp+=2*dim, powert=t*t; i<order; i++, powert*=t, cp +=dim)
       {
 	 bincoeff *= order-i;
 	 bincoeff *= (GLuint) inv_tab[i];

 	 for(k=0; k<dim; k++)
 	    out[k] = s*out[k] + bincoeff*powert*cp[k];
       }
    }
    else /* order=1 -> constant curve */
    {
       for(k=0; k<dim; k++)
 	 out[k] = cp[k];
    }
 }

 /*
  * Tensor product Bezier surfaces
  *
  * Again the Horner scheme is used to compute a point on a
  * TP Bezier surface. First a control polygon for a curve
  * on the surface in one parameter direction is computed,
  * then the point on the curve for the other parameter
  * direction is evaluated.
  *
  * To store the curve control polygon additional storage
  * for max(uorder,vorder) points is needed in the
  * control net cn.
  */

 void
 _math_horner_bezier_surf(GLfloat *cn, GLfloat *out, GLfloat u, GLfloat v,
 			 GLuint dim, GLuint uorder, GLuint vorder)
 {
    GLfloat *cp = cn + uorder*vorder*dim;
    GLuint i, uinc = vorder*dim;

    if(vorder > uorder)
    {
       if(uorder >= 2)
       {
 	 GLfloat s, poweru;
 	 GLuint j, k, bincoeff;

 	 /* Compute the control polygon for the surface-curve in u-direction */
 	 for(j=0; j<vorder; j++)
 	 {
 	    GLfloat *ucp = &cn[j*dim];

 	    /* Each control point is the point for parameter u on a */
 	    /* curve defined by the control polygons in u-direction */
 	    bincoeff = uorder-1;
 	    s = 1.0-u;

 	    for(k=0; k<dim; k++)
 	       cp[j*dim+k] = s*ucp[k] + bincoeff*u*ucp[uinc+k];

 	    for(i=2, ucp+=2*uinc, poweru=u*u; i<uorder;
 		i++, poweru*=u, ucp +=uinc)
 	    {
 	       bincoeff *= uorder-i;
 	       bincoeff *= (GLuint) inv_tab[i];

 	       for(k=0; k<dim; k++)
 		  cp[j*dim+k] = s*cp[j*dim+k] + bincoeff*poweru*ucp[k];
 	    }
 	 }

 	 /* Evaluate curve point in v */
 	 _math_horner_bezier_curve(cp, out, v, dim, vorder);
       }
       else /* uorder=1 -> cn defines a curve in v */
 	 _math_horner_bezier_curve(cn, out, v, dim, vorder);
    }
    else /* vorder <= uorder */
    {
       if(vorder > 1)
       {
 	 GLuint i;

 	 /* Compute the control polygon for the surface-curve in u-direction */
 	 for(i=0; i<uorder; i++, cn += uinc)
 	 {
 	    /* For constant i all cn[i][j] (j=0..vorder) are located */
 	    /* on consecutive memory locations, so we can use        */
 	    /* horner_bezier_curve to compute the control points     */

 	    _math_horner_bezier_curve(cn, &cp[i*dim], v, dim, vorder);
 	 }

 	 /* Evaluate curve point in u */
 	 _math_horner_bezier_curve(cp, out, u, dim, uorder);
       }
       else  /* vorder=1 -> cn defines a curve in u */
 	 _math_horner_bezier_curve(cn, out, u, dim, uorder);
    }
 }

 /*
  * The direct de Casteljau algorithm is used when a point on the
  * surface and the tangent directions spanning the tangent plane
  * should be computed (this is needed to compute normals to the
  * surface). In this case the de Casteljau algorithm approach is
  * nicer because a point and the partial derivatives can be computed
  * at the same time. To get the correct tangent length du and dv
  * must be multiplied with the (u2-u1)/uorder-1 and (v2-v1)/vorder-1.
  * Since only the directions are needed, this scaling step is omitted.
  *
  * De Casteljau needs additional storage for uorder*vorder
  * values in the control net cn.
  */

 void
 _math_de_casteljau_surf(GLfloat *cn, GLfloat *out, GLfloat *du, GLfloat *dv,
 			GLfloat u, GLfloat v, GLuint dim,
 			GLuint uorder, GLuint vorder)
 {
    GLfloat *dcn = cn + uorder*vorder*dim;
    GLfloat us = 1.0-u, vs = 1.0-v;
    GLuint h, i, j, k;
    GLuint minorder = uorder < vorder ? uorder : vorder;
    GLuint uinc = vorder*dim;
    GLuint dcuinc = vorder;

    /* Each component is evaluated separately to save buffer space  */
    /* This does not drasticaly decrease the performance of the     */
    /* algorithm. If additional storage for (uorder-1)*(vorder-1)   */
    /* points would be available, the components could be accessed  */
    /* in the innermost loop which could lead to less cache misses. */

 #define CN(I,J,K) cn[(I)*uinc+(J)*dim+(K)]
 #define DCN(I, J) dcn[(I)*dcuinc+(J)]
    if(minorder < 3)
    {
       if(uorder==vorder)
       {
 	 for(k=0; k<dim; k++)
 	 {
 	    /* Derivative direction in u */
 	    du[k] = vs*(CN(1,0,k) - CN(0,0,k)) +
 	       v*(CN(1,1,k) - CN(0,1,k));

 	    /* Derivative direction in v */
 	    dv[k] = us*(CN(0,1,k) - CN(0,0,k)) +
 	       u*(CN(1,1,k) - CN(1,0,k));

 	    /* bilinear de Casteljau step */
 	    out[k] =  us*(vs*CN(0,0,k) + v*CN(0,1,k)) +
 	       u*(vs*CN(1,0,k) + v*CN(1,1,k));
 	 }
       }
       else if(minorder == uorder)
       {
 	 for(k=0; k<dim; k++)
 	 {
 	    /* bilinear de Casteljau step */
 	    DCN(1,0) =    CN(1,0,k) -   CN(0,0,k);
 	    DCN(0,0) = us*CN(0,0,k) + u*CN(1,0,k);

 	    for(j=0; j<vorder-1; j++)
 	    {
 	       /* for the derivative in u */
 	       DCN(1,j+1) =    CN(1,j+1,k) -   CN(0,j+1,k);
 	       DCN(1,j)   = vs*DCN(1,j)    + v*DCN(1,j+1);

 	       /* for the `point' */
 	       DCN(0,j+1) = us*CN(0,j+1,k) + u*CN(1,j+1,k);
 	       DCN(0,j)   = vs*DCN(0,j)    + v*DCN(0,j+1);
 	    }

 	    /* remaining linear de Casteljau steps until the second last step */
 	    for(h=minorder; h<vorder-1; h++)
 	       for(j=0; j<vorder-h; j++)
 	       {
 		  /* for the derivative in u */
 		  DCN(1,j) = vs*DCN(1,j) + v*DCN(1,j+1);

 		  /* for the `point' */
 		  DCN(0,j) = vs*DCN(0,j) + v*DCN(0,j+1);
 	       }

 	    /* derivative direction in v */
 	    dv[k] = DCN(0,1) - DCN(0,0);

 	    /* derivative direction in u */
 	    du[k] =   vs*DCN(1,0) + v*DCN(1,1);

 	    /* last linear de Casteljau step */
 	    out[k] =  vs*DCN(0,0) + v*DCN(0,1);
 	 }
       }
       else /* minorder == vorder */
       {
 	 for(k=0; k<dim; k++)
 	 {
 	    /* bilinear de Casteljau step */
 	    DCN(0,1) =    CN(0,1,k) -   CN(0,0,k);
 	    DCN(0,0) = vs*CN(0,0,k) + v*CN(0,1,k);
 	    for(i=0; i<uorder-1; i++)
 	    {
 	       /* for the derivative in v */
 	       DCN(i+1,1) =    CN(i+1,1,k) -   CN(i+1,0,k);
 	       DCN(i,1)   = us*DCN(i,1)    + u*DCN(i+1,1);

 	       /* for the `point' */
 	       DCN(i+1,0) = vs*CN(i+1,0,k) + v*CN(i+1,1,k);
 	       DCN(i,0)   = us*DCN(i,0)    + u*DCN(i+1,0);
 	    }

 	    /* remaining linear de Casteljau steps until the second last step */
 	    for(h=minorder; h<uorder-1; h++)
 	       for(i=0; i<uorder-h; i++)
 	       {
 		  /* for the derivative in v */
 		  DCN(i,1) = us*DCN(i,1) + u*DCN(i+1,1);

 		  /* for the `point' */
 		  DCN(i,0) = us*DCN(i,0) + u*DCN(i+1,0);
 	       }

 	    /* derivative direction in u */
 	    du[k] = DCN(1,0) - DCN(0,0);

 	    /* derivative direction in v */
 	    dv[k] =   us*DCN(0,1) + u*DCN(1,1);

 	    /* last linear de Casteljau step */
 	    out[k] =  us*DCN(0,0) + u*DCN(1,0);
 	 }
       }
    }
    else if(uorder == vorder)
    {
       for(k=0; k<dim; k++)
       {
 	 /* first bilinear de Casteljau step */
 	 for(i=0; i<uorder-1; i++)
 	 {
 	    DCN(i,0) = us*CN(i,0,k) + u*CN(i+1,0,k);
 	    for(j=0; j<vorder-1; j++)
 	    {
 	       DCN(i,j+1) = us*CN(i,j+1,k) + u*CN(i+1,j+1,k);
 	       DCN(i,j)   = vs*DCN(i,j)    + v*DCN(i,j+1);
 	    }
 	 }

 	 /* remaining bilinear de Casteljau steps until the second last step */
 	 for(h=2; h<minorder-1; h++)
 	    for(i=0; i<uorder-h; i++)
 	    {
 	       DCN(i,0) = us*DCN(i,0) + u*DCN(i+1,0);
 	       for(j=0; j<vorder-h; j++)
 	       {
 		  DCN(i,j+1) = us*DCN(i,j+1) + u*DCN(i+1,j+1);
 		  DCN(i,j)   = vs*DCN(i,j)   + v*DCN(i,j+1);
 	       }
 	    }

 	 /* derivative direction in u */
 	 du[k] = vs*(DCN(1,0) - DCN(0,0)) +
 	    v*(DCN(1,1) - DCN(0,1));

 	 /* derivative direction in v */
 	 dv[k] = us*(DCN(0,1) - DCN(0,0)) +
 	    u*(DCN(1,1) - DCN(1,0));

 	 /* last bilinear de Casteljau step */
 	 out[k] =  us*(vs*DCN(0,0) + v*DCN(0,1)) +
 	    u*(vs*DCN(1,0) + v*DCN(1,1));
       }
    }
    else if(minorder == uorder)
    {
       for(k=0; k<dim; k++)
       {
 	 /* first bilinear de Casteljau step */
 	 for(i=0; i<uorder-1; i++)
 	 {
 	    DCN(i,0) = us*CN(i,0,k) + u*CN(i+1,0,k);
 	    for(j=0; j<vorder-1; j++)
 	    {
 	       DCN(i,j+1) = us*CN(i,j+1,k) + u*CN(i+1,j+1,k);
 	       DCN(i,j)   = vs*DCN(i,j)    + v*DCN(i,j+1);
 	    }
 	 }

 	 /* remaining bilinear de Casteljau steps until the second last step */
 	 for(h=2; h<minorder-1; h++)
 	    for(i=0; i<uorder-h; i++)
 	    {
 	       DCN(i,0) = us*DCN(i,0) + u*DCN(i+1,0);
 	       for(j=0; j<vorder-h; j++)
 	       {
 		  DCN(i,j+1) = us*DCN(i,j+1) + u*DCN(i+1,j+1);
 		  DCN(i,j)   = vs*DCN(i,j)   + v*DCN(i,j+1);
 	       }
 	    }

 	 /* last bilinear de Casteljau step */
 	 DCN(2,0) =    DCN(1,0) -   DCN(0,0);
 	 DCN(0,0) = us*DCN(0,0) + u*DCN(1,0);
 	 for(j=0; j<vorder-1; j++)
 	 {
 	    /* for the derivative in u */
 	    DCN(2,j+1) =    DCN(1,j+1) -    DCN(0,j+1);
 	    DCN(2,j)   = vs*DCN(2,j)    + v*DCN(2,j+1);

 	    /* for the `point' */
 	    DCN(0,j+1) = us*DCN(0,j+1 ) + u*DCN(1,j+1);
 	    DCN(0,j)   = vs*DCN(0,j)    + v*DCN(0,j+1);
 	 }

 	 /* remaining linear de Casteljau steps until the second last step */
 	 for(h=minorder; h<vorder-1; h++)
 	    for(j=0; j<vorder-h; j++)
 	    {
 	       /* for the derivative in u */
 	       DCN(2,j) = vs*DCN(2,j) + v*DCN(2,j+1);

 	       /* for the `point' */
 	       DCN(0,j) = vs*DCN(0,j) + v*DCN(0,j+1);
 	    }

 	 /* derivative direction in v */
 	 dv[k] = DCN(0,1) - DCN(0,0);

 	 /* derivative direction in u */
 	 du[k] =   vs*DCN(2,0) + v*DCN(2,1);

 	 /* last linear de Casteljau step */
 	 out[k] =  vs*DCN(0,0) + v*DCN(0,1);
       }
    }
    else /* minorder == vorder */
    {
       for(k=0; k<dim; k++)
       {
 	 /* first bilinear de Casteljau step */
 	 for(i=0; i<uorder-1; i++)
 	 {
 	    DCN(i,0) = us*CN(i,0,k) + u*CN(i+1,0,k);
 	    for(j=0; j<vorder-1; j++)
 	    {
 	       DCN(i,j+1) = us*CN(i,j+1,k) + u*CN(i+1,j+1,k);
 	       DCN(i,j)   = vs*DCN(i,j)    + v*DCN(i,j+1);
 	    }
 	 }

 	 /* remaining bilinear de Casteljau steps until the second last step */
 	 for(h=2; h<minorder-1; h++)
 	    for(i=0; i<uorder-h; i++)
 	    {
 	       DCN(i,0) = us*DCN(i,0) + u*DCN(i+1,0);
 	       for(j=0; j<vorder-h; j++)
 	       {
 		  DCN(i,j+1) = us*DCN(i,j+1) + u*DCN(i+1,j+1);
 		  DCN(i,j)   = vs*DCN(i,j)   + v*DCN(i,j+1);
 	       }
 	    }

 	 /* last bilinear de Casteljau step */
 	 DCN(0,2) =    DCN(0,1) -   DCN(0,0);
 	 DCN(0,0) = vs*DCN(0,0) + v*DCN(0,1);
 	 for(i=0; i<uorder-1; i++)
 	 {
 	    /* for the derivative in v */
 	    DCN(i+1,2) =    DCN(i+1,1)  -   DCN(i+1,0);
 	    DCN(i,2)   = us*DCN(i,2)    + u*DCN(i+1,2);

 	    /* for the `point' */
 	    DCN(i+1,0) = vs*DCN(i+1,0)  + v*DCN(i+1,1);
 	    DCN(i,0)   = us*DCN(i,0)    + u*DCN(i+1,0);
 	 }

 	 /* remaining linear de Casteljau steps until the second last step */
 	 for(h=minorder; h<uorder-1; h++)
 	    for(i=0; i<uorder-h; i++)
 	    {
 	       /* for the derivative in v */
 	       DCN(i,2) = us*DCN(i,2) + u*DCN(i+1,2);

 	       /* for the `point' */
 	       DCN(i,0) = us*DCN(i,0) + u*DCN(i+1,0);
 	    }

 	 /* derivative direction in u */
 	 du[k] = DCN(1,0) - DCN(0,0);

 	 /* derivative direction in v */
 	 dv[k] =   us*DCN(0,2) + u*DCN(1,2);

 	 /* last linear de Casteljau step */
 	 out[k] =  us*DCN(0,0) + u*DCN(1,0);
       }
    }
 #undef DCN
 #undef CN
 }


 /*
  * Do one-time initialization for evaluators.
  */
 void _math_init_eval( void )
 {
    GLuint i;

    /* KW: precompute 1/x for useful x.
     */
    for (i = 1 ; i < MAX_EVAL_ORDER ; i++)
       inv_tab[i] = 1.0 / i;
 }
	/* $Id: m_eval.c,v 1.2 2001/03/07 05:06:12 brianp Exp $ */

	/*
	* Mesa 3-D graphics library
	* Version: 3.5
	*
	* Copyright (C) 1999-2000 Brian Paul All Rights Reserved.
	*
	* Permission is hereby granted, free of charge, to any person obtaining a
	* copy of this software and associated documentation files (the "Software"),
	* to deal in the Software without restriction, including without limitation
	* the rights to use, copy, modify, merge, publish, distribute, sublicense,
	* and/or sell copies of the Software, and to permit persons to whom the
	* Software is furnished to do so, subject to the following conditions:
	*
	* The above copyright notice and this permission notice shall be included
	* in all copies or substantial portions of the Software.
	*
	* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
	* BRIAN PAUL BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN
	* AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
	* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
	*/


	/*
	* eval.c was written by
	* Bernd Barsuhn (bdbarsuh@cip.informatik.uni-erlangen.de) and
	* Volker Weiss (vrweiss@cip.informatik.uni-erlangen.de).
	*
	* My original implementation of evaluators was simplistic and didn't
	* compute surface normal vectors properly. Bernd and Volker applied
	* used more sophisticated methods to get better results.
	*
	* Thanks guys!
	*/


	#include "glheader.h"
	#include "config.h"
	#include "m_eval.h"

	static GLfloat inv_tab[MAX_EVAL_ORDER];



	/*
	* Horner scheme for Bezier curves
	*
	* Bezier curves can be computed via a Horner scheme.
	* Horner is numerically less stable than the de Casteljau
	* algorithm, but it is faster. For curves of degree n
	* the complexity of Horner is O(n) and de Casteljau is O(n^2).
	* Since stability is not important for displaying curve
	* points I decided to use the Horner scheme.
	*
	* A cubic Bezier curve with control points b0, b1, b2, b3 can be
	* written as
	*
	* (([3] [3] ) [3] ) [3]
	* c(t) = (([0]sb0 + [1]tb1)s + [2]t^2b2)s + [3]t^2b3
	*
	* [n]
	* where s=1-t and the binomial coefficients [i]. These can
	* be computed iteratively using the identity:
	*
	* [n] [n ] [n]
	* [i] = (n-i+1)/i * [i-1] and [0] = 1
	*/


	void
	_math_horner_bezier_curve(const GLfloat cp, GLfloat out, GLfloat t,
	GLuint dim, GLuint order)
	{
	GLfloat s, powert;
	GLuint i, k, bincoeff;

	if(order >= 2)
	{
	bincoeff = order-1;
	s = 1.0-t;

	for(k=0; k<dim; k++)
	out[k] = scp[k] + bincoefft*cp[dim+k];

	for(i=2, cp+=2dim, powert=tt; i<order; i++, powert*=t, cp +=dim)
	{
	bincoeff *= order-i;
	bincoeff *= (GLuint) inv_tab[i];

	for(k=0; k<dim; k++)
	out[k] = sout[k] + bincoeffpowert*cp[k];
	}
	}
	else /* order=1 -> constant curve */
	{
	for(k=0; k<dim; k++)
	out[k] = cp[k];
	}
	}

	/*
	* Tensor product Bezier surfaces
	*
	* Again the Horner scheme is used to compute a point on a
	* TP Bezier surface. First a control polygon for a curve
	* on the surface in one parameter direction is computed,
	* then the point on the curve for the other parameter
	* direction is evaluated.
	*
	* To store the curve control polygon additional storage
	* for max(uorder,vorder) points is needed in the
	* control net cn.
	*/

	void
	_math_horner_bezier_surf(GLfloat cn, GLfloat out, GLfloat u, GLfloat v,
	GLuint dim, GLuint uorder, GLuint vorder)
	{
	GLfloat cp = cn + uordervorder*dim;
	GLuint i, uinc = vorder*dim;

	if(vorder > uorder)
	{
	if(uorder >= 2)
	{
	GLfloat s, poweru;
	GLuint j, k, bincoeff;

	/* Compute the control polygon for the surface-curve in u-direction */
	for(j=0; j<vorder; j++)
	{
	GLfloat ucp = &cn[jdim];

	/* Each control point is the point for parameter u on a */
	/* curve defined by the control polygons in u-direction */
	bincoeff = uorder-1;
	s = 1.0-u;

	for(k=0; k<dim; k++)
	cp[jdim+k] = sucp[k] + bincoeffuucp[uinc+k];

	for(i=2, ucp+=2uinc, poweru=uu; i<uorder;
	i++, poweru*=u, ucp +=uinc)
	{
	bincoeff *= uorder-i;
	bincoeff *= (GLuint) inv_tab[i];

	for(k=0; k<dim; k++)
	cp[jdim+k] = scp[jdim+k] + bincoeffpoweru*ucp[k];
	}
	}

	/* Evaluate curve point in v */
	_math_horner_bezier_curve(cp, out, v, dim, vorder);
	}
	else /* uorder=1 -> cn defines a curve in v */
	_math_horner_bezier_curve(cn, out, v, dim, vorder);
	}
	else /* vorder <= uorder */
	{
	if(vorder > 1)
	{
	GLuint i;

	/* Compute the control polygon for the surface-curve in u-direction */
	for(i=0; i<uorder; i++, cn += uinc)
	{
	/* For constant i all cn[i][j] (j=0..vorder) are located */
	/* on consecutive memory locations, so we can use */
	/* horner_bezier_curve to compute the control points */

	_math_horner_bezier_curve(cn, &cp[i*dim], v, dim, vorder);
	}

	/* Evaluate curve point in u */
	_math_horner_bezier_curve(cp, out, u, dim, uorder);
	}
	else /* vorder=1 -> cn defines a curve in u */
	_math_horner_bezier_curve(cn, out, u, dim, uorder);
	}
	}

	/*
	* The direct de Casteljau algorithm is used when a point on the
	* surface and the tangent directions spanning the tangent plane
	* should be computed (this is needed to compute normals to the
	* surface). In this case the de Casteljau algorithm approach is
	* nicer because a point and the partial derivatives can be computed
	* at the same time. To get the correct tangent length du and dv
	* must be multiplied with the (u2-u1)/uorder-1 and (v2-v1)/vorder-1.
	* Since only the directions are needed, this scaling step is omitted.
	*
	* De Casteljau needs additional storage for uorder*vorder
	* values in the control net cn.
	*/

	void
	_math_de_casteljau_surf(GLfloat cn, GLfloat out, GLfloat du, GLfloat dv,
	GLfloat u, GLfloat v, GLuint dim,
	GLuint uorder, GLuint vorder)
	{
	GLfloat dcn = cn + uordervorder*dim;
	GLfloat us = 1.0-u, vs = 1.0-v;
	GLuint h, i, j, k;
	GLuint minorder = uorder < vorder ? uorder : vorder;
	GLuint uinc = vorder*dim;
	GLuint dcuinc = vorder;

	/* Each component is evaluated separately to save buffer space */
	/* This does not drasticaly decrease the performance of the */
	/* algorithm. If additional storage for (uorder-1)(vorder-1) /
	/* points would be available, the components could be accessed */
	/* in the innermost loop which could lead to less cache misses. */

	#define CN(I,J,K) cn[(I)uinc+(J)dim+(K)]
	#define DCN(I, J) dcn[(I)*dcuinc+(J)]
	if(minorder < 3)
	{
	if(uorder==vorder)
	{
	for(k=0; k<dim; k++)
	{
	/* Derivative direction in u */
	du[k] = vs*(CN(1,0,k) - CN(0,0,k)) +
	v*(CN(1,1,k) - CN(0,1,k));

	/* Derivative direction in v */
	dv[k] = us*(CN(0,1,k) - CN(0,0,k)) +
	u*(CN(1,1,k) - CN(1,0,k));

	/* bilinear de Casteljau step */
	out[k] = us(vsCN(0,0,k) + v*CN(0,1,k)) +
	u(vsCN(1,0,k) + v*CN(1,1,k));
	}
	}
	else if(minorder == uorder)
	{
	for(k=0; k<dim; k++)
	{
	/* bilinear de Casteljau step */
	DCN(1,0) = CN(1,0,k) - CN(0,0,k);
	DCN(0,0) = usCN(0,0,k) + uCN(1,0,k);

	for(j=0; j<vorder-1; j++)
	{
	/* for the derivative in u */
	DCN(1,j+1) = CN(1,j+1,k) - CN(0,j+1,k);
	DCN(1,j) = vsDCN(1,j) + vDCN(1,j+1);

	/* for the `point' */
	DCN(0,j+1) = usCN(0,j+1,k) + uCN(1,j+1,k);
	DCN(0,j) = vsDCN(0,j) + vDCN(0,j+1);
	}

	/* remaining linear de Casteljau steps until the second last step */
	for(h=minorder; h<vorder-1; h++)
	for(j=0; j<vorder-h; j++)
	{
	/* for the derivative in u */
	DCN(1,j) = vsDCN(1,j) + vDCN(1,j+1);

	/* for the `point' */
	DCN(0,j) = vsDCN(0,j) + vDCN(0,j+1);
	}

	/* derivative direction in v */
	dv[k] = DCN(0,1) - DCN(0,0);

	/* derivative direction in u */
	du[k] = vsDCN(1,0) + vDCN(1,1);

	/* last linear de Casteljau step */
	out[k] = vsDCN(0,0) + vDCN(0,1);
	}
	}
	else /* minorder == vorder */
	{
	for(k=0; k<dim; k++)
	{
	/* bilinear de Casteljau step */
	DCN(0,1) = CN(0,1,k) - CN(0,0,k);
	DCN(0,0) = vsCN(0,0,k) + vCN(0,1,k);
	for(i=0; i<uorder-1; i++)
	{
	/* for the derivative in v */
	DCN(i+1,1) = CN(i+1,1,k) - CN(i+1,0,k);
	DCN(i,1) = usDCN(i,1) + uDCN(i+1,1);

	/* for the `point' */
	DCN(i+1,0) = vsCN(i+1,0,k) + vCN(i+1,1,k);
	DCN(i,0) = usDCN(i,0) + uDCN(i+1,0);
	}

	/* remaining linear de Casteljau steps until the second last step */
	for(h=minorder; h<uorder-1; h++)
	for(i=0; i<uorder-h; i++)
	{
	/* for the derivative in v */
	DCN(i,1) = usDCN(i,1) + uDCN(i+1,1);

	/* for the `point' */
	DCN(i,0) = usDCN(i,0) + uDCN(i+1,0);
	}

	/* derivative direction in u */
	du[k] = DCN(1,0) - DCN(0,0);

	/* derivative direction in v */
	dv[k] = usDCN(0,1) + uDCN(1,1);

	/* last linear de Casteljau step */
	out[k] = usDCN(0,0) + uDCN(1,0);
	}
	}
	}
	else if(uorder == vorder)
	{
	for(k=0; k<dim; k++)
	{
	/* first bilinear de Casteljau step */
	for(i=0; i<uorder-1; i++)
	{
	DCN(i,0) = usCN(i,0,k) + uCN(i+1,0,k);
	for(j=0; j<vorder-1; j++)
	{
	DCN(i,j+1) = usCN(i,j+1,k) + uCN(i+1,j+1,k);
	DCN(i,j) = vsDCN(i,j) + vDCN(i,j+1);
	}
	}

	/* remaining bilinear de Casteljau steps until the second last step */
	for(h=2; h<minorder-1; h++)
	for(i=0; i<uorder-h; i++)
	{
	DCN(i,0) = usDCN(i,0) + uDCN(i+1,0);
	for(j=0; j<vorder-h; j++)
	{
	DCN(i,j+1) = usDCN(i,j+1) + uDCN(i+1,j+1);
	DCN(i,j) = vsDCN(i,j) + vDCN(i,j+1);
	}
	}

	/* derivative direction in u */
	du[k] = vs*(DCN(1,0) - DCN(0,0)) +
	v*(DCN(1,1) - DCN(0,1));

	/* derivative direction in v */
	dv[k] = us*(DCN(0,1) - DCN(0,0)) +
	u*(DCN(1,1) - DCN(1,0));

	/* last bilinear de Casteljau step */
	out[k] = us(vsDCN(0,0) + v*DCN(0,1)) +
	u(vsDCN(1,0) + v*DCN(1,1));
	}
	}
	else if(minorder == uorder)
	{
	for(k=0; k<dim; k++)
	{
	/* first bilinear de Casteljau step */
	for(i=0; i<uorder-1; i++)
	{
	DCN(i,0) = usCN(i,0,k) + uCN(i+1,0,k);
	for(j=0; j<vorder-1; j++)
	{
	DCN(i,j+1) = usCN(i,j+1,k) + uCN(i+1,j+1,k);
	DCN(i,j) = vsDCN(i,j) + vDCN(i,j+1);
	}
	}

	/* remaining bilinear de Casteljau steps until the second last step */
	for(h=2; h<minorder-1; h++)
	for(i=0; i<uorder-h; i++)
	{
	DCN(i,0) = usDCN(i,0) + uDCN(i+1,0);
	for(j=0; j<vorder-h; j++)
	{
	DCN(i,j+1) = usDCN(i,j+1) + uDCN(i+1,j+1);
	DCN(i,j) = vsDCN(i,j) + vDCN(i,j+1);
	}
	}

	/* last bilinear de Casteljau step */
	DCN(2,0) = DCN(1,0) - DCN(0,0);
	DCN(0,0) = usDCN(0,0) + uDCN(1,0);
	for(j=0; j<vorder-1; j++)
	{
	/* for the derivative in u */
	DCN(2,j+1) = DCN(1,j+1) - DCN(0,j+1);
	DCN(2,j) = vsDCN(2,j) + vDCN(2,j+1);

	/* for the `point' */
	DCN(0,j+1) = usDCN(0,j+1 ) + uDCN(1,j+1);
	DCN(0,j) = vsDCN(0,j) + vDCN(0,j+1);
	}

	/* remaining linear de Casteljau steps until the second last step */
	for(h=minorder; h<vorder-1; h++)
	for(j=0; j<vorder-h; j++)
	{
	/* for the derivative in u */
	DCN(2,j) = vsDCN(2,j) + vDCN(2,j+1);

	/* for the `point' */
	DCN(0,j) = vsDCN(0,j) + vDCN(0,j+1);
	}

	/* derivative direction in v */
	dv[k] = DCN(0,1) - DCN(0,0);

	/* derivative direction in u */
	du[k] = vsDCN(2,0) + vDCN(2,1);

	/* last linear de Casteljau step */
	out[k] = vsDCN(0,0) + vDCN(0,1);
	}
	}
	else /* minorder == vorder */
	{
	for(k=0; k<dim; k++)
	{
	/* first bilinear de Casteljau step */
	for(i=0; i<uorder-1; i++)
	{
	DCN(i,0) = usCN(i,0,k) + uCN(i+1,0,k);
	for(j=0; j<vorder-1; j++)
	{
	DCN(i,j+1) = usCN(i,j+1,k) + uCN(i+1,j+1,k);
	DCN(i,j) = vsDCN(i,j) + vDCN(i,j+1);
	}
	}

	/* remaining bilinear de Casteljau steps until the second last step */
	for(h=2; h<minorder-1; h++)
	for(i=0; i<uorder-h; i++)
	{
	DCN(i,0) = usDCN(i,0) + uDCN(i+1,0);
	for(j=0; j<vorder-h; j++)
	{
	DCN(i,j+1) = usDCN(i,j+1) + uDCN(i+1,j+1);
	DCN(i,j) = vsDCN(i,j) + vDCN(i,j+1);
	}
	}

	/* last bilinear de Casteljau step */
	DCN(0,2) = DCN(0,1) - DCN(0,0);
	DCN(0,0) = vsDCN(0,0) + vDCN(0,1);
	for(i=0; i<uorder-1; i++)
	{
	/* for the derivative in v */
	DCN(i+1,2) = DCN(i+1,1) - DCN(i+1,0);
	DCN(i,2) = usDCN(i,2) + uDCN(i+1,2);

	/* for the `point' */
	DCN(i+1,0) = vsDCN(i+1,0) + vDCN(i+1,1);
	DCN(i,0) = usDCN(i,0) + uDCN(i+1,0);
	}

	/* remaining linear de Casteljau steps until the second last step */
	for(h=minorder; h<uorder-1; h++)
	for(i=0; i<uorder-h; i++)
	{
	/* for the derivative in v */
	DCN(i,2) = usDCN(i,2) + uDCN(i+1,2);

	/* for the `point' */
	DCN(i,0) = usDCN(i,0) + uDCN(i+1,0);
	}

	/* derivative direction in u */
	du[k] = DCN(1,0) - DCN(0,0);

	/* derivative direction in v */
	dv[k] = usDCN(0,2) + uDCN(1,2);

	/* last linear de Casteljau step */
	out[k] = usDCN(0,0) + uDCN(1,0);
	}
	}
	#undef DCN
	#undef CN
	}


	/*
	* Do one-time initialization for evaluators.
	*/
	void _math_init_eval( void )
	{
	GLuint i;

	/* KW: precompute 1/x for useful x.
	*/
	for (i = 1 ; i < MAX_EVAL_ORDER ; i++)
	inv_tab[i] = 1.0 / i;
	}