%Look below for the word ``PCTeX'' for a special remark.
%
%                    Anders Thorup
%                 PSPICTURE (for DVIPS)
%
%Syntax:
%
%    \PSpicture<xdimen><ydimen><picture material>\endPSpicture
%
%

\chardef\atcode=\catcode`\@
\catcode`\@=11

\def\PSpicture#1#2{\vbox to #2\bgroup\gdef\PSstr@ng{"}%
  \offinterlineskip
  \hrule height0pt width #1\relax\vfill\PSpictsetup}
\def\endPSpicture{\special{\PSstr@ng}\egroup}

%\PSpicture starts the "picture box", a \vbox of width xdimen and 
%height ydimen. Most commands of the PICTURE MATERIAL will result in
%commands for the PostScript interpreter. These commands are appended
%to the \PSstr@ng. 
%The \endPSpicture ships the \PSstr@ng with the \special
%command, and ends the box.

%Inside the picture box, coordinates and lengths are dimensionless.
%They are interpreted as multiples of \PSunit. 
\newdimen\PSunit \PSunit=1bp

%The lower left vertex of the picture box has coordinates (0,0). 
%To change the coordinates of the lower left vertex, 
%give the command \llcoords{xcoord ycoord} as the very first
%command in PSPICTURE.
\def\llcoords#1{\begingroup\c@l=0 \st@ck#1 \end
  \ifnum\c@l=2 \else
      \errmessage{Two coordinates needed for \string\llcoords}\fi
   \dimen@x=\ll@x\dimen@y=\ll@y
   \global\ll@y=-\Acol \advance\c@l-1 \global\ll@x=-\Acol
   \advance\ll@x-\dimen@x \advance\ll@y-\dimen@y
   \PS{\The\ll@x\space\The\ll@y\space translate}\endgroup}

%Default line width for lines and curves drawn in PSPICTURE:
\newdimen\PSlinewidth %\PSlinewidth=0.5pt
\PSlinewidth=1pt

%The \PSlinewidth can be changed inside PSPICTURE. The change
%affects the following paths. Note that, contrary to the other
%commands inside PSpicture, the argument of \setPSlinewidth is a
%(usual) TeX dimension.
\def\setPSlinewidth#1{\PSlinewidth=#1 %
   \PS{\The\PSlinewidth\space\The\PSunit\space div setlinewidth}}

%The picture material to go inside the box could be direct commands
%to the PostScript interpreter, set with \PS. Example:
%          \PS{0 0 moveto 1.2 2.0 lineto stroke}
%will produce a line from (0,0) to (1.2,2.0).
\def\PS#1{\global\let\tmp=\PSstr@ng
   \global\edef\PSstr@ng{\tmp\space#1}}

%The picture setup macro allocates a couple of local registers,
%it sets the origin of the default coordinate system,
%and passes the global scale factor and the line width to PS
\def\PSpictsetup{\newlocregs \ll@x=0pt \ll@y=0pt %
  \count@=\PSunit\PS{\the\count@\space 65781.76 div dup scale}%
  \setPSlinewidth\PSlinewidth
  \defPSarrowh
  \setPSarrowheadlength{10\PSlinewidth}}

%The picture material also includes (horizontal) TeX material to be
%set in an \hbox OF WIDTH 0 at specified coordinates with the command
%\TX. The box is centered vertically. Example:
%         \TX{$A$\hss}{0 1}
%will produce a math A at the point with coordinates (0,1). 
%The trailing \hss will set the A right to the point.
\def\TX#1#2{\begingroup\c@l=0 \st@ck#2 \end
  \ifnum\c@l=2 \else
     \errmessage{Two coordinates needed for \string\TX}\fi
  \advance\ll@y\Acol \advance\c@l-1 \advance\ll@x\Acol
  \ll@y=\The\ll@y\PSunit \ll@x=\The\ll@x\PSunit
  \vbox to 0pt{\kern-\ll@y\vss
       \hbox to 0pt{\kern\ll@x#1\kern-\ll@x}
       \vss\kern\ll@y}\endgroup}

%The following construct:
%         \TX\markbox{1 2}
%will produce a small black box centered at the point (1,2).
\def\markbox{\hss\vrule height 3\PSlinewidth 
      width 3\PSlinewidth \hss}

%The main purpose of this macro package is to define the macros
%\spline and \cspline. The command inside PSpicture:
%    \spline{1 1.2 4 5.5 6.7 4.5 2.2 1}
%will produce a natural (third degree) spline through the points 
%         (1 , 1.2) (4 , 5.5) (6.7 , 4.5) (2.2 , 1). 
%and \cspline (with the same arguments) will produce a closed spline.
%The number of points has to be at least 3 (and the number of space
%separated coordinates has to be even and at least equal to 6).
%The best result for the natural spline is obtained when the points
%are evenly distributed along the curve.
%
%The spline is broken into Bezier curves drawn by PostScript. 
%To do so, we let TeX compute the coordinates of the control points,
%using the algorithms described in the text following "endinput" below.
%
%Input to the algorithm is n points A1,...,An
%where each point is a pair of coordinates Ai=Bi,Ci. They are stored
%in 2n registers Aj; the odd Aj's are the Bi and the even Aj's are the
%Ci.
%
%Output is the set of n control points Zi=Xi,Yi. They are stored in
%2n registers Zj; the odd Zj's are the Xi and the even are the Yi.
%
%TeX is bad at arithmetics. In the code, each factor fi is between
%2-\sqrt3 and 1. We store fi in an integer register as Fi=fi*base,
%where base=2^15, approx the square root of maxint.

\newcount\base \base=32768 
%Set instead \base=10000 for test, or for (old versions of) PCTeX!!
%A TeX count register should be able to hold a number less than the
%number 2^{31}, and in particular \base*\base. This is not the case
%for the TeX version produced by PCTeX, so for use with PCTeX you
%have to decrease the \base, resulting in loss of accuracy.

%The multiplication by a factor, Z:=Z*f in the algorithm, is
%interpreted as Z:=Z*F/B, implemented as follows: 
%Split Z=Zh*B+Zl. Then Z*F/B=Zh*F+Zl*F/B, that is:
%  Zh:=Z/B, Zl:=Z-Zh*B, Z:=Zh*F+Zl*F/B
\def\m@ltf#1#2{\dimen@y=#1\divide\dimen@y\base
  \dimen@x=-\dimen@y \multiply\dimen@x\base\advance\dimen@x#1%
  \multiply\dimen@y#2\multiply\dimen@x#2\divide\dimen@x\base
  \advance\dimen@y\dimen@x#1=\dimen@y}

%The inversion f:=(m-f)^{-1} (where m is 2 or 4) is interpreted 
%as F=B^2/(m*B-F)*B, implemented as follows:
\def\inv@rt#1#2{\count@@=#1\multiply\count@@\base 
  \advance\count@@-#2#2=\base \multiply#2\base \divide#2\count@@}

%We want to allocate (locally) registers:
\def\newlocdimen#1{\advance\count11 by \@ne
  \ifnum\count11 <\insc@unt
     \else\errmessage{No room for the spline}\fi
  \dimendef#1=\count11 }
\def\newloccount#1{\advance\count10 by \@ne
  \ifnum\count10 <\insc@unt
     \else\errmessage{No room for the splines}\fi
  \countdef#1=\count10 }
\def\newdimenvar#1{\newlocdimen\currvar
   \expandafter\let\csname#1\the\c@l\endcsname=\currvar
   \wlog{\expandafter\string\csname#1\the\c@l\endcsname
      =\string\dimen\the\count11 }}
\def\newcountvar#1{\newloccount\currvar
   \expandafter\let\csname#1\the\c@l\endcsname=\currvar
  \wlog{\expandafter\string\csname#1\the\c@l\endcsname
   =\string\count\the\count10 }}

\def\newlocregs{\newloccount\c@l\newloccount\numc@ls
  \newloccount\prevF\newloccount\prevQ
  \newloccount\count@\newloccount\count@@
  \newlocdimen\dimen@\newlocdimen\dimen@x\newlocdimen\dimen@y
  \newlocdimen\ll@x\newlocdimen\ll@y}

\def\Acol{\csname A\the\c@l\endcsname}
\def\Fcol{\csname F\the\c@l\endcsname}
\def\Qcol{\csname Q\the\c@l\endcsname}
\def\Zcol{\csname Z\the\c@l\endcsname}

