\chapter{A GP tutorial}

This chapter is intended to be a guided tour and a tutorial to the GP calculator.
Many examples will be given, but each time a new function is used, the reader
should look at the appropriate section in the manual for detailed explanations.
Hence although this chapter can be read independently of the rest (for example
to get rapidly acquainted with the possibilities of GP without having to read
the whole reference manual), the reader will profit most from it by reading it
in conjunction with the reference manual.
\medskip
\section{Greetings!}
\medskip
So you are sitting in front of your workstation (or terminal, or PC,\dots), and you
type {\bf gp} to get the program started (remember to always hit the {\tt enter} key
and not the {\tt return} key on a Macintosh computer). 

It says hello in its particular manner, and then waits for you after its {\tt prompt},
initially {\bf ?}.

Type $2+2$. What happens? Maybe not what you expect. First of all,
of course, you should tell GP that your input is finished, and this is done by hitting
the {\tt return} (or {\tt newline}) key, or the {\tt enter} key on the Mac. If you
do exactly this, you will get the expected answer. However some of you may be used
to other systems like Macsyma or Maple. In this case, you will have subconsciously
ended the line with a semicolon {\bf ;} before hitting {\tt return}, since this is
how it is done on those systems. In that case, you will simply see GP answering
you with a smug expression, i.e. a new prompt and no answer!
This is because a semicolon at the
end of a line in GP tells it to keep the result, but not to print it (you will
certainly want to use this feature if the output is several pages long).

Try $27*37$. Wow! even multiplication works.

Now this session is getting lengthy, so the second thing one needs to learn is to
quit. Each system has its quit signal (to name a few: {\tt quit}, {\tt quit();},
{\tt exit}, {\tt system},\dots). In GP it is {\tt \bs q} (backslash), the
{\tt q} being of course for quit. Try it.

Now you've done it! You're out of GP, so how do you want to continue studying this
tutorial? Get back in please (see above).

Let's get to more serious stuff. Let's see, I seem to remember that the decimal
expansion of $1/7$ has some interesting properties. Let's see what GP has to say
about this. Type $1/7$. What? This computer is making fun of me, it just spits back
to me my own input, that's not what I want! 

Now stop complaining, and think a little. This system has been written mainly for
pure mathematicians, and not for physicists (although they are welcome to use it :-).
And mathematically, $1/7$ is an element of the field $\Bbb Q$ of rational numbers,
so how else but $1/7$ can the computer give the answer to you? (well maybe $2/14$, but
why complicate matters?). Seriously, the basic point here is that PARI (hence GP)
will almost always try to give you a result which is as precise as possible
(we will see why ``almost'' later), hence since here you gave an operation whose
result can be represented exactly, that's what it gives you.

OK, but I still want the decimal expansion of $1/7$. No problem. Type one of the
following: $1./7$, $1/7.$, $1./7.$, $0.+1/7$,\dots. Immediately 28 decimals of this
fraction appear, and the repeating pattern is $142857$. The reason is that you have
included in the operations numbers like $0.$, $1.$ or $7.$ which are {\sl imprecise\/}
real numbers, hence GP cannot give you an exact result.

Why 28 decimals by the way? well, it is the default initial precision, as indicated
when you launch GP. This has been chosen so that the computations are very fast,
and give already 12 decimals more accuracy than conventional double precision
floating point operations. Only large mainframes or supercomputers have 28 digits
of precision in their standard libraries, and that is the absolute limit.

Not here of course. You can extend the precision (almost) as much as you like as
we will see in a moment.

I'm getting bored, why don't we get on with some more exciting stuff?
Well, try {\tt exp(1)}. Presto, comes out the value of $e$ to 28 digits.
Try {\tt log(exp(1))}. Well, it's not exactly equal to 1, but pretty close!
That's what you lose by working numerically.

Now try {\tt exp(pi$*$sqrt(163))}. Hmmm, since from our last example we suspect
that the last digit may be wrong, can this really be an integer?
This is the time to change precision. Type {\tt \bs precision=50},
then try {\tt exp(pi$*$sqrt(163))} again. We were right to suspect that the last
decimal was incorrect, since we get even more nines than before, but it is now
convincingly clear that this is not an integer. Maybe it's a bug in PARI, and
the result is really an integer? Type {\tt sqr(log(\%)/pi)} immediately after
the preceding computation
({\tt \%} means the result of the last computed expression. More generally, the
results are numbered {\tt \%1, \%2, \dots} {\bf including} the results that you
do not want to see printed by putting a semicolon at the end of the line, and
you can evidently use all these quantities in any further computations. {\tt sqr}
is the square function ({\tt sqr(x)=x$*$x}), not to be confused with {\tt sqrt}
which is the square root function). The result seems to be indistinguishable
from 163, hence it does not seem to be a bug.

In fact it is a known fact that {\tt exp(pi$*$sqrt(n))} not only is not an integer
or a rational number, but is even a transcendental number when $n$ is a positive
rational number.

So GP is just a fancy calculator, able to give me more decimals than I will ever
need? Not so, GP is incredibly more powerful than an ordinary calculator, even
independently of its arbitrary precision possibilities.
\medskip
{\bf Additional comments}
\medskip
You are supposed to skip this at first, and come back later.
\smallskip
1) It seems that the text implicitly says that as soon as an imprecise number
is entered, the result will be imprecise. Is this always true? There is a unique
exception: when you multiply an imprecise number by the exact number 0, you will
get the exact 0. Try: $0*1.4$, and compare with $0.*1.4$.
\smallskip
2) Not only can the number of decimal places of real numbers be large, but the
number of digits of integers also. Try $100!$. It is never necessary to tell
in advance to GP the size of the integers that it will encounter, since this is
adjusted automatically. On the other hand, for many computations with
real numbers, it is necessary to specify a default precision (initially 28 digits).
\smallskip
3) Come back to 28 digits of precision ({\tt \bs precision=28}), and type
{\tt exp(24$*$pi)}. What is this $*$ that you get at the end of the result?

Well, GP never wants you to believe that a result is correct when it is not.
We are working with 28 digits of precision, but the integer part of {\tt exp(24$*$pi)}
has 33 decimal digits. Hence GP dutifully prints out 33 digits, but it prints a
$*$ to warn you that the last few digits are wrong (without telling you how
many). The $*$ can thus be interpreted as a decimal point with a warning message.
\smallskip
4) To avoid this $*$, there are two ways. One is of course to increase the precision
to more than 33 decimals. Let's try it. To give it a wide margin, we set the
precision to 40 decimals. Then we recall our last result ({\tt \%} or {\tt \%n}
where {\tt n} is the number of the result). What? We still have the $*$! Do you
understand why?

Again let's try to see what's happening. The number you recalled had been computed only
to 28 decimals, and even if you set the precision to 1000 decimals, GP knows that
your number has only 28 digits of accuracy but an integral part with 33 digits. So
you haven't improved things by increasing the precision. Or have you? What if we
retype {\tt exp(24$*$pi)} now that we have 40 digits? Try it. Now we do not have
any warning, and we see the beginning of what's after the decimal point. We also
see that at 28 decimals the last 6 digits were wrong.
\smallskip
5) There is still another method to avoid the warning: put yourself back into
precision 28. Now type {\tt \bs format=e0.28}. The $e$ instead of the $g$
which is printed when GP greets you, means that real numbers (by opposition to
integers or rational numbers) will be printed in {\sl exponential\/} format, i.e.
with a E xxx denoting powers of 10. Try again {\tt exp(24$*$pi)} in precision 28.
Now we don't have any warning, and the reason is clear: GP has to print only a
single digit for the integer part, and then it prints out as many decimals as
was asked by the format (here 28).
\smallskip
6) {\sl Warning\/}. Try the following: starting in precision 28, type 
successively {\tt \bs format=e0.50}, then {\tt exp(24$*$pi)}. Do you
understand why the result is so bad, and why there are lots of zeros at the end?
Convince yourself by typing {\tt log(exp(1))}. The moral is that the {\tt format}
command changes only the output format, but {\sl not\/} the default precision.
On the other hand, the {\tt precision} command changes both the precision and
the output format.
\medskip
\section{Warming up}
\medskip
Another thing you better get used to pretty fast is error messages. Try typing
$1/0$. Couldn't be clearer. Taking again our universal example in precision 28,
type {\tt floor(exp(24$*$pi))} ({\tt floor} is the mathematicians' integer part, not
to be confused with {\tt trunc}, which is the computer scientists': 
{\tt floor$(-3.4)$} is equal to $-4$ while {\tt trunc$(-3.4)$} is equal to $-3$).
You get a more cryptic error message, which you would immediately understand if you
had read the additional comments of the preceding section. Since I told you not
to read them, the explanation is simply that GP is unable to compute the integer
part of {\tt exp(24$*$pi)} given only 28 decimals of accuracy, since it has 33 digits.

Some error messages are even much more cryptic than that and are sometimes not so
easy to understand (well, it's nothing compared to \TeX's error messages! :-).

For instance, try {\tt log(x)}. Not really clear, is it? But no matter, it simply
tells you that GP simply does not understand what {\tt log(x)} is (although it
does know the {\tt log} function).

Now let's try {\tt sqrt($-$1)} to see what error message we get now. Haha! GP even
knows about complex numbers, so impossible to trick it that way. Similarly, try
typing {\tt log($-$2)}, {\tt exp(i$*$pi)}, {\tt i\^{}i},\dots. So we have a lot of real and complex
analysis at our disposal (note that there is always a specific branch of
multivalued complex transcendental functions which is taken, specified in the manual).

Just for fun, let's try {\tt 6$*$zeta(2)/pi\^{}2}. Pretty good? (If you
don't know what the Riemann zeta function is, don't worry, we won't
use it in this tutorial).
\medskip
Now GP didn't seem to know what {\tt log(x)} was, although it did know how to compute
numerical values of {\tt log}. This is annoying. Maybe it knows the exponential
function? Let's give it a try. Type {\tt exp(x)}. What's this? If you had had any experience
with other systems, the answer should have simply been {\tt exp(x)} again. But here the
answer is the Taylor expansion of the function around {\tt x=0}, to 16 terms (16 is the
default {\tt serieslength}, which can be changed by typing {\tt \bs serieslength=n},
where {\tt n} is the number of terms that you want in your power series (note the {\tt O(x\^{}16)}
which ends the series, and which is trademark of this type of object in GP. It is the familiar
``bigoh'' notation of analysis).

You will thus automatically get the Taylor expansion of any function that can be expanded
around {\tt x=0}, and incidentally this explains
why we weren't able to do anything with {\tt log(x)} which is not defined at 0. By if we try
{\tt log(1$+$x)}, then it works. But what if we wanted the expansion around a point different
from 0? Well, you're able to change {\tt x} into {\tt x$-$a}, aren't you? So for instance
you can type {\tt log(x$+$2)} to have the expansion of {\tt log} around {\tt x=2}. 

As exercises, try {\tt cos(x)}, {\tt cos(x)\^{}2$+$sin(x)\^{}2}, {\tt exp(cos(x))},
{\tt exp(exp(x)$-$1)}, {\tt gamma(1$+$x)}, with different values of {\tt serieslength} if you like.

Let's try something else: type {\tt (1$+$x)\^{}3}. No bigoh here, since the result is a polynomial.
Haha, but I have learnt that if you do not take exponents which are integers greater or equal
to 0, you obtain a power series with an infinite number of nonzero terms. Let's try.
Type {\tt (1$+$x)\^{}($-3$)} (the parentheses around $-3$ are not necessary but make things 
easier to read). Surprise! Contrary to what we expected, we don't get a power series but a
rational function. Again this is for the same reason that $1/7$ just gave you $1/7$: the result
being exact, PARI doesn't see any reason to make it nonexact.

But I still want that power series. To obtain it, just do as in the $1/7$ 
example: type 

{\tt (1$+$x)\^{}($-$3)$+$O(x\^{}{16})}, {\tt (1$+$O(x\^{}{16}))$*$(1$+$x)\^{}($-$3)},
{\tt (1$+$x$+$O(x\^{}{16}))\^{}($-$3)},etc\dots

Now try {\tt (1$+$x)\^{}(1/2)}. Now we obtain a power series, since the result is an object
which PARI does not know how to represent exactly (we could teach PARI about algebraic functions,
but then take {\tt (1$+$x)\^{}pi} as another example). This gives us still another solution to
our preceding exercise: we can type {\tt (1$+$x)\^{}($-3.$)}. Since $-3.$ is not an exact
quantity, PARI has no means to know that we are dealing with a rational function, and
will instead give you the power series.

Finally, if you want to be really fancy, you can type {\tt taylor((1$+$x)\^{}($-3$),x)}
(look at the entry for {\tt taylor} for the description of the syntax), but this is in fact
almost never used.
\smallskip
To summarize, in this section we have seen that in addition to integers,
real numbers and rational numbers, PARI can handle complex numbers,
polynomials, power series, rational functions. A large number of functions 
exist which handle these types, but in this tutorial we will only look at
a few. \medskip
{\bf Additional comments}
\medskip
As before, you are supposed to skip this at first, and come back later.
\smallskip
1) A complex number has a real part and an imaginary part (who would have guessed?).
However, beware that when the imaginary part is the exact integer zero, it is not printed,
but the complex number is not converted to a real number. Hence it may {\sl look\/} like
a real number without being one, and this may cause some confusion in programs which expect
real numbers. For example, type {\tt 3$+$i$-$i}. The answer is, as expected 3. But it is still
a complex number. Hence if you then type {\tt (1$+$x)\^{}\%}, instead of getting the polynomial
{\tt 1$+$3$*$x$+$3$*$x\^{}2$+$x\^{}3} as expected, you obtain a power series. Worse, when you try to
apply an arithmetic function, say the Euler {\tt phi} function, you get the error message
which says that ``arithmetic functions want integer arguments'', which you would have
guessed yourself, but the message is difficult to understand since 3 looks like a genuine
integer!!! (Please read again if this is not clear. It is a common source of incomprehension).

