> restart: with(plots):
>  with(linalg):
Warning, new definition for norm
Warning, new definition for trace
> subdivision:=proc(degree, points,anzahl, it)
> # Plotten einer B-Spline-Kurve mit Kontrollpunkten points
> # degree: polynomialer Grad der B-Spline-Kurve
> # anzahl: Anzahl der Kontrollpunkte
> # it : Anzahl der Iterationen 
> # option trace;
> local aa,j,k,l,dd1,N,bsplinekurve,cpoints,anz,pt,deg1,deg2 ;
> global curve;
> #Berechnung der Maskenkoeffizienten
> aa:=vector(degree+2);
> for j from 1 to degree+2 do aa[j]:=binomial(degree+1, j-1)/2^(degree) od;
> # y-Komponente
> N:=2^(it)*(anzahl+degree)-degree;
> curve:=array(1..N,1..2);
> dd1:=vector(N,0);
> for j from 1 to anzahl do curve[j,2]:=points[j] od;
> anz:=anzahl;
> for k from 1 to it do 
>    deg1:=floor(degree/2); deg2:=floor((degree-1)/2);
>    for j from 0 to anz+deg2 do
>      dd1[2*j+1]:=0;
>      for l from 0 to deg2+1 do if (j-l+1 >0) and (j-l<anz) then dd1[2*j+1]:=dd1[2*j+1]+curve[j-l+1,2]*aa[2*l+1] fi od
>    od;
>    for j from 0 to anz+deg1-1 do
>      dd1[2*j+2]:=0;
>      for l from 0 to deg1 do if (j-l+1 >0) and (j-l<anz) then dd1[2*j+2]:=dd1[2*j+2]+curve[j-l+1,2]*aa[2*l+2] fi od
>    od;
>    for j from 1 to 2*anz+degree do curve[j,2]:=dd1[j] od;
>    anz:=2*anz+degree;
> od;
> # x-Komponente
> for j from 1 to N do curve[j,1]:=(j+(degree+1)/2-1)/2^(it) od;
> bsplinekurve:=plot(convert(curve,listlist)):
> pt:=array(1..anzahl+2,1..2);
> for j from 1 to anzahl do pt[j+1,1]:=evalf(j+(degree-1)/2); pt[j+1,2]:=points[j] od;
> pt[1,2]:=0; pt[anzahl+2,2]:=0; pt[1,1]:=(degree-1)/2; pt[anzahl+2,1]:=anzahl+1+(degree-1)/2;
> cpoints:=listplot([seq([pt[j,1],pt[j,2]],j=1..anzahl+2)],color=blue):
> RETURN (display([bsplinekurve,cpoints])):
> end:
> P:=[-1,3,3,-2,-3];
> subdivision(3,P,5,3);
> 
> 
> 

                       P := [-1, 3, 3, -2, -3]


> P:=[1];
> subdivision(4,P,1,4);

                               P := [1]


>