%Def of the two spline macros:
\def\spline{\begingroup\spl@ne}
\def\cspline{\begingroup\let\d@spline=\d@cspline
   \let\PSw@finish=\PSw@cfinish\spl@ne}
\def\spl@ne#1{\maybem@rks\c@l=0 \st@ck#1 \end
   \ch@ckeven\ch@ckthree
   \d@spline\d@PSw\PSw@finish\stroke\endgroup}
\def\ch@ckthree{\ifnum\numc@ls<6 %
      \errmessage{At lest three points needed for SPLINE}\fi}
\def\ch@cktwo{\ifnum\numc@ls<4 %
      \errmessage{At lest two points needed}\fi}
\def\ch@ckeven{\ifodd\numc@ls 
       \errmessage{Even number of coordinates needed}\fi}
\def\d@spline{\fwd@init\loop\fwd@body\ifnum\c@l<\numc@ls\repeat
  \fwd@finish\s@lve}
\def\d@cspline{\fwd@cinit\loop\fwd@cbody\ifnum\c@l<\numc@ls\repeat
  \fwd@cfinish\cs@lve}


%The coordinates are read and put into the dimen registers \Acol
%by \st@ck ... \end. Each coordinate is a real in decimal
%notation, and they are separated by one or more spaces. 
%The number of coordinates read is stored in \numc@ls
{\def\\{\global\let\sptoken= }\\ } %now \sptoken is a space token
\def\st@ck{\futurelet\nexttok\testst@ck}
\def\testst@ck{\ifx\nexttok\sptoken
   \let\next\eatspacest@ck
   \else\ifx\nexttok\end\let\next\stackfinish
   \else\let\next\st@ckit\fi\fi\next}
\def\eatspacest@ck#1 {\st@ck}
\def\stackfinish#1{\numc@ls=\c@l}
\def\st@ckit#1 {\advance\c@l1 \newdimenvar A\Acol=#1pt %
  \markoptions\st@ck}

%After having read the arguments, \numc@ls=\c@l is equal to 2n.
%The linear equations are solved by Gauss elimination, see the text
%after "endinput".
%Initialization of the forward loops:
\def\fwd@init{%        %set F2=\prevF:=base/2, Z1:=A1+2*A3, Z2:=A2+2*A4
  \c@l=1 \newdimenvar Z\currvar=\Acol 
    \advance\c@l2 \advance\currvar 2\Acol
  \c@l=2 \prevF=\base \divide\prevF2 \newcountvar F\Fcol=\prevF
    \newdimenvar Z\currvar=\Acol 
    \advance\c@l2 \advance\currvar 2\Acol
  \c@l=4 }
\def\fwd@cinit{%       %Q1=\prevQ:=base, F2=\prevF:=base/2, 
                   %Z1:=4A1+2A3, Z2:=4A2+2A4
                   %lastF=4*base, lastX=4A{2n-1}+2A1, lastY=4A{2n}+2A2
  \c@l=\numc@ls
   \newcountvar F\let\lastF=\currvar \lastF\base \multiply\lastF4 %
   \newdimenvar Z\let\lastY=\currvar \lastY=4\Acol 
  \advance\c@l-1 \newdimenvar Z\let\lastX=\currvar\lastX=4\Acol 
  \c@l=1 \advance\lastX 2\Acol 
   \prevQ=\base\newcountvar Q\Qcol=\prevQ
   \newdimenvar Z\Zcol=4\Acol\advance\c@l2 \advance\currvar 2\Acol
  \c@l=2 \advance\lastY 2\Acol
   \prevF=\base \divide\prevF4 \newcountvar F\Fcol=\prevF
   \newdimenvar Z\Zcol=4\Acol\advance\c@l2 \advance\currvar 2\Acol
  \c@l=4 }

%The bodies of the two forward loops:
\def\fwd@body{%      %Zi=:-Z{i-2}*\prevF+4Ai+2A{i+2}. If i is even,
                    %Fi=\prevF:=(4-\prevF)^{-1}
                    %At entrance, \c@l=i+1, at exit \c@l=i+2
  \advance\c@l-1 \newdimenvar Z%
  \advance\c@l-2 \currvar=-\Zcol \m@ltf\currvar\prevF 
  \advance\c@l2 \advance\currvar 4\Acol
  \advance\c@l2 \advance\currvar 2\Acol 
  \advance\c@l-2 \ifodd\c@l\else
     \inv@rt4\prevF\newcountvar F\Fcol=\prevF\fi
  \advance\c@l2 }
\def\fwd@cbody{%     %Zi=:-Z{i-2}*\prevF+4Ai+2A{i+2}. If i is even,
                    %\tmpQ=\prevQ
                    %Q{i-1}=\prevQ:=-\prevQ*\prevF
                    %\lastF:=\lastF+\prevQ*\tmpQ
                    %\lastX:=\lastX+\prevQ*Z{i-1}
                    %\lastY:=\lastY+\prevQ*Zi
                    %\prevQ:=\newQ
                    %Fi=\prevF:=(4-\prevF)^{-1}
  \advance\c@l-1 \newdimenvar Z%
  \advance\c@l-2\currvar=-\Zcol \m@ltf\currvar\prevF 
  \advance\c@l2 \advance\currvar 4\Acol 
  \advance\c@l2 \advance\currvar 2\Acol
  \advance\c@l-2 \ifodd\c@l\else
     \count@=\prevQ 
     \prevQ=-\prevQ \multiply\prevQ\prevF \divide\prevQ\base
     \multiply\count@\prevQ \divide\count@\base
     \advance\lastF\count@
     \advance\c@l-2 \dimen@=\Zcol \m@ltf\dimen@\prevQ
       \advance\lastY\dimen@
     \advance\c@l-1 \dimen@=\Zcol \m@ltf\dimen@\prevQ
       \advance\lastX\dimen@
     \advance\c@l2 \newcountvar Q\Qcol=\prevQ
     \advance\c@l1 \inv@rt4\prevF\newcountvar F\Fcol=\prevF\fi
  \advance\c@l2 }

%Finishing the forward loops:
\def\fwd@finish{%     %Z{2n-1}:=3A{2n-1}-Z{2n-3}*\prevF
                     %Z{2n}:=3A{2n}-Z{2n-2}*\prevF
                     %F{2n}=\prevF:=(2-\prevF)^{-1}
  \c@l=\numc@ls
  \advance\c@l-1 \newdimenvar Z%
  \advance\c@l-2 \currvar=-\Zcol\m@ltf\currvar\prevF
  \advance\c@l2 \advance\currvar 3\Acol
  \advance\c@l1 \newdimenvar Z%
  \advance\c@l-2 \currvar=-\Zcol\m@ltf\currvar\prevF
  \advance\c@l2 \advance\currvar 3\Acol
  \inv@rt2\prevF \newcountvar F\Fcol=\prevF}
\def\fwd@cfinish{%    %\prevQ:=1+\prevQ
                     %tmpQ=-\prevQ*\prevF
                     %\lastX:=\lastX+tmpQ*Z{2n-3}
                     %\lastY:=\lastY+tmpQ*Z{2n-2}
                     %\lastF:=\lastF+tmpQ*\prevQ
                     %\lastF:=1/\lastF
  \c@l=\numc@ls \advance\prevQ\base 
    \count@=-\prevQ \multiply\count@\prevF\divide\count@\base 
  \advance\c@l-2 \dimen@=\Zcol
    \m@ltf\dimen@\count@\advance\lastY\dimen@
  \advance\c@l-1 \dimen@=\Zcol
    \m@ltf\dimen@\count@\advance\lastX\dimen@
    \multiply\count@\prevQ \divide\count@\base\advance\lastF\count@ 
  \prevF=\base \multiply\prevF\base \divide\prevF\lastF \lastF=\prevF
  \c@l=\numc@ls}

%Now the two backwards substitutions to solve the equations:
\def\s@lve{\c@l=\numc@ls \m@ltf\Zcol\prevF
  \advance\c@l-1 \m@ltf\Zcol\prevF
  \loop\advance\c@l-1 %
    \ifodd\c@l \else \prevF=\Fcol\fi
    \advance\c@l2 \dimen@=-\Zcol 
    \advance\c@l-2 \advance\Zcol\dimen@\m@ltf\Zcol\prevF
    \ifnum\c@l>1 \repeat}
