%To print this file, make the line with \endinput below a comment
%line (with a leading `%'). To use it as a sty-file, remove the `%'.
%
%Alternative: tex and print a file prt.tex with the following contents
%      \let\endinput=\relax
%      \input spline.sty
%      \bye
% 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    Anders Thorup
%                 PSPICTURE (for DVIPS)
%
%Syntax:
%
%    \PSpicture<xdimen><ydimen><picture material>\endPSpicture
%
%

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

\def\wl@g#1{\relax}
%FOR TESTS:
%\let\wl@g=\wlog

\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

%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}


\def\LTX#1{\TX{\hss#1}}
\def\CTX#1{\TX{\hss#1\hss}}
\def\RTX#1{\TX{#1\hss}}

%\rotateRTX is like \RTX with an extra rotation parameter determining the
%the amount of rotation. Thus
%   \rotateRTX{30}{text}{0 1}
%   will set the hbox containing the word text with its lower left
%   corner at (0,1), and rotate it 30 deg around the point.
%   
%   \rotateRTX{0}{tag}{0 1} is esentially equivalent to
%   \RTX{tag}{0 1}
%   except that the latter sets the hbox with tag centered vertically.
\def\rotateRTX#1#2#3{\begingroup\c@l=0 \st@ck#3 \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
    \special{ps: gsave currentpoint currentpoint translate #1
    rotate neg exch neg exch translate}#2\special{ps: grestore}%
    \hss\kern-\ll@x}\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, \cspline, \dspline. 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
   \wl@g{\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
  \wl@g{\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} %If \c@l=3, then \Acol=\A3
\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 three spline macros.  Structure
% def SPLINE#1
%             maybemarks         %Make \st@ck put marks at points if wanted
%            \c@l=0              %\c@l is number of arguments read so far
%            \st@ck #1           %Read arguments, set their number,
%                                %put them in dimen-registers Ai,
%                                %and (maybe) set marks on the points. 
%             checkarguments     %Even number and at least ....
%             Gauss              %Solve the equations: determine Z_i
%             PSwSPL@NE          %Write the Bezier cubics in special PS
%             PSwCfinish         %Finish with an extra cubic for \cspline
%             stroke             %stroke the spline
 
\def\spline@N#1{\begingroup
  \maybem@rks
  \c@l=0 \st@ck#1 \end
  \ch@ckeven\ch@ckthree
  \gauss@N
  \PSwSPL@NE \stroke
  \endgroup}
\def\spline@C#1{\begingroup
  \maybem@rks
  \c@l=0 \st@ck#1 \end
  \ch@ckeven\ch@ckthree
  \gauss@C   
  \PSwSPL@NE \PSw@Cfinish \stroke
  \endgroup}
\def\spline@D#1{\begingroup
  \maybem@rks
  \c@l=-2 \st@ck#1 \end
  \advance\numc@ls2 \ch@ckeven\ch@ckfour 
  \advance\numc@ls-4
         \ifnum\numc@ls>4 \let\next\gauss@D
                     \else\let\next\d@BEZIER\fi \next
  \PSwSPL@NE \stroke
  \endgroup}

%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\def\next{\expandafter\st@ckit}\fi\fi\next}%If arg=\cs
\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}

%The three spline algorithms are implemented in Gauss:
\def\gauss@N{\fwd@Ninit
      \loop\fwd@Nbody\ifnum\c@l<\numc@ls\repeat
      \fwd@Nfinish\solve@N}
\def\gauss@C{\fwd@Cinit
      \loop\fwd@Cbody\ifnum\c@l<\numc@ls\repeat
      \fwd@Cfinish\solve@C}
\def\gauss@D{\fwd@Dinit
      \loop\fwd@Nbody\ifnum\c@l<\numc@ls\repeat
      \fwd@Dfinish\solve@D}
\def\d@BEZIER{%Special case with only two points for a \dspline
  %Z1:=2A1-A{-1}, Z2:=2A2-A0 (Z_1)
  \c@l=1 \newdimenvar Z\currvar=2\Acol
    \advance\c@l-2 \advance\currvar -\Acol
  \c@l=2 \newdimenvar Z\currvar=2\Acol
    \advance\c@l-2 \advance\currvar -\Acol
  \fwd@Dfinish}