Similarly, {\tt 3$+$x$-$x} is not the integer 3 but the constant polynomial equal
to 3.

If you want the final expression to be in the simplest form possible (for
example before applying an arithmetic function, or simply because things
will go faster afterwards), apply the function {\tt simplify} to the 
result.

2) As already stated, power series expansions are always implicitly around {\tt x=0}. If you
want them around {\tt x=a}, replace {\tt x} by {\tt z$+$a} in the function that you want to
expand. For example, to have the expansion of {\tt exp(x)} around {\tt x=1}, we simply type
{\tt exp(z$+$1)}, where {\tt z} stands for {\tt x$-$a}. For complicated functions, it may be
simpler to use the substitution function {\tt subst}. For example, if 
{\tt p=1/(x\^{}4$+$3$*$x\^{}3$+$5$*$x\^{}2$-$6$*$x$+$7)}, you may not want to retype this, replacing {\tt x} by {\tt z$+$a},
so you can write {\tt subst(p,x,z$+$a)} (look up the exact description of the {\tt subst} function).

Now try typing {\tt p=1$+$x$+$x\^{}2$+$O(x\^{}10)}, then {\tt subst(p,x,z$+$1)}. Do you 
understand why you get an error message?
\medskip
\section{The Rest of the PARI Types}
\medskip
Let's handle first some more the basic PARI types.

Type {\tt p=x$*$exp($-$x)}. As expected, you get the power series expansion
to 16 terms (if you have not changed the default). Now type
{\tt pr=reverse(p)}. You are asking here for the {\sl reversion\/} of the
power series $p$, in other words the inverse function. This is possible only
for power series whose first non-zero coefficient is that of $x^1$.
To check the correctness of the result, you can type {\tt subst(p,x,pr)}
or {\tt subst(pr,x,p)} and you should get back {\tt x$+$O(x\^{}17)}.

Now the coefficients of {\tt pr} obey a very simple formula. First, we would
like to multiply the coefficient of {\tt x\^{}n} by {\tt n!} (in the case of the
exponential function, this would simplify things considerably!). The PARI
function {\tt laplace} does just that. So type {\tt ps=laplace(pr)}. The
coefficients now become integers, which can be immediately recognized by
inspection. The coefficient of {\tt x\^{}n} is now equal to $\text{n}^{\text{n-1}}$. In other words, we have 
\def\n{{\rm n}}
$$\text{pr}=\sum_{\n\ge1}\dfrac{\n^{\n-1}}{\n!}X^{\n}.$$Do you know how to
prove this? (If you have never seen this, the proof is difficult.)
\smallskip
Of course PARI knows about vectors (rows and columns are distinguished, even
though mathematically there is no difference) and matrices. Type for example
{\tt [1,2,3,4]}. This gives the row vector whose coordinates are 1, 2, 3 and 4.
If you want a column vector, type {\tt [1,2,3,4]\~{}}, the tilde meaning of
course transpose. You don't see much difference in the output, except for the
tilde at the end. However, now type {\tt \bs b}: lo and behold, the vector
does become a column. This command is used mainly for this purpose.

Type {\tt m=[a,b,c;d,e,f]}. You have just entered a matrix with 2 rows and 3
columns. Note that the matrix is entered by {\sl rows\/} and the rows are
separated by semicolons ``;''. The matrix is printed naturally in a rectangle
shape. If you want it printed horizontally just as you typed it, type
{\tt \bs r}, or if you want this type of printing to be the permanent default
type {\tt \bs p} twice (the first time puts GP in a mode which in the present
version you do not want to use). Type {\tt \bs p} a third time if you want
to come back to the original default.

Now type {\tt m[1,2]}, {\tt m[1,]}, {\tt m[,2]}. Are explanations necessary?
(In an expression such as {\tt m[j,k]}, the {\tt j} always refers to the
row number, and the {\tt k} to the column number, and the first index is
always 1, never 0. This default cannot be changed.)

Even better, type {\tt m[1,2]=5;m} (the semicolon allows us to put several
instructions on the same line. It is an absolute necessity when we program
in GP). Now type {\tt m[1,]=[15,$-$17]}. No problem. Finally type
{\tt m[,2]=[j,k]}. You have an error message since you have typed a row
vector, while {\tt m[,2]} is a column vector. If you type instead
{\tt m[,2]=[j,k]\~{}} it works.
\smallskip
Type now {\tt h=hilbert(20)}. You get the so-called ``Hilbert matrix'' whose 
coefficient of row {\tt i} and column {\tt j} is equal to
{\tt 1$/$(i$+$j$-$1)}. Incidentally, the matrix {\tt h} takes too
much room. If you don't want to see it, simply type a semi-colon ``;'' at the
end of the line ({\tt h=hilbert(20);}). This is an example of a ``precomputed''
matrix, built-in PARI. There are only a few. We will see later an example
of a much more general construction.

What is interesting about Hilbert matrices is that first their inverses and
determinants can be computed explicitly (and the inverse has integer
coefficients), and second they are numerically very unstable, which make
them a severe test for linear algebra packages in numerical analysis.
Of course with PARI, no such problem can occur: since the coefficients are
given as rational numbers, the computation will be done exactly, so there
cannot be any numerical error. Try it. Type {\tt d=det(h)} (you have to be
a little patient, this is quite a complicated computation). The result is
a rational number (of course) of numerator equal to 1 and denominator having
226 decimal digits. How do I know, by the way? I did not count! Instead,
simply type {\tt 1.$*$d}. The result is now in exponential format, and the
exponent gives us the answer.

Now type {\tt hr=1.$*$h;} (do not forget the semicolon, we don't want to see
all the junk!), then {\tt dr=det(hr)}. You notice two things. First the
computation, although not instantaneous, is much faster than in the rational
case. The reason for this is that PARI is handling real numbers with 28 digits
of accuracy, while in the rational case it is handling integers having up to
226 decimal digits.

The second more important fact is that the result is terribly wrong. If you
compare with {\tt 1.$*$d} computed earlier, which is correct, you will see that
only 2 decimals agree! This catastrophic instability is as already 
mentioned one of the characteristics of Hilbert matrices. In fact, the 
situation is much worse than that. Type {\tt norml2(1/h$-$1/(h$*$1.))}
(the function {\tt norml2} gives the square of the $L^2$ norm, i.e. the sum of
the squares of the coefficients). Again be patient since computing {\tt 1/h}
will take even more time (not much) than computing {\tt det(h)}. The result
is larger than $10^{50}$, showing that some coefficients of {\tt 1/(h$*$1.)}
are wrong by as much as $10^{24}$ (the largest error is in fact equal
to $7.644 10^{24}$ for the coefficient of row 15 and column 14, which is a
28 digit integer).

To obtain the correct result after rounding for the inverse, we have to use
a default precision of 56 digits (try it).
\smallskip
Although vectors and matrices can be entered manually, by typing explicitly
their elements, very often the elements satisfy a simple law and one uses a
different syntax. For example, assume that you want a vector whose {\tt k}-th
coordinate is equal to {\tt k\^{}2}. No problem, type for example
{\tt vector(10,k,k\^{}2)} if you want a vector of length 10. Similarly, if
you type 

\centerline{{\tt matrix(5,5,j,k,1$/$(j$+$k$-$1))}}

\noindent you will get the Hilbert
matrix of order 5 (hence the {\tt hilbert} function is unnecessary).
The {\tt j} and {\tt k} represent dummy variables which are used to number
the rows and columns respectively (in the case of a vector only one is present
of course). You must not forget, in addition to the dimensions of the vector
or matrix, to indicate explicitly the names of these variables.

{\bf Important warning.} The letter {\tt i} is reserved for the complex number
equal to the square root of $-1$. Hence it is absolutely forbidden to use it
as a variable. Try typing {\tt vector(10,i,i\^{}2)}, the error message that you
get clearly indicates that GP does not consider {\tt i} as a variable.
Also remember that (for now) GP is case insensitive so the variable {\tt I} is
just as forbidden.

When creating vectors or matrices, it is often useful to use boolean operators
and the {\tt if()} statement (see the section on programming for details on
using this statement). Indeed, an {\tt if} expression has a value, which is of
course equal to the evaluated part of the {\tt if}. So for example you can
type

\centerline{{\tt matrix(8,8,j,k,if((j$-$k)$\%$2,x,0))} }

to get a checkerboard matrix of {\tt x} and {\tt 0}. Note however that a
vector or matrix must be {\it created} first before being used. For example,
it is forbidden to write 

\centerline{{\tt for(j=1,5,v[j]=1/j)}}

if the vector {\tt v} has not been created beforehand (for example by a
{\tt v=vector(5,j,0)} command).
\medskip
{\bf Important remark}
\medskip
A drawback of the present GP implementation is that each time a vector or
matrix is modified, even in a single component, the whole vector is copied
(this is of course not the case in library mode). This implies that the 
necessary time for filling up a vector recursively is not proportional to
its length but to the the square of it length, which can be painful if the
length is at all large. The only exception to this, and an important one of
course, is the direct use of the {\tt vector} or {\tt matrix} functions which
do {\it not} fill the vector or matrix recursively but directly, hence take
linear time.

Consider the following example. We want to store the first 400 values of the
sequence defined by $u_0=0$, $u_1=5$, $u_{n+1}=\sqrt{u_n^2+3}-u_{n-1}$.
A reasonable way to do this is through the following GP program:

{\tt a=0;b=5;v=[a,b];for(j=3,400,c=sqrt(b$*$b$+$3)$-$a;a=b;b=c;v=concat(v,c))}

If you time the execution of this program and then try the same thing with
400 replaced by 800 you will see that the time has almost been multiplied by
4 although the computations are not harder for the last 400 iterations than for
the first.

One method to bypass this drawback of the GP implementation is to try to limit
to short vectors. In the example above, we would for example store each hunk of
100 consecutive values in {\it different} vectors, and start with the last 
values reached to fill in the next hunk. At the end we can concatenate all the
vectors together to get a single vector as desired.
\medskip
The last PARI types which we have not yet played with are types which are
closely linked to number theory (hence people not interested in number theory
can skip).

The first is the type ``integermod''. Let us see an example. Type
{\tt n=10\^{}15$+$3}. We want to know whether {\tt n} is prime or not. Of course
we could make use of the builtin facilities of PARI, but let us do otherwise.
We first trial divide by the builtin table of primes, by default up to
500000. We slightly cheat here and use the function {\tt smallfact} which
does exactly this. So type {\tt smallfact(n)}. The result is a 2 column
matrix (as all factoring programs), the first column giving the primes and
the second their exponents. Here we get a single row, telling us that {\tt n}
is not divisible by any prime up to the default (500000). We could now
trial divide further, or even cheat completely and call the PARI function
{\tt factor}, but before we do this let us see how to get an answer ourselves.

By Fermat's little theorem, if {\tt n} is prime we must have 
$\text{a}^{\text{n-1}}\equiv1\pmod{\text{n}}$ for all {\tt a} not divisible by
{\tt n}. Hence we could try this with {\tt a=2} for example. But 
$2^{\text{n-1}}$ is a number with approximately $3\cdot10^{14}$ digits, hence
impossible to write down, even less to compute. But instead type
{\tt a=mod(2,n)}. This creates the number 2 considered now as an element of
the ring $R=\Z/\text{n}\Z$. The elements of $R$, called integermods, can
always be represented by numbers smaller than {\tt n}, hence very small.
Fermat's theorem can be rewritten {\tt $\text{a}^{\text{n-1}}=$mod(1,n)} in
the ring $R$, and this can be computed very efficiently. Type
{\tt a\^{}(n$-$1)}. The result is definitely {\sl not\/} equal to {\tt mod(1,n)},
thus {\sl proving\/} that {\tt n} is not a prime (if we had obtained 
{\tt mod(1,n)} on the other hand, it would have given us a hint that {\tt n}
is maybe prime, but never a proof).

To find the factors is another story. One must use less naive techniques than
trial division (or be very patient). To finish this example, type 
{\tt factor(n)} to see the factors. Since the smallest factor is 14902357,
you would have had to be very patient with trial division!
\smallskip
The second specifically number-theoretic type is the $p$-adic numbers. I have
no room for definitions, so please skip if you have no use for such beasts.
A $p$-adic number is entered as a rational or integer valued expression to which
is added {\tt O(p\^{}n)} (or simply {\tt O(p)} if {\tt n=1}) where {\tt p} is
the prime and {\tt n} the $p$-adic precision. Apart from the usual arithmetic
operations, you can apply a number of transcendental functions.
For example, type {\tt n=569$+$O(7\^{}8)}. Then type {\tt s=sqrt(n)}, you
obtain one of the square roots of {\tt n} (if you want to check, type 
{\tt s$*$s$-$n}). Type now {\tt l=log(n)}, then {\tt e=exp(l)}. If you
know about $p$-adic logarithms, you will not be surprised that {\tt e} is not
equal to {\tt n}. Type {\tt (n/e)\^{}6}: {\tt e} is in fact equal to {\tt n}
times a (p$-$1)-st root of unity.