\def\cs@lve{\c@l=\numc@ls \m@ltf\Zcol\prevF
  \advance\c@l-1 \m@ltf\Zcol\prevF
  \loop\advance\c@l-1 %
    \ifodd\c@l \dimen@=-\lastX 
         \else \dimen@=-\lastY \prevF=\Fcol
               \advance\c@l-1 \prevQ=\Qcol\advance\c@l1 \fi
    \m@ltf\dimen@\prevQ
    \advance\c@l2 \advance\dimen@-\Zcol
    \advance\c@l-2 \advance\dimen@\Zcol \m@ltf\dimen@\prevF 
    \Zcol=\dimen@ \ifnum\c@l>1\repeat}


%Next we come to the PostScript writing:
%Dirty tricks, the TeX book page 375:
{\catcode`p=12 \catcode`t=12 \gdef\STRIPPT#1pt{#1}}
\def\The#1{\expandafter\STRIPPT\the#1}

%The start of the writing is identical for both splines:
\def\PSw@#1{\advance\c@l-1 \PS{\The#1}\advance\c@l1 \PS{\The#1}}
\def\PSw@Wcol{\advance\c@l1 \dimen@x=2\Acol\advance\dimen@x-\Zcol
  \advance\c@l1 \dimen@y=2\Acol\advance\dimen@y-\Zcol
  \PS{\The\dimen@x\space\The\dimen@y}}
\def\PSw@init{\c@l=2 \newpath\PSw@\Acol\moveto}
\def\PSw@ZWA{%write: Zi   Wi=2A{i+1}-Z{i+1}   A{i+1}   curveto
             %at entrance, \c@l=2i, at exit \c@l=\c@l+2
  \PSw@\Zcol\PSw@Wcol\PSw@\Acol\PS{curveto}}
%The finish is different:
\def\PSw@finish{\relax}
\def\PSw@cfinish{\PSw@\Zcol
  \c@l=0 \PSw@Wcol\PSw@\Acol\PS{curveto}}
\def\d@PSw{\PSw@init
  \loop\PSw@ZWA\ifnum\c@l<\numc@ls\repeat}
\def\normalpaths{\def\stroke{\PS{stroke}}%
  \def\moveto{\PS{moveto}}\def\newpath{\PS{newpath}}}
%The default is normal paths:
\normalpaths


\def\PSline{\begingroup\PSl@ne}
\def\PSl@ne#1{\maybem@rks\c@l=0 \st@ck#1 \end
  \ch@ckeven\ch@cktwo
  \c@l=2 \newpath\PSw@\Acol\moveto
  \loop\advance\c@l2 \PSw@\Acol\PS{lineto}%
  \ifnum\c@l<\numc@ls\repeat\PSl@nefinish\stroke\endgroup}
\let\PSl@nefinish=\relax
\def\PScline{\begingroup\def\PSl@nefinish{\PS{closepath}}\PSl@ne}

\def\PSarc{\begingroup\PS@rc}
\def\PSarcn{\begingroup\let\PS@arcfinish=\PS@arcnfinish\PS@rc}
\def\PS@rc#1#2{\c@l=0 \st@ck#2 #1 \end
  \ifnum\c@l=5 \else
     \errmessage{Wrong number of arguments for circles/arcs}\fi
  \newpath\c@l=0 \loop
       \advance\c@l1 \PS{\The\Acol}\ifnum\c@l<5 \repeat
  \PS@arcfinish\stroke\endgroup}
\def\PS@arcfinish{\PS{arc}}
\def\PS@arcnfinish{\PS{arcn}}
\def\PScircle{\PSarc{0 360}}

%Various marks:
\let\markoptions=\relax
\def\putTeXmarks{\def\maybem@rks{\let\markoptions=\TeXm@rks}}
\def\putPSmarks{\def\maybem@rks{\let\markoptions=\PSm@rks}}
\def\nomarks{\let\maybem@rks=\relax}
 %make nomarks the global default
\nomarks
\def\TeXm@rks{\ifodd\c@l\else
  \begingroup \advance\ll@y\Acol \advance\c@l-1 \advance\ll@x\Acol
  \ll@y=\The\ll@y\PSunit \ll@x=\The\ll@x\PSunit
  \vbox to 0pt{\kern-\ll@y\vss
       \hbox to 0pt{\kern\ll@x\markbox\kern-\ll@x}
       \vss\kern\ll@y}\endgroup\fi}
\let\PSmarkradius=\PSlinewidth
\def\PSm@rks{\ifodd\c@l\else\PS{newpath}\PSw@\Acol
   \PS{\The\PSmarkradius\space\The\PSunit\space div 
         0 360 arc stroke}\fi}
\def\TeXmark{\begingroup\let\markoptions=\TeXm@rks\m@rk}
\def\PSmark{\begingroup\let\markoptions=\PSm@rks\m@rk}
\def\m@rk#1{\c@l=0 \st@ck#1 \end\ch@ckeven\endgroup}

\def\setdash#1{\dimen@=#1\relax
  \PS{[\The\dimen@\space\The\PSunit\space div] 0 setdash}}
\def\nodashes{\PS{[] 0 setdash}}

\def\longpath#1{\c@l=0 \st@ck#1 \end
  \ifnum\c@l=2 \else
     \errmessage{Wrong number of arguments for \string\longpath}\fi
  \PSw@\Acol\PS{moveto}%
  \let\newpath\relax\let\stroke\relax\def\moveto{\PS{lineto}}}
\def\shortpaths{\def\moveto{\PS{moveto}}}
\def\subpaths{\PS{newpath}\let\newpath\relax
   \let\stroke\relax\def\moveto{\PS{moveto}}}

%arrows:
\def\defPSarrowh{\PS{/arrowh {
 /curmtx matrix currentmatrix def
 currentpoint translate
 arrlen dup scale
 rotate
 -1  0
 -2.5  0.5
 -3  1
 curveto
 -2.5  0.5
 -2.5  -0.5
 -3  -1
 curveto
 -2.5  -0.5
 -1  0
 0 0
 curveto closepath fill
 curmtx setmatrix} def}}
\def\setPSarrowheadlength#1{\dimen@=#1\relax
  \PS{/arrlen \The\dimen@\space\The\PSunit\space div 3 div def}}

\def\PSarrow#1{\begingroup
  \c@l=0 \st@ck#1 \end
  \ifnum\c@l=4 \else
  \errmessage{4 coordinates needed for \string\PSarrow}\fi
  \PS{newpath}\c@l=2 \PSw@\Acol\PS{moveto}\c@l=4 \PSw@\Acol
  \PS{lineto}\stroke\PS{newpath}\PSw@\Acol\PS{moveto}
  \dimen@y=\Acol\advance\c@l-1 \dimen@x=\Acol
  \advance\c@l-1 \advance\dimen@y-\Acol
  \advance\c@l-1 \advance\dimen@x-\Acol
  \ifdim\dimen@x=0pt\ifdim\dimen@y=0pt
      \errmessage{equal endpoints of \string\PSarrow}\fi\fi
  \PS{\The\dimen@y\space\The\dimen@x\space atan arrowh}\endgroup}

\catcode`\@=\atcode

%END
\PSunit=1truecm
\def\LTX#1{\TX{\hss #1}}
\def\RTX#1{\TX{#1\hss}}
\def\CTX#1{\TX{\hss #1\hss}}
\endinput

\def\uncatcodespecials{\def\do##1{\catcode`##1=12 }\dospecials}
\def\setupverbatim{\tt\def\par{\leavevmode\endgraf}\obeylines
    \uncatcodespecials \obeyspaces}
\def\verbatim{\begingroup\setupverbatim\doverbatim}
{\obeyspaces\global\let =\ }
\def\doverbatim#1{\def\next##1#1{##1\endgroup}\next}
\def\listverbatim#1{\par\begingroup\setupverbatim\input#1 \endgroup}

\magnification=\magstep1
\parindent=0pt \parskip=\smallskipamount
\hsize=15,5truecm \vsize=24truecm
\def\head#1{\bigskip\noindent{\bf#1.}\smallskip}
\def\bull{\medskip\noindent\hbox{$\bullet$\ \ }}


\centerline{ANDERS THORUP: PSpictures for DVIPS, 5 April 1995}


