3D グラフィックスのための数学入門 157 ページをコードに．
コードは，kodatuno を参考にさせてもらった．
double	BS_basis	(const double t,c_v_double& knot,const long order,const long i)
{
	if (order <= 0)                  			{	return	0. ;	}
	if (order == 1) {	//	Ni,0(t)
		if (t == knot[knot.size()-1])
		{
			if (knot[i]<=t && t<=knot[i+1])		{	return	1. ;	}
			else                        		{	return	0. ;	}
			}
		else {
			if (knot[i]<=t && t< knot[i+1])		{	return	1. ;	}
			else                        		{	return	0. ;	}
			}
		}
	else
	{
		double	n1 = 0 ;                                	//	(*/0)  -->  n = 0
		double	n2 = 0 ;
		double	d1 = knot[i+order-1] - knot[i  ] ;      	//	T(i+k)  -T(i)
		double	d2 = knot[i+order-0] - knot[i+1] ;      	//	T(i+k+1)-T(i+1)
		if (!::V1_is_near<double<(d1,0))	{
		    	n1 = t - knot[i] ;                      	//	t - T(i)
		    	n1 = n1 * BS_basis(t,knot,order-1,i  ) ;	//	N(i,  k-1)(t)
		    	n1 = n1 / d1 ;
			}
		if (!::V1_is_near<double<(d2,0))	{
		    	n2 = knot[i+order-0] - t ;              	//	T(i+k+1)-t
		    	n2 = n2 * BS_basis(t,knot,order-1,i+1) ;	//	N(i+1,k-1)(t)
		    	n2 = n2 / d2 ;
			}
		return	(n1 + n2) ;
		}
	return	-1. ;
	}
double	BS_basis	(const double t,const long order,const long i)
{
	if (order <= 0)                  			{	return	0. ;	}
	if (order == 1) {	//	Ni,0(t)
		if (i<=t && t< i+1)           			{	return	1. ;	}
		else                            		{	return	0. ;	}
		}
	else
	{
		double	n1 = 0 ;
		double	n2 = 0 ;
		double	d1 = order - 1 ;                    		//	k
		double	d2 = order - 1 ;                    		//	k
		if (!::V1_is_near<double<(d1,0))	{
		    	n1 = t - i ;                        		//	t-i
		    	n1 = n1 * BS_basis(t,order-1,i  ) ; 		//	N(i  ,k-1)
		    	n1 = n1 / d1 ;
			}
		if (!::V1_is_near<double<(d2,0))	{
		    	n2 = (i+order-0) - t ;              		//	i+k+1 - t
		    	n2 = n2 * BS_basis(t,order-1,i+1) ; 		//	N(i+1,k-1)
		    	n2 = n2 / d2 ;
			}
		return	(n1 + n2) ;
		}
	return	-1. ;
	}
	long	pt_count = 7 ;
	v_vd_n.resize(pt_count) ;
	for (size_t index=0 ; index<10 ; index++) {
		for (size_t index_d=0 ;  index_d<div_c ; index_d++) {
			double	t = double(index*div_c + index_d)/div_c ;
			for (size_t i=0 ; i<pt_count ; i++) {
				double	r = ::BS_basis(t,order,i) ;
				Vd3A	v3a = v_vd_n[i] ;
				Vd3	v(v3a.size(),r*100,0) ;
				v3a.push_back(v) ;
				v_vd_n[i] = v3a ;
				}
			}
		}

Vd4 として w に設定してみたが，これで良いのか？

Vd3A	b_spline	(const Vd4& q0,const Vd4& q1,const Vd4& q2,const Vd4& q3,const long div_c)
{
	Vd3A	vd3a ;
	for (long index=0 ; index<=div_c ; index++) {
		double	t  = 1./div_c*index ;
		double	n0 = 1./6 * (1.-t)*(1.-t)*(1.-t) ;
		double	n1 = 1./2 * t*t*t - t*t + 2./3 ;
		double	n2 =-1./2 * t*t*t + 1./2 * t*t + 1./2 *t + 1./6 ;
		double	n3 = 1./6 * t*t*t ;
		double	px = ( (n0*q0.x)*q0.w + (n1*q1.x)*q1.w + (n2*q2.x)*q2.w + (n3*q3.x)*q3.w ) / ((n0)*q0.w + (n1)*q1.w + (n2)*q2.w + (n3)*q3.w) ;
		double	py = ( (n0*q0.y)*q0.w + (n1*q1.y)*q1.w + (n2*q2.y)*q2.w + (n3*q3.y)*q3.w ) / ((n0)*q0.w + (n1)*q1.w + (n2)*q2.w + (n3)*q3.w) ;
		double	pz = ( (n0*q0.z)*q0.w + (n1*q1.z)*q1.w + (n2*q2.z)*q2.w + (n3*q3.z)*q3.w ) / ((n0)*q0.w + (n1)*q1.w + (n2)*q2.w + (n3)*q3.w) ;
		vd3a.push_back(Vd3(px,py,pz)) ;
		}
	return	vd3a ;
	}