Incidentally, if you want to get back the integer 569 from the $p$-adic number
{\tt n}, type {\tt trunc(n)}.
\smallskip
The third number-theoretic type is the type ``quadratic number''. This type
is specially tailored so that we can easily work in a quadratic extension
of a base field (usually $\Bbb Q$). It is a generalization of the type
``complex''. To start, we must specify which quadratic field we want to work
in. For this, we use the function {\tt quadgen} applied to the 
{\sl discriminant\/} {\tt d} (as opposed to the radicand) of the quadratic
field. This returns a number (always printed as {\tt w}) equal to
{\tt (d$+$a)$/$2} where {\tt a} is equal to 0 or 1 according to whether {\tt d}
is even or odd. The behavior of {\tt quadgen} is a little special: although
its result is always printed as {\tt w}, the variable {\tt w} itself is not
set to that value. Hence it is necessary to write systematically
{\tt w=quadgen(d)} using the variable name {\tt w} (or {\tt w1} etc.. if you
have several quadratic fields), otherwise you will be confused.

So type {\tt w=quadgen($-163$)}, then {\tt char(w,x)} which asks for the
characteristic polynomial of {\tt w} expressed with the variable {\tt x}.
The result shows what {\tt w} will represent. You can also ask for {\tt 1.$*$w}
to see which root of the quadratic has been taken, but this is rarely 
necessary. We can now play in the field $\Bbb Q(\sqrt{-163})$. Type for example
{\tt w\^{}10, norm(3$+$4$*$w), 1$/$(4$+$w)}. More interesting, type
{\tt a=mod(1,23)$*$w} then {\tt b=a\^{}264}. This is a generalization of Fermat's
theorem to quadratic fields. If you do not want to see the 23 all the time,
\def\x{{\tt x}}
type {\tt lift(b)}. Another example, type 
{\tt p$=\x^2+\text{w}*\x+5*\text{w}+7$}, then {\tt norm(p)} (or even
better {\tt simplify(norm(p))}, see above). We thus
obtain the quartic equation over $\Q$ corresponding to the relative quadratic
extension over $\Q(\text{w})$ defined by {\tt p}.

On the other hand, if you type {\tt wr=sqrt(w\^{}2)}, do not expect to get back
{\tt w}. Instead, you get the numerical value, the function {\tt sqrt} being
considered as a ``transcendental'' function, even though it is algebraic.
Type {\tt algdep(wr,2)}. This is one way to get {\tt w} back. Similarly,
type {\tt algdep(sqrt(3$*$w$+$5),4)}. See the users manual for the function
{\tt algdep}.
\smallskip
The fourth number-theoretic type is the type ``polymod'', i.e. polynomial
modulo another polynomial. This type is used to work in algebraic extensions,
for example elements of number fields (if the base field is $\Bbb Q$), or
elements of finite fields (if the base field is $\Z/p\Z$ for a prime $p$,
defined by an integermod). In a sense it is a generalization of the type
quadratic number. The syntax used is the same as for integermods.
For example, instead of typing {\tt w=quadgen($-$163)}, you can type
{\tt w=mod(x,x\^{}2$-$x$+$41)}. Then exactly as in the quadratic case you can
type {\tt w\^{}10, norm(3$+$4$*$w), 1$/$(4$+$w), a=mod(1,23)$*$w, b=a\^{}264},
obtaining of course the same results (type {\tt lift(...)} if you don't want
to see the polynomial {\tt x\^{}2$-$x$+$41} repeated all the time). Of course,
the basic interest is that you can work in any degree, not only quadratic.

There is however a slight difference in behavior. Keeping our polymod {\tt w},
type {\tt 1.$*$w}. As you can see, the result is not the same. Type
{\tt sqrt(w)}. Here, we obtain a vector with 2 components, the two components
being the principal branch of the square root of all the possible embeddings
of {\tt w} in $\Bbb C$ (NOT the two square roots). More generally, if {\tt w}
was of degree n, we would get an n-component vector, and similarly for
transcendental functions.

We have at our disposal the usual arithmetic functions, plus a few others.
For example, type {\tt a=mod($\x,\x^3-\x-1$)} defining a cubic extension.
We can for example ask for {\tt b=a\^{}7}. But now assume we want to express
{\tt a} as a polynomial in {\tt b}. This is possible since {\tt b} is also a
generator of the same field. No problem, type {\tt modreverse(b)}. This gives
a new defining polynomial for the same field (i.e. the characteristic 
polynomial of {\tt b}), and expresses {\tt a} in terms of this new polymod, i.e.
in terms of {\tt a}. We will see this in much more detail in the number field
section.
\medskip
\section{Elementary Arithmetic Functions}
\medskip
Since PARI is aimed at number theorists, it is not surprising that there
exists a large number of arithmetic functions (see the list in the 
corresponding section of the users manual). We have already seen several, such
as {\tt factor}. Note that {\tt factor} factors not only integers, but also
polynomials. Type for example {\tt factor(x\^{}15$-$1)}. You can also ask to
factor a polynomial modulo a prime p ({\tt factmod}) and even in a finite
field which is not the prime field ({\tt factfq}).

Evidently, you have functions for computing GCD's, extended GCD's, solving
the Chinese remainder theorem and so on.

In addition to the factoring facilities, you have a few functions related to
primality testing such as {\tt isprime}, {\tt ispsp} and {\tt nextprime}.
Note that in all these cases, only strong pseudoprimes are produced, there
is no sophisticated primality test.