\head{Introduction}

 Let $a_1,\dots,a_n$ be $n$ points in the plane. A {\it spline\/}
through the $n$ points is a sequence of curves parameterized by
polynomials $a_i(t)$ for $i=1,\dots,n-1$, such that the curve
$a_i(t)$ begins at $a_i$ and ends at $a_{i+1}$ and such that the
total curve from $a_1$ to $a_n$ is smooth.  

It would be desirable to be able to include a spline in a \TeX{}
manuscript with a command like the following:

\hskip2truecm\verbatim|\spline{|$\,a_1\ \dots\ a_n\,$\verbatim|}|

 where the $n$ points $a_i$ should be replaced by $n$ pairs of
coordinates.  Naturally, the command should be packed into an
environment giving an interpretation of the coordinates as multiples
of a given unit length.  The result of the command should be a smooth
curve through the points on the paper determined by the $n$ pairs of
coordinates. 

\head{The Package}

This manual is the printout of a \TeX{} file containing splines. The
first part of the file is a package defining the necessary macros.
It uses Bezier cubics as the polynomials $a_i(t)$, and it leaves the
work of drawing the curves to the PostScript printer. To do so,
PostScript needs, in addition to the coordinates of the $n$ points,
the coordinates of the so called control points. Each Bezier cubic is
determined by two control points. For a spline through $n$ points,
the coordinates of the control points are determined by two systems
of linear equations with $n$ unknowns, see the section on splines at
the end of the manual. The main innovation of the package is to let
\TeX\ solve the two systems of linear equations.

In addition, the package defines several other commands for curve
drawing. To use the package as an input file in other \TeX{}
documents, find the following line in the file:

\hskip2truecm\verbatim|%\endinput|

and erase the leading \verbatim|%|.  Then save the file in your
favorite input directory. The macros of the package have been tested
with plain \TeX, with LaTeX, and with AmSTeX.  It should be
emphasized that the code resulting from the package is {\it not\/}
portable. It assumes that DVIPS is used for transforming the dvi file
into a PostScript printable file.

\head{The Macros}
\vskip-\bigskipamount

\bull 
 The picture command sets up the overall environment in which all
curve drawing commands have to be included. It has the following form:

 \smallskip
\hskip2truecm\verbatim|\PSpicture{|%
 {\it xdimen}\verbatim|}{|{\it ydimen}\verbatim|}|
  {\it picture-material} \verbatim|\endPSpicture|
 \smallskip
 
The result of the command is a vertical (picture) box of width {\it
xdimen\/} and height {\it ydimen\/} containing the {\it
picture-material}.  The {\it picture-material\/} consists of various
picture commands setting objects in a coordinate system; by default,
the origin of the coordinate system is the lower left vertex of the
picture box.

\bull
 Coordinates (and lengths) inside the picture environment are
interpreted as multiples of \verbatim|\PSunit|. The default value
is 1bp (1in=72bp),
 \smallskip
\hskip2truecm\verbatim|\PSunit=1bp|
 \smallskip
 It can be reset to any desired value outside the PSpicture
environment. Inside the picture environment, coordinates and lengths
are real numbers in decimal notation. Numerically, the maximal
acceptable coordinate is $2^{15}-1$. It is easy to let the two
arguments, {\it xdimen\/} and {\it ydimen\/}, determining the
dimensions of the picture box, depend on the \verbatim|\PSunit|.  For
instance, \verbatim|\PSpicture{3.5\PSunit}{2\PSunit}| is a valid
start of the PSpicture command.

\bull
Splines are drawn with \verbatim|\spline| and \verbatim|\cspline|:

 \smallskip
\hskip2truecm\verbatim|\spline{|{\it point-list}\verbatim|}|
 \hskip1truecm\verbatim|\cspline{|{\it point-list}\verbatim|}|
 \smallskip
 The {\it point-list\/} is a space separated list of the coordinates:
{\it xcoord ycoord xcoord ycoord \dots ycoord}. In particular, the
number of coordinates has to be even. In addition, it has to be at
least $6$, corresponding to the coordinates of at least $3$ points.
The effect of \verbatim|\spline| is to draw a natural spline through
the given points; \verbatim|\cspline| draws a closed spline.

The spline commands make heavy use of the \TeX{} registers, and
\TeX{} will spend quite some time on the computations involved when
the number of points is large.  For a spline through $n$ points, the
number of dimension registers used is $4n+5$. \TeX{} has only 256
dimension registers, and some are in use when the spline command is
started. So, depending on the \TeX{} dialect in use, the maximal
number of points for each spline command is less than 50.
The spline commands will tell you if your number of points exceeds the
maximum. 

\bull
Text may be inserted in the picture by the command \verbatim|\TX|:
 \smallskip
\hskip2truecm\verbatim|\TX{|{\it TeX-stuff-of-width-zero}%
 \verbatim|}{|{\it xcoord ycoord}\verbatim|}|
 \smallskip
 The command will put the {\it TeX-stuff\/} at the point with
coordinates ({\it xcoord,ycoord\/}). The {\it TeX-stuff\/} is any
horizontal material of width zero; so is has to include some
shrinkability.  For instance, the command %
\verbatim|\TX{ $a_i$\hss}{1 2}|
 will place $a_i$ (preceded by a space) right to the point with
coordinates $(1,2)$. With a leading \verbatim|\hss| instead, the
$a_i$ would be placed to the left, and with both a leading and a
trailing \verbatim|\hss| it would be centered. The box containing the
{\it TeX-stuff\/} will be centered vertically around the specified
point. 

\bull
 Circles may be drawn with \verbatim|\PScircle|:
 \smallskip
\hskip2truecm\verbatim|\PScircle{|{\it xcoord ycoord radius}\verbatim|}|
 \smallskip
 The center of the circle will be ({\it xcoord,ycoord}). Us usual,
the {\it radius\/} is interpreted as a multiple of
\verbatim|\PSunit|. 

\bull
 More generally, arcs may be drawn with \verbatim|\PSarc|:
 \smallskip
\hskip2truecm\verbatim|\PSarc{|{\it degree$_1$ degree$_2$}%
 \verbatim|}{|{\it xcoord ycoord radius}\verbatim|}|
 \smallskip
 Thus \verbatim|\PSarc{300 60}{0 0 1}| will draw the arc of the
circle (counter clockwise oriented) from the point with angle $-\pi
/3$ to the point with angle $\pi /3$. The command
\verbatim|\PScircle| is equivalent to \verbatim|\PSarc{0 360}|.
\verbatim|\PSarcn| draws clockwise.

\bull
Piece wise linear curves may be  drawn by
\verbatim|\PSline| and \verbatim|\PScline|:  
 \smallskip
\hskip2truecm\verbatim|\PSline{|{\it point-list}\verbatim|}|
 \hskip1truecm\verbatim|\PScline{|{\it point-list}\verbatim|}|
 \smallskip
 The {\it point-list\/} has to consists of the coordinates of at
least two points. For instance, \verbatim|\PSline{1 0 0 1 0 0 1 1}|,
in the picture environment and with \verbatim|\PSunit=8pt|, produces
this symbol:
 {\PSunit=8pt\PSpicture\PSunit\PSunit 
 \PSline{1 0 0 1 0 0 1 1}\endPSpicture}. 
 With only two points in the list, that is, with four coordinates as
argument, \verbatim|\PSline| draws a straight line.  It is an error
to give only one point in the list. A closed polygon is drawn with
\verbatim|\PScline|. 

\bull
 An arrow from one point to another may be drawn with
\verbatim|\PSarrow|:
 \smallskip
\hskip2truecm\verbatim|\PSarrow{|{\it xcoord$_1$ ycoord$_1$
 xcoord$_2$ ycoord$_2$}\verbatim|}|
 \smallskip
 The result is an arrow with its tip at the point ({\it
xcoord$_2$,ycoord$_2$}). The two points have to be different, since
otherwise the direction of the arrowhead is undefined.

\bull
 Marks can be set at the points supporting the spline commands and
the line commands.
 \smallskip
\hskip2truecm\verbatim|\putPSmarks|\qquad
     \verbatim|\putTeXmarks|
 \smallskip

 After \verbatim|\putPSmarks|, a small circle is drawn around the