%The Bezier cubics are written to PostScript by \PSwSPL@NE:
\def\PSwSPL@NE{\PSw@init
  \loop\PSw@ZWA\ifnum\c@l<\numc@ls\repeat}
%Extra cubic for the closed spline
\def\PSw@Cfinish{\PSw@\Zcol
  \c@l=0 \PSw@Wcol\PSw@\Acol\PS{curveto}}


%\def\spline{\begingroup\spl@ne}
%\def\cspline{\begingroup\let\d@spline=\d@cspline
%   \let\PSw@finish=\PSw@Cfinish\spl@ne}
%\def\dspline{\begingroup\dspl@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\dspl@ne#1{\maybem@rks\c@l=0 \dst@ck#1 \end
%%AT\def\dspl@ne#1{\maybem@rks\c@l=0 \st@ck#1 \end
%\def\dspl@ne#1{\maybem@rks\c@l=-2 \st@ck#1 \end
%%AT   \ch@ckeven\ch@ckfour\advance\numc@ls-2 %
%   \ch@ckeven\ch@ckthree\advance\numc@ls-2 %
%%AT   \ifnum\numc@ls>6 %
%   \ifnum\numc@ls>4 %
%     \let\next\d@dspline 
%   \else
%     \let\next\d@dspecial
%   \fi\next\let\PSw@init=\PSw@dinit\d@PSw\PSw@finish\stroke\endgroup}

% We need special stack handling that doesn't
% put marks at control points for a directed spline
% NO, we will also mark the control points for a directed spline.
% So following code not needed.
%\def\dst@ck{\futurelet\nexttok\testdst@ck}
%\def\testdst@ck{\ifx\nexttok\sptoken
%   \let\next\eatspacedst@ck
%   \else\ifx\nexttok\end\let\next\dstackfinish
%   \else\def\next{\expandafter\dst@ckit}\fi\fi\next}
%\def\eatspacedst@ck#1 {\dst@ck}
%\def\dstackfinish#1{\numc@ls=\c@l}
%\def\dst@ckit#1 {\ifnum\c@l>3 \markoptions\fi %
% \advance\c@l1 \newdimenvar A\Acol=#1pt \dst@ck}


%Check arguments:
\def\ch@ckfour{\ifnum\numc@ls<8 %
      \errmessage{At least four points needed for directed spline}\fi}
\def\ch@ckthree{\ifnum\numc@ls<6 %
      \errmessage{At least three points needed for SPLINE}\fi}
\def\ch@cktwo{\ifnum\numc@ls<4 %
      \errmessage{At least two points needed}\fi}
\def\ch@ckeven{\ifodd\numc@ls
       \errmessage{Even number of coordinates needed}\fi}

%\def\d@spline{\fwd@Ninit\loop\fwd@Nbody\ifnum\c@l<\numc@ls\repeat
%  \fwd@Nfinish\solve@N}
%\def\d@cspline{\fwd@Cinit\loop\fwd@Cbody\ifnum\c@l<\numc@ls\repeat
%  \fwd@Cfinish\solve@C}
%\def\d@dspline{\fwd@Dinit%
%  \loop\fwd@Nbody\ifnum\c@l<\numc@ls\repeat\fwd@Dfinish\solve@D}
%\def\d@dspecial{%Special case with only two points for a \dspline
%  %Z3:=2A3-A1, Z4:=2A4-A2 (Z_1)
%%AT  \c@l=3 \newdimenvar Z\currvar=2\Acol
%  \c@l=1 \newdimenvar Z\currvar=2\Acol
%    \advance\c@l-2 \advance\currvar -\Acol
%%AT  \c@l=4 \newdimenvar Z\currvar=2\Acol
%  \c@l=2 \newdimenvar Z\currvar=2\Acol
%    \advance\c@l-2 \advance\currvar -\Acol
%  \fwd@Dfinish}
%