The usual multiplicative arithmetic functions exist (the M\"obius $\mu$ 
function ({\tt mu}) the Euler $\phi$ function ({\tt phi}), the $\omega$
and $\Omega$ functions ({\tt omega} and {\tt bigomega}), the sum of divisors
functions ({\tt sigma} and {\tt sigmak}), etc\dots

You can compute continued fractions. For example, if you type 
{\tt\bs precision=1000} then {\tt cf(exp(1))} you obtain the continued fraction
of the base of natural logarithms, which as you can see obeys a very simple
pattern (can you prove it?).

In many cases, one wants to perform some task only when an arithmetic condition
is satisfied. GP gives you the following functions: {\tt isprime} as mentioned
above, {\tt issquare}, {\tt isfund} to test whether an integer is a fundamental
discriminant (i.e. the discriminant of the ring of integers of a quadratic 
field), and the {\tt forprime} and {\tt divsum} loops. Assume for example that
we want to compute the product of all the divisors of a positive integer $n$.
One way is to write:

\centerline{\tt round(exp(divsum(n,d,log(d))))}

(there is a very simple formula for this product: find and prove it).

If we want to know the list of primes $p$ less than 1000 such that 2 is a
primitive root modulo $p$, one way would be to write:

\centerline{\tt forprime(p=3,1000,if(primroot(p)$==$2,print(p),))}

Note that this assumes that {\tt primroot} returns the smallest primitive root,
and this is indeed the case.

Functions related to quadratic fields, binary quadratic forms and general
number fields will be seen in the next sections.
\medskip
\section{Working in Quadratic Number Fields}
\medskip
The simplest of all number fields outside $\Q$ are quadratic fields. Such 
fields are characterized by their discriminant, and even better, any nonsquare
integer $D$ congruent to 0 or 1 modulo 4 is the discriminant of a specific
order in a quadratic field. We can check whether this order is maximal by
using the command {\tt isfund(D)}.
Elements of a quadratic field are of the form $a+b\omega$, where $\omega$ is
chosen as $\sqrt{D}/2$ if $D$ is even and $(1+\sqrt{D})/2$ if $D$ is odd,
and are represented in Pari by quadratic numbers. To initialize working in
a quadratic order, one should start by the command 

\centerline{\tt w$=$quadgen(D)}

This sets {\tt w} equal to $\omega$ as above, and is printed {\tt w}. Note
however that if several different quadratic orders are used, a printed {\tt w}
may have several different meanings. For example if you type

\centerline{\tt a=quadgen($-$23);b=quadgen($-$15);}

Then ask for the value of {\tt a} and {\tt b}, both will be equal to {\tt w},
but of course they are not equal. Hence beware when dealing with several 
quadratic orders at once. 
\smallskip
In addition to elements of a quadratic order, we also want to be able to handle
ideals of such orders. In the quadratic case, it is equivalent to handle
binary quadratic forms, and this has been chosen in Pari. For negative
discriminants, quadratic forms are triples $(a,b,c)$ representing the form
$ax^2+bxy+cy^2$. Such a form will be printed as (and can be created by)

\centerline{\tt qfi(a,b,c)}

({\tt i} stands for imaginary).

Such forms can be multiplied, divided, powered as many Pari objects using
the usual operations, and they can also be reduced using the function
{\tt redimag} (it is not the purpose of this tutorial to explain what all these
things mean). In addition, Shanks's NUCOMP algorithm has been implemented
(functions {\tt nucomp} and {\tt nupow}), and this is usually a little faster.

Finally, you have at your disposal the functions {\tt classno} which gives 
(usually) the class number, the function {\tt hclassno} which gives the
Hurwitz class number, and the much more sophisticated {\tt buchimag} function
which gives the class number and class group structure, even for large 
discriminants.

Let us see examples of all this at work.

Type {\tt classno($-$10007)}. GP tells us that the result is 77. However, you may
have noticed in the explanation above that the result is only usually correct.
This is because the implementers of the algorithm have been lazy and have not
put the complete Shanks algorithm in Pari, causing it to fail in certain very
rare cases. In practice, it is almost always correct, and the much more 
powerful {\tt buchimag} program, which {\it is} complete, can give
confirmation.

So we may be a little suspicious of this class number. Let us check it.
First, we need to find a quadratic form of discriminant $-10007$. Since this
discriminant is congruent to 1 modulo 8, we know that there is an ideal of 
norm equal to 2, i.e. a binary quadratic form $(a,b,c)$ with $a=2$. To compute it
we type {\tt f=pf($-$10007,2)} ({\tt pf} stands for primeform, i.e. binary forms
corresponding to prime ideals). OK, now we have a form. If the class
number is correct, the very least is that this formed raised to the
power 77 should equal the identity. Let's check this.
Type {\tt f\^{}77}. We get a form starting
with 1, i.e. the identity, so this test is OK. Raising f to the power 11 and 7
does not give the identity, hence we now know that the order of f is exactly
77, hence the class number is a multiple of 77. But how can we be sure that
it is exactly 77 and not a multiple? Well, type

\centerline{\tt prodeuler(p=2,500,1./(1$-$kro($-$10007,p)/p))$*$sqrt(10007)/pi}

This is nothing else than an approximation to the Dirichlet class number
formula. The function {\tt kro} is the Kronecker symbol, in this case simply
the Legendre symbol. Note also that we have written {\tt 1./(1$-$...)} with
a dot after the 1. Otherwise, Pari may want to compute the whole thing as
a rational number, which would be terribly long and useless. In fact Pari
does no such thing in this particular case, but you never know. Better safe
than sorry!

We find 77.77, pretty close to 77, so things seem in order. Explicit bounds
on the prime limit to be used in the Euler product can be given which make
the above reasoning rigorous.

Let us try the same thing with $D=-3299$. {\tt classno} and the Euler product
convince us that the class number must be 27. However, we get stuck when we try
to prove this in the simple-minded way above. Indeed, we type
{\tt f=pf($-$3299,3)} (here 2 is not the norm of a prime ideal but 3 is), and
we see that f raised to the power 9 is equal to the identity. This is the
case for any other quadratic form we choose. So we suspect that the class
group is not cyclic. Indeed, if we list all 9 distinct powers of f, we see
that the primeform {\tt pf($-$3299,5)} is not on the list (although its cube
is as it must). This implies that the class group is probably equal to a
product of a cyclic group of order 9 by a cyclic group of order 3. The
Euler product plus explicit bounds prove this.

Another way to check it is to use the {\tt buchimag} function by typing
for example

\centerline{\tt buchimag($-$3299,1,6)}

Indeed the class group is as we thought (see Chapter 3 for the description of
the different components of the result).
\medskip
Working in real quadratic fields instead of complex ones, i.e. with $D>0$, is
not very different.

The same {\tt quadgen} function is used to create elements. Ideals are again
represented by binary quadratic forms $(a,b,c)$, this time indefinite. However,
the Archimedean valuations of the number field start to come into play (as
is clear if one considers ideles instead of ideals), hence in fact quadratic
forms with positive discriminant will be represented as a quadruplet
$(a,b,c,d)$ where the quadratic form itself is $ax^2+bxy+cy^2$ with $a$,
$b$ and $c$ integral, and $d$ is the Archimedean component, a real number.
For people familiar with the notion, $d$ represents a ``distance'' as defined
by Shanks and Lenstra.

To create such a form we write

\centerline{\tt qfr(a,b,c,d)}

({\tt r} stands for real).

Again these forms can be multiplied, divided, powered, and they can be reduced
using the function {\tt redreal}. This function is in fact a succession of
elementary reduction steps corresponding essentially to a continued fraction
expansion, and a single one of these steps can be achieved using the function
{\tt rhoreal}. Since handling the fourth component $d$ usually
involves computing logarithms, the functions {\tt redrealnod} and 
{\tt rhorealnod} are also provided ({\tt nod} meaning no distance) so that the
fourth component is ignored. Finally, it is sometimes useful to operate on
forms of positive discriminant without performing any reduction (this is
useless in the negative case), the functions {\tt comprealraw} and 
{\tt powrealraw} are available which avoid any reduction.

Again, the function {\tt pf} gives a prime form, but the form which is given
corresponds to an ideal of prime norm which is usually not reduced. If
desired, it should be reduced with one of the function {\tt redreal} or
{\tt redrealnod}.

Finally, you have at your disposal the function {\tt classno} which gives
the class number (this time guaranteed correct), {\tt regula} which gives
the regulator, and the much more sophisticated {\tt buchreal} giving both
the class group structure and its generators, and the regulator, even for
large discriminants.
The {\tt classno} and {\tt regula} program use a $O(\sqrt D)$ algorithm
hence become very slow for discriminants having more than 10 digits. 
{\tt buchreal} can be used on a much larger range.

Let us see examples of all this at work and learn some little known
number theory at the same time. First of all, type {\tt d=3$*$3299;classno(d)}.
We see that the class number is 3 (we know in advance that it must be
divisible by 3 from the {\tt d=$-$3299} case above and Sholtz's theorem).
Let us create a form
by typing {\tt f=redrealnod(pf(d,2),isqrt(d))}. This gives us a prime ideal
of norm equal to 2. Is this ideal principal ? Well, one way to check this,
which is not the most efficient but will suffice for now, is to look at
the complete cycle of reduced forms equivalent to {\tt f}, using the
function {\tt rhoreal}. Type {\tt g=f;for(j=1,20,g=rhoreal(g);print(g))}.
We see that we come back to the form {\tt f} without having the principal
form (starting with $\pm1$) in the cycle, so the ideal corresponding to
{\tt f} is not principal. 

Since the class number is equal to 3, we know however that {\tt f\^{}3} will be
a principal ideal $\alpha\Z_K$. How do we find $\alpha$? For this, type
{\tt f3=powrealraw(f,3)}. This computes the cube of {\tt f}, without reducing
it. Hence it corresponds to an ideal of norm equal to $8=2^3$, so we already
know that the norm of $\alpha$ is equal to $\pm8$. We need more information,
and this will be given by the fourth component of the form. Apply {\tt rhoreal}
to {\tt f3} until you reach the unit form (you must do this 6 times).

Extract the Archimedean component by typing {\tt c=compo(\%,4)}. By
definition of this component, we know that 
$${\alpha\over{\sigma(\alpha)}}=\pm e^{2c},$$
where $\sigma$ denotes real conjugation in our quadratic field. Thus, if we type 

\centerline{\tt a=sqrt(8$*$exp(2$*$c))}

and then {\tt sa=8/a}, we know that up to sign, {\tt a} and {\tt sa} are
numerical approximations of $\alpha$ and $\sigma(\alpha)$. Of course, $\alpha$
can always be chosen to be positive, and a quick numerical check shows that
the difference of {\tt a} and {\tt sa} is close to an integer, and not the sum,
so that in fact the norm of $\alpha$ is equal to $-8$ and the numerical 
approximation to $\sigma(\alpha)$ is {\tt $-$sa}. Thus we type

\centerline{\tt p=x\^{}2$-$round(a$-$sa)$*$x$-$8}

and this is the characteristic polynomial of $\alpha$. We can check that
the discriminant of this polynomial is a square multiple of {\tt d}, so
$\alpha$ is indeed in our field. More precisely, solving in $\alpha$
and using the numerical approximation that we have to solve the sign 
ambiguity in the square root, we get explicitly $\alpha=(15221+153\sqrt d)/2$.
Note that this can also be done automatically using the functions {\tt polred2}
and {\tt modreverse}, as we will see later in the general number field case,
or by solving a system of 2 linear equations in 2 variables.

Exercise: now that we have $\alpha$ explicitly, check that it is indeed a 
generator of the ideal corresponding to the form {\tt f3}.
\medskip
Let us now play a little with cycles. Type

\centerline{\tt buchreal(10000001,0,1,6)}

We get as a result a 5-component vector, which tells us that the class number
is equal to 1, and the regulator is equal to $2641.5\dots$.
Now just as an example, let's assume that we want the regulator to 500
decimals, say. I claim that by simply knowing the crude approximation above,
this can be computed with no effort (of course we could redo the whole 
computation of {\tt buchreal} after having changed the precision, but it
would be a waste of time).

This time, we want to start with the unit form. Since $D=10000001$ is odd,
after having set $D$ equal to that value we can type

\centerline{\tt u=redrealnod(qfr(1,1,(1$-$D)/4,0.),isqrt(D))}

We use the function with no distance since we want the initial distance to be
equal to 0. The function {\tt redrealnod} asks for a second argument which is
the integer part of the square root of the discriminant, obtained using the 
function {\tt isqrt}.

Now we type  {\tt f=rhoreal(u)}. This is the first form encountered along the
principal cycle. For the moment, keep the precision low, for example the
initial default precision. The distance from the identity of {\tt f} is
around 4.253. Very crudely, since we want a distance of $2641.5$, this should
be encountered approximately at $2641.5/4.253=621$ times the distance of 
{\tt f}. Hence, as a first try, we type {\tt f\^{}621}. Oops, we overshot,
since the distance is now $3173.02$. Now we can refine our initial estimate
and believe that we should be close to the correct distance if we raise {\tt f}
to the power $621*2641.5/3173$ which is close to $517$. Now if we compute
{\tt f\^{}517} we hit the principal form right on the dot. Note that this is
not a lucky accident: we will always land extremely close to the correct target
using this method, and usually at most one {\tt rhoreal} correction step is
necessary. Of course, only the distance component can tell us where we are
along the cycle.

Up to now, we have only worked to low precision. The goal was to obtain this
unknown integer $517$. Note that this number has absolutely no mathematical
significance: indeed the notion of reduction of a form with positive 
discriminant is not well defined since there are usually many reduced forms
equivalent to a given form. However, when Pari makes its computations, the
specific order and reductions that it performs is dictated entirely by the
coefficients of the quadratic form itself, and not by the distance component,
hence the precision used has no effect.

Hence we now start again by setting the precision to (for example) 500,
we retype the definition of {\tt u} (why is this necessary?), and then
{\tt f=rhoreal(u)} and finally {\tt f\^{}517}. Of course we know in advance
that we land on the unit form, and the fourth component gives us the regulator
to 500 decimal places with no effort at all.

In a similar way, we could obtain the so-called {\it compact representation}
of the fundamental unit itself, or $p$-adic regulators. I leave this as
exercises for the interested reader.
\medskip
\section{Working in General Number Fields}
\medskip
({\sl Note: the sequence of instructions described in this section can be found
as a GP-readable file in the directory {\tt examples} under the name
{\tt tutnf.gp}.})
\smallskip
The situation here is of course more difficult. First of all, we need to find a
reasonable equation defining our field. This is done using a variant of the
{\tt polred} function which initializes a number of things which will be
needed. For example assume that we want to work in the number field $K$ defined
by one of the roots of the equation $x^4+24x^2+585x+1791=0$. We can type

\centerline{\tt T=x\^{}4$+$24$*$x\^{}2$+$585$*$x$+$1791;polred(T)}

and see from 
the third component that the polynomial $x^4-x^3-21x^2+17x+133$ defines the
same field with smaller coefficients. The {\tt polred} function gives a 
(usually) simpler polynomial, and also sometimes some information on 
the existence of subfields. For example in this case, the second component of
{\tt polred} tells us that the field defined by $x^2-x+1=0$, i.e. the field
of cube roots of unity, is a subfield of our number field $K$. 

However to start working in a number field it is preferable to use a variant of
{\tt polred} which gives a single polynomial and stores some useful
information which will be needed later. This is the function {\tt initalg},
which of course stands for initialization of algebraic number field.

Type {\tt nf=initalg(T)}. You get a vector containing 9 components.
The first is the new polynomial which Pari chooses
to use for working in your number field. The second contains the signature
$(r_1,r_2)$ of the field, the third the field discriminant, the seventh
contains an integral basis, etc\dots See Chapter 3 for the complete list.

Call $\theta$ a root of your initial polynomial, and $\alpha$ a root of the
one that {\tt initalg} has found. These are algebraic numbers, and as already
mentioned are represented as polymods. For example, in our special case
$\theta$ is equal to the polymod 

\centerline{\tt mod(x,x\^{}4$+$24$*$x\^{}2$+$585$*$x$+$1791)}

while $\alpha$ is equal to the polymod

\centerline{\tt mod(x,x\^{}4$-$x\^{}3$-$21$*$x\^{}2$+$17$*$x$+$133)}

Here of course we are considering only the algebraic aspect, and hence ignore
completely {\it which} root $\theta$ or $\alpha$ is chosen.

Now probably you may have a number of elements of your number field which
are expressed as polymods with respect to your old polynomial, i.e. as
polynomials in $\theta$. Since we are now going to work with $\alpha$ instead,
it is necessary to convert these numbers to a representation using $\alpha$.
This is done using the function {\tt initalg2}. For example if you type
{\tt initalg2(T)} you obtain as a result the polymod

\centerline{\tt mod(x\^{}2$+$x$-$2,x\^{}4$-$x\^{}3$-$21$*$x\^{}2$+$17$*$x$+$133)}

meaning that $\theta=\alpha^2+\alpha-2$, and hence the conversion from a 
polynomial in $\theta$ to one in $\alpha$ is easy. If we want to do the
reverse, i.e. go back from a representation in $\alpha$ to a representation
in $\theta$, we use the function {\tt modreverse} on the polymod obtained
by {\tt initalg2}. Try it. The result has a big denominator (147), essentially
because our initial polynomial $T$ was too large.
\smallskip
From now on, we completely forget about our initial polynomial $T$, and use
only the polynomial $x^4-x^3-21x^2+17x+133$ defining $\alpha$, and
the components of the
vector {\tt nf} which gives information on our number field $K$. Note that 
even if the polynomial has small coefficients,
it is necessary to use the {\tt initalg} function before being able to do
{\bf ideal} operations in the number field. Of course {\tt initalg} is not
really necessary if only {\bf element} operations are performed.
\smallskip
Let us first perform some operations on elements. Type {\tt A=nf[1]} to get
the new polynomial, and type 

\centerline{\tt u=mod(x\^{}3$-$5$*$x\^{}2$-$8$*$x$+$56,A)$/$7} 

We may for instance compute {\tt u\^{}3}. Try it. Or, we can type {\tt 1/u}.
Better yet, if we want to know the norm from $K$ to $\Q$ of {\tt u},
we type {\tt norm(u)} (what else?). Still better, if we want all the
symmetric functions of {\tt u} and not only the norm, we type
{\tt char(u,x)} (do not forget the {\tt x} variable which tells GP which
variable to use for the characteristic polynomial. Note that this gives
the characteristic polynomial of {\tt u}, and not in general the minimal
polynomial. Exercises: how does one (easily) find the minimal polynomial
from this? Find a simpler expression for {\tt u}.)
\smallskip
OK, now let's work on the field itself. The {\tt initalg} command already gave
us some information. The field is totally complex (signature $[0,2]$), its 
discriminant is $D=18981$ and 
$(1,\alpha,\alpha^2,(\alpha^3+2\alpha^2+6\alpha)/7)$ is an integral basis.
The Galois group of its Galois closure can be obtained by typing
{\tt galois(A)}. The answer ($[8,-1,1]$) shows that it is equal to $D_4$, the 
dihedral group with 8 elements, i.e. the group of symmetries of a square. 
\smallskip
This implies that the field is ``partially Galois'', i.e. that there exists at
least one non-trivial field isomorphism which fixes $K$. To find out which it
is, we use the function {\tt galoisconj}. This uses the LLL algorithm to find 
linear relations. Unfortunately, in the present version \vers{} of Pari this is
quite unstable so the results are not as nice as we would like. Anyway, let's
type {\tt galoisconj(A)}. Oops, we get an error message (in certain versions
you may not, in which case ignore the next 2 sentences). This is due to the 
bad LLL implementation and should disappear someday. In this particular 
instance, the problem disappears if we {\it lower}
the precision to 18 decimals, say. The result tells us that, apart from the
trivial automorphism, the map $\alpha\mapsto(-\alpha^3+5\alpha^2+\alpha-49)/7$
is a field automorphism. Indeed, if we type

\centerline{\tt char(mod(($-$x\^{}3$+$5$*$x\^{}2$+$x$-$49)$/$7,A),x)}

we obtain the polynomial {\tt A} once again.
\smallskip
Before working with ideals, let us note one more thing. The main part of the
work of the {\tt polred} or {\tt initalg} function is to compute an integral 
basis, i.e. a $\Z$-basis of the maximal order $\Z_K$ of $K$. For a large
polynomial, this implies factoring the discriminant of the polynomial, which
is very often out of the question. There are two ways in which the situation 
may be improved. First, it is often the case that the polynomial that one 
considers is of quite a special type, giving some information on the
discriminant. For example, one may know in advance that the discriminant is a
square. Hence we can ``help'' Pari by giving it that information. More
precisely, using the extra information that we have we may be able to factor
the discriminant of the polynomial. We can then use the function 
{\tt factoredbase} or {\tt factoredpolred} (see Chapter 3)
to inform Pari of this factorization. The function {\tt factoredinitalg} 
has not been written, but it would of course be a trivial task.

The second way in which the situation may be improved is that often we do not
need the complete information on the maximal order, but only require that the
order be $p$-maximal for a certain number of primes $p$. This can again be 
done using the function {\tt factoredbase} by cheating, and giving on purpose
an incomplete factorization of the discriminant involving only the primes
which we want. For example coming back to our initial polynomial $T$, 
the discriminant of the polynomial is $3^7\cdot7^6\cdot19\cdot37$. If we only
want the $7$-maximal order, we simply type

\centerline{\tt factoredbase(T,[7,6;1537461,1])}

and the factors of 1537461 will not be looked at.
\medskip
We now would like to work with ideals (and even with ideles) and not only
with elements. An ideal will be represented by a square matrix giving the
Hermite Normal Form of a $\Z$-basis of the ideal expressed on the
integral basis. Note that, except when mentioned otherwise, we always
represent elements on the integral basis and not on a power basis like
the powers of $\alpha$. An idele will be represented as a 2-element
vector, the first element being the HNF of the corresponding ideal,
which summarizes the non-archimedean information, and the second
element is a vector of real and complex numbers with $r_1+r_2$ components,
the first $r_1$ being real, the second complex. In fact, to avoid
certain ambiguities, the first $r_1$ components are allowed to
have an imaginary part which is a multiple of $\pi$. These $r_1+r_2$
components correspond to the Archimedean places of the number field
$K$.

Finally, a prime ideal may either be represented as an ordinary ideal
as explained above, but usually will be represented as a 5-element
vector {\tt [p,a,e,f,b]} as explained in the {\tt primedec} heading of
Chapter 3.
\medskip
Let us keep our number field $K$ above, and hence the vector {\tt nf}.
Type 

\centerline{\tt pd=primedec(nf,7)}

This gives the decomposition of the prime
number 7 into prime ideals. We have chosen 7 because it is the index
of $\Z[\theta]$ in $\Z_K$ (which is given as {\tt nf[4]}) hence is the
most difficult case to treat.

The result is a vector with 4 components, showing that 7 is totally
split in the field $K$ into prime ideals of norm 7. Let us take
one of these ideals, say the first, so type {\tt pr1=pd[1]}.
To work further, we would like to have the Hermite Normal Form of this
ideal. No problem: since ideal multiplication always gives the result
in HNF, we will simply multiply our prime ideal by $\Z_K$, which is
represented by the identity matrix. So type

\centerline{\tt hp=idealmul(nf,idmat(4),pr1)}

and we have the desired HNF. 

We can now perform ideal operations. For example type

\centerline{\tt hp3=idealmul(nf,hp,idealmul(nf,hp,hp))}

or more simply

\centerline{\tt hp3=idealpow(nf,hp,3)}

to get the cube of the ideal {\tt pr1}.
Since the norm of this ideal is equal to $343=7^3$,
to check that it is really the cube of {\tt pr1} and not of other ideals
above 7, type

\centerline{\tt for(j=1,4,print(idealval(nf,hp3,pd[j])))}

and we see that the valuation at {\tt pr1} is equal to 3, while
the others are equal to zero.

Let us now ``idelize'' hp3 by typing {\tt hpi3=[hp3,[0.,0.]]}.
(we need $r_1+r_2=2$ components for the second vector).
Then type {\tt hr1=ideallllred(nf,hpi3,0)}. We get a new ideal which
is equivalent to the first, and the Archimedean component is
non-trivial. Now, just for fun type

\centerline{\tt hr=hr1;for(j=1,3,hr=ideallllred(nf,hr,[1,5]);print(hr))}

We see that the third {\tt hr} is equal to the initial {\tt hr1}. This
means that we have found a unit in our field, and it is easy to
extract this unit given the Archimedean information (clearly it
would be impossible without). We first form the difference of the
Archimedean contributions of {\tt hr} and {\tt hr1} by typing

\centerline{\tt arch=hr[2]$-$hr1[2];l1=arch[1];l2=arch[2];}

From this, we obtain the logarithmic embedding of the unit by typing

\centerline{\tt s=real(l1$+$l2)$/$4;v1=[l1,l2,conj(l1),conj(l2)]\~{}$/$2$-$[s,s,s,s]\~{};}

This 4-component vector contains by definition the logarithms of the
four complex embeddings of the unit. Since the matrix {\tt nf[5][1]}
contains the values of the $r1+r2$ embeddings of the elements of the
integral basis, we can obtain the representation of the unit on the
integral basis by typing

{\tt m1=nf[5][1];m=matrix(4,4,j,k,if(j<=2,m1[j,k],conj(m1[j$-$2,k])));}

{\tt v=exp(v1);au=gauss(m,v);}

{\tt vu=round(real(au))}

Then {\tt vu} is the representation of the unit on the integral basis.
The closeness of the approximation of {\tt au} to {\tt vu} gives us
confidence that we have made no numerical mistake. To be sure that {\tt vu}
represents a unit, type {\tt u=mod(nf[7]$*$vu,A)}, then typing
{\tt norm(u)} we see that it is equal to 1 so {\tt u} is a unit.

There is of course no reason for {\tt u} to be a fundamental unit.
Let us see if it is a square. Type {\tt f1=factor(subst(char(u,x),x,x\^{}2))}.
We see that the characteristic polynomial of {\tt u} where {\tt x}
is replaced by {\tt x\^{}2} stays irreducible, hence {\tt u} is not a square.
However, type 

{\tt f1=factor(subst(char($-$u,x),x,x\^{}2))}. 

The characteristic
polynomial of {\tt $-$u} factors into 2 polynomials of degree 4 
(the second being the same as the first with $x$ replaced by $-x$),
showing that {\tt $-$u} is the square of a unit, the characteristic polynomial
of which is given as the component {\tt f1[1,1]} (or {\tt f1[2,1]} which is
the characteristic polynomial of the opposite unit). 

We now want to find the square root of {\tt $-$u}. We can again use the 
{\tt gauss} function as above. For this we need to take the square
root of each element of the vector {\tt $-$v}, and hence there are
sign ambiguities. Let's do it anyway. Type {\tt v=sqrt($-$v)}. We see that
{\tt v[1]} and {\tt v[3]} are conjugates, as well as {\tt v[2]} and {\tt v[4]},
so for the moment the signs seem OK. Now try {\tt au=gauss(m,v)}. The numbers
obtained are clearly not integers, hence the last remaining sign change must
be performed. Type {\tt v[1]=$-$v[1];v[3]=$-$v[3]} (they must stay conjugate)
and then again {\tt au=gauss(m,v)}. This time the components are close to 
integers, so we are done (after typing {\tt vu=round(real(au))} as before).

Anyway, we find that a square root {\tt u2} of {\tt $-$u} is represented by the 
vector {\tt vu=[$-$4,1,1,$-$1]\~{}} on the integral basis, and this is in fact a
fundamental unit. If we want it as an algebraic number, we simply type
{\tt u2=mod(nf[7]$*$vu,A)}.
\medskip
The use of the function {\tt polred2} gives us another method to find
{\tt u2} as follows. We type {\tt q=polred2(f1[1,1])}. We
recognize the polynomial {\tt A} as the component {\tt q[3,2]}. To 
obtain the square root of our unit we then simply type
{\tt up2=modreverse(mod(q[3,1],f1[1,1]))} (Exercise: why?). We find that
{\tt up2} is represented by the vector {\tt [$-$3,$-$1,0,0]\~{}}
on the integral basis, which is not the result that we have found before
nor its opposite. Where is the error ? (Please think about this before
reading on. There is a mathematical subtlety hidden here.)

Have you solved the problem? Good for you! Yes, the problem occurs because
as mentioned before (but you may not have noticed since it is not stressed
in standard textbooks) although the number field $K$ is not Galois over
$\Q$, there does exist non-trivial automorphisms, and we have found it
above by using the function {\tt galoisconj}. Indeed, if we apply this
automorphism to {\tt up2} (by typing {\tt mod(subst(lift(up2),x,aut),A)}
where {\tt aut} is the nontrivial component of {\tt galoisconj(A)}
computed above), we find the opposite of {\tt u2}, which is OK.
\smallskip
Still another method which avoids all sign ambiguities and automorphism
problems is as follows. Type {\tt r=f1[1,1]\%(x\^{}2$+$u)} to find the remainder
of the characteristic polynomial of {\tt u2} divided by {\tt x\^{}2$+$u}.
This will be a polynomial of degree 1 in {\tt x} (with polymod coefficients)
and we know that {\tt u2}, being a root of both polynomials will be the
root of {\tt r}, hence can be obtained by typing 
{\tt u2=$-$coeff(r,0)$/$coeff(r,1)}. Indeed, we immediately find the
correct result with no trial and error.
\smallskip
Now ideals can be used in a wide variety of formats. We have already seen
HNF-representations, ideles and prime ideals. We can also use algebraic
numbers. For example type {\tt al=mod(x\^{}2$-$9,A)}, then 
{\tt principalidele(nf,al)}. We obtain the idele corresponding to {\tt al}
(see the manual for the exact description). However it is usually not necessary
to compute this explicitly since this is done automatically inside Pari.
For example, you can type

\centerline{\tt for(j=1,4,print(j,": ",idealval(nf,al,pd[j])))}

We see that the valuation is non-zero (equal to 1) at the prime ideals 
{\tt pd[1]} and {\tt pd[4]}. In addition, typing {\tt norm(al)} shows that
{\tt al} is of norm $49=7^2$. Let's check this differently. Type 
{\tt pd14=idealmul(nf,pd[1],pd[4])} and then {\tt idealmul(nf,al,idmat(4))}.
We see that the results are the same, hence the product of the two prime
ideals {\tt pd[1]} and {\tt pd[4]} is equal to the principal ideal generated
by {\tt al}. There is still something to be done with this examples as we
shall see below with the function {\tt isprincipal}.
Essentially all functions that you would want on ideals are available:

{\tt idealadd}, {\tt idealdiv}, {\tt idealfactor}, {\tt idealintersect},
{\tt idealinv}, {\tt ideallllred}, {\tt idealmul}, {\tt idealmulred},
{\tt idealpow}, {\tt idealpowred}, {\tt idealval}, {\tt principalideal}, 
{\tt principalidele}. We will also see {\tt isprincipal} below (don't try
it now because you need a much larger {\tt initalg}).

We suggest you play with some of these to get the feel of the algebraic
number theory package. Remember simply that when a matrix (usually in Hermite
normal form) is output, it is always a $\Z$-basis of the result expressed on
the {\it integral basis} {\tt nf[7]} of the number field, {\it not} on
the power basis, except of course if it is an integral basis.
\medskip
Apart from the above functions you have at your disposal the very powerful
function {\tt buchgenfu} which is of the same type as {\tt buchimag}
and {\tt buchreal} seen above, but for general number fields, and hence
much slower. See Chapter 3 for a detailed explanation of its use.

Type for example {\tt m=buchgenfu(A,0.2,0.2)}, where {\tt A} is the same 
polynomial as before. After some work, we get a matrix with one column, which
we extract by typing {\tt v=m[,1];} Then $v$ is a vector with 10 components.
We immediately recognize the first component as being the polynomial {\tt A}
once again. The really new things start as the fifth component. We see that the
class number $v[5][1]$ is equal to 4, that the class group $v[5][2]$ is cyclic,
(and $v[5][3]$ gives us a generator in HNF), that the regulator $v[6]$ is 
equal to $3.794\dots$. Then $v[7]$ which is close to 1 tells us that the 
results should be correct under a lot of assumptions (i.e. the class group or
the regulator are not smaller). $v[8][1]=6$ tells us that the roots of unity in
the number field are exactly the sixth roots of 1 and $v[8][2]$ gives us a 
primitive sixth root of 1 as a polynomial in $\alpha$, and finally the 
component $v[9][1]$ gives us a fundamental unit as a polynomial in $\alpha$.
This unit must be linked to the unit {\tt u2} found above since the unit
rank is equal to 1. 

To find the relation, let us use the logarithmic embeddings. Type
{\tt uf=mod(v[9][1],A)}, {\tt uu=mod(v[8][2],A)} to get the generators of the 
unit group, then 

{\tt cu2=log(conjvec(u2));cuf=log(conjvec(uf));cuu=log(conjvec(uu));}

to get the (complex) logarithmic embeddings. Then type

{\tt lindep(real([cu2[1],cuf[1],cuu[1]]))} 

to find a linear dependence. Unfortunately, the result $[0,0,1]$ is to be
expected since {\tt uu} is a root of unity, hence the components of {\tt cuu}
are pure imaginary. Hence we must not take the real part. However, in that
case the logs are defined only up to a multiple of $2i\pi$. Hence we type

{\tt lindep([cu2[1],cuf[1],cuu[1],2$*$i$*$pi])} 

The result $[1,-1,-3,0]$ shows that {\tt u2$/$uf} is probably equal to 
{\tt uu\^{3}}, i.e. to $-1$ (since {\tt uu} is of order 6).
Indeed, typing {\tt u2$/$uf} results in {\tt $-1$} (modulo {\tt A}).

Note that this works, but is not really a good way to obtain the result.
A better way will be seen below using the function {\tt isunit}.

Warning: since the fundamental unit obtained depends on the
random seed, you may obtain a different unit from the one given here. To be
sure to obtain the same result, type {\tt setrand(1)} before the 
{\tt buchgenfu} instruction.

You may ignore the last component of {\tt v}.

You will notice that the manual explains that the units are given as 
polynomials in $\alpha$. You may prefer them as vectors of coordinates on
the integral basis. You can go back and forth between those two representations
using the components {\tt nf[7]} and {\tt nf[8]} which were given by
{\tt initalg} above. For example, to obtain the representation on the integral
basis of the primitive sixth root of unity $v[8][2]$, type
{\tt ru=nf[8]$*$vvector(4,j,coeff(v[8][2],j$-$1))}. Conversely, to come back to
the initial representation type {\tt nf[7]$*$ru}.
\smallskip
The {\tt buchgenfu} program is in fact even more powerful. Instead of using
{\tt buchgenfu}, we are going to use a variant, which is going to give us
even more information. Type 

{\tt setrand(1);bnf=buchinitfu(A,0.2,0.2);}.

This performs exactly the same computations as {\tt buchgenfu}, but keeps much
more information. In particular, you don't {\it want} to see the result, 
whence the semicolon.

To get most of the {\tt buchgenfu} components, type {\tt bnf[8]}. This gives
you the class group and generators, the regulator, the check component
and the unit group. If you don't want the unit group explicitly (it is 
implicitly contained in the rest of {\tt bnf}), instead of using
{\tt buchinitfu} use {\tt buchinit}. Similarly, if you use {\tt buchgenfu}
and don't want the units, use instead {\tt buchgen}.

To initialize the number field, no need to use {\tt initalg}: simply type
{\tt nf=bnf[7]}. Of course in our case this had already been done.
We are now ready to perform more sophisticated operations in the class group.
For example, let us take again our prime ideal {\tt hp} above 7. Since we
know that the class group is of order 4, we deduce that {\tt hp} raised to
the fourth power must be principal. Type {\tt hp4=idealpow(nf,hp,4)} then
{\tt vis=isprincipal(bnf,hp4)}. The function {\tt isprincipal} now uses all the
information contained in {\tt buchinit} and tells us that indeed {\tt hp4}
is principal. More precisely, the first component, $[0]$, gives us the position
of the ideal in the class group. Here, $[0]$ means that it is up to equivalence
equal to the 0-th power of the generator given by {\tt buchinit} or 
{\tt buchgen}, in other words that it is a principal ideal.

The second component gives us the algebraic number $\alpha$ such that
{\tt hp4$=\alpha\Z_K$}, where $\Z_K$ is the ring of integers of our number
field, $\alpha$ being expressed on the integral basis nf[7]. To get $\alpha$
as an algebraic number, we thus type {\tt alpha=mod(nf[7]$*$vis[2],A)}.

Let us check that the result is correct: first, type {\tt norm(alpha)}.
It is indeed equal to $7^4=2401$, which is the norm of hp4 (it could also
have been equal to $-2401$). This is only a first check. The complete check
is obtained by computing the HNF of the principal ideal generated by 
{\tt alpha}. To do this, type {\tt idealmul(nf,idmat(4),mat(vis[2]))}.
Here, {\tt idmat(4)} represents the ideal $\Z_K$ itself, and we must convert
{\tt vis[2]} into a one-column matrix using the function {\tt mat}, otherwise
{\tt idealmul} will complain. The result that we obtain is identical to
{\tt hp4}, thus showing that {\tt alpha} is correct (not that there was any
doubt!). You may ignore the third component of {\tt vis} which just tells us
that the accuracy to which {\tt alpha} was computed was largely sufficient.

But {\tt isprincipal} also gives us information for nonprincipal ideals.
For example, type 

{\tt vit=isprincipal(bnf,hp)}.

The component {\tt vit[1]}
is now equal to $[1]$, and tells us that {\tt hp} is ideal-equivalent to
the first power of the generator {\tt g} given by {\tt buchinit} or {\tt buchgen}.
The second component {\tt vit[2]} gives us $\alpha$ on the integral basis
such that {\tt hp$=\alpha g^1$}.
\smallskip
Let us now take the example of the principal ideal {\tt pd14} that we have 
seen above. We know that {\tt pd14} is principal, but of course we have
forgotten the generator that we had found, so we want one. For this,
we type {\tt pp=isprincipal(bnf,pd14)}. The first component of the result is
{\tt [0]}, telling us that the ideal is indeed principal, as we knew already. 
But the second component gives us the components of a generator of {\tt pd14}
on the integral basis. To get it on the power basis, type
{\tt al2=mod(nf[7]$*$pp[2],A)}. Now we remember suddenly that we already
had a generator {\tt al} for the same ideal. Hence type {\tt u3=al2/al}.
This must be a unit. To check this, computing its norm would not be
sufficient, we need the complete characteristic polynomial, so we type
{\tt char(u3,x)}, and indeed it has integral coefficients with constant term 1
so {\tt u3} is a unit. You can again find out what unit it is as we did above.
However, as we mentioned, this is not really the best method.

To find the unit, we use explicitly the third component of the vector 
{\tt bnf} given by {\tt buchinit}. This contains an $(r+1)\times r$ complex
matrix whose columns represent the complex logarithmic embedding of the
fundamental units. Here $r=r_1+r_2-1$ is the unit rank. We first compute
the component of {\tt u3} on the torsion-free part of the group of units
by proceeding as follows. Type

{\tt me=concat(bnf[3],[2,2]~)}. Indeed, this is a variant of the regulator
matrix and is much more practical to use since it is more symmetric and
avoids suppressing one row arbitrarily. Now type

{\tt cu3=principalidele(nf,u3)[2]} 

to get the complex logarithmic embedding of {\tt u3}.
We could of course also have computed this logarithmic embeddings directly 
using {\tt conjvec} as we did above, but then we must take care of the factors
of 1 and 2 occurring.

Then type {\tt xc=gauss(real(me),real(cu3))}

Whatever field we are in, if {\tt u3} is a unit this {\it must} end with a 0
because of the ``spurious'' vector $[2,2]$, and the other components (here
only one) give the exponents on the fundamental units. Here the only other
component is the first, with a coefficient of $1$, so we know that {\tt u}
is equal to {\tt uf} multiplied by a root of unity.

To find this root of unity, we type {\tt xd=cu3$-$me$*$xc} then
{\tt xu=principalidele(nf,uu)[2]} and finally {\tt xd[1]/xu[1]}. We find
$-3$ as a result, so finally our unit {\tt u3} must be equal to 
{\tt $uu$\^{}3*$uf=$-$uf}, which is clear.

Although the above operations are easily performed, they are a little tedious,
hence Pari has incorporated the above in a function called {\tt isunit}.
So instead of doing the above, you simply type {\tt isunit(bnf,u3)}
and the result $[1,\mod(3,6)]$ tells us that {\tt u3} is equal to
{\tt uf\^{}(1).uu\^{}mod(3,6)}, i.e. to {\tt $uu$\^{}3*$uf} as before, since
{\tt uu} is of order 6.
\medskip
\section{Using Transcendental Functions}
\medskip
All the elementary transcendental functions and several higher transcendental
functions (gam\-ma function, incomplete gamma function, error function, 
exponential integral, K-bessel functions, confluent hypergeometric functions, 
Riemann zeta function, polylogarithms, Weber functions, theta functions) are
provided. More could be written if the need arose.

In this type of functions, the default precision plays an essential role.
In almost all cases transcendental functions work in the following way.
If the argument is exact, the result will be computed using the current
default precision. If the argument is not exact, the precision of the
argument is used for the computation. A note of warning however: even in this
case the {\it printed\/} value will be the current real format (usually the
same as the default precision). 

Let's assume that we have 28 decimals of default precision (this is what
we get automatically at the start of a GP session). Type {\tt e=exp(1)}.
We get the number $e=2.718\dots$ to 28 decimals. Let us check how many
correct decimals we really have. The hard (but reasonable) way is to proceed
as follows. Change the precision to a substantially higher value, for example
by typing {\tt\bs precision=50}. Then type {\tt e}, then {\tt exp(1)} once
again. This last value is the correct value of the mathematical constant $e$
to 50 decimals, while the variable {\tt e} shows the value that was
computed to 28 decimals. Clearly they coincide to exactly 29 significant
digits.