points in the list of the following spline and line commands.
After \verbatim|\putTeXmarks|, a small black square is drawn. The
default is to set no marks, and it can be obtained by
\verbatim|\nomarks|. 

\bull
A list of points, and in particular a single point, can be
marked using
 \smallskip
\hskip2truecm\verbatim|\TeXmark{|{\it point-list}\verbatim|}|\qquad
    \verbatim|\PSmark{|{\it point-list}\verbatim|}|
 \smallskip
 The result is a mark set at each point in the list.

\bull
 By default, the lower left point of the picture box is the origin of
the coordinate system. It can be reset with the command,
\smallskip
\hskip2truecm\verbatim|\llcoords{|{\it xcoord ycoord}\verbatim|}|
\smallskip
 Thus, after \verbatim|\llcords{2 -0.5}|, the lower left point of the
picture box is the point with coordinates $(0,-1/2)$. The change
influences the following picture commands. Normally, the
\verbatim|\llcoords| command, if present, should be one of the first
commands of the picture environment. However, if parts of the picture
commands refer to different coordinate systems, it may be convenient
to use several \verbatim|\llcoords| commands.


\bull
Curves resulting from picture commands have a thickness equal to the
value of \verbatim|\PSlinewidth|. By default,
 \smallskip
\hskip2truecm\verbatim|\PSlinewidth=0.5pt|
 \smallskip
 It can be given a different value globally. Inside the picture
environment, it can be reset with \verbatim|\setPSlinewidth|:
 \smallskip
\hskip2truecm\verbatim|\setPSlinewidth{|{\it dimen}\verbatim|}|
 \smallskip
 The result is a change in the value of \verbatim|\PSlinewidth|, and
it influences the following picture commands in the same picture
environment. Note that the argument to \verbatim|\setPSlinewidth| is
a usual \TeX{} dimension. 


\bull
 The curves generated by the picture commands can be dashed.
The command,
 \smallskip
\hskip2truecm
 \verbatim|\setdash{|{\it dimen}\verbatim|}|
 \smallskip
 will make the following curves in the picture environment appear as
dashed curves with dashes of length equal to {\it dimen}. For
instance, after \verbatim|\setdash{5\PSlinewidth}|, the dashes will
have a length equal to $5$ times the linewidth. The default is
reestablished with
 \smallskip
\hskip2truecm \verbatim|\nodashes|
 \smallskip

\bull
 The size of the arrowhead produced with \verbatim|\PSarrow| can be
changed: 
 \smallskip
\hskip2truecm
 \verbatim|\setPSarrowheadlength{|{\it dimen}\verbatim|}|
 \smallskip
 As a default, the length of the arrow head is equal to $10$ times the
\verbatim|\PSlinewidth|. 

\bull
 With \verbatim|\putTeXmarks| or \verbatim|\TeXmark|, the points are
actually marked with the content of \verbatim|\markbox| defined as
 \smallskip
\hskip2truecm
 \verbatim|\def\markbox{\hss\vrule height 3\PSlinewidth 
                                   width 3\PSlinewidth\hss}|
 \smallskip
 It is easy to redefine \verbatim|\markbox| to produce different
marks. With \verbatim|\putPSmarks|, each point is marked with a small
circle of radius equal to \verbatim|\PSmarkradius|. By default,
 \smallskip
\hskip2truecm
 \verbatim|\let\PSmarkradius=\PSlinewidth|
 \smallskip
 In other words, as a dimension register, \verbatim|\PSmarkradius| is
equal to the register \verbatim|\PSlinewidth|.  As a result, the
radius will vary with the line width. Hence, if you change the
\verbatim|\PSlindewidth| with \verbatim|\setPSlinewidth|, you will
get a proportional change of the marks. If you want a different
radius, you have either to let \verbatim|\PSmarkradius| be equal to a
different dimension register, or you have to declare
\verbatim|\PSradius| as a dimension register (with
\verbatim|\newdimen|) and give it the appropriate value. 

\bull
 The PostScript interpreter can be addressed directly with the command
\verbatim|\PS|:
 \smallskip
\hskip2truecm\verbatim|\PS{|{\it PostScript commands}\verbatim|}|
 \smallskip
 For instance, the command

 \hskip2truecm
 \verbatim|\PS{newpath 0 1 moveto 2 3 4 5 6 7 curveto stroke}|

 will draw a Bezier cubic from $(0,1)$ to $(6,7)$. It is easy to make
syntax errors addressing PostScript directly. The result of a syntax
error will most probably be that the printing process is stopped.

\bull
 Actually, each curve generating command creates a path in the
PostScript memory. At the beginning of the command, the control
sequence \verbatim|\newpath| is executed and at the end
\verbatim|\stroke| is called. The default effects of
\verbatim|\newpath| and \verbatim|\stroke| are, respectively, to
start a new path and to stroke the current path. In fact,
 \smallskip
\hskip2truecm\verbatim|\def\newpath{\PS{newpath}}|\qquad
   \verbatim|\def\stroke{\PS{stroke}}|
 \smallskip
 However, with some knowledge of the PostScript language, you may
obtain different effects.  For instance, after

\hskip2truecm
 \verbatim|\def\stroke{\PS{gsave 0.8 setgray fill grestore stroke}}|

each of the following curves generated by the picture commands will
have their enclosed area filled with gray before the curve is
stroked. Also, with temporary settings like
\verbatim|\let\newpath=\relax|, \verbatim|\let\stroke=\relax|, you
may combine several curves into a single path object. Three commands
deal directly with combinations of paths:
 \smallskip
\hskip2truecm\verbatim|\subpaths|\qquad
  \verbatim|\longpath{|{\it xcoord ycoord}\verbatim|}|\qquad
  \verbatim|\shortpaths|
 \smallskip
 After \verbatim|\subpaths|, each of the following curve drawing
commands will add a subpath to the current path. After
\verbatim|\subpaths|, the command \verbatim|\longpath| starts a new
subpath at the point ({\it xcoord, ycoord}), and the following curve
drawing commands will append to the subpath: the last point of one
curve drawing command will be connected by a straight line to the
first point of the next. Again, after \verbatim|\shortpaths|, each
command will contribute with a single subpath. At the end, you have
to instruct PostScript explicitly on what to do with the path, for
instance to fill it with \verbatim|\PS{fill}|, or to stroke it with
\verbatim|\PS{stroke}|.  In addition, you have to restore the
defaults with
 \smallskip
\hskip2truecm\verbatim|\normalpaths|
 \smallskip

\bull 
 The picture commands are meant to be used inside the picture environment.
 If used outside, most of them will provoke a syntax error, typically
of the form 
  
\hskip2truecm\verbatim|! Undefined control sequence.  
                   ... \c@l|

\vfill\eject

\head{Examples}

\bull
A closed spline through $4$ points on a circle, and through $6$
points:

\medskip
\centerline{\PSunit=1truecm \PSpicture{6\PSunit}{2\PSunit}
\putPSmarks
\llcoords{-1 -1}
\cspline{
1 0
0 1
-1 0
0 -1
}
\llcoords{-5 -1}
\putTeXmarks
\cspline{
1    0
0.5  0.8660
-0.5 0.8660
-1   0
-0.5 -0.8660
0.5  -0.8660}
\endPSpicture}
 \medskip
The drawing was produced with the following code:

{\parskip=0pt
\verbatim|\centerline{\PSunit=1truecm \PSpicture{6\PSunit}{2\PSunit}
\putPSmarks
\llcoords{-1 -1}
\cspline{
1 0
0 1
-1 0
0 -1}
\llcoords{-5 -1}
\putTeXmarks
\cspline{
1    0
0.5  0.8660
-0.5 0.8660
-1   0
-0.5 -0.8660
0.5  -0.8660}
\endPSpicture}
|}

The \verbatim|\centerline| centers the picture on the paper.  Note
the double use of the command \verbatim|\llcoords|. With respect to
the second circle, the lower left point of the picture box has
coordinates $(-5,-1)$. The second spline is quite close to a circle,
the first is not. The same code, with no marks, and an
additional ``true circle'' drawn with \verbatim|\PScircle{0 0 1}| for
each of the two parts,
produces the following: 