%After having read the arguments, \numc@ls=\c@l is equal to 2n.
%The linear equations are solved by Gauss elimination in Gauss,
%using the algorithms described in the manual. . 
%Macros for the individual parts of the algorithm:
%
%Initialization of the forward loops:
\def\fwd@Ninit{%        %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/4,
                   %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 }
\def\fwd@Dinit{%       %set F6=\prevF:=base/4
                       %    Z3:=2A3-A1, Z4:=2A4-A2 (Z_1)
                       %    Z5:=4A5+2A7-Z3, Z6:=4A6+2A8-Z4 (Z_2)
%AT  \c@l=3 \newdimenvar Z\currvar=2\Acol
  \c@l=1 \newdimenvar Z\currvar=2\Acol
    \advance\c@l-2 \advance\currvar -\Acol
%AT  \c@l=4 \newdimenvar Z\currvar=2\Acol
  \c@l=2 \newdimenvar Z\currvar=2\Acol
    \advance\c@l-2 \advance\currvar -\Acol
%AT  \c@l=5 \newdimenvar Z\currvar=4\Acol
  \c@l=3 \newdimenvar Z\currvar=4\Acol
    \advance\c@l2 \advance\currvar 2\Acol
    \advance\c@l-4 \advance\currvar -\Zcol
%AT  \c@l=6 \prevF=\base \divide\prevF4 \newcountvar F\Fcol=\prevF
  \c@l=4 \prevF=\base \divide\prevF4 \newcountvar F\Fcol=\prevF
    \newdimenvar Z\currvar=4\Acol
    \advance\c@l2 \advance\currvar 2\Acol
    \advance\c@l-4 \advance\currvar -\Zcol
%AT  \c@l=8 }
  \c@l=6 }

%The bodies of the forward loops (same bodies of spline@N and spline@D):
\def\fwd@Nbody{%      %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 }
%\fwd@Dbody = \fwd@Nbody

%Finishing the forward loops:
\def\fwd@Nfinish{%     %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}
\def\fwd@Dfinish{%   %Z{2n-1}:=A{2n+1}, Z{2n}:=A{2n+2}
  \c@l=\numc@ls
  \advance\c@l-1 \newdimenvar Z%
  \advance\c@l2 \currvar=\Acol
  \advance\c@l-1 \newdimenvar Z%
  \advance\c@l2 \currvar=\Acol}

%Now the backwards substitutions to solve the equations:
\def\solve@N{\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\solve@C{\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}
\def\solve@D{\c@l=\numc@ls \advance\c@l-1
  \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
%AT    \ifnum\c@l>5 \repeat}
    \ifnum\c@l>3 \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}
%AT\def\PSw@dinit{\c@l=4 \newpath\PSw@\Acol\moveto}
%\def\PSw@dinit{\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}
\newdimen\PSmarkradius \PSmarkradius=2\PSlinewidth
%Old version had: \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 fill 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{Two points 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}

\let\spline=\spline@N
\let\cspline=\spline@C
\let\dspline=\spline@D
\catcode`\@=\atcode

%END OF spline.sty
\endinput

\def\uncatcodespecials{\def\do##1{\catcode`##1=12 }%
     \dospecials\do\@} %AT added \do @ 7/2-02 for print with AMSTeX
\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$\ \ }}
\def\Indent{\par\noindent\hskip1.33truecm}

\centerline{ANDERS THORUP: PSpictures for DVIPS, 10 September 2004}


\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 in a \TeX{}
manuscript a spline with commands like the following:

\noindent \qquad
 \verbatim|\spline{|$\,a_1\ \dots\ a_n\,$\verbatim|}|, \quad