A simpler way to see how many decimals we have is to type {\tt\bs x} after
the initial {\tt e=exp(1)}. The {\tt 02010005} which is the first codeword
of the real number indicates that we are dealing with a real number having
exactly $3=5-2$ mantissa words, and since $3\ln(2^32)/\ln(10\approx28.9$
we see that we have 28 or 29 significant digits.
\smallskip
Come back to 28 decimals ({\tt\bs precision=28}). If we type {\tt exp(1.)}
you can check that we also obtain 28 decimals. However, type 
{\tt f=exp(1$+$10.\^{}($-$30))}. Although the default precision is still 28,
you can check using the one of the two methods above that we have in fact 59
significant digits! The reason is that {\tt 1$+$10.\^{}($-$30)} is computed
according to the Pari philosophy, i.e.~to the best possible precision. Since
$10.$ has 29 significant digits and 1 has ``infinite'' precision, the number
{\tt 1$+$10.\^{}($-$30)} will have $59=29+30$ significant digits, hence {\tt f}
also.

Now type {\tt cos(10.\^{}($-$15))}. The result is printed as $1.0000\dots$,
but is of course not exactly equal to 1. Using {\tt\bs x}, we see that the
result has 7 mantissa words, giving us the possibility of having 67 
correct significant digits. In fact (look in precision 100), only 60 are 
correct. Pari gives you as much as it can, and since 6 mantissa words
would have given you only 57 digits, it uses 7. But why does it give so
precise a result? Well, it is the same reason as before. When $x$ is close
to 1, $\cos(x)$ is close to $1-x^2/2$, hence the precision is going to be
approximately the same as this quantity, which will be $1-0.5*10^{-30}$ where
$0.5*10^{-30}$ is considered with 28 significant digit accuracy, hence the
result will have approximately $28+30=58$ significant digits.

Unfortunately, this philosophy cannot go too far. For example, when you
type {\tt cos(0)}, GP should give you exactly 1. Since it is reasonable for
a program to assume that a transcendental function never gives you an exact
result, GP gives you $1.000\dots$ to one more mantissa word than the current
precision.
\medskip
OK, now let's see some more transcendental functions at work. Type
{\tt gamma(10)}. No problem (type {\tt 9!} to check). Type {\tt gamma(100)}.
The number now ends with a $*$, indicating that the printed integral part
may be incorrect (in fact it is, type {\tt 99!} to get the correct answer).
To get the correct value, there are two solutions. The first and most natural
one is to increase the precision. Since {\tt gamma(100)} has 157 decimal
digits, type {\tt\bs precision=170} (to be on the safe side), then
{\tt gamma(100)} once again. After some work, the printed result is this time
perfectly correct.

However, this is probably not the best way to proceed. Come back first to
the default precision (type {\tt\bs precision=28}). The gamma function 
increasing very rapidly, one usually uses its logarithm. Type 
{\tt lngamma(100)}. This time the result has a reasonable size, and is
exactly equal to {\tt log(99!)}.

Try {\tt gamma(1/2$+$10$*$i)}. No problem, we have the complex gamma function.
Now type 

{\tt t=1000;z=gamma(1$+$i$*$t)$*$t\^{}($-$1/2)$*$exp(pi/2$*$t)/sqrt(2$*$pi)},

then {\tt norm(z)}. We see that {\tt norm(z)} is very close to 1, in accordance
with the complex Stirling formula.
\smallskip
Let's play now with the Riemann zeta function. First put the timer on (type
{\tt \#}). Type {\tt zeta(2)$*$6/pi\^{}2}. This seems correct. Type 
{\tt zeta(3)}. All this takes essentially no time at all. However, type
{\tt zeta(3.)}. Although the result is the same, you will notice that the
time is substantially larger (if your machine is too fast to see the 
difference, increase the precision!). This is because Pari uses special 
formulas to compute {\tt zeta(n)} when {\tt n} is a (positive or negative)
integer, even when {\tt n} is odd and positive.

Type {\tt zeta(1$+$i)}. This also works. Now for fun, let us compute in a very
na\"\i ve way the first complex zero of {\tt zeta}. We know that it is of the
form {\tt 1/2$+$i$*$t} with between 14 and 15. Thus, we can use the following
series of instructions. Type:

\centerline{\tt t1=1/2$+$14$*$i;t2=1/2$+$15$*$i;z1=zeta(t1);z2=zeta(t2);}

\centerline{\tt if(norm(z2)<norm(z1),t3=t1;t1=t2;t2=t3;z1=z2,);}

\centerline{\tt while(norm(z2)>10\^{}$-$50,t2=(t1$+$t2)/2.;z2=zeta(t2);\bs}

\centerline{\tt if(norm(z2)<norm(z1),t3=t1;t1=t2;t2=t3;z1=z2,);print(t1,": ",z1))}

We thus obtain the first zero to 25 significant digits. 
\medskip
As mentioned at the beginning of this tutorial, some transcendental functions
can also be applied to $p$-adic numbers. This is as good a time as any to
familiarize yourself with them. Type {\tt a=exp(7$+$O(7\^{}10))}, then
{\tt log(a)}. All seems in order. Now type {\tt b=log(5$+$O(7\^{}10))}, then
{\tt exp(b)}. Is something wrong, since we don't recover the number we
started with? Absolutely not. Type {\tt exp(b)$*$teich(5$+$O(7\^{}10))}, and
we indeed recover our initial number {\tt 5$+$O(7\^{}10)}. The function
{\tt teich(x)} is the Teichm\"uller character, and is characterized by the
fact that it is the unique $p-1$-st root of unity (here with $p=7$) which
is congruent to {\tt x} modulo $p$, assuming that {\tt x} is a $p$-adic
unit. 
\smallskip
Let us come back to real numbers for the moment. Type {\tt agm(1,sqrt(2))}. 
This gives the arithmetic-geometric mean of 1 
and $\sqrt2$, and is the basic method for computing (complete) elliptic
integrals. In fact, type 

{\tt pi/2/intnum(t=0,pi/2,1/sqrt(1$+$sin(t)\^{}2))},

and the result is the same. The elementary transformation {\tt x=sin(t)}
gives the mathematical equality
$$\int_0^1\dfrac{dx}{\sqrt{1-x^4}}=\dfrac{\pi}{2\text{agm}(1,\sqrt2)}\enspace,$$
which was one of Gauss's remarkable discoveries in his youth.

Now type {\tt 2$*$agm(1,i)/(1$+$i)}. As you see, the complex AGM also works, 
although one must be careful with its definition. The result found is
almost identical to the previous one. Exercise: do you see why?

Finally, type {\tt agm(1,1$+$7$+$O(7\^{}10))}. As you see, $p$-adic AGM also
work. Note however that since the square root of a $p$-adic number is not
represented in Pari, only certain $p$-adic AGM can be computed. In addition,
when $p=2$, the congruence restriction is that {\tt agm(a,b)} can be computed
only when {\tt a/b} is congruent to 1 modulo $16$ (not 8 as could be expected).
\smallskip
Now type {\tt ?3}. This gives you the list of all transcendental functions.
Instead of continuing with more examples, we suggest that you experiment
yourself with the list of functions. In each case, try integer, real, complex
and $p$-adic arguments. You will notice that some have not been implemented
(or do not have a reasonable definition).
\medskip
\section{Using Numerical Tools}
\medskip
Although not written to be a numerical analysis package, Pari can nonetheless
perform some numerical computations. We leave for a subsequent section
linear algebra and polynomial computations. A few of the functions are even
quite fun and surprising.

You of course know the formula $\pi=4(1-\dfrac13+\dfrac15-\dfrac17+\cdots)$
which is deduced from the power series expansion of {\tt atan(x)}. You also
know that $\pi$ cannot be computed from this formula, since the convergence
is so slow. Right? Wrong! Type {\tt\bs precision=100} (just to make it
more astonishing), then {\tt 4$*$sumalt(k=0,($-$1)\^{}k/(2$*$k$+$1))}. In a
few seconds (admittedly more than simply typing {\tt pi}), we get $\pi$ to
100 significant digits (type {\tt pi} to check). This is due to the combination
of a method due to Euler for accelerating alternating sums, and a programming
trick due to the Dutch mathematician van Wijngarten (see Numerical Recipes
for an explanation).

Similarly, type {\tt\bs precision=28} (otherwise the time will be too long)
and {\tt sumpos(k=1,1/k\^{}2)}. Although once again the convergence is
slow, a similar trick allows to compute the sum when the terms are positive
(compare with the exact result {\tt pi\^{}2/6}. This is much less impressive
because it is much slower, but is still very useful.

Even better, {\tt sumalt} can be used to sum divergent series! Type

\centerline{\tt f(s)=sumalt(k=1,($-$1)\^{}(k$-$1)/k\^{}s)/(1$-$2\^{}(1$-$s))}

(this example has already been given in Chapter 3).

Then for positive values of {\tt s} different from 1, {\tt f(s)} is equal
to {\tt zeta(s)} and the series cvonverges, albeit slowly (sumalt doesn't
care however). For negative {\tt s}, the series diverges, but {\tt f(s)}
still gives the correct result! Try {\tt f($-$1)}, {\tt f($-$2)}, {\tt f($-$1.5)},
and compare with the corresponding values of {\tt zeta}.
You should not push the game too far: {\tt f($-$2.5)} will probably not work.

Try {\tt f(i)}, and compare with {\tt zeta(i)}. Even (some) complex values 
work, even though the sum is not alternating any more!

Similarly, try {\tt sumalt(n=1,($-1$)\^{}n/(n+i))}.
\medskip
More traditional functions are the numerical integration functions.
For example, type {\tt intnum(t=1,2,1/t)} and presto! you get 26 decimals
of {\tt log(2)}. Look at Chapter 3 to see the different integration functions
which are available.

With Pari, however, you can go further since complex types are allowed.
For example, assume that we want to know the location of the zeros of the
function $h(z)=e^z-z$. We use Cauchy's theorem, which tells us that the
number of zeros in a disk of radius $r$ centered around the origin is
equal to
$$\dfrac{1}{2i\pi}\int_{C_r}\dfrac{h'(z)}{h(z)}\,dz\enspace,$$
where $C_r$ is the circle of radius $r$ centered at the origin.
Hence type 

\centerline{\tt f(z,u)=u=exp(z);(u$-$1)/(u$-$z)}

\centerline{\tt g(r)=intnum(t=0,2$*$pi,f(r$*$exp(i$*$t))$*$exp(i$*$t))$*$r/(2$*$pi)}

The function {\tt g(r)} will then count the number of zeros (we have simply
made the change of variable {\tt z=r$*$exp(i$*$t)}).

Now type {\tt\bs precision=9} (otherwise the computation would take too long,
and anyway we know that the result is an integer), then {\tt g(1)},
{\tt g(1.5)}. The results tells us that there are no zeros inside the unit
disk, but that there are two (necessarily complex conjugate) whose modulus
is between $1$ and $1.5$. For the sake of completeness, let us compute them.
Let $z$ be such a zero, and write $z=x+iy$ with $x$ and $y$ real. Then
the equation $e^z-z=0$ implies, after elementary transformations, that
$e^{2x}=x^2+y^2$ and that $e^x\cos(y)=x$. Hence $y=\pm\sqrt{e^{2x}-x^2}$
and hence $e^x\cos(\sqrt{e^{2x}-x^2})=x$. Therefore, type

\centerline{\tt s(x,u)=u=exp(x);u$*$cos(sqrt(u$*$u$-$x$*$x))$-$x}

Then {\tt s(0)} is positive while {\tt s(1)} is negative. Come back to
precision 28 ({\tt\bs precision=28}) and type

\centerline{\tt x0=solve(x=0,1,s(x))}

This quickly gives us the value of {\tt x}, and we can then type
{\tt z=x0$+$i$*$sqrt(exp(2$*$x0)$-$x0\^{}2)}. This (together with its complex
conjugate) is the required zero. As a check, type {\tt exp(z)$-$z}, and also
{\tt abs(z)}.

Of course you can integrate over contours which are more complicated than
circles, but you must perform yourself the changes of variable as we have
done above to reduce the integral to a number of integrals on line segments.
\smallskip
The example above also shows the use of the {\tt solve} function. To use
{\tt solve} on functions of a complex variable, it is necessary to reduce
the problem to a real one. For example, to find the first complex zero
of the Riemann zeta function as above, we could try typing

{\tt solve(t=14,15,real(zeta(1/2$+$i$*$t)))},

but this would not work because
the real part is positive for {\tt t=14} and {\tt t=15}. As it happens, the
imaginary part works. Type 

{\tt solve(t=14,15,imag(zeta(1/2$+$i$*$t)))},

and this now works. We could also narrow the search interval and type for
instance 

{\tt solve(t=14,14.2,real(zeta(1/2$+$i$*$t)))} 

which would also work.
\medskip
\section{Working with Elliptic Curves}
\medskip
We have seen that the first thing to do to work in a number field is to use
the function {\tt initalg}. Well, similarly, the first thing to do to work
with an elliptic curve is to use the function {\tt initell}. So type
{\tt e=initell([6,$-3$,9,$-16$,$-14$])}. This computes a number of things
about the elliptic curve defined by the affine equation
$$y^2+6xy+9y=x^3-3x^2-16x-14\enspace.$$ We see that the discriminant
e[12] is equal to 37, hence the conductor of the curve is 37. Of course in 
general it is not so trivial. In fact, the equation of the curve is clearly
not minimal, so type {\tt r=globalred(e)}. The first component r[1] tells
us that the conductor is 37 as we already knew. The second component is a
4-component vector which will allow us to get the minimal equation:
simply type {\tt e=chell(e,r[2])} and the new {\tt e} is now our minimal
equation with corresponding data. You can for the moment ignore the third
component r[3].

The new {\tt e} tells us that the minimal equation is $y^2+y=x^3-x$.
Let us now play a little with points on {\tt e}. Type {\tt q=[0,0]}, which is
clearly on the curve (type {\tt isoncurve(e,q)} to check). Well, {\tt q} may
be a torsion point. Type {\tt hell(e,q)}, which computes the canonical
Neron-Tate height of q. This is non-zero, hence {\tt q} is not torsion. To see
this even better, type 

{\tt for(k=1,20,print(powell(e,q,k)))},

and we see the
characteristic parabolic explosion of the size of the points. As a further
check, type {\tt hell(e,powell(e,q,20))/hell(e,q)}. We indeed find $400=20^2$
as it should be. You can also type {\tt orderell(e,q)} which returns 0,
telling you that {\tt q} is non-torsion.
\smallskip
Let's try another curve. Type {\tt e=initell([0,$-$1,1,0,0])}. This corresponds
to the equation $y^2+y=x^3-x^2$. Again from the discriminant e[12] we see that
the conductor is equal to 11, and if you type {\tt globalred(e)} you will see
that the equation  for {\tt e} is minimal. Type {\tt q=[0,0]} which is clearly
a point on {\tt e}, and {\tt hell(e,q)}. This time we obtain a value which is
very close to zero, hence {\tt q} must be a torsion point. Indeed, typing
{\tt for(k=1,5,print(powell(e,q,k)))} we see that {\tt q} is a point of order
5 (note that the point at infinity is represented as [0]). More simply,
you can type {\tt orderell(e,q)}.
\smallskip
Let's try still another curve. Type {\tt e=initell([0,0,1,$-7$,6])} to get
the curve $y^2+y=x^3-7x+6$. Typing {\tt globalred(e)} shows that this is a
minimal equation and that the conductor, equal to the discriminant, is 5077.
There are some trivial integral points on this curve, but let's try to be
more systematic. 

First, let's study the torsion points. Typing {\tt torsell(e)} shows that
the torsion subgroup is trivial, so we won't have to worry about torsion
points. Next, the component e[14] gives us the 3 roots of the minimal
equation over $\Bbb C$, i.e. $Y^2=X^3-7X+25/4$ (set $(X,Y)=(x,y+1/2)$)
so if $(x,y)$ is a real point on the curve, $x$ must be at least equal to the 
smallest root, i.e. $x\ge-3$. Finally, if $(x,y)$ is on the curve, its
opposite is clearly $(x,-y-1)$. So we are going to use the {\tt ordell} 
function ({\it not} the orderell function, don't confuse the two) and type

\centerline{\tt v=[];for(x=$-$3,1000,s=ordell(e,x);if(length(s),v=concat(v,[[x,s[1]]]),));v}

We thus get a large number (18) of integral points. Together with their 
opposites and the point at infinity, this makes a total of 37 integral points,
which is large for a curve having such a small conductor. So we suspect
(if we don't know already, since this curve is quite famous!) that the rank
of this curve must be high. Let's try and put some order in this (note that
we work only with the integral points, but in general rational points should
also be considered).

Type {\tt hv=hell(e,v)}. This gives the vector of canonical heights. Let us 
order the points according to their height. For this, type 

{\tt iv=indsort(hv)}, then {\tt hv=extract(hv,iv)} and {\tt v=extract(v,iv)}.

It seems reasonable to take the numbers with smallest height as generators
of the Mordell-Weil group. Let's try the first 4: type

{\tt m=mathell(e,extract(v,[1,2,3,4]));det(m)}

Since the curve has no torsion, The fact that the determinant is close to 
zero shows that the first four points are dependent. To find the dependency,
it is enough to find the kernel of the matrix {\tt m}. Since its determinant
is not exactly equal to zero, it will be empty. We need an approximate kernel
function, and this is the Pari function {\tt kerr}. So type {\tt kerr(m)}, 
and we indeed get a non-trivial kernel, and the coefficients are (close to)
integers as they should. Typing {\tt addell(e,v[1],v[3])} does indeed show
that it is equal to {\tt v[4]}.

Taking any other four points, we would in fact always find a dependency. Let's
find them all. Type {\tt vp=[v[1],v[2],v[3]]\~;m=mathell(e,vp);det(m)}. This is
now clearly non-zero so the first 3 points are linearly independent, showing
that the rank of the curve is at least equal to 3 (it is in fact equal to 3,
and is the curve of smallest conductor having rank 3). We would like to see
whether the other points are dependent. For this, we use the function
{\tt bilhell}. Indeed, if {\tt Q} is some point which is dependent on
{\tt v[1],v[2]} and {\tt v[3]}, then {\tt gauss(m,bilhell(e,vp,Q))} will by
definition give the coefficients of the dependence relation. If these 
coefficients are close to integers, then there is a dependency, otherwise not.
This is much safer than using the {\tt kerr} function.

Thus, type {\tt w=vector(18,k,gauss(m,bilhell(e,vp,v[k])))}. We ``see'' that
the coefficients are all very close to integers, and we can prove it by
typing {\tt wr=round(w);sqrt(norml2(w-wr))}, which gives an upper bound
on the maximum distance to an integer. Thus {\tt wr} is the vector expressing
all the components of {\tt v} on its first 3. We are thus led to strongly
believe that the curve has rank exactly 3, and this can be proved to be
the case.
\smallskip
Let's now explore a few more elliptic curve related functions. Let us keep
our curve {\tt e} of rank 3, and type 

{\tt v1=[1,0];v2=[2,0];} then {\tt z1=zell(e,v1)} and {\tt z2=zell(e,v2)}.

We thus get the complex parametrization of the curve. To add the points
{\tt v1} and {\tt v2}, we should of course type {\tt addell(e,v1,v2)},
but we can also type {\tt pointell(e,z1$+$z2)} which of course has the
disadvantage of giving complex numbers, but illustrates how the group law
on {\tt e} is obtained from the addition law on $\Bbb C$.

Type {\tt f=x$*$series(anell(e,30),x)}. This gives a power series which,
is the Fourier expansion of a modular form of weight 2 on $\Gamma_0(5077)$
(this has been proved directly, but also follows from Wiles's result since
{\tt e} is semi-stable). In fact, to find the modular parametrization of
the curve, type {\tt modul=taniyama(e)}, then {\tt u=modul[1];v=modul[2];}.
Type {\tt (v\^{}2$+$v)$-$(u\^{}3$-$7$*$u$+$6)} to see that this parametrizes
Now type {\tt x$*$u'$/$(2$*$v$+$1)}, and we see that this is equal to the
modular form {\tt f} found above. The functions {\tt u} and {\tt v}, considered
on the upper half plane (with $x=e^{2i\pi\tau}$) are in fact modular 
{\it functions} on $\Gamma_0(5077)$.
\smallskip
Finally, let us come back once more to the curve defined by typing
{\tt e=initell([0,$-$1,1,0,0])} where we had seen that the point {\tt q=[0,0]}
was of order 5. We know that the conductor of this curve is equal to 11
(type {\tt globalred(e)}). We want the sign of the functional equation.
Type {\tt lseriesell(e,1,-11,1.1)}, then {\tt lseriesell(e,1,-11,1)}. Since
the values are clearly different, the sign cannot be $-$ (look at Chapter 3
for the explanation of {\tt lseriesell}). In fact there is an algebraic 
algorithm which would allow to compute this sign, but it has not yet been
completely written, although in case of conductors prime to 6 it is very
simple.

Now type {\tt ls=lseriesell(e,1,11,1)}, and just as a check
{\tt lseriesell(e,1,11,1.1)}. The values agree (approximately) as they should,
and give the value of the L-function of {\tt e} at 1. Now according to
the Birch and Swinnerton-Dyer conjecture (which is proved for this curve),
{\tt ls} is given by the following formula (in this case):
$$L(E,1)=\dfrac{\Omega\cdot c\cdot\text{Sha}}{|E_{\text{tors}}|^2}\enspace,$$
where $\Omega$ is the real period of $E$, $c$ is the global Tamagawa number,
produuct of the local $c_p$ for primes $p$ dividing the conductor, Sha is the
order of the Tate-Shafarevitch group, and $E_{\text{tors}}$ is the torsion
group of $E$.

Now we know many of these quantities: $\Omega$ is equal to {\tt e[15]}
(if there had been 3 real roots instead of 1 in {\tt e[14]}, $\Omega$ would
be equal to {\tt 2$*$e[15]}). The Tamagawa number $c$ is given as the last
component of {\tt globalred(e)}, and here is equal to 1. We already know
that the torsion subgroup of $E$ contains a point of order 5, and 
typing {\tt torsell(e)} shows that it is of order exactly 5. Hence type 
{\tt ls$*$25/e[15]}. We find $1.00000\dots$ as a result, which shows that 
Sha must be equal to 1.
\medskip
\section{Functions Related to Polynomials and Power Series}
\medskip
First a word of warning to the unwary: it is essential to understand the
crucial difference between exact and inexact objects (see Chapter 1). This
is especially important in the case of polynomials. Let's immediately take
a plunge into these problems. Type {\tt gcd$(\x^2-1,\x^2-3*\x+2)$}. No problem,
the result is {\tt x$-$1} as expected. But now type 
{\tt gcd$(\x^2-1.,\x^2-3.*\x+2.)$}. You are lucky, the result is almost
correct except from a bizzare factor of 3 which come from the way Pari does
the computation. In any cas, it is still essentially a reasonable result.
But now type {\tt gcd(x$-$pi,x\^{}2$-$6$*$zeta(2))}. Although this should be
equal to {\tt x$-$pi}, Pari finds 1 as a result. This is because the notion
of GCD of nonexact polynomials doesn't make much sense. However, type
{\tt resultant(x$-$pi,x\^{}2$-$6$*$zeta(2))}. The result is extremely close
to zero, showing that indeed the GCD is non-trivial, without telling us
what it is. This being said, we will usually use polynomials (and power series)
with exact coefficients in our examples.
\smallskip
Type {\tt pol=cyclo(15)}. This gives the $15$-th cyclotomic polynomial.
Type {\tt r=roots(pol)}. You have the 8 complex roots of pol given to 28
significant digits. To see them better, type {\tt\bs b}. As you see, they are
given as pairs of complex conjugate roots, in a random order. The only ordering
done by the function {\tt roots} concern the real roots, which are given
first, and in increasing order.

The roots of {\tt pol} are by definition the primitive $15$-th roots of unity.
To check this, simply type {\tt rc=r\^{}15}. Why, we get an error message! 
Well, fair enough, vectors cannot be multiplied (even less raised to a power) 
that easily. However, type {\tt rc=r\^{}15.} with a $.$ at the end. Now it 
works, because powering to a non-integer (i.e. here real) power is a
transcendental function and hence is applied termwise. Note that the fact that
$15.$ is a real number which is represented exactly as an integer has nothing 
to do with the problem.

We see that the components of the result are very close to 1. It is however
tedious to look at all these real and imaginary part. It would be impossible
if we had many more. Let's do it automatically. Type {\tt rr=round(real(rc))},
then {\tt sqrt(norml2(rc-rr))}. We see that {\tt rr} is indeed all 1's, and
that the L2-norm of {\tt rc-rr} is around $10^-27$, reasonable enough when we
work with 28 significant digits. Note that the function {\tt norml2}, contrary
to what it may seem, does not give the L2 norm but its square, hence we must
take the square root (well, not really in this case!).
\smallskip
However, it is necessary to stress the fact that the {\tt roots} function of
Pari is {\it not} very robust, and may sometime give wrong answers. This
may be corrected in some future version, but the user should be aware of this.
It has been desigened essentially for speed. A slower and sometimes better
version is {\tt rootslong}, but even that sometimes fails. For example. type
{\tt pol=(x$-$10000)$*$(x$-$10001)}, then {\tt roots(pol)}. The result is
clearly wrong. If however we type {\tt rootslong(pol)} the results are now
correct. But now type {\tt pol2=subst(pol,x,x$*$10000)$/$10\^{}8}, then
{\tt rootslong(pol2)}. The result is wrong, even with {\tt rootslong}. So be
careful. Essentially, as is well known, the problem occurs when the polynomial
has two roots which are very close. Still, {\tt roots} is very useful.
\smallskip
Now type {\tt pol=x\^{}5$+$x\^{}4$+$2$*$x\^{}3$-$2$*$x\^{}2$-$4$*$x$-$3},
then {\tt factor(pol)}. The polynomial {\tt pol} factors over $\Q$ (or $\Z$)
as a product of two factors. Now type {\tt f(p)=factorpadic(pol,p,10)}.
This creates a function {\tt f(p)} which factors {\tt pol} over $\Q_p$ to 
$p$-adic precision 10. Type {\tt factor(disc(pol))}. We see that the primes
dividing the discriminant are $11$, $23$ and $37$. Type {\tt f(5)}, 
{\tt f(11)}, {\tt f(23)}, and {\tt f(37)} to see different splittings.

Similarly, we can type {\tt g(p)=lift(factmod(pol,p))}, and {\tt g(2)},
{\tt g(11)}, {\tt g(23)} and {\tt g(37)} which show the different 
factorizations, this time over $\F_p$.

In fact, even better: type {\tt pol2=x\^{}3$+$x\^{}2$+$x$-$1}, then

\centerline{\tt centerlift(lift(factfq(x\^{}3$+$x\^{}2$+$x$-$1,3,t\^{}3$+$t\^{}2$+$t$-$1)))}.

This factors the polynomial {\tt pol2} over the field $\F_3$ adjoined
with a root of the polynomial {\tt t\^{}3$+$t\^{}2$+$t$-$1}, which is a form
of the field $\F_{27}$. We know that Gal$(\F_{27}/\F_3)$ is cyclic of order
3 generated by the Frobenius homomorphism $u\mapsto u^3$, and the roots that
we have found give the action of the powers of the Frobenius on {\tt t}
(if you do not know what I am talking about, please learn it with the help
of explicit examples given by Pari).

Similarly, type {\tt pol3=x\^{}4$-$4$*$x\^{}2$+$16} and
{\tt fn=factornf(pol3,t\^{}2$+$1)}, and we get the factorization of the 
polynomial {\tt pol3} over the number field defined by {\tt t\^{}2$+$1},
i.e. over $\Q(i)$. To see the result even better, type {\tt lift(fn)},
remembering that {\tt t} stands for the generator of the number field
(here equal to $i=\sqrt{-1}$). 

Note that it is possible, although ill advised, to use the same variable
for the polynomial and the number field. You may for example type
{\tt fn2=factornf(pol3,x\^{}2$+$1)}, and the result is correct. However,
the Pari object thus created may give unreasonable results. For example,
if you type {\tt lift(fn2)} in the example above, you will get a strange
object, with a symbol such as {\tt x$*$x} typed. This is because
Pari knows that the dummy variable {\tt x} is not the same as the explicit
variable {\tt x}, but since it must print it when you lift, it has to
do something.
\smallskip
To summarize, in addition to being able to factor integers, you can
factor polynomials over $\C$ and $\R$ (this is the function {\tt roots}),
over $\F_p$ (the function {\tt factmod}, over $\F_{p^k}$ (the function
{\tt factfq}), over $\Q_p$ (the function {\tt factorpadic}), over $\Q$ or
$\Z$ (the function {\tt factor}), and over number fields (the function
{\tt factnf}). In the present version \vers{}, it is {\it not} possible to
factor over other rings, for example it cannot factor multivariate polynomials.
\smallskip
Functions related to factoring are {\tt apprpadic}, {\tt rootmod}, 
{\tt rootpadic}, {\tt sturm} and {\tt sturmpart}. Play with them a little.

Now let's type {\tt polsym(pol3,20)}, where {\tt pol3} is the same polynomial
as above. This gives the sum of the $k$-th powers of the roots of {\tt pol3}
up to $k=20$, of course computed using Newton's formula and not using
{\tt roots}. You notice that every odd sum is zero (this is trivial since the
polynomial is even), but also that the signs follow a regular pattern and
that the  (nonzero) absolute values are powers of 2. This is true: prove it,
and more precisely find an explicit formula for the $k$-th symmetric power
not involving (non-rational) algeberaic numbers.
\smallskip
You can also play a little with the function {\tt karamul} (for Karatsuba
multiplication). Put the timer on ({\tt \#}), and then set 
{\tt pol4=(x$+$1)\^{}200;} (do not forget the semicolon, you don't want to
see the result, just the timing. Similarly for the next few commands.
Then type {\tt pol5=pol4$*$pol4;}, and look at the timing. Now type instead
{\tt pol6=karamul(pol4,pol4,5);}. This should be about twice as fast (it is
implementation dependent so the gain may vary). Try with other values than 5.
Also type {\tt pol5$-$pol6} just to be sure that the result is correct.
So if you have polynomials which are that large, it may be worthwhile to use
the function {\tt karamul}.
\medskip
Now let's play a little with power series. We have already done so a little
at the beginning.  Type
\centerline{\tt 8$*$x$+$prod(1$+$O(x\^{}40),n=1,39,if(n\%4,1$-$x\^{}n,1))\^{}8}
You obtain a power series which has apparently only even powers of {\tt x}
appearing. This is surprising, but can be proved using the theory of modular
forms. Note that we have initialized the product to {\tt 1$+$O(x\^{}40)} and
not simply to 1 otherwise the whole computation would have been done with
polynomials, and this would first have been extremely long (try if you want,
but don't forget to type {\tt Ctrl-C} when you get tired of waiting) and also
totally useless since the coefficients of {\tt x\^{}40} and above are 
irrelevant anyhow if we stop the product at {\tt n=39}.

While we are on the subject of modular forms (which, together with taylor
series expansions of usual functions is another great source of power series),
type {\tt\bs serieslength=122} then {\tt d=x$*$eta(x)\^24}. This gives the
first 122 terms of the (Fourier) series expansion of the modular discriminant
function $\Delta$ of Ramanujan, its coefficients being by definition the
Ramanujan $\tau$ function which has a number of marvelous properties (look at
any book on modular forms for explanations). We would like to see its 
properties modulo 2. Type {\tt d\%2}. Hmm, apparently Pari doesn't like to
reduce coefficients of power series (or polynomials for that matter) directly.
Can we do it without writing a little program? No problem. Type instead
{\tt lift(mod(1,2)$*$d)} and now this works like a charm. 

The pattern in the result is clear. Of course, it now remains to prove it
(see Antwerp III). Similarly, type {\tt centerlift(mod(1,3)$*$d)}. This time
the pattern is less clear, but nonetheless there is one. Refer to Anwerp III
again.
\medskip
\section{Performing Linear Algebra}
\medskip
All the standard linear algebra programs available of course, and many more.
In addition, linear algebra over $\Z$, i.e. work on lattices, can also be
performed. Let us see how this works.
First of all, a vector space (more generally a module) will be given by a
generating set of vectors (often a basis) which will be representend as
{\it column} vectors. This set of vectors will in turn be represented as
a matrix: in Pari, we have chosen to consider matrices as row vectors of
column vectors. The base field (or ring) can be any ring type Pari supports.
However, certain operations are specifically written for a real or complex
base field, while others are written for $\Z$ as the base ring.





%Table of contents of the tutorial:
%\section{Greetings!}
%\section{Warming up}
%\section{The Rest of the PARI Types}
%\section{Elementary Arithmetic Functions}
%\section{Working in Quadratic Number Fields}
%\section{Working in General Number Fields}
%\section{Using Transcendental Functions}
%\section{Using Numerical Tools}
%\section{Working with Elliptic Curves}
%\section{Functions Related to Polynomials and Power Series}
%\section{Performing Linear Algebra}
%\section{GP Programming}


\vfill\eject