\medskip
\centerline{\PSunit=1truecm \PSpicture{6\PSunit}{2\PSunit}
\llcoords{-1 -1}
\PScircle{0 0 1}
\cspline{
1 0
0 1
-1 0
0 -1
}
\llcoords{-5 -1}
\PScircle{0 0 1}
\cspline{
1    0
0.5  0.8660
-0.5 0.8660
-1   0
-0.5 -0.8660
0.5  -0.8660}
\endPSpicture}
 \medskip

\bull The sine function from $0$ to $2\pi$, drawn as a natural spline
through $10$ points of distance $2\pi /9$. The coordinates of the
$10$ points were extracted from a table of the sine function.

\medskip
\centerline{\PSunit=1truecm
\PSpicture{7.2\PSunit}{2.4\PSunit}
\llcoords{0 -1.2}
\PSarrow{0 0 7 0}
\PSarrow{0 -1 0 1}
\putPSmarks
\spline{
0   0
%0.1745   0.1736
%0.3491   0.3420
%0.5236   0.5000
0.6981   0.6428
%0.8727   0.7660
%1.0472   0.8660
%1.2217   0.9397
1.3963   0.9848
%1.5708   1.0000
%1.7453   0.9848
%1.9199   0.9397
2.0944   0.8660
%2.2689   0.7660
%2.4435   0.6428
%2.6180   0.5000
2.7925   0.3420
%2.9671   0.1736
%3.1415   0.0000
%3.3160   -0.1736
3.4906   -0.3420
%3.6651   -0.5000
%3.8396   -0.6428
%4.0142   -0.7660
4.1887   -0.8660
%4.3632   -0.9397
%4.5378   -0.9848
%4.7123   -1.0000
4.8868   -0.9848
%5.0614   -0.9397
%5.2359   -0.8660
%5.4104   -0.7660
5.5850   -0.6428
%5.7595   -0.5000
%5.9340   -0.3420
%6.1086   -0.1736
6.2830   -0.0000
}
\TX{   $x$\hss}{0   1}
\TX{   $t$\hss}{7   0}
\TX{   $x=\sin   t$\hss}{5   0.5}
\endPSpicture}

\medskip
The spline for the sine function was produced by the following code:

{\parskip=0pt
\verbatim|\centerline{\PSunit=1truecm
\PSpicture{7.2\PSunit}{2.4\PSunit}
\llcoords{0 -1.2}
\PSarrow{0 0 7 0}
\PSarrow{0 -1 0 1}
\putPSmarks
\spline{
0      0
0.6981 0.6428
|\noindent\hskip2cm {\it plus 7 more pairs of coordinates: $t$ $\sin t$}

\verbatim|6.2830 0.0000
}
\TX{ $x$\hss}{0 1}
\TX{ $t$\hss}{7 0}
\TX{ $x=\sin t$\hss}{5 0.5}
\endPSpicture}
|
}

\bull
 The corresponding spline for the cosine is rather poor near
the end points:

 \medskip
\centerline{\PSunit=1truecm
\PSpicture{7.2\PSunit}{2.4\PSunit}
\llcoords{0 -1.2}
\PSline{0 0 7 0}
\PSline{0 -1 0 1}
\TX{ $x=\cos t$\hss}{5 -0.5}
\putTeXmarks
\spline{
0    1.0000
%0.1745  0.9848
%0.3491  0.9397
%0.5236 0.8660
0.6981   0.7660
%0.8727  0.6428
%1.0472  0.5000
%1.2217 0.3420
1.3963   0.1736
%1.5708  0.0000
%1.7453  -0.1736
%1.9199 -0.3420
2.0944   -0.5000
%2.2689  -0.6428
%2.4435  -0.7660
%2.6180 -0.8660
2.7925   -0.9397
%2.9671  -0.9848
%3.1415  -1.0000
%3.3160 -0.9848
3.4906   -0.9397
%3.6651  -0.8660
%3.8396  -0.7660
%4.0142 -0.6428
4.1887   -0.5000
%4.3632  -0.3420
%4.5378  -0.1736
%4.7123  0 
4.8868   0.1736
%5.0614  0.3420
%5.2359  0.5000
%5.4104  0.6428
5.5850   0.7660
%5.7595  0.8660
%5.9340  0.9397
%6.1086 0.9848
6.2830   1 
}
\endPSpicture
}
 \medskip
 The misbehavior is caused by the fact that the cosine
function has non vanishing second order derivatives at $0$ and $2\pi
$. So the natural spline gives a poor approximation at the end
points. With more supporting points the result is better:

 \medskip
\centerline{\PSunit=1truecm
\PSpicture{7.2\PSunit}{2.4\PSunit}
\llcoords{0 -1.2}
\PSline{0 0 7 0}
\PSline{0 -1 0 1}
\putPSmarks
\spline{
0    1.0000
0.1745  0.9848
0.3491  0.9397
%0.5236 0.8660
0.6981   0.7660
%0.8727  0.6428
1.0472  0.5000
%1.2217 0.3420
1.3963   0.1736
%1.5708  0.0000
1.7453  -0.1736
%1.9199 -0.3420
2.0944   -0.5000
%2.2689  -0.6428
2.4435  -0.7660
%2.6180 -0.8660
2.7925   -0.9397
%2.9671  -0.9848
3.1415  -1.0000
%3.3160 -0.9848
3.4906   -0.9397
%3.6651  -0.8660
3.8396  -0.7660
%4.0142 -0.6428
4.1887   -0.5000
%4.3632  -0.3420
4.5378  -0.1736
%4.7123  0 
4.8868   0.1736
%5.0614  0.3420
5.2359  0.5000
%5.4104  0.6428
5.5850   0.7660
%5.7595  0.8660
5.9340  0.9397
6.1086 0.9848
6.2830   1 
}
\TX{ $x=\cos t$\hss}{5 -0.5}
\endPSpicture
}

\bull
The $6$ points on the circle in a different order, with various
strokes: 



\bigskip
\centerline{\PSunit=1truecm
\PSpicture{10\PSunit}{2\PSunit}
\llcoords{-1 -1}
\putPSmarks
\cspline{
1    0
-1   0
-0.5 0.8660
0.5 -0.8660
-0.5 -0.8660
0.5  0.8660
}
\nomarks
\llcoords{-5 -1}
\setdash{5\PSlinewidth}
\cspline{
1    0
-1   0
-0.5 0.8660
0.5 -0.8660
-0.5 -0.8660
0.5  0.8660
}
\nodashes
\llcoords{-9 -1}
\def\stroke{\PS{gsave 0.8 setgray eofill grestore stroke}}
\cspline{
1    0
-1   0
-0.5 0.8660
0.5 -0.8660
-0.5 -0.8660
0.5  0.8660
}
\endPSpicture}

\bigskip\bigskip\bigskip
\centerline{\PSunit=5mm
\PSpicture{4.5\PSunit}{4\PSunit}
   \setPSlinewidth{2mm}
   \PScline{0 0 4 0 4.5 4 0.5 4}
   \TX{\hss \sl Square\hss}{2.25 2}
\endPSpicture}


\bigskip

The square above was produced by the following code:

{\parskip=0pt
\verbatim|\centerline{\PSunit=5mm
\PSpicture{4.5\PSunit}{4\PSunit}
   \setPSlinewidth{2mm}
   \PScline{0 0 4 0 4.5 4 0.5 4}
   \TX{\hss \sl Square\hss}{2.25 2}
\endPSpicture}|
}

\head{Splines}

A possible choice of the curves $a_i(t)$ in the general spline
problem is to let each be a {\it Bezier cubic}. Assume given, in
addition to the two points $a_i$ and $a_{i+1}$, two control points
$z_i$ and $w_i$. Then, by definition, the Bezier cubic from $a_i$ to
$a_{i+1}$ determined by the two control points is the cubic curve
$a_i(t)$ for $0\le t\le 1$ given by the polynomial,
 $$a_i(t)=a_i(1-t)^3+3z_it(1-t)^2+3w_it^2(1-t)+a_{i+1}t^3.
 $$
 At the extremal points, the values are $a_i(0)=a_i$ and
$a_i(1)=a_{i+1}$, and the derivatives are $a_i'(0)=3(z_i-a_i)$ and
$a_i'(1)=3(-w_i+a_{i+1})$. Thus the two control points determine the
tangents at the extremal points. More precisely, the vector from
$a_i$ to $z_i$ is one third of the derivative $a_i'(0)$ and the
vector from $w_i$ to $a_{i+1}$ is one third of the derivative
$a_i'(1)$.

