Additions to tutorial

[Ilya Zakharevich <ilya@math.ohio-state.edu>]

\section{Plotting}

Plotting functions use a multitude of flags, which are mostly power-of-2. To
simplify understanding give these flags symbolic names.  The following
code also decreases the default precision, so that the calculations
are performed much quickier (this is very important, since plotting
involves calculation of functions at many points).

\bprog
parametric=1; recursive = 2; norescale=4;
no\_x\_axis=8; no\_y\_axis=16; noframe=32;
points=64; points\_lines = 128; splines = 256; 
nw=0; sw=2; se=4; ne=6; left = 0; center = 1; right = 2; hgap = 16;
bottom = 0; vcenter = 4; top = 8; vgap = 32; relative=1;
\b{p} 9
\eprog

Start with something really simple:

\kbd{ploth(X=-2,2,sin(X\pow 7))}

You can see the limitations of the ``straightforward'' mode of plotting:
while the first several cycles of \kbd{sin} reach $-1$ and $1$, the cycles
which are closer to the left and right border do not. This is understandable,
since PARI is calculating $\sin X^7$ at many points, but these points
have no direct relationship to the ``interesting'' point on the graph
of this function.  No value close enough to the maxima and minima are
calculated, which leads to wrong turning points of the graph.

There is a way to fix this: one can ask PARI to use variable step which
smaller at the points where the graph of the function is more curved:

\kbd{ploth(X=-2,2,sin(X\pow 7),recursive)}

\noindent The precision near the edges of the graph is much better now.
However, the recursive plotting (named so since PARI subdivides intervals
until the graph becomes almost straight) can its own pitfalls.  Try

\kbd{ploth(X=-2,2,sin(X*7),recursive)}

\noindent Note that the graph looks correct far away, but it has a straight
interval near the origin, and some sharp corners as well.  This happens
because the graph is symmetric with respect to the origin, thus the middle 3
points calculated during the initial subdivision of $[-2,2]$ are exactly on
the same line.  To PARI this indicates that no further subdivision is needed,
and it plots the graph on this subinterval as a straight line.

There are many ways to circumvent this.  Say, one can make the right limit
2.1.  Or one can ask PARI for an initial subdivision into 16 points instead
of default 15:

\kbd{ploth(X=-2,2,sin(X*7),recursive,16)}

All these arrangements break the symmetry of the initial subdivision, thus
make the problem go away.  Eventually PARI will be able to better detect such
pathological cases, but currently some manual intervention may be required.

Function \kbd{ploth} has some additional enhancements which allow graphing
in situations when the calculation of the function takes a lot of time.  Let
us plot $\zeta(0.5+it)$:

\kbd{ploth(t=100,110,real(zeta(0.5+I*t)),,1000)}

\noindent This can take quite some time.  (1000 is close to the default for
many plotting devices, we want to specify it explicitely so that the result
do not depend on the output device.)  Try the recursive plot:

\kbd{ploth(t=100,110,real(zeta(0.5+I*t)),recursive)}

It takes approximately the same time.  Now try specifying fewer points,
but make PARI approximate the data by a smooth curve:

\kbd{ploth(t=100,110,real(zeta(0.5+I*t)),splines,100)}

This takes much less time, and the output is practically the same.  How to
compare these two outputs?  We will see it shortly.  Right now let us
plot both real and complex parts of $\zeta$ on the same graph:

\bprog
f(t) = z=zeta(0.5+I*t);[real(z),imag(z)]
ploth(t=100,110,f(t),,1000)
\eprog

Note how one half of the roots of the real and imaginary parts coincide.  Why
did we define a function \kbd{f(t)}?  To avoid calculation of $\zeta(0.5+it)$
twice for the same value of t.  Similarly, we can plot parametric graph:

\kbd{ploth(t=100,110,f(t),parametric,1000)}

Again, one can speed up the calculation with

\kbd{ploth(t=100,110,f(t),parametric+splines,100)}

If your plotting device supports it, you may ask PARI to show the points
in which it calculated your function:

\kbd{ploth(t=100,110,f(t),parametric+splines+points\_lines,100)}

As you can see, the points are very dense on the graph.  To see some crude
graph, one can even decrease the number of points to 30.  However, if you
decrease the number of points to 20, you can see that the approximation to
the graph now misses zero.  Using splines, one can create reasonable graphs
for larger values of t, say with