\verbatim|\cspline{|$\,a_1\ \dots\ a_n\,$\verbatim|}|, \quad
\verbatim|\dspline{|$\,a_0\ a_1\ \dots\ a_{n+1}\,$\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 first command should be a smooth
curve through the points on the paper determined by the $n$ pairs of
coordinates. The second command should result in a closed spline, and
the last in a directed spline: the points $a_0$ and $a_{n+1}$ are
control points determining the tangents at $a_1$ and $a_n$. 

\centerline{\PSunit=0.7truecm\PSpicture{20\PSunit}{4.5\PSunit}
\putPSmarks
\llcoords{-1 -1}
\spline{0 0 3 1 2.5 2 0 2} \RTX{\it spline}{0 3}
\CTX{$_{a_n}$}{0.1 -0.25} %Let tex choose the subscript font!
\LTX{$_{a_1}$}{-0.1 2}
\RTX{$_{a_2}$}{2.5 2.2}
\llcoords{-9 -1}
\cspline{0 0 3 1 2.5 2 0 2} \RTX{\it cspline}{0 3}
\CTX{$_{a_n}$}{0.1 -0.25}\LTX{$_{a_1}$}{-0.1 2}
\RTX{$_{a_2}$}{2.5 2.2}
\llcoords{-17 -1}
\dspline{-1 -1 0 0 3 1 2.5 2 0 2 -1 3} \RTX{\it dspline}{0 3}
\CTX{$_{a_n}$}{0.1 -0.25}\LTX{$_{a_1}$}{-0.1 2}
\RTX{$_{a_2}$}{2.5 2.2}
\RTX{$_{a_{n+1}}$}{-0.85 -1.1}\CTX{$_{a_{0}}$}{-1 2.7}
\endPSpicture}



\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 macros for curve
drawing. To use the package as an input file in other \TeX{}
documents, find the following line in the file:

\Indent\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.

\vfill\eject

\head{The Macros}


\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
\Indent\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
\Indent\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|,  \verbatim|\cspline|, and
\verbatim|\dspline|: 

 \smallskip
\Indent\verbatim|\spline{|{\it point-list}\verbatim|}|,  \quad
\verbatim|\cspline{|{\it point-list}\verbatim|}|, \quad
\verbatim|\dspline{|{\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 points $a_1,\dots,a_n$ given by the $2n$ coordinates in
the list; \verbatim|\cspline| draws a closed spline. Finally,
\verbatim|\dspline| draws a directed spline, see the section on
splines. 

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 with \verbatim|\LTX|,
\verbatim|\RTX|, \verbatim|\CTX|:
 \smallskip
\Indent\verbatim|\LTX{|{\it TeX-stuff}%
 \verbatim|}{|{\it xcoord ycoord}\verbatim|}|
 \smallskip
The {\it TeX-stuff} is any horizontal \TeX\ material. 
\verbatim|\LTX| will set the {\it TeX-stuff\/} in an hbox left to the
point with coordinates ({\it xcoord,ycoord\/}), and centered
vertically.  Similarly, 

\Indent\verbatim|\RTX{|{\it TeX-stuff}%
 \verbatim|}{|{\it xcoord ycoord}\verbatim|}|, and 
 \verbatim|\CTX{|{\it TeX-stuff}%
 \verbatim|}{|{\it xcoord ycoord}\verbatim|}|
\smallskip

 will, respectively, set the {\it TeX-stuff} to the right of the point
 and centered 
horizontally around the point.  Occasionally, the command
\verbatim|\TX{|{\it TeX-stuff-of-width-zero}%
 \verbatim|}{|{\it xcoord ycoord}\verbatim|}| is useful.
In addition, the {\it TeX-stuff} may be rotated with
\verbatim|\rotateRTX|:

\smallskip
\Indent\verbatim|\rotateRTX{|{\it angle}\verbatim|}{|{\it TeX-stuff}\verbatim|}{|{\it xcoord ycoord}% 
       \verbatim|}|
\smallskip
 rotates the {TeX-stuff}, placed in an hbox with its lower left point at
the given coordinates, clock-wise the specified angle around the
point.
 
\bull
 Circles may be drawn with \verbatim|\PScircle|:
 \smallskip
\Indent\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
\Indent\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
\Indent\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
\Indent\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 latter point has to be different from
the first point in the list, 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
\Indent\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
\Indent\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
\Indent\verbatim|\llcoords{|{\it xcoord ycoord}\verbatim|}|
\smallskip

 Thus, after \verbatim|\llcoords{2 -0.5}|, the lower left point of the
picture box is the point with coordinates $(2,-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 the various picture commands have a thickness
equal to the value of \verbatim|\PSlinewidth|. By default,
 \smallskip
\Indent\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
\Indent\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
\Indent
 \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 line width. The default is
reestablished with
 \smallskip
\Indent \verbatim|\nodashes|
 \smallskip

\bull
 The size of the arrowhead produced with \verbatim|\PSarrow| can be
changed:
 \smallskip
\Indent
 \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
\Indent
 \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
filled circle of radius equal to \verbatim|\PSmarkradius|. The
value of this dimension register is called whenever a PSmark is set,
so you may change the value of \verbatim|\PSmarkradius| at any time.
The default value is set as follows:
 \smallskip
\Indent
 \verbatim|\PSmarkradius=2\PSlinewidth|,
 \smallskip
 but if you want this coupling between the size of the marks and the
width of the curves, you have to repeat the above setting after each use of
\verbatim|\setPSlinewidth|.


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

 \Indent
 \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 of the curve macros (\verbatim|\spline|,
\verbatim|\PSline|, \verbatim|\PSarc|, etc) described in this package
results in PostScript code defining a path in the PostScript memory.
The beginning of the code is the expansion of the control sequence
\verbatim|\newpath|, and the end of the code is the expansion of
\verbatim|\stroke|. The default effects of the expansions of these
two control sequences are,
 %\verbatim|\newpath| and \verbatim|\stroke| are,
 respectively, to start a new path and to stroke the current path. In
fact,
 \smallskip
\Indent\verbatim|\def\newpath{\PS{newpath}}|\qquad
   \verbatim|\def\stroke{\PS{stroke}}|
 \smallskip
 However, with some knowledge of the PostScript language, you may
redefine the two control sequences to obtain different effects.  For
instance, after

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

each of the following curves generated by the curve macros will have
their enclosed area filled with gray before the curve is stroked.
Also, with temporary settings like \verbatim|\let\newpath=|
\verbatim|\relax|,
\verbatim|\let\stroke=\relax|, you may combine several curves into a
single path object. Three commands deal directly with combinations of
paths:
 \smallskip
\Indent\verbatim|\subpaths|\qquad
  \verbatim|\longpath{|{\it xcoord ycoord}\verbatim|}|\qquad
  \verbatim|\shortpaths|
 \smallskip
 After the first macro, \verbatim|\subpaths|, each of the following
curve macros will add a subpath to the current path. These curve macros
are typically mixed with occurrences of the two other macros:
After \verbatim|\shortpaths|, each
curve macro will contribute with a single subpath, and after
\verbatim|\longpath| each curve macro will append to the subpath
started at the point ({\it xcoord, ycoord}), with the beginning of
the curve connected with a straight line to the end of the previous
curve.  (Note that the arc defining commands need a little special
care: even after the command \verbatim|\shortpaths|, you have to
specify where the arc begins: use
  \verbatim|\PS{0 2 moveto}\PSarc{90 180}{0 0 2}}| to make sure that
the arc starts at the point $(0,2)$ (as intended).  Also, as the
PSmarks are drawn with small circles, don't use them after
\verbatim|\shortpaths|.)

When the path has been created, you have to instruct PostScript
explicitly on what to do with the path, for instance to fill it with
\verbatim|\PS{fill}| (or \verbatim|\PS{eofill}|), or to stroke it
with \verbatim|\PS{stroke}|.  In addition, you may want to restore
the defaults with

\Indent\verbatim|\normalpaths|
 

For instance, in the following code 

\Indent\verbatim|\def\stroke{\PS{gsave 1 setgray fill grestore}}|
\Indent\verbatim|\PScircle{0 0 1}\PScline{0 0 1 0 1 2 0 2}\normalpaths|
\smallskip
the \verbatim|\PScircle| produces a white disk and the
\verbatim|\PScline| a white rectangle.  Figures like this may be
used to cover parts of previously set black material. 


\bull
 The picture commands are meant to be used inside the picture environment.
 If used outside, most of them will provoke mysterious syntax errors,
typically of the form

\Indent\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 6.2830 0}
\PSline{0 -1 0 1}
\TX{ $x=\cos t$\hss}{4.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. And with a
directed spline, giving the horizontal tangents at the end points via
two extra control points, fewer points are needed:

 \medskip
\centerline{\PSunit=1truecm
\PSpicture{15.2\PSunit}{2.4\PSunit}
\llcoords{0 -1.2}
\PSline{0 0 6.2830 0}
\PSline{0 -1 0 1}
\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
}
\TX{ $x=\cos t$\hss}{4.5 -0.5}
\llcoords{-8.5 -1.2}\nomarks
\PSline{0 0 6.2830 0}
\PSline{0 -1 0 1}
\putTeXmarks
\dspline{
-0.45 1
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
6.7330   1
}
\TX{ $x=\cos t$\hss}{4.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}|
}

\bull
A dirty trick: The drawing

\centerline{\PSunit=3truemm
\PSpicture{10\PSunit}{5\PSunit}
%
%newpath 0.0 1.0 moveto
% 1.26195 1.006 2.52391 1.012 3.0 2.0 curveto
% 3.47609 2.988 3.16669 4.95839 4.0 5.0 curveto
% 4.83331 5.04161 6.80963 3.15477 8.0 2.0 curveto
% 9.19037 0.84523 9.59546 0.42259 10.0 0.0 curveto
% stroke
%newpath 10.0 3.0 moveto
% 9.64455 2.60002 9.28911 2.20006 8.0 2.0 curveto
% 6.71089 1.79994 4.48886 1.8 3.0 2.0 curveto
% 1.51114 2.2 0.75554 2.6001 0.0 3.0 curveto
% stroke
\PS{newpath 3.0 2.0 moveto
 3.47609 2.988 3.16669 4.95839 4.0 5.0 curveto
 4.83331 5.04161 6.80963 3.15477 8.0 2.0 curveto
 6.71089 1.79994 4.48886 1.8 3.0 2.0 curveto
 gsave 0.8 setgray fill grestore}
\spline{0 1 3 2 4 5 8 2 10 0}\spline{10 3 8 2 3 2 0 3}
\endPSpicture}

was made as follows: First, the two curves of the drawing were produced
by two splines:

\smallskip \hskip1.5truecm
 \verbatim|\spline{0 1 3 2 4 5 8 2 10 0} \spline{10 3 8 2 3 2 0 3}|
\smallskip

 Note that the two points of intersection, $(3,2)$ and $(8,2)$,
appear in both point lists. When compiled, the two splines resulted
in PostScript code, easily found in a special command in the dvi file:

{\parskip=0pt\parindent=1truecm
\noindent\verbatim|
newpath 0.0 1.0 moveto
 1.26195 1.006 2.52391 1.012 3.0 2.0 curveto
 3.47609 2.988 3.16669 4.95839 4.0 5.0 curveto
 4.83331 5.04161 6.80963 3.15477 8.0 2.0 curveto
 9.19037 0.84523 9.59546 0.42259 10.0 0.0 curveto
 stroke
newpath 10.0 3.0 moveto
 9.64455 2.60002 9.28911 2.20006 8.0 2.0 curveto
 6.71089 1.79994 4.48886 1.8 3.0 2.0 curveto
 1.51114 2.2 0.75554 2.6001 0.0 3.0 curveto
 stroke|
}

It is easy to identify from this code the $4+3$ Bezier cubics forming
the two curves.  In particular, it was easy to extract the PostScript
code of the 3 Bezier cubics forming the boundary of the shaded region.
In turn, this extract was inserted with shading commands directly in
the tex file, before the two spline commands, as follows:

{\parskip=0pt\parindent=1truecm
\noindent\verbatim|
\PS{newpath 3.0 2.0 moveto
 3.47609 2.988 3.16669 4.95839 4.0 5.0 curveto
 4.83331 5.04161 6.80963 3.15477 8.0 2.0 curveto
 6.71089 1.79994 4.48886 1.8 3.0 2.0 curveto
 gsave 0.8 setgray fill grestore}
|
}
%\vfill\eject
\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
}

\def\Imath#1$#2${\smallskip
   \Indent\hbox to 1truecm{\hss#1\qquad}$#2$
   \smallskip\noindent\ignorespaces }


For cubic splines, the condition of smoothness is that the total
curve is two times continuously differentiable.  For the Bezier cubic
above, the derivatives are given by the equations,
\Imath$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,$
  \vskip-\medskipamount
 \Imath$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 continuity 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,
 \Imath$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
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$,
 \Imath$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 continuity for the second order derivatives,
$a_i''(1)=a_{i+1}''(0)$, is the following set of
equations,

\Imath{(*)}
 $z_i+4z_{i+1}+z_{i+2}=4a_{i+1}+2a_{i+2},\quad i=1,\dots,n-2.$
 % \smallskip\noindent
% $$\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$:

\Imath{(**)}
 $2z_1+z_2=a_1+2a_2, \qquad z_{n-1}+2z_n=3a_n.
 $

\noindent
 The spline through $a_1,\dots,a_n$ determined by the equations (*)
and (**) is called the {\it natural spline}. It should be emphasized
that there are other natural conditions to impose on a spline. For
instance, the spline is determined by giving the tangents at the two
extremal points, or, equivalently, by giving the first control point
$z_1$ and the last control point $w_{n-1}=2a_n-z_n$.  We will call
this spline the {\it directed spline}. In the package, the two
control points are given by adding two extra points, $a_0$ and
$a_{n+1}$, to the list of points $a_i$, with the interpretation that
$z_1-a_1=a_1-a_0$ and $z_n=a_{n+1}$; so we add the following
equations: 

 \Imath{(***)}
 $z_1=2a_1-a_0, \quad  z_{n}=a_{n+1}.$

  A third 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 equations:

\Imath
 $z_1+z_2=(a_1+5a_2)/3,\quad z_{n-1}+z_n=(a_{n-1}+5a_n)/3.
 $
   This third possibility is not implemented in the package.

 A different problem is to look for a {\it closed spline}, that is,
to look for an additional 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 set of linear equations for the control points $z_i$, for the
natural spline, for the directed spline, and for the closed spline,
have a unique solution.  For the natural spline, the equations are
(*) and (**), and 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 equations (*) and (**) are 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.45cm{\hfill\rm#1}}
\smallskip
\begingroup\parskip0pt\bf
\cmt{SPLINE:\hskip7.5truemm} 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

For the closed spline, the equations (*) modulo $n$ give a
more complicated matrix: 
 $$ \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 a more complicated algorithm (the factor $f_n$ is inverted at the
end of the forward loop):
\smallskip
\begingroup\parskip0pt \bf
\cmt{CSPLINE:\hskip7.5truemm} 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}, \  
 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{(loop)} 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

 \medskip
 Finally, for the directed spline, 
 the equations (*) and (***) correspond to the following
matrix:
 $$
\left[\ \vcenter{\halign{$\strut#$\hfill&&\quad\hfill$#$\hfill\cr
1&&&&\cr
1&4&1\cr
 &\ddots&\ddots&\ddots\cr
&&1&4&1\cr
&&&&1\cr}}\ \right]
$$
 The equations (*) and (***) are solved by the following algorithm:

\def\cmt#1{\par\noindent\hbox to 2.2cm{\hfill\rm#1}}
 \smallskip
\begingroup\parskip0pt\bf
\cmt{DSPLINE:\hskip6.5truemm} input $n,\,a_0,a_1,\dots,a_n, a_{n+1}$
\cmt{} registers $i,\,z_1,\dots,z_n,\,f_2,\dots,f_{n-1}$
\cmt{(initialize)} $f_2\gets 1/4$, \ $z_1\gets 2a_1-a_0$, \ %
     $z_2\gets 4a_2+2a_3-z_1$
\cmt{(forward loop)} for $i=3,\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 a_{n+1}$
\cmt{(loop)} for $i=n-1,\dots,3,2$ 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




\head{Implementation}
 The three 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 other
algorithms 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
 The main part of the \TeX-code was written in 1995. The code
implementing the DSPLINE algorithm was written by Ulrik Buchholtz in
September 2003.  The code for determining the PSmarkradius was
changed as described above in
October 2003. Compared to the previous version, the default circles produced
by \verbatim|\PSmark| and \verbatim|\putPSmarks| have doubled their
sizes. 


\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