\medskip
\centerline{\PSunit=0,5cm
\PSpicture{5.1\PSunit}{5.1\PSunit}
\PS{newpath 0 0 moveto 2 5 4 4 5 2 curveto stroke}
\setPSlinewidth{0.1pt}
\PS{newpath 0 0 moveto 2 5 lineto stroke}
\PS{newpath 5 2 moveto 4 4 lineto stroke}
\TX{ $a_i$\hss}{0 0.05}
\TX{ $z_i$\hss}{2 5}
\TX{ $w_i$\hss}{4 4}
\TX{ $a_{i+1}$\hss}{5 2}
\PScircle{0 0.0 0.05}
\PScircle{2 5 0.05}
\PScircle{4 4 0.05}
\PScircle{5 2 0.05}
\endPSpicture
}


For cubic splines, the condition of smoothness is that the total
curve is two times differentiable.  For the Bezier cubic above, the
derivatives are given by the equations, 
 $$a_i'(t)/3=[z_i-a_i](1-t)^2+2[w_i-z_i]t(1-t)
+[a_{i+1}-w_i]t^2,
 $$
 $$a_i''(t)/6=[-2z_i+w_i+a_i](1-t)+[z_i-2w_i+a_{i+1}]t.
 $$
In particular, the condition of smoothness for the first order
derivatives, $a_i'(1)=a_{i+1}'(0)$ for $i=1,\dots,n-2$, is the
following set of $n-2$ equations,
 $$w_i+z_{i+1}=2a_{i+1}.
 $$
 It is natural to define $z_n:=2a_n-w_{n-1}$. Then the $n-1$ control
points $w_i$ are determined by the $n$ points $z_i$ through the
displayed equations above. So, the sequences of Bezier cubics
$a_i(t)$ forming a differentiable curve from $a_1$ to $a_n$ are
parameterized by the sets of $n$ points $z_i$.

In terms of the $z_i$, the second order derivatives at the extremal
points are given as follows, for $i=1,\dots,n-1$,
 $$a_i''(0)/6=-2z_i-z_{i+1}+a_i+2a_{i+1},\quad
a_i''(1)/6=z_i+2z_{i+1}-3a_{i+1}.
 $$
 Hence, the condition of smoothness for the second order derivatives,
$a_i''(1)=a_{i+1}''(0)$, is the following set of
equations,
 $$\llap{(*)\hskip1cm}
  z_i+4z_{i+1}+z_{i+2}=4a_{i+1}+2a_{i+2},\quad i=1,\dots,n-2.
 $$
 To determine the spline completely, we need two more independent
conditions on the $n$ control points $z_i$. One natural choice is to
require that the second order derivatives vanish at the extremal
points $a_1$ and $a_n$, that is, $a_1''(0)=0$ and $a_{n-1}''(1)=0$.
It imposes the following two conditions on the $z_i$:
 $$\llap{(**)\hskip1cm}
2z_1+z_2=a_1+2a_2, \qquad z_{n-1}+2z_n=3a_n.
 $$

The spline determined by the equations (*) and (**) is called the
{\it natural spline\/} through $a_1,\dots,a_n$. It should be
emphasized that there are other natural conditions to impose on a
spline. One possibility is to look for a spline with given tangents
at the two extremal points. In other words, assume that the first
control point $z_1$ and the last control point $w_{n-1}=2a_n-z_n$ are
given.  A second possibility is to determine the spline by the
condition that the second order derivatives at the extremal points
$a_1$ and $a_n$ are, respectively, equal to the second order
derivatives at the (nearby) points $a_2$ and $a_{n-2}$. The latter
condition replaces (**) by the following:
 $$\llap{(***)\hskip.7cm} 
z_1+z_2=(a_1+5a_2)/3,\quad z_{n-1}+z_n=(a_{n-1}+5a_n)/3
 $$
 Finally, the spline could be determined by giving the values of the
derivatives at the two end points. 
  
 A different problem is to look for a {\it closed spline}, that is,
to look in addition for a further cubic $a_n(t)$ connecting $a_n$ to
$a_1$. The further condition is then that the total closed curve is
smooth.  Reading the indices modulo $n$, the further condition simply
adds the two equations (*) for $i=n-1$ and $i=n$.

\smallskip 
There are several reasons for choosing (Bezier) cubics to
define the spline. The most important is simply that the PostScript
printer is able to draw a Bezier cubic given the coordinates of
the four points. This package implements the natural spline and the
closed spline.


\head{Solution of the equations}

The two sets of linear equations, for the natural spline and for the
closed spline, have a unique solution $(z_1,\dots,z_n)$. For the
natural spline, the matrix of coefficients is tridiagonal:
$$ \left[\ \vcenter{\halign{$\strut#$\hfill&&\quad\hfill$#$\hfill\cr
2&1\cr 1&4&1\cr
 &\ddots&\ddots&\ddots\cr
&&1&4&1\cr
&&&1&2\cr}}\ \right]
$$

The equation is solved by Gaussian elimination resulting in the
following algorithm (where the reciprocals of the diagonal
elements are stored in the factors $f_i$):

\def\cmt#1{\par\noindent\hbox to 2.2cm{\hfill\rm#1}}
\smallskip
\begingroup\parskip0pt\bf
\cmt{} input $n,\,a_1,\dots,a_n$
\cmt{} registers $i,\,z_1,\dots,z_n,\,f_1,\dots,f_n$
\cmt{(initialize)} $f_1\gets 1/2\ ,\ \ z_1\gets a_1+2{}a_2$ 
\cmt{(forward loop)} for $i=2,\dots,n-1$ do  
\cmt{}\hskip0.8cm $z_i\gets 4a_i+2a_{i+1}-z_{i-1}{}f_{i-1}$
\cmt{}\hskip0.8cm $f_i\gets 1/(4-f_{i-1})$
\cmt{} end
\cmt{(finish)} $z_n\gets 3a_n-z_{n-1}{}f_{n-1}\ ,\ \ %
       f_n\gets 1/(2-f_{n-1})$ 
\cmt{(solve)}  $z_n\gets z_n{}f_n$  
\cmt{(loop)} for $i=n-1,\dots,2,1$ do
\cmt{}\hskip0.8cm $z_i\gets (z_i-z_{i+1})f_i$
\cmt{} end
\cmt{} output $z_1,\dots,z_n$
\endgroup
\smallskip

In the loop, the factor $f_i$ converges to $2-\sqrt 3\approx 0.3$. So
the multiplications by $f_i$ and the inversions of $4-f_i$ cause no
problems of arithmetic overflow. 

\medskip

The matrix for the closed spline problem is slightly more complicated:

$$
\left[\ \vcenter{\halign{$\strut#$\hfill&&\quad\hfill$#$\hfill\cr
4&1&&&1\cr
1&4&1\cr
 &\ddots&\ddots&\ddots\cr
&&1&4&1\cr
1&&&1&4\cr}}\ \right]
$$

and so is the algorithm (the factor $f_n$ is inverted at the end of
the forward loop):
\smallskip
\begingroup\parskip0pt \bf
\cmt{} input $n,\,a_1,\dots,a_n$
\cmt{} registers $i,\,z_1,\dots,z_n,\,f_1,\dots,f_n,\,%
      q_1,\dots,q_{n}$
\cmt{(initialize)} $f_1\gets 1/2\ ,\ \ q_1\gets1\ ,\ \ %
    z_1\gets 4a_1+2{}a_2\ ,\ \ f_n\gets 4\ ,\ \ %
    z_n\gets 4a_n+2a_1$
\cmt{(forward loop)} for $i=2,\dots,n-1$ do  
\cmt{}\hskip0.8cm $z_i\gets 4a_i+2a_{i+1}-z_{i-1}{}f_{i-1}$
\cmt{}\hskip0.8cm $q_i\gets -q_{i-1}{}f_{i-1}\ ,\ \ 
         f_i\gets 1/(4-f_{i-1})$
\cmt{}\hskip0.8cm $f_n\gets f_n+q_i{}q_{i-1}\ ,\ \ %
         z_n\gets z_n+q_i{}z_{i-1}$