\kbd{ploth(t=10000,10004,f(t),parametric+splines+points\_lines,50)}

How can we compare two graphs of the same function plotted by different
methods?  Documentation shows that \kbd{ploth} does not provide any direct
method to do so.  However, it is possible, and even not very complicated.

The solution comes from the other direction.  PARI has a power mix of high
level plotting function with low level plotting functions, and these functions
can be combined together to obtain many different effects.  Return back
to the graph of $\sin X^7$.  Suppose we want to create an additional 
rectangular frame around our graph.  No problem!

First create a ``rectangle'' (this is a PARI term for a target of low-level
drawing primitives), name it rectangle 2 (any number between 0 andn 15 will
go):
\kbd{plotinit(2)}, \kbd{plotscale(2, 0, 1, 0, 1)}.  This makes a rectangle
which the size equal to the size of the output device, and makes the coordinate
system inside this rectangle go from 0 to 1 (for both $x$ and $y$).
Create a rectangular ``frame'' along the boundary of this rectangle:
\kbd{plotmove(2,0,0); plotbox(2,1,1)}.  Suppose we want to draw the graph
inside a subrectangle of this with the gaps on the bottom and on the left
of 0.1 of the size of the larger rectangle, and gaps 0.02 on the right
and on the top.

Create another rectangle (call it 3) of linear size 0.88 of the size of the
initial rectangle: \kbd{plotinit(3,0.88,0.88,relative)}.  Now graph the
function inside rectangle 3: \kbd{plotrecth(3,X=-2,2,sin(X\pow 7),recursive)}
and output rectangles 2 and 3 on the same plot:

\kbd{plotdraw([2,0,0,3,0.1,0.02],relative)}

\noindent (Funny offsets of the
rectangle 3 appear since the offsets are counted from upper-left corner.)
The output misses something comparing to the output of \kbd{ploth}: there
are no coordinates of the corners of the internal rectangle.  If your output
device supports mouse operations, you can find coordinates of particular
points of the graph, but it is nice to have something printed on
a hardcopy too.

However, it is easy to put $x$- and $y$-limits on the graph.  In the coordinate
system of the rectangle 2 the corners are $(0.1,0.1)$, $(0.1,0.98)$,
$(0.98,0.1)$, $(0.98,0.98)$.  We can mark lower $x$-limit
by doing \kbd{plotmove(2,0.1,0.1); plotstring(2,"-2.000",left+top+vgap)}.
$y$-coordinate is a little bit trickier, since in principle we do not know
it in advance (though for $\sin X^7$ it is very easy to guess).  Fortunately,
\kbd{plotrecth} returns the $x$- and $y$-limits.

Here is the complete program (rounding could have been done better,
is it possible in PARI?):

\bprog
plotinit(3,0.88,0.88,relative)
lims=plotrecth(3,X=-2,2,sin(X\pow 7),recursive)
limits=vector(4,i,round(1e4*lims[i])/1e4)
plotinit(2);          plotscale(2, 0, 1, 0, 1)
plotmove(2,0,0);      plotbox(2,1,1)
plotmove(2,0.1,0.1);  plotstring(2,limits[1],left+top+vgap)
plotstring(2,limits[3],bottom+vgap+right+hgap)
plotmove(2,0.98,0.1); plotstring(2,limits[2],right+top+vgap)
plotmove(2,0.1,0.98); plotstring(2,limits[4],right+hgap+top)
plotdraw([2,0,0,3,0.1,0.02],relative)
\eprog

We started with a trivial requirement: have an additional frame around
the graph, and it took some effort to do so.  But at least it was possible,
and PARI did the hardest part: creating the actual graph.
Now do a different thing: plot together the ``exact'' graph of $\zeta(0.5+it)$
together with one obtained from splines approximation.  We can emit
these graphs into two rectangles, say 0 and 1, then put these two
rectangles together on one plot.  Or we can emit these graphs into one
rectangle 0.

However, a problem arises: note how we
introduced a coordinate system in rectangle 2 of the above example, but we
did not introduce a coordinate system in rectangle 3.  Plotting a
graph into rectangle 3 automatically created a coordinate system
inside this rectangle (you could see this coordinate system in action
if your output device supports mouse operations).  If we use two different
methods of graphing, the bounding boxes of the graphs will not be exactly
the same, thus outputting the rectangles may be tricky.  Thus during
the second plotting we ask \kbd{plotrecth} to use the coordinate system of
the first plotting.  Let us add another plotting with fewer
points too:

\bprog
plotinit(0,0.9,0.9,relative)
plotrecth(0,t=100,110,f(t),parametric,300)
plotrecth(0,t=100,110,f(t),parametric+splines+points\_lines+norescale,30);
plotrecth(0,t=100,110,f(t),parametric+splines+points\_lines+norescale,20);
plotdraw([0,0.05,0.05],relative)
\eprog

This achieves what we wanted: we may compare different ways to plot a graph,
but the picture is confusing: which graph is what, and why there are multiple
boxes around the graph?  At least with some output devices one can control
how the output curves look like, so we can use this to distinguish different
graphs.  And the mystery of multiple boxes is also not that hard to solve:
they are bounding boxes for calculated points on each graph.  We can disable
output of bounding boxes with appropriate options for \kbd{plotrect}.
With these frills the script becomes:

\bprog
plotinit(0,0.9,0.9,relative)
plotpointtype(-1,0)                     \q\q\bs\bs set color of graph points
plotpointsize(0,0.4)                    \q\q\bs\bs use tiny markers (if available)
plotrecth(0,t=100,110,f(t),parametric+points,300)
plotpointsize(0,1)                      \q\q\bs\bs normal-size markers
plotlinetype(-1,1)                      \q\q\bs\bs set color of graph lines
plotpointtype(-1,1)                     \q\q\bs\bs set color of graph points
curve\_only = norescale + noframe + no\_x\_axis + no\_y\_axis
plotrecth(0,t=100,110,f(t),parametric+splines+points\_lines+curve\_only,30);
plotlinetype(-1,2)                      \q\q\bs\bs set color of graph lines
plotpointtype(-1,2)                     \q\q\bs\bs set color of graph points
plotrecth(0,t=100,110,f(t),parametric+splines+points\_lines+curve\_only,20);
plotdraw([0,0.05,0.05],relative)
\eprog

\noindent Plotting axes on the second and third graph would not hurt, but 
is not needed too, this is why we omit axes.  One can see that the
discrepancies between the exact graph and one based on 30 points exist,
but are pretty small.  Decreasing the number of points to 20 creates quite
noticable discrepancies.

Keep in mind that \kbd{plotlinetype}, \kbd{plotpointtype},
\kbd{plotpointsize} may do nothing on some terminals.  

What if we
want to create a high-resolution hardcopy of the plot?  There may be several
possible solutions.  First, the display output device may allow a
high-resolution hardcopy itself.  Say, PM display (with gnuplot output on
OS/2) pretends that its resolution is $19500\times 12500$, thus the data
PARI sends to it are already high-resolution, and printing is available
through the menubar.  Alternatively, with gnuplot output one can switch
the output plotting device to many different hardcopy devices:
\kbd{plotfile("plot.tex")}, \kbd{plotterm("texdraw")}.
After this all the plotting will go into file {\tt plot.tex} with whatever
output conventions gnuplot format {\tt texdraw} provides.  To switch output
back to normal, one needs to restore the initial plotting terminal, and
restore the initial output file by doing \kbd{plotfile("-")}.

One can combine PARI programming capabilities to produce multiple plots:

\bprog
plotfile("manypl1.gif")       \q\q\bs\bs Avoid switching STDOUT to binary mode
plotterm("gif=300,200")
wpoints = plothsizes()[1]       \q\q\bs\bs 300 x 200 is advice only
$\{$  for( k=1,6,
\q\q	plotfile("manypl"k".gif");
\q\q	ploth(x=-1,3,sin(x\pow k),0,wpoints)
\q   )
$\}$
\eprog

\noindent This plots 6 graphs of $\sin x^k$, $k=1\dots 6$ into 
$300\times 200$ GIF files {\tt manypl1.gif}\dots {\tt manypl6.gif}.

Additionally, one can ask PARI to output a plot into a PS file: just
use the command \kbd{psdraw} instead of \kbd{plotdraw} in the above examples
(or \kbd{psploth} instead of \kbd{ploth}).  See \kbd{psfile} argument
to \kbd{default} for how to change the output file for this operation.  Keep
in mind that the precision of PARI PS output will be the same as for the
current output device.