\cmt{} end
\cmt{(finish)} $q_n\gets -(1+q_{n-1})f_{n-1}$
\cmt{(loop)} $z_n\gets z_n+q_n{}z_{n-1}$
\cmt{} $f_n\gets f_n+q_n{}(1+q_{n-1})\ ,\ \ f_n\gets 1/f_n$
\cmt{(solve)}  $z_n\gets z_n{}f_n$  
\cmt{} for $i=n-1,\dots,2,1$ do
\cmt{}\hskip0.8cm $z_i\gets (z_i-z_{i+1}-q_i{}z_n)f_i$
\cmt{} end
\cmt{} output $z_1,\dots,z_n$
\endgroup


\head{Implementation}
 The two algorithms are implemented in the package. First, while
reading the arguments to \verbatim|\spline|, the arguments are stored
in dimension registers $A$ allocated dynamically. The factors $f$
appearing in the algorithms are stored in count registers $F$ as the
integers $F:=f*B$, where the base is $B=2^{15}$.  This choice of base
allows the square $B^2$ to be a valid \TeX{} integer, and it allows
the following procedure for the multiplications $z*f$ in the
algorithm: split $z$ into $z=zh*B+zl$ where the two (digits) $zh$ and
$zl$ are less than $B$. Then $zf=zh*F+zl*F/B$.  Similarly, the
inversion $f\gets(m-f)^{-1}$ (where $m=2$ or $m=4$) is obtained as
$F\gets B^2/(mB-F)$. The factors $q$ appearing in the algorithm for the
closed spline are treated similarly. 
 
The count registers $F$ and the dimension registers $Z$ are allocated
during the forward loop.  Finally, when the equations are solved, the
contents of the relevant registers are passed to the PostScript
interpreter in a \verbatim|\special|.

\bigskip
\bigskip\bigskip
\centerline{\PSunit=1truecm
\PSpicture{1.2\PSunit}{1.2\PSunit}
 \subpaths
 \cspline{0 0 1 0 1 0 1 1 1 1 0 1 0 1 0 0}
 \PS{gsave 0.8 setgray fill grestore stroke}
\endPSpicture
 %New picture
\hskip1truecm\PSunit=0.5truecm\PSpicture{2\PSunit}{2\PSunit}
 \llcoords{-1 -1}
%the black yin:
 \subpaths    %\longpath is not needed, since arcs are automatically
              %appended. 
 \PSarc{270 90}{0 0 1}
 \PSarc{90 270}{0 0.5 0.5}
 \PSarcn{90 270}{0 -0.5 0.5}
 \PS{fill}\normalpaths
%the black circle bounding yang:
 \PScircle{0 0 1}
%the black spot in yang:
 \def\stroke{\PS{fill}}             %fill fills with black
 \PScircle{0 -0.5 0.07}
%the white spot in yin:
 \def\stroke{\PS{1 setgray fill}}   %Now fill fills with white
 \PScircle{0 0.5 0.07}
\endPSpicture
%newpicture:
\hskip1truecm\PSunit=0.17truecm \PSpicture{10\PSunit}{6\PSunit}
 \llcoords{-5 -3}
 \subpaths
 \cspline{5 0 0 3 -5 0 0 -3}
 \longpath{2 0}
 \spline{2 0 0 0.5 -2 0}
 \spline{-2 0 0 -0.4 2 0}
 \PS{closepath gsave 0.8 setgray eofill grestore stroke}
 \normalpaths
 \PSline{2 0 2.5 0.25}
 \PSline{-2 0 -2.5 0.25}
\endPSpicture
%newpicture
\hskip1truecm \PSunit=0.3333truecm \PSpicture{3\PSunit}{3\PSunit}
 \llcoords{-1.5 -1.5}
 \def\stroke{\PS{fill}}
 \PScline{-1.5 -1.5 1.5 -1.5 1.5 1.5 -1.5 1.5}
 \def\stroke{\PS{1 setgray stroke}}
 \cspline{0 0 1 1 0 0 1 -1 0 0 -1 -1 0 0 -1 1}
\endPSpicture}

%      DEF of ANDERS THORUP signature.
%      Data were obtained by TeX-ing splsign.tex and then cutting the
%      resulting PostScript code out of the DVI-file. Then 
%      all decimal points and all trailing decimals were removed.
%      The result assumes \magnification=\magstep1 in source file.

\def\signature{Anders Thorup\special{" /c {curveto} def
%The overall scalefactor:
15000 65781.76 div dup scale
%Thickness:
0.8 0.2845 div setlinewidth
%Position:
-300 50 translate
%Anders:
newpath 40 49 moveto 38 46 37 44 40 44 c 42 43 47 44 50 49 c 52 53
53 62 50 64 c 46 65 37 61 30 54 c 22 46 14 37 11 30 c 7 22 6 16 8 13
c 9 9 11 6 17 10 c 22 13 31 22 37 29 c 42 35 44 38 47 41 c stroke
newpath 47 41 moveto 43 37 39 33 37 29 c 34 24 32 19 34 16 c 35 12 39
8 44 10 c 48 11 53 16 58 22 c stroke
newpath 58 22 moveto 56 17 54 13 55 10 c 55 6 58 4 62 6 c 65 7 70 12
73 15 c 75 17 76 16 76 15 c 75 13 75 11 76 9 c 76 6 79 5 81 5 c 82 4
83 6 88 10 c 92 13 101 18 107 20 c 112 21 115 19 118 18 c stroke
newpath 118 18 moveto 113 19 109 20 107 19 c 104 17 102 14 102 10 c
101 5 101 1 110 8 c 118 14 134 32 141 42 c 147 51 144 50 139 44 c 133
37 124 23 122 15 c 119 6 123 3 130 6 c 136 8 144 16 146 20 c 148 23
145 23 143 20 c 141 16 139 8 140 5 c 140 1 144 1 150 4 c 155 6 161 11
165 15 c 168 18 168 19 168 20 c stroke
newpath 168 20 moveto 168 19 168 19 168 17 c 167 14 164 8 166 7 c 167
5 171 7 175 10 c 178 12 181 15 185 18 c stroke
newpath 185 18 moveto 185 12 185 6 184 5 c 182 3 180 5 178 8 c stroke
%Thorup 
195 0 translate
newpath 6 15 moveto 6 13 6 12 6 10 c 5 7 4 4 5 4 c 5 3 7 4 11 10 c 14
15 20 24 24 30 c 27 35 29 35 32 36 c stroke
newpath 10 38 moveto 10 35 10 33 9 35 c 7 36 5 40 6 43 c 6 45 10 45
18 46 c 25 46 37 46 45 47 c 52 47 54 48 56 49 c stroke
newpath 55 48 moveto 54 48 53 48 49 42 c 44 35 34 21 30 14 c 25 6 25
5 30 8 c 34 10 41 17 44 18 c 46 18 43 13 43 10 c 42 6 43 3 50 6 c 56
8 68 15 74 18 c 79 20 79 19 79 19 c stroke
newpath 79 19 moveto 78 19 77 19 75 18 c 72 16 66 12 66 10 c 65 7 69
5 73 6 c 76 6 79 10 81 13 c 82 15 81 17 81 19 c stroke
newpath 81 19 moveto 83 17 85 16 90 17 c 94 17 99 20 103 21 c 106 21
107 19 106 17 c 104 14 100 12 100 10 c 99 7 101 5 105 6 c 108 6 112 9
117 12 c 121 14 125 17 129 20 c stroke
newpath 129 20 moveto 126 16 123 13 123 10 c 122 6 124 4 128 5 c 131
5 135 9 138 12 c 140 14 142 16 144 18 c stroke
newpath 144 18 moveto 142 14 140 10 140 8 c 139 5 141 4 146 6 c 150 7
157 12 165 17 c stroke
newpath 165 17 moveto 162 11 160 6 157 0 c 153 -6 150 -15 147 -20 c
143 -24 141 -24 144 -20 c 146 -15 153 -5 160 0 c 166 5 172 8 177 10 c
181 11 185 11 187 10 c 188 8 186 3 185 2 c 183 0 182 0 183 2 c 183 3
184 4 186 10 c 187 15 189 24 190 30 c 190 35 189 37 183 40 c 176 42
165 45 160 45 c 154 44 155 41 160 40 c 164 38 173 38 183 41 c 192 43
202 47 205 50 c 207 52 203 54 200 56 c stroke}}

\bigskip\bigskip\bigskip
\hbox to 14truecm{\hss\signature}



\bye