Now suppose we want to join many different small graphs into one picture.
We cannot use one rectangle for all the output as we did in the example
with $\zeta(0.5+it)$, since the graphs should go into different places.
Rectangles are a scarce commodity in PARI, there are only
16 user-accessible rectangles.  Does it mean that we cannot have more than
16 graphs on one picture?  Thanks to an additional operation of PARI
plotting engine, there is no such restrictions.  This operation is 
\kbd{plotcopy}.

The following script puts 4 different graphs on one plot using 2 rectangles
only, \kbd{2} and \kbd{3}:

\bprog
plotinit(2);         plotscale(2, 0, 1, 0, 1)
\vskip 0.15truecm
plotinit(3, 0.42, 0.42, relative);
plotrecth(3, x= -5, 5, sin(x), recursive)
plotcopy(3, 2, 0.05, 0.05, relative + nw)
\vskip 0.15truecm
plotmove(2, 0.05 + 0.42/2, 1 - 0.05/2)
plotstring(2,"Graph", center + vcenter)
\vskip 0.15truecm
plotinit(3, 0.42, 0.42, relative);
plotrecth(3, x= -5, 5, [sin(x),cos(2*x)], 0)
plotcopy(3, 2, 0.05, 0.05, relative + ne)
\vskip 0.15truecm
plotmove(2, 1 - 0.05 - 0.42/2, 1 - 0.05/2)
plotstring(2,"Multigraph", center + vcenter)
\vskip 0.15truecm
plotinit(3, 0.42, 0.42, relative);
plotrecth(3, x= -5, 5, [sin(3*x), cos(2*x)], parametric)
plotcopy(3, 2, 0.05, 0.05, relative + sw)
\vskip 0.15truecm
plotmove(2, 0.05 + 0.42/2, 0.5)
plotstring(2,"Parametric", center + vcenter)
\vskip 0.15truecm
plotinit(3, 0.42, 0.42, relative);
plotrecth(3, x= -5, 5, [sin(x), cos(x), sin(3*x),cos(2*x)], parametric)
plotcopy(3, 2, 0.05, 0.05, relative + se)
\vskip 0.15truecm
plotmove(2, 1 - 0.05 - 0.42/2, 0.5)
plotstring(2,"Multiparametric", center + vcenter)
\vskip 0.15truecm
plotlinetype(2,3)
plotmove(2, 0, 0)
plotbox(2, 1, 1)
\vskip 0.15truecm
plotdraw([2, 0, 0])
\bs\bs psdraw([2, 0, 0], relative)          \bs\bs if hardcopy is needed
\eprog

The rectangle \kbd{2} plays the role of accumulator, rectangle \kbd{3} is
used as a target to \kbd{plotrecth} only.  Immediately after plotting into 
rectangle \kbd{3} the contents is copied to accumulator.  Let us explain
numbers which appear in this example: we want to create 4 internal rectangles
with a gap 0.06 between them, and the outside gap 0.05 (of the size of the
plot).  This leads to the size 0.42 for each rectangle.  We also
put captions above each graph, centered in the middle of each gap.  There
is no need to have a special rectangle for captions, they go into accumulator
too.

To simplify positioning of the rectangles, the above example uses relative
"geographic" notation: the last argument of \kbd{plotcopy} specifies the
corner of the graph (say, northwest) one counts offset from.
(Positive offsets go into the rectangle.)

To demonstrate yet another useful plotting function, design a program to
plot Taylor polynomials for a $\sin x$ about 0.  For simplicity, make the
program good for any function, but assume that a function is odd, so only
odd-numbered Taylor series about 0 should be plotted.  Start with defining
some useful shortcuts

\bprog
xlim=13;    ordlim=25;     f(x)=sin(x);
default(seriesprecision,ordlim)
$\{$ farray(x)=
\q\q    vector((ordlim+1)/2,k,
\q\q\q\q\q\q            subst(Pol(f(1.*tt+tt\pow (2*k)*O(tt)),tt),
\q\q\q\q\q\q\q\q\q                  tt,x))
$\}$
\eprog

\noindent Function \kbd{farray(x)} returns a vector of Taylor polynomials for
$f(x)$.  We want to plot $f(x)$ into a rectangle, then make the rectangle
which is 1.2 times higher, and plot Taylor polynomials into the larger
rectangle.  Assume that the larger rectangle takes 0.9 of the final plot.

First of all, we need to measure the height of the smaller rectangle:

\bprog
plotinit(3, 0.9, 0.9/1.2, 1);
curve\_only=no\_x\_axis+no\_y\_axis+noframe;
lims = plotrecth(3, x= -xlim, xlim, f(x), recursive+curve\_only,16);
h = lims[4] - lims[3];
\eprog

\noindent Next step is to create a larger rectangle, and plot the Taylor
polynomials into the larger rectangle:

\bprog
plotinit(4, 0.9, 0.9,     1);
plotscale(4, lims[1], lims[2], lims[3]- h/10, lims[4] + h/10)
plotrecth(4, x= -xlim, xlim, farray(x), norescale);
\eprog

Here comes the central command of this example:

\kbd{plotclip(4)}

\noindent What does it do?  The command \kbd{plotrecth(\dots, norescale)}
scales the graphs according to coordinate system in the 
rectangle, but it does not pay any other attension to the size of
the rectangle.  Since \kbd{xlim} is 13, the Taylor polynomials take 
very large values in the interval \kbd{-xlim...xlim}.  In particular,
significant part of the graphs is going to be {\it outside} of the rectangle.
\kbd{plotclip} removes the parts of the
plottings which fall off the boundary of the rectangle.

\bprog
plotinit(2)
plotscale(2, 0.0, 1.0, 0.0, 1.0)
plotmove(2,0.5,0.975)
plotstring(2,"Multiple Taylor Approximations",center+vcenter)

plotdraw([2, 0, 0, 3, 0.05, 0.05 + 0.9/12, 4, 0.05, 0.05], relative)
\eprog

These commands draw a caption, and combine 3 rectangles (one with the
caption, one with the graph of the function, and one with graph of Taylor
polynomials) together.

This finishes our short survey of PARI plotting functions, but let us add
some remarks.  First of all, for a typical output device the picture is
composed of small colored squares (pixels) (as a very large checkerboard).
Each output rectangle is in fact a union of such squares.  Each drop
of paint in the rectangle will color a whole square it is in.  Since the
rectangle has a coordinate system, sometimes it is important to know how
this coordinate system is positioned with respect to the boundaries of
these squares.

The command \kbd{plotscale} describes a range of $x$ and $y$ in the
rectangle.  The limit values of $x$ and $y$ in the coordinate system are
coordinates {\it of the centers} of corner squares.  In particular,
if ranges of $x$ and $y$ are $[0,1]$, then the segment which connects (0,0)
with (0,1) goes along the {\it middle} of the left column of the rectangle.
In particular, if we made tiny errors in calculation of endpoints of this
segment, this will not change which squares the segment intersects, thus
the resulting picture will be the same.  (It is important to know such details
since many calculations in PARI are approximate.)

Another consideration is that all examples we did in this section were
using relative sizes and positions for the rectangles.  This is nice, since
different output devices will have very similar pictures, while we
did not need to care about particular resolution of the output device.
On the other hand,
using relative positions does not guarantie that the pictures will be
similar.  Why?  Even if two output devices have the same resolution,
the picture may be different.  The devices may use fonts of different
size, or may have a different ``unit of length''.

The information about the device in PARI is encoded in 6 numbers: resolution,
size of a character cell of the font, and unit of length, all separately
for horizontal and vertical direction.  These sizes are expressed as
numbers of pixels.  To inspect these numbers one may use the function
\kbd{plothsizes}.  The ``units of length'' are currently used to calculate
right and top gaps near graph rectangle of \kbd{ploth}, and gaps for
\kbd{plotstring}.  Left and bottom gaps near graph rectangle are calculate
using both units of length, and sizes of character boxes (so that there
is enough place to print limits of the graphs).

What does it show?  Using relative sizes during plotting produces {\tt
approximately} the same plotting on different devices, but does not
ensure that the plottings ``look the same''.  Moreover, ``looking the same''
is not a desirable target, ``looking tuned for the environment'' will
be much better.  If you want to produce such fine-tuned plottings, you need
to abandone a relative-size model, and do your plottings in units of
pixels.  To do this one removes flag \kbd{relative} from the above examples,
which will make size and offset arguments interpreted in units of pixels.
After querying sizes with \kbd{plothsizes} one can fine-tune sizes
and locations of subrectangles to the details of an arbitrary plotting device.

To check how good your fine-tuning is, you may test your graphs with a
medium-resolution plotting (as many display output devices are), and
with a low-resolution plotting (say, with \kbd{plotterm("dumb")} of gnuplot).