Annotation of OpenXM_contrib/pari-2.2/doc/tutorial.tex, Revision 1.1.1.1
1.1 noro 1: % $Id: tutorial.tex,v 1.7 2000/11/06 18:59:00 karim Exp $
2: % Copyright (c) 2000 The PARI Group
3: %
4: % This file is part of the PARI/GP documentation
5: %
6: % Permission is granted to copy, distribute and/or modify this document
7: % under the terms of the GNU Free Documentation License
8:
9: % This should be compiled with plain TeX
10: \def\TITLE{A Tutorial for Pari-GP}
11: \input parimacro.tex
12: \ifPDF
13: \input pdfmacs.tex
14: \fi
15: \def\maketitle#1{\currentlabel.\ #1}
16: \begintitle
17: \vskip2.5truecm
18: \centerline{\mine A Tutorial}
19: \vskip1.truecm
20: \centerline{\mine for}
21: \vskip1.truecm
22: \centerline{\mine PARI / GP}
23: \vskip1.truecm
24: \authors
25: \centerline{last updated November 5, 2000}
26: \centerline{(this document distributed with version \vers)}
27: \endtitle
28:
29: \copyrightpage
30:
31: This booklet is intended to be a guided tour and a tutorial to the GP
32: calculator. Many examples will be given, but each time a new function is
33: used, the reader should look at the appropriate section in the User's Manual
34: for detailed explanations. Hence although this chapter can be read
35: independently (for example to get rapidly acquainted with the possibilities
36: of GP without having to read the whole manual), the reader will profit most
37: from it by reading it in conjunction with the reference manual.
38:
39: \section{Greetings!}
40:
41: So you are sitting in front of your workstation (or terminal, or PC,\dots),
42: and you type \kbd{gp} to get the program started (remember to always hit the
43: \kbd{Enter} key and not the \kbd{Return} key on a Macintosh computer).
44:
45: It says hello in its particular manner, and then waits for you after its
46: \kbd{prompt}, initially \kbd{?} (or something like {\bf gp}~\kbd{>}).
47:
48: Type \kbd{2 + 2}. What happens? Maybe not what you expect. First of all, of
49: course, you should tell GP that your input is finished, and this is done by
50: hitting the \kbd{Return} (or \kbd{Newline}) key, or the \kbd{Enter} key on
51: the Mac. If you do exactly this, you will get the expected answer. However
52: some of you may be used to other systems like Macsyma or Maple. In this case,
53: you will have subconsciously ended the line with a semicolon ``\kbd{;}''
54: before hitting \kbd{Return}, since this is how it is done on those systems.
55: In that case, you will simply see GP answering you with a smug expression,
56: i.e.~a new prompt and no answer! This is because a semicolon at the end of a
57: line in GP tells it to keep the result, but not to print it (you will
58: certainly want to use this feature if the output is several pages long).
59:
60: Try \kbd{27 * 37}. Wow! even multiplication works. Actually, maybe those
61: spaces are not necessary after all. Let's try \kbd{27*37}. Yup, seems to be
62: ok. We'll still insert them in this document since it makes things easier to
63: read, but as GP does not care about them, you don't have to type them all!
64:
65: Now this session is getting lengthy, so the second thing one needs to learn
66: is to quit. Each system has its quit signal (to name a few: \kbd{quit},
67: \kbd{exit}, \kbd{system},\dots). In GP, you can use \kbd{quit} or \b{q}
68: (backslash q), the \kbd{q} being of course for quit. Try it.
69:
70: Now you've done it! You're out of GP, so how do you want to continue studying
71: this tutorial? Get back in please (see above).
72:
73: Let's get to more serious stuff. Let's see, I seem to remember that the
74: decimal expansion of $1/7$ has some interesting properties. Let's see what GP
75: has to say about this. Type \kbd{1 / 7}. What? This computer is making fun of
76: me, it just spits back to me my own input, that's not what I want!
77:
78: Now stop complaining, and think a little. This system has been written mainly
79: for pure mathematicians, and not for physicists (although they are welcome to
80: use it). And mathematically, $1/7$ is an element of the field $\Q$ of
81: rational numbers, so how else but $1/7$ can the computer give the answer to
82: you? (well maybe $2/14$, but why complicate matters?). Seriously, the basic
83: point here is that PARI (hence GP) will almost always try to give you a
84: result which is as precise as possible (we will see why ``almost'' later),
85: hence since here you gave an operation whose result can be represented
86: exactly, that's what it gives you.
87:
88: OK, but I still want the decimal expansion of $1/7$. No problem. Type one of
89: the following:
90: \bprog
91: 1./ 7
92: 1 / 7.
93: 1./ 7.
94: 1 / 7 + 0.@com etc \dots
95: @eprog
96: Immediately a number of decimals of this fraction appear (28 on most systems,
97: 38 on the others), and the repeating pattern is $142857$. The reason is that
98: you have included in the operations numbers like \kbd{0.}, \kbd{1.} or \kbd{7.}
99: which are {\it imprecise\/} real numbers, hence GP cannot give you an exact
100: result.
101:
102: Why 28 / 38 decimals by the way? Well, it is the default initial precision,
103: as indicated when you start GP. This has been chosen so that the
104: computations are very fast, and gives already 12 decimals more accuracy than
105: conventional double precision floating point operations. The precise value
106: depends on a technical reason: if your machine supports 64-bit integers (the
107: standard library can handle integers up to $2^{64}$), the default precision
108: will be 38 decimals, and 28 otherwise. As the latter is most probably the
109: case (if your machine is not a \kbd{DEC alpha}, that is), we'll assume it
110: henceforth.
111:
112: Only large mainframes or supercomputers have 28 digits of precision in their
113: standard libraries, and that is their absolute limit. Not here of course. You
114: can extend the precision (almost) as much as you like as we will see in a
115: moment.
116:
117: I'm getting bored, why don't we get on with some more exciting stuff? Well,
118: try \kbd{exp(1)}. Presto, comes out the value of $e$ to 28 digits. Try
119: \kbd{log(exp(1))}. Well, it's not exactly equal to 1, but pretty close! That's
120: what you lose by working numerically.
121:
122: Let's see, what could we try now. Hum, \kbd{pi}? The answer is not that
123: enlightening. \kbd{Pi}? Ah. This works better. But let's remember that GP
124: distinguishes between uppercase and lowercase letters. \kbd{pi} was as
125: meaningless to it as \kbd{stupid garbage} would have been: in both cases GP
126: will just create a variable with that funny unknown name you just used. Try
127: it! Note that it is actually equivalent to type \kbd{stupidgarbage}: all
128: spaces are suppressed from the input. In the \kbd{27~*~37} example it was
129: not so conspicuous as we had an operator to separate the two operands. This
130: has important consequences for the writing of GP scripts (more about this
131: later).
132:
133: By the way, you can ask GP about any identifier you think it might know
134: about: just type it, prepending a question mark ``\kbd{?}''. Try \kbd{?Pi}
135: and \kbd{?pi} for instance. On some systems, an extended online help might
136: be available: try doubling the question mark to check whether it's the case on
137: yours: \kbd{??Pi}. In fact the GP header already gave you that information if
138: it was the case (just before the copyright message). As well, if it says
139: something like ``\kbd{readline enabled}'' then you should have a look at
140: \secref{se:readline} in the User's Manual before you go on: it will be much
141: easier to type in examples and correct typos after you've done that.
142:
143: Now try \kbd{exp(Pi * sqrt(163))}. Hmmm, since from our last example we
144: suspect that the last digit may be wrong, can this really be an integer?
145: This is the time to change precision. Type \kbd{\b{p} 50}, then try
146: \kbd{exp(Pi * sqrt(163))} again. We were right to suspect that the last
147: decimal was incorrect, since we get even more nines than before, but it is
148: now convincingly clear that this is not an integer. Maybe it's a bug in PARI,
149: and the result is really an integer? Type \kbd{sqr(log(\%) / Pi)} immediately
150: after the preceding computation (\kbd{\%} means the result of the last
151: computed expression). More generally, the results are numbered
152: \kbd{\%1, \%2, \dots} {\it including} the results that you do not want to see
153: printed by putting a semicolon at the end of the line, and you can evidently
154: use all these quantities in any further computations. \kbd{sqr} is the square
155: function (\kbd{sqr(x) = x * x}), not to be confused with \kbd{sqrt} which is
156: the square root function). The result seems to be indistinguishable from $163$,
157: hence it does not seem to be a bug.
158:
159: In fact, it is known that $\exp(\pi*\sqrt{n})$ not only is not an integer
160: or a rational number, but is even a transcendental number when $n$ is a
161: non-zero rational number.
162:
163: So GP is just a fancy calculator, able to give me more decimals than I will
164: ever need? Not so, GP is incredibly more powerful than an ordinary
165: calculator, even independently of its arbitrary precision possibilities.
166:
167: \noindent {\bf Additional comments} (you are supposed to skip this at first,
168: and come back later)
169:
170: 1) If you are a PARI old timer, you will notice pretty soon (or maybe you
171: have already?) that many many things changed between the older 1.39.xx
172: versions and this one. And conspicuously, most function names have been
173: changed. We sincerely think it's for the best since they are much more
174: logical now and the systematic prefixing is very convenient coupled with the
175: automatic completion mechanism: it's now very easy to know what functions are
176: available to deal with, say, elliptic curves since they all share the prefix
177: \kbd{ell}.
178:
179: Of course, this is going to break all your nice old scripts. Well, you can
180: either change the compatibility level (typing \kbd{default(compatible, 3)}
181: will send you back to the stone-age behaviour of good ol' version 1.39.15),
182: or rewrite the scripts. We really advise you to do the latter if they are not
183: too long, since they can now be written much more cleanly than before
184: (especially with the new control statements: \kbd{break}, \kbd{next},
185: \kbd{return}), and besides it'll be as good a way as any to get used to the new
186: names. We {\it might} provide an automatic transcriptor with future versions.
187:
188: To know how a specific function was changed, just type
189: \kbd{whatnow({\rm function})}.
190:
191: 2) It seems that the text implicitly says that as soon as an imprecise number
192: is entered, the result will be imprecise. Is this always true? There is a
193: unique exception: when you multiply an imprecise number by the exact number 0,
194: you will get the exact 0. Compare \kbd{0 * 1.4} and \kbd{0.~*~1.4}. \smallskip
195: %
196: 3) Not only can the number of decimal places of real numbers be large, but the
197: number of digits of integers also. Try \kbd{100!}. It is never necessary to
198: tell GP in advance the size of the integers that it will encounter, since this
199: is adjusted automatically. On the other hand, for many computations with real
200: numbers, it is necessary to specify a default precision (initially 28 digits).
201: \smallskip
202: %
203: 4) Come back to 28 digits of precision (\kbd{\b{p} 28}), and type
204: \kbd{exp(24 * Pi)}. As you can see the result is printed in exponential format.
205: This is because GP never wants you to believe that a result is correct when it
206: is not.
207:
208: We are working with 28 digits of precision, but the integer part of
209: $\exp(24*\pi)$ has 33 decimal digits. Hence if GP had dutifully printed out 33
210: digits, the last few digits would have been wrong. Hence GP wants to print only
211: 28 significant digits, but to do so it has to print in exponential format.
212: \smallskip
213: %
214: 5) There are two ways to avoid this. One is of course to increase the precision
215: to more than 33 decimals. Let's try it. To give it a wide margin, we set the
216: precision to 40 decimals. Then we recall our last result (\kbd{\%}
217: or \kbd{\%n} where \kbd{n} is the number of the result). What? We still have
218: an exponential format! Do you understand why?
219:
220: Again let's try to see what's happening. The number you recalled had been
221: computed only to 28 decimals, and even if you set the precision to 1000
222: decimals, GP knows that your number has only 28 digits of accuracy but an
223: integral part with 33 digits. So you haven't improved things by increasing
224: the precision. Or have you? What if we retype \kbd{exp(24 * Pi)} now that we
225: have 40 digits? Try it. Now we do not have an exponential format, and we see
226: that at 28 decimals the last 6 digits would have been wrong if they had been
227: printed in fixed-point format.
228: \smallskip
229: %
230: 6) {\it Warning\/}. Try the following: starting in precision 28, type first
231: \kbd{default(format, "e0.50")}, then \kbd{exp(24 * Pi)}. Do you understand why
232: the result is so bad, and why there are lots of zeros at the end? Convince
233: yourself by typing \kbd{log(exp(1))}. The moral is that the
234: \kbd{default(format,)} command changes only the output format, but {\it not\/}
235: the default precision. On the other hand, the \b{p} command changes both the
236: precision and the output format.
237:
238: 7) What if I forget what the current precision is and I don't feel like
239: counting all the decimals? Well, you can learn about GP internal variables
240: (and change them!) using \kbd{default}. Type \kbd{default(realprecision)},
241: then \kbd{default(realprecision, 38)}. Huh? In fact this last command is
242: strictly equivalent to \kbd{\b{p} 38}! (admittedly more cumbersome to type).
243: There are more ``defaults'' than just \kbd{format} and \kbd{realprecision}:
244: type \kbd{default} by itself now, they are all there.
245:
246: 8) Note that the \kbd{default} command reacts differently according to the
247: number of input arguments. This is not an uncommon behaviour for GP functions.
248: You can see this from the online help (or the complete description in
249: Chapter~3): any argument surrounded by braces \kbd{\obr\cbr} in the function
250: prototype is optional (which really means that a {\it default} argument will be
251: supplied by GP). You can then check out from the text what effect a given value
252: will have, and in particular the default one.
253:
254: \section{Warming up}
255:
256: Another thing you better get used to pretty fast is error messages. Try
257: typing \kbd{1/0}. Couldn't be clearer. Taking again our universal example in
258: precision 28, type \kbd{floor(exp(24 * Pi))} (\kbd{floor} is the
259: mathematician's integer part, not to be confused with \kbd{truncate}, which is
260: the computer scientist's: \kbd{floor(-3.4)} is equal to $-4$ whereas
261: \kbd{truncate(-3.4)} is equal to $-3$). You get a more cryptic error message,
262: which you would immediately understand if you had read the additional
263: comments of the preceding section. Since I told you not to read them, the
264: explanation is simply that GP is unable to compute the integer part of
265: \kbd{exp(24 * Pi)} given only 28 decimals of accuracy, since it has 33 digits.
266:
267: Some error messages are even much more cryptic than that and are sometimes
268: not so easy to understand (well, it's nothing compared to \TeX's error
269: messages!).
270:
271: For instance, try \kbd{log(x)}. Not really clear, is it? But no matter, it
272: simply tells you that GP simply does not understand what \kbd{log(x)} is
273: (although it does know the \kbd{log} function, as \kbd{?log} will readily
274: tell us).
275:
276: Now let's try \kbd{sqrt(-1)} to see what error message we get now. Haha! GP
277: even knows about complex numbers, so impossible to trick it that way.
278: Similarly, try typing \kbd{log(-2)}, \kbd{exp(I*Pi)}, \kbd{I\pow I},\dots
279: So we have a lot of real and complex analysis at our disposal (note that
280: there is always a specific branch of multivalued complex transcendental
281: functions which is taken, specified in the manual). Again, beware that
282: \kbd{I} and \kbd{i} are not the same thing. Compare \kbd{I\pow2} with
283: \kbd{i\pow2} for instance.
284:
285: Just for fun, let's try \kbd{6*zeta(2) / Pi\pow2}. Pretty close, no?
286:
287: \medskip
288: Now GP didn't seem to know what \kbd{log(x)} was, although it did know how to
289: compute numerical values of \kbd{log}. This is annoying. Maybe it knows the
290: exponential function? Let's give it a try. Type \kbd{exp(x)}. What's this? If
291: you had had any experience with other systems, the answer should have simply
292: been \kbd{exp(x)} again. But here the answer is the Taylor expansion of the
293: function around $\kbd{x}=0$, to 16 terms. 16 is the default
294: \kbd{seriesprecision}, which can be changed by typing \kbd{\b{ps} $n$} or
295: \kbd{default(seriesprecision, $n$)} where $n$ is the number of terms that you
296: want in your power series (note the \kbd{O(x\pow16)} which ends the series,
297: and which is trademark of this type of object in GP. It is the familiar
298: ``big--oh'' notation of analysis).
299:
300: You will thus automatically get the Taylor expansion of any function that can
301: be expanded around $0$, and incidentally this explains why we weren't
302: able to do anything with \kbd{log(x)} which is not defined at 0. But if we
303: try \kbd{log(1+x)}, then it works. But what if we wanted the expansion
304: around a point different from 0? Well, you're able to change $x$ into
305: $x-a$, aren't you? So for instance you can type \kbd{log(x+2)} to
306: have the expansion of \kbd{log} around $\kbd{x}=2$. As exercises you can
307: try
308: \bprog
309: cos(x) cos(x)^2 + sin(x)^2 exp(cos(x))
310: gamma(1 + x) exp(exp(x) - 1) 1 / tan(x)
311: @eprog
312: for different values of \kbd{serieslength}.
313:
314: Let's try something else: type \kbd{(1 + x)\pow 3}. No \kbd{O(x)} here, since
315: the result is a polynomial. Haha, but I have learnt that if you do not take
316: exponents which are integers greater or equal to 0, you obtain a power series
317: with an infinite number of non-zero terms. Let's try. Type
318: \kbd{(1 + x)\pow (-3)} (the parentheses around \kbd{-3} are not necessary but
319: make things easier to read). Surprise! Contrary to what we expected, we don't
320: get a power series but a rational function. Again this is for the same reason
321: that \kbd{1 / 7} just gave you $1/7$: the result being exact, PARI doesn't see
322: any reason to make it non-exact.
323:
324: But I still want that power series. To obtain it, you can do as in the $1/7$
325: example and type
326: \bprog
327: (1 + x)^(-3) + O(x^16)
328: (1 + O(x^16)) * (1 + x)^(-3)
329: (1 + x + O(x^16))^(-3)@com, etc \dots
330: @eprog
331: You can also use the series constructor which transforms any object into a
332: power series, using the current \kbd{seriesprecision}, and simply type
333: \bprog
334: Ser( (1 + x)^(-3) )
335: @eprog
336:
337: Now try \kbd{(1 + x)\pow (1/2)}: we obtain a power series, since the
338: result is an object which PARI does not know how to represent exactly (we
339: could teach PARI about algebraic functions, but then take \kbd{(1 + x)\pow Pi}
340: as another example). This gives us still another solution to our preceding
341: exercise: we can type \kbd{(1 + x)\pow (-3.)}. Since \kbd{-3.} is not an exact
342: quantity, PARI has no means to know that we are dealing with a rational
343: function, and will instead give you the power series, this time with real
344: instead of integer coefficients.
345:
346: Finally, if you want to be really fancy, you can type
347: \kbd{taylor((1 + x)\pow (-3), x)} (look at the entry for \kbd{taylor} for the
348: description of the syntax), but this is in fact almost never used.
349: \smallskip
350:
351: To summarize, in this section we have seen that in addition to integers, real
352: numbers and rational numbers, PARI can handle complex numbers, polynomials,
353: power series, rational functions. A large number of functions exist which
354: handle these types, but in this tutorial we will only look at a few.
355:
356: \noindent {\bf Additional comments} (as before, you are supposed to skip this
357: at first reading)
358:
359: 1) To be able to duplicate the following example, first type \b{y} to
360: suppress automatic simplification.
361:
362: A complex number has a real part and an imaginary part (who would have
363: guessed?). However, beware that when the imaginary part is the exact integer
364: zero, it is not printed, but the complex number is not converted to a real
365: number. Hence it may {\it look\/} like a real number without being one, and
366: this may cause some confusion in programs which expect real numbers. For
367: example, type \kbd{n = 3 + I - I}. The answer reads \kbd{3}, as expected. But
368: it is still a complex number. Check it with \kbd{type(n)}. Hence if you then
369: type \kbd{(1+x)\pow n}, instead of getting the polynomial
370: \kbd{1 + 3*x + 3*x\pow 2 + x\pow 3} as expected, you obtain a power series.
371: Worse, when you try to apply an arithmetic function, say the Euler totient
372: function (known as \kbd{eulerphi} to GP), you get a sententious error message
373: recalling you that ``arithmetic functions want integer arguments''. You would
374: have guessed yourself, but the message is difficult to understand since 3 looks
375: like a genuine integer!!! (Please read again if this is not clear. It is a
376: common source of incomprehension).
377:
378: Similarly, \kbd{3 + x - x} is not the integer 3 but a constant polynomial
379: (in the variable \kbd{x}), equal to $3 = 3x^0$.
380:
381: If you want the final expression to be in the simplest form possible (for
382: example before applying an arithmetic function, or simply because things will
383: go faster afterwards), apply the function \kbd{simplify} to the result. In
384: fact, by default GP does this automatically at the end of a GP command. Note
385: that it does {\it not\/} simplify in intermediate expressions. This default
386: can be switched off and on by the command \b{y}. This is why I asked you to
387: type this before starting.
388:
389: 2) As already stated, power series expansions are always implicitly around
390: $\kbd{x}=0$. When we wanted them around $\kbd{x}=\kbd{a}$, we replaced
391: \kbd{x} by \kbd{z + a} in the function we wanted to expand. For complicated
392: functions, it may be simpler to use the substitution function \kbd{subst}.
393: For example, if \kbd{p~= 1 / (x\pow 4 + 3*x\pow 3 + 5*x\pow 2 - 6*x + 7)},
394: you may not want to retype this, replacing \kbd{x} by \kbd{z~+ a}, so you can
395: write \kbd{subst(p, x, z+a)} (look up the exact description of the
396: \kbd{subst} function).
397:
398: Now try typing \kbd{p = 1 + x + x\pow 2 + O(x\pow 10)}, then
399: \kbd{subst(p, x, z+1)}. Do you understand why you get an error message?
400:
401: \section{The Remaining PARI Types}
402: Let's talk some more about the basic PARI types.
403:
404: Type \kbd{p = x * exp(-x)}. As expected, you get the power series expansion
405: to 16 terms (if you have not changed the default). Now type
406: \kbd{pr = serreverse(p)}. You are asking here for the {\it reversion\/} of the
407: power series \kbd{p}, in other words the inverse function. This is possible
408: only for power series whose first non-zero coefficient is that of $x^1$. To
409: check the correctness of the result, you can type \kbd{subst(p, x, pr)} or
410: \kbd{ subst(pr, x, p)} and you should get back \kbd{x + O(x\pow 17)}.
411:
412: Now the coefficients of \kbd{pr} obey a very simple formula. First, we would
413: like to multiply the coefficient of \kbd{x\pow n} by \kbd{n!} (in the case of
414: the exponential function, this would simplify things considerably!). The PARI
415: function \kbd{serlaplace} does just that. So type \kbd{ps = serlaplace(pr)}.
416: The coefficients now become integers, which can be immediately recognized by
417: inspection. The coefficient of $x^n$ is now equal to
418: $n^{n-1}$. In other words, we have
419: %
420: $$\kbd{pr} = \sum_{n\ge1}\dfrac{n^{n-1}}{n!} X^{n}.$$
421: %
422: Do you know how to prove this? (If you have never seen this, the proof is
423: difficult.)
424: \smallskip
425: %
426: Of course PARI knows about vectors (rows and columns are distinguished, even
427: though mathematically there is no difference) and matrices. Type for example
428: \kbd{[1,2,3,4]}. This gives the row vector whose coordinates are 1, 2, 3 and
429: 4. If you want a column vector, type \kbd{[1,2,3,4]\til}, the tilde meaning
430: of course transpose. You don't see much difference in the output, except for
431: the tilde at the end. However, now type \b{b}: lo and behold, the vector does
432: become a column. The \b{b} command is used mainly for this purpose.
433:
434: Type \kbd{m = [a,b,c; d,e,f]}. You have just entered a matrix with 2 rows and
435: 3 columns. Note that the matrix is entered by {\it rows\/} and the rows are
436: separated by semicolons ``\kbd{;}''. The matrix is printed naturally in a
437: rectangle shape. If you want it printed horizontally just as you typed it,
438: type \b{a}, or if you want this type of printing to be the permanent default
439: type \kbd{default(output, 1)}. Type \kbd{default(output, 0)} if you want to
440: come back to the original default.
441:
442: Now type \kbd{m[1,2]}, \kbd{m[1,]}, \kbd{m[,2]}. Are explanations necessary?
443: (In an expression such as \kbd{m[j,k]}, the \kbd{j} always refers to the
444: row number, and the \kbd{k} to the column number, and the first index is
445: always 1, never 0. This default cannot be changed.)
446:
447: Even better, type \kbd{m[1,2] = 5; m} (the semicolon also allows us to put
448: several instructions on the same line. The final result will be the output of
449: the last statement on the line). Now type \kbd{m[1,] = [15,-17,8]}. No
450: problem. Finally type \kbd{m[,2] = [j,k]}. You have an error message since you
451: have typed a row vector, while \kbd{m[,2]} is a column vector. If you type
452: instead \kbd{m[,2] = [j,k]\til} it works. \smallskip
453: %
454: Type now \kbd{h = mathilbert(20)}. You get the so-called ``Hilbert matrix''
455: whose coefficient of row $i$ and column $j$ is equal to $(i+j-1)^{-1}$.
456: Incidentally, the matrix \kbd{h} takes too much room. If you don't want to
457: see it, simply type a semi-colon ``\kbd{;}'' at the end of the line
458: (\kbd{h = mathilbert(20);}). This is an example of a ``precomputed'' matrix,
459: built into PARI. There are only a few. We will see later an example of a much
460: more general construction.
461:
462: What is interesting about Hilbert matrices is that first their inverses and
463: determinants can be computed explicitly (and the inverse has integer
464: coefficients), and second they are numerically very unstable, which make them
465: a severe test for linear algebra packages in numerical analysis. Of course
466: with PARI, no such problem can occur: since the coefficients are given as
467: rational numbers, the computation will be done exactly, so there cannot be
468: any numerical error. Try it. Type \kbd{d~=~matdet(h)} (you have to be a
469: little patient, this is quite a complicated computation). The result is a
470: rational number (of course) of numerator equal to 1 and denominator having
471: 226 decimal digits. How do I know, by the way? I did not count! Instead,
472: simply type \kbd{1.* d}. The result is now in exponential format, and the
473: exponent gives us the answer. Another, more natural, way would be to simply
474: type \kbd{sizedigit(1/d)}.
475:
476: Now type \kbd{hr = 1.* h;} (do not forget the semicolon, we don't want to see
477: all the junk!), then \kbd{dr = matdet(hr)}. You notice two things. First the
478: computation, although not instantaneous, is much faster than in the rational
479: case. The reason for this is that PARI is handling real numbers with 28
480: digits of accuracy, while in the rational case it is handling integers having
481: up to 226 decimal digits.
482:
483: The second more important fact is that the result is terribly wrong. If you
484: compare with \kbd{1.$*$d} computed earlier, which is correct, you will see
485: that only 2 decimals agree! This catastrophic instability is as already
486: mentioned one of the characteristics of Hilbert matrices. In fact, the
487: situation is much worse than that. Type \kbd{norml2(1/h - 1/hr)} (the
488: function \kbd{norml2} gives the square of the $L^2$ norm, i.e.~the sum of the
489: squares of the coefficients). Again be patient since computing \kbd{1/h} will
490: take even more time (not much) than computing \kbd{matdet(h)}. The result is
491: larger than $10^{50}$, showing that some coefficients of \kbd{1/hr} are
492: wrong by as much as $10^{24}$ (the largest error is in fact equal to $7.644
493: \cdot 10^{24}$ for the coefficient of row 15 and column 14, which is a 28
494: digit integer).
495:
496: To obtain the correct result after rounding for the inverse, we have to use a
497: default precision of 56 digits (try it).
498: \smallskip
499:
500: Although vectors and matrices can be entered manually, by typing explicitly
501: their elements, very often the elements satisfy a simple law and one uses a
502: different syntax. For example, assume that you want a vector whose $i$-th
503: coordinate is equal to $i^2$. No problem, type for example
504: \kbd{vector(10,i, i\pow 2)} if you want a vector of length 10. Similarly, if
505: you type
506:
507: \centerline{\kbd{matrix(5,5,i,j, 1/(i+j-1))}}
508:
509: \noindent you will get the Hilbert matrix of order 5 (hence the
510: \kbd{mathilbert} function is redundant). The \kbd{i} and \kbd{j} represent
511: dummy variables which are used to number the rows and columns respectively
512: (in the case of a vector only one is present of course). You must not forget,
513: in addition to the dimensions of the vector or matrix, to indicate explicitly
514: the names of these variables.
515:
516: \misctitle{Warning.} The letter \kbd{I} is reserved for the complex number
517: equal to the square root of $-1$. Hence it is absolutely forbidden to use
518: it as a variable. Try typing \kbd{vector(10,I, I\pow 2)}, the error message
519: that you get clearly indicates that GP does not consider \kbd{I} as a
520: variable. There are two other reserved variable names: \kbd{Pi} and
521: \kbd{Euler}. All function names are forbidden as well. On the other hand
522: there is nothing special about \kbd{i}, \kbd{pi} or \kbd{euler}.
523:
524: When creating vectors or matrices, it is often useful to use boolean
525: operators and the \kbd{if()} statement (see the section on programming for
526: details on using this statement). Indeed, an \kbd{if} expression has a value,
527: which is of course equal to the evaluated part of the \kbd{if}. So for
528: example you can type
529: \bprog
530: matrix(8,8, i,j, if ((i-j)%2, x, 0))
531: @eprog
532: \noindent to get a checkerboard matrix of \kbd{x} and \kbd{0}. Note however
533: that a vector or matrix must be {\it created} first before being used. For
534: example, it is forbidden to write
535: \bprog
536: for (i = 1, 5, v[i] = 1/i)
537: @eprog
538: \noindent if the vector \kbd{v} has not been created beforehand (for example
539: by a \kbd{v = vector(5,j,0)} command).
540:
541: \medskip The last PARI types which we have not yet played with are types
542: which are closely linked to number theory (hence people not interested in
543: number theory can skip ahead).
544:
545: The first is the type ``integer--modulo''. Let us see an example. Type
546: \kbd{n = 10\pow 15 + 3}. We want to know whether this number is prime or not. Of
547: course we could make use of the built-in facilities of PARI, but let us do
548: otherwise. (Note that PARI does not actually prove that a number is prime: any
549: strong pseudoprime for a number of bases is declared to be prime.) We first
550: trial divide by the built-in table of primes. We slightly cheat here and use a
551: variant of the function \kbd{factor} which does exactly this. So type
552: \kbd{factor(n, 200000)} (the last argument tells \kbd{factor} to trial divide
553: up to the given bound and stop at this point. You can set it to 0 to trial
554: divide by the full set of built-in primes, which goes up to $500000$ by
555: default).
556:
557: The result is a 2 column matrix (as for all factoring functions), the first
558: column giving the primes and the second their exponents. Here we get a single
559: row, telling us that \kbd{n} is not divisible by any prime up to $200000$. We
560: could now trial divide further, or even cheat completely and call the PARI
561: function \kbd{factor} without the optional second argument, but before we do
562: this let us see how to get an answer ourselves.
563:
564: By Fermat's little theorem, if $n$ is prime we must have $a^{n-1}\equiv 1
565: \pmod{n}$ for all $a$ not divisible by $n$. Hence we could try this with $a=2$
566: for example. But $2^{n-1}$ is a number with approximately $3\cdot10^{14}$
567: digits, hence impossible to write down, let alone to compute. But instead type
568: \kbd{a = Mod(2,n)}. This creates the number $2$ considered now as an element
569: of the ring $R = \Z/\kbd{n}\Z$. The elements of $R$, called integermods, can
570: always be represented by numbers smaller than \kbd{n}, hence very small.
571: Fermat's theorem can be rewritten
572: %
573: $\kbd{a}^{n-1} = \kbd{Mod(1,n)}$
574: %
575: in the ring $R$, and this can be computed very efficiently. Type
576: \kbd{a\pow (n-1)}. The result is definitely {\it not\/} equal to
577: \kbd{Mod(1,n)}, thus {\it proving\/} that \kbd{n} is not a prime (if we had
578: obtained \kbd{Mod(1,n)} on the other hand, it would have given us a hint that
579: \kbd{n} is maybe prime, but never a proof).
580:
581: To find the factors is another story. One must use less naive techniques than
582: trial division (or be very patient). To finish this example, type
583: \kbd{factor(n)} to see the factors. Since the smallest factor is 14902357,
584: you would have had to be very patient with trial division! Note that none of
585: the factors in the decomposition are {\it proven} primes.
586: \smallskip
587:
588: The second specifically number-theoretic type is the $p$-adic numbers. I have
589: no room for definitions, so please skip ahead if you have no use for such
590: beasts. A $p$-adic number is entered as a rational or integer valued
591: expression to which is added \kbd{O(p\pow n)} (or simply \kbd{O(p)} if
592: $\kbd{n}=1$) where \kbd{p} is the prime and \kbd{n} the $p$-adic precision.
593: Apart from the usual arithmetic operations, you can apply a number of
594: transcendental functions. For example, type \kbd{n = 569 + O(7\pow 8)}, then
595: \kbd{s~=~sqrt(n)}, you obtain one of the square roots of \kbd{n} (if you want
596: to check, type \kbd{s\pow 2 - n}). Type now \kbd{l~=~log(n)}, then \kbd{e =
597: exp(l)}. If you know about $p$-adic logarithms, you will not be surprised that
598: \kbd{e} is not equal to \kbd{n}. Type \kbd{(n/e)\pow 6}: \kbd{e} is in fact
599: equal to \kbd{n} times a $(p-1)$-st root of unity.
600:
601: Incidentally, if you want to get back the integer 569 from the $p$-adic
602: number \kbd{n}, type \kbd{truncate(n)}.
603: \smallskip
604:
605: The third number-theoretic type is the type ``quadratic number''. This type
606: is specially tailored so that we can easily work in a quadratic extension of
607: a base field (usually $\Q$). It is a generalization of the type
608: ``complex''. To start, we must specify which quadratic field we want to work
609: in. For this, we use the function \kbd{quadgen} applied to the
610: {\it discriminant\/} \kbd{d} (as opposed to the radicand) of the quadratic
611: field. This returns a number (always printed as \kbd{w}) equal to
612: $(\kbd{d}+a) / 2$ where $a$ is equal to 0 or 1 according to whether \kbd{d} is
613: even or odd. The behavior of \kbd{quadgen} is a little special: although its
614: result is always printed as \kbd{w}, the variable \kbd{w} itself is not set
615: to that value. Hence it is necessary to write systematically
616: \kbd{w = quadgen(d)} using the variable name \kbd{w} (or \kbd{w1} etc. if you
617: have several quadratic fields), otherwise things will get confusing.
618:
619: So type \kbd{w = quadgen(-163)}, then \kbd{charpoly(w)} which asks for the
620: characteristic polynomial of \kbd{w} (in the variable \kbd{x};
621: you can type \kbd{charpoly(w, y)} to directly express it in terms of some
622: other variable without resorting to the \kbd{subst} function). The result
623: shows what \kbd{w} will represent. You can also ask for \kbd{1.*w} to see
624: which root of the quadratic has been taken, but this is rarely necessary. We
625: can now play in the field $\Q(\sqrt{-163})$. Type for example
626: \kbd{w\pow 10}, \kbd{norm(3 + 4*w)}, \kbd{1 / (4+w)}. More interesting, type
627: \kbd{a = Mod(1,23) * w} then \kbd{b = a\pow 264}. This is a generalization of
628: Fermat's theorem to quadratic fields. If you do not want to see the modulus 23
629: all the time, type \kbd{lift(b)}.
630:
631: Another example: type \kbd{p = x\pow 2 + w*x + 5*w + 7}, then \kbd{norm(p)}. We
632: thus obtain the quartic equation over $\Q$ corresponding to the relative
633: quadratic extension over $\Q(\kbd{w})$ defined by \kbd{p}.
634:
635: On the other hand, if you type \kbd{wr = sqrt(w\pow 2)}, do not expect to get
636: back \kbd{w}. Instead, you get the numerical value, the function \kbd{sqrt}
637: being considered as a ``transcendental'' function, even though it is
638: algebraic. Type \kbd{algdep(wr,2)} (this looks for algebraic relations
639: involving the powers of \kbd{w} up to degree 2). This is one way to get
640: \kbd{w} back. Similarly, type \kbd{algdep(sqrt(3*w + 5), 4)}. See the user's
641: manual for the function \kbd{algdep}.\smallskip
642:
643: The fourth number-theoretic type is the type ``polynomial--modulo'', i.e.
644: polynomial modulo another polynomial. This type is used to work in general
645: algebraic extensions, for example elements of number fields (if the base
646: field is $\Q$), or elements of finite fields (if the base field is
647: $\Z/p\Z$ for a prime $p$, defined by an integermod). In a sense it is a
648: generalization of the type quadratic number. The syntax used is the same as
649: for integermods. For example, instead of typing \kbd{w = quadgen(-163)}, you
650: can type \kbd{w = Mod(x, x\pow 2 - x + 41)}. Then, exactly as in the
651: quadratic case, you can type \kbd{w\pow 10}, \kbd{norm(3 + 4*w)},
652: \kbd{1 / (4+w)}, \kbd{a = Mod(1,23)*w}, \kbd{b = a\pow 264}, obtaining of
653: course the same results (type \kbd{lift(\dots)} if you don't want to see the
654: polynomial \kbd{x\pow 2 - x + 41} repeated all the time). Of course, the basic
655: interest is that you can work in any degree, not only quadratic (of course,
656: even for quadratic moduli, the corresponding elementary operations will be
657: slower than with quadratic numbers).
658:
659: There is however a slight difference in behavior. Keeping our polmod \kbd{w},
660: type \kbd{1.*w}. As you can see, the result is not the same. Type
661: \kbd{sqrt(w)}. Here, we obtain a vector with 2 components, the two components
662: being the principal branch of the square root of all the possible embeddings
663: of \kbd{w} in $\C$ ({\it not\/} the two square roots). More generally, if
664: \kbd{w} was of degree $n$, we would get an $n$-component vector, and similarly
665: for transcendental functions.
666:
667: We have at our disposal the usual arithmetic functions, plus a few others.
668: Type \kbd{a = Mod(x, x\pow 3 - x - 1)} defining a cubic extension. We can for
669: example ask for \kbd{b = a\pow 5}. Now assume we want to express \kbd{a}
670: as a polynomial in \kbd{b}. This is possible since \kbd{b} is also a
671: generator of the same field. No problem, type \kbd{modreverse(b)}. This gives
672: a new defining polynomial for the same field (i.e.~the characteristic
673: polynomial of \kbd{b}), and expresses \kbd{a} in terms of this new polmod,
674: i.e.~in terms of \kbd{a}. We will see this in much more detail in the number
675: field section.
676:
677: \section{Elementary Arithmetic Functions}
678:
679: Since PARI is aimed at number theorists, it is not surprising that there
680: exists a large number of arithmetic functions (see the list in the
681: corresponding section of the users manual). We have already seen several,
682: such as \kbd{factor}. Note that \kbd{factor} handles not only integers, but
683: also (univariate) polynomials. Type for example \kbd{factor(x\pow 15 - 1)}.
684: You can also ask to factor a polynomial modulo a prime $p$ (\kbd{factormod})
685: and even in a finite field which is not a prime field (\kbd{factorff}).
686:
687: Evidently, you have functions for computing GCD's (\kbd{gcd}), extended GCD's
688: (\kbd{bezout}), solving the Chinese remainder theorem (\kbd{chinese}) and so
689: on.
690:
691: In addition to the factoring facilities, you have a few functions related to
692: primality testing such as \kbd{isprime}, \kbd{ispseudoprime},
693: \kbd{precprime}, and \kbd{nextprime}. As previously mentioned, only strong
694: pseudoprimes are produced; there is no sophisticated primality test.
695:
696: We also have the usual multiplicative arithmetic functions: the M\"obius $\mu$
697: function (\kbd{moebius}), the Euler $\phi$ function (\kbd{eulerphi}), the
698: $\omega$ and $\Omega$ functions (\kbd{omega} and \kbd{bigomega}), the
699: $\sigma_k$ functions (\kbd{sigma}), which compute sums of $k$-th powers of the
700: positive divisors of a given integer, etc\dots
701:
702: You can compute continued fractions. For example, type \kbd{\b{p} 1000}, then
703: \kbd{contfrac(exp(1))}: you obtain the continued fraction of the base of
704: natural logarithms, which as you can see obeys a very simple pattern (can you
705: prove it?).
706:
707: In many cases, one wants to perform some task only when an arithmetic
708: condition is satisfied. GP gives you the following functions: \kbd{isprime}
709: as mentioned above, \kbd{issquare}, \kbd{isfundamental} to test whether an
710: integer is a fundamental discriminant (i.e.~$1$ or the discriminant of the
711: ring of integers of a quadratic field), and the \kbd{forprime}, \kbd{fordiv}
712: and \kbd{sumdiv} loops. Assume for example that we want to compute the
713: product of all the divisors of a positive integer \kbd{n}. The easiest way is
714: to write
715:
716: \centerline{\tt p=1; fordiv(n,d, p *= d); p }
717:
718: \noindent
719: (there is a very simple formula for this product: find and prove it). The
720: notation \kbd{p *= d} (inherited from the C programming language) is just a
721: shorthand for \kbd{p = p * d}.
722:
723: If we want to know the list of primes $p$ less than 1000 such that 2 is a
724: primitive root modulo $p$, one way would be to write:
725: \bprog
726: forprime(p=3,1000, if (znprimroot(p) == 2, print(p)))
727: @eprog\noindent
728: %
729: Note that this assumes that \kbd{znprimroot} returns the smallest primitive
730: root, and this is indeed the case. Had we not known about this, we could
731: have written
732: \bprog
733: forprime(p=3,1000, if (znorder(Mod(2,p)) == p-1, print(p)))
734: @eprog\noindent
735: %
736: (which is actually faster since we only compute the order of $2$ in $\Z/p\Z$,
737: instead of looking for a generator by trying successive elements whose orders
738: have to be computed as well.)
739:
740: Functions related to quadratic fields, binary quadratic forms and general
741: number fields will be seen in the next sections.
742:
743: \section{Performing Linear Algebra}
744: All the standard linear algebra programs are available of course, and many
745: more. In addition, linear algebra over $\Z$, i.e.~work on lattices, can also
746: be performed. Let us see how this works. First of all, a vector space (more
747: generally a module) will be given by a generating set of vectors (often a
748: basis) which will be representend as {\it column} vectors. This set of vectors
749: will in turn be represented as a matrix: in PARI, we have chosen to consider
750: matrices as row vectors of column vectors. The base field (or ring) can be any
751: ring type PARI supports (except $p$-adics which are currently not correctly
752: handled by the linear algebra package). However, certain operations are
753: specifically written for a real or complex base field, while others are
754: written for $\Z$ as the base ring.
755:
756: ----- TO BE COMPLETED -----
757:
758:
759:
760: \section{Using Transcendental Functions}
761:
762: All the elementary transcendental functions and several higher transcendental
763: functions (gam\-ma function, incomplete gamma function, error function,
764: exponential integral, $K$-bessel functions, confluent hypergeometric functions,
765: Riemann $\zeta$ function, polylogarithms, Weber functions, theta functions)
766: are provided. More may be written if the need arises.
767:
768: In this type of functions, the default precision plays an essential role.
769: In almost all cases transcendental functions work in the following way.
770: If the argument is exact, the result will be computed using the current
771: default precision. If the argument is not exact, the precision of the
772: argument is used for the computation. A note of warning however: even in this
773: case the {\it printed\/} value will be the current real format (usually the
774: same as the default precision). In the present chapter we assume that your
775: machine works with 32-bit long integers. If it is not the case, we leave it
776: to you as a very good exercise to make the necessary modifications.
777:
778: Let's assume that we have 28 decimals of default precision (this is what we
779: get automatically at the start of a GP session on 32-bit machines). Type
780: \kbd{e = exp(1)}. We get the number $e=2.718\dots$ to 28 decimals. Let us check
781: how many correct decimals we really have. The hard (but reasonable) way is to
782: proceed as follows. Change the precision to a substantially higher value, for
783: example by typing \kbd{\b{p} 50}. Then type \kbd{e}, then \kbd{exp(1)} once
784: again. This last value is the correct value of the mathematical constant $e$ to
785: 50 decimals, while the variable \kbd{e} shows the value that was computed to 28
786: decimals. Clearly they coincide to exactly 29 significant digits.
787:
788: A simpler way to see how many decimals we have, is to ask for \kbd{length(e)}
789: which indicates we have exactly $3$ mantissa words. Since
790: $3\ln(2^{32}) / \ln(10)\approx28.9$ we see that we have 28 or 29 significant
791: digits (on 32-bit machines).
792:
793: \smallskip
794: Come back to 28 decimals (\kbd{\b{p} 28}). If we type \kbd{exp(1.)}
795: you can check that we also obtain 28 decimals. However, type
796: \kbd{f = exp(1 + 10.\pow(-30))}. Although the default precision is still 28,
797: you can check using one of the two methods above that we have in fact 59
798: significant digits! The reason is that \kbd{1 + 10.\pow(-30)} is computed
799: according to the PARI philosophy, i.e.~to the best possible precision. Since
800: \kbd{10.} has 29 significant digits and 1 has ``infinite'' precision, the
801: number \kbd{1 + 10.\pow(-30)} will have $59=29+30$ significant digits,
802: hence \kbd{f} also.
803:
804: Now type \kbd{cos(10.\pow(-15))}. The result is printed as $1.0000\dots$, but
805: is of course not exactly equal to 1. Using \kbd{length(\%)}, we see that the
806: result has 7 mantissa words, giving us the possibility of having 67
807: correct significant digits. In fact (look in precision 100), only 60 are
808: correct. PARI gives you as much as it can, and since 6 mantissa words
809: would have given you only 57 digits, it uses 7. But why does it give so
810: precise a result? Well, it is the same reason as before. When $x$ is close
811: to 1, $\cos(x)$ is close to $1-x^2/2$, hence the precision is going to be
812: approximately the same as this quantity, which will be $1-0.5*10^{-30}$ where
813: $0.5*10^{-30}$ is considered with 28 significant digit accuracy, hence the
814: result will have approximately $28+30=58$ significant digits.
815:
816: Unfortunately, this philosophy cannot go too far. For example, when you
817: type \kbd{cos(0)}, GP should give you exactly 1. Since it is reasonable for
818: a program to assume that a transcendental function never gives you an exact
819: result, GP gives you $1.000\dots$ to one more mantissa word than the current
820: precision.
821: \medskip
822: OK, now let's see some more transcendental functions at work. Type
823: \kbd{gamma(10)}. No problem (type \kbd{9!} to check). Type \kbd{gamma(100)}.
824: The number is now written in exponential format because the default
825: accuracy is too small to give the correct result (type \kbd{99!} to get the
826: complete answer).
827: To get the complete value, there are two solutions. The first and most natural
828: one is to increase the precision. Since \kbd{gamma(100)} has 156 decimal
829: digits, type \kbd{\b{p} 170} (to be on the safe side), then \kbd{gamma(100)}
830: once again. After some work, the printed result is this time perfectly
831: correct.
832:
833: However, this is probably not the best way to proceed. Come back first to the
834: default precision (type \kbd{\b{p} 28}). As the gamma function increases
835: very rapidly, one usually uses its logarithm. Type \kbd{lngamma(100)}. This
836: time the result has a reasonable size, and is exactly equal to \kbd{log(99!)}.
837:
838: Try \kbd{gamma(1/2 + 10*I)}. No problem, we have the complex gamma function.
839: Now type
840:
841: \kbd{t = 1000; z = gamma(1 + I*t) * t\pow(-1/2) * exp(Pi/2*t)/sqrt(2*Pi)},
842:
843: \noindent then \kbd{norm(z)}. We see that \kbd{norm(z)} is very close to 1,
844: in accordance with the complex Stirling formula. \smallskip
845: %
846: Let's play now with the Riemann zeta function. First turn on the timer (type
847: \kbd{\#}). Type \kbd{zeta(2)}, then \kbd{Pi\pow 2/6}. This seems correct. Type
848: \kbd{zeta(3)}. All this takes essentially no time at all. However, type
849: \kbd{zeta(3.)}. Although the result is the same, you will notice that the
850: time is substantially larger (if your machine is too fast to see the
851: difference, increase the precision!). This is because PARI uses special
852: formulas to compute \kbd{zeta(n)} when \kbd{n} is a (positive or negative)
853: integer.
854:
855: Type \kbd{zeta(1 + I)}. This also works. Now for fun, let us compute in a
856: very naive way the first complex zero of \kbd{zeta}. We know that it is
857: of the form $1/2 + i*t$ with $t$ between 14 and 15. Thus, we can use the
858: following series of instructions. But instead of typing them directly, write
859: them into a file, say \kbd{zeta.gp}, then type \kbd{\b{r} zeta.gp} under GP to
860: read it in:
861: \bprog
862: {
863: t1 = 1/2 + 14*I;
864: t2 = 1/2 + 15*I; eps = 10^(-50);
865: z1 = zeta(t1);
866: until (norm(z2) < eps,
867: z2 = zeta(t2);
868: if (norm(z2) < norm(z1),
869: t3 = t1; t1 = t2; t2 = t3; z1 = z2
870: );
871: t2 = (t1+t2) / 2.;
872: print(t1 ": " z1)
873: )
874: }
875: @eprog
876:
877: Don't forget the braces: they tell GP that a sequence of instructions is going
878: to span many lines (another, less convenient, way would be to type \b{} at the
879: end of each line to tell the parser that the input was not yet finished).
880: By the way, you don't need to type in the suffix~\kbd{.gp} it will be
881: supplied by GP, if you forget it (the suffix is not mandatory either, it is
882: just more convenient to have all your GP scripts labeled in the same
883: distinctive way).
884:
885: We thus obtain the first zero to 25 significant digits.
886: \medskip
887: %
888: As mentioned at the beginning of this tutorial, some transcendental functions
889: can also be applied to $p$-adic numbers. This is as good a time as any to
890: familiarize yourself with them. Type \kbd{a~=~exp(7 + O(7\pow 10))}, then
891: \kbd{log(a)}. All seems in order. Now type \kbd{b = log(5 + O(7\pow 10))},
892: then \kbd{exp(b)}. Is something wrong, since we don't recover the number we
893: started with? Absolutely not. Type
894: \kbd{exp(b) * teichmuller(5 + O(7\pow 10))},
895: and we indeed recover our initial number. The function \kbd{teichmuller(x)}
896: is the Teichm\"uller character, and is characterized by the fact that it is
897: the unique \hbox{$(p-1)$-st} root of unity (here with $p=7$) which is
898: congruent to \kbd{x} modulo $p$, assuming that \kbd{x} is a $p$-adic
899: unit.\smallskip
900: %
901: Let us come back to real numbers for the moment. Type \kbd{agm(1,sqrt(2))}.
902: This gives the arithmetic-geometric mean of 1 and $\sqrt2$, and is the basic
903: method for computing (complete) elliptic integrals. In fact, type
904:
905: \kbd{Pi/2 / intnum(t=0,Pi/2, 1 / sqrt(1 + sin(t)\pow 2))},
906:
907: \noindent and the result is the same. The elementary transformation
908: \kbd{x = sin(t)} gives the mathematical equality
909: $$\int_0^1 \dfrac{dx}{\sqrt{1-x^4}} = \dfrac{\pi}{2\text{agm}(1,\sqrt2)}
910: \enspace,$$
911: which was one of Gauss's remarkable discoveries in his youth.
912:
913: Now type \kbd{2 * agm(1,I) / (1+I)}. As you see, the complex AGM also works,
914: although one must be careful with its definition. The result found is
915: almost identical to the previous one. Exercise: do you see why?
916:
917: Finally, type \kbd{agm(1, 1 + 7 + O(7\pow 10))}. So we also have $p$-adic
918: AGM. Note however that since the square root of a $p$-adic number is not
919: in general an element of the same $p$-adic field,
920: only certain $p$-adic AGMs can be computed. In addition,
921: when $p=2$, the congruence restriction is that \kbd{agm(a,b)} can be computed
922: only when \kbd{a/b} is congruent to 1 modulo $16$ (not 8 as could be
923: expected).\smallskip
924: %
925: Now type \kbd{?3}. This gives you the list of all transcendental functions.
926: Instead of continuing with more examples, we suggest that you experiment
927: yourself with the list of functions. In each case, try integer, real, complex
928: and $p$-adic arguments. You will notice that some have not been implemented
929: (or do not have a reasonable definition).
930:
931: \section{Using Numerical Tools}
932:
933: Although not written to be a numerical analysis package, PARI can
934: nonetheless perform some numerical computations. We leave for a subsequent
935: section linear algebra and polynomial computations.
936:
937: You of course know the formula $\pi = 4(1-\dfrac13+\dfrac15-\dfrac17+\cdots)$
938: which is deduced from the power series expansion of \kbd{atan(x)}. You also
939: know that $\pi$ cannot be computed from this formula, since the convergence
940: is so slow. Right? Wrong! Type \kbd{\b{p} 100} (just to make it more
941: interesting), then \kbd{4~*~sumalt(k=0, (-1)\pow k/(2*k + 1))}. In a split
942: second (admittedly more than simply typing \kbd{Pi}), we get $\pi$ to 100
943: significant digits (type \kbd{Pi} to check). In version 1.38, the method used
944: was a combination of a method due to Euler for accelerating alternating sums,
945: and a programming trick due to the Dutch mathematician van Wijngaarden (see
946: one of the Numerical Recipes books for an explanation). The method which we
947: presently use is considerably better, and is based on a combination of ideas of
948: F.~Villegas, D.~Zagier and H.~Cohen.
949:
950: Similarly, type \kbd{\b{p} 28} (otherwise the time will be too long) and
951: \kbd{sumpos(k=1, 1 / k\pow 2)}. Although once again the convergence is slow, a
952: similar trick allows to compute the sum when the terms are positive (compare
953: with the exact result \kbd{Pi\pow 2/6}). This is much less impressive because
954: it is much slower, but is still useful.
955:
956: Even better, \kbd{sumalt} can be used to sum divergent series! Type
957:
958: \centerline{\tt zet(s)= sumalt(k=1, (-1)\pow(k-1) / k\pow s) / (1 - 2\pow(1-s))}
959:
960: Then for positive values of \kbd{s} different from 1, \kbd{zet(s)} is equal
961: to \kbd{zeta(s)} and the series cvonverges, albeit slowly (sumalt doesn't
962: care however). For negative \kbd{s}, the series diverges, but \kbd{zet(s)}
963: still gives the correct result! Try \kbd{zet(-1)}, \kbd{zet(-2)},
964: \kbd{zet(-1.5)}, and compare with the corresponding values of \kbd{zeta}.
965: You should not push the game too far: \kbd{zet(-14.5)}, for example,
966: gives a completely wrong answer.
967:
968: Try \kbd{zet(I)}, and compare with \kbd{zeta(I)}. Even (some) complex values
969: work, although the sum is not alternating any more!
970:
971: Similarly, try \kbd{sumalt(n=1, (-1)\pow n/(n+I))}.
972: \medskip
973: %
974: More traditional functions are the numerical integration functions.
975: For example, type \kbd{intnum(t=1,2, 1/t)} and presto! you get 26 decimals
976: of $\log(2)$. Look at Chapter 3 to see the different integration functions
977: which are available.
978:
979: With PARI, however, you can go further since complex types are allowed.
980: For example, assume that we want to know the location of the zeros of the
981: function $h(z)=e^z-z$. We use Cauchy's theorem, which tells us that the
982: number of zeros in a disk of radius $r$ centered around the origin is
983: equal to
984: $$\dfrac{1}{2i\pi}\int_{C_r}\dfrac{h'(z)}{h(z)}\,dz\enspace,$$
985: where $C_r$ is the circle of radius $r$ centered at the origin.
986: Hence type
987: \bprog
988: fun(z) =
989: {
990: local(u);
991: u = exp(z);
992: (u-1) / (u-z)
993: }
994: zero(r) = r/(2*Pi) * intnum(t=0, 2*Pi, fun(r*exp(I*t)) * exp(I*t))
995: @eprog
996: \noindent (Here \kbd{u} is a local variable to the function \kbd{f}: whenever
997: a function is called, GP fills its argument list with the actual arguments
998: given, and initializes the other declared parameters and local variables to
999: 0. It will then restore their former values upon exit. If, however, we had
1000: not declared \kbd{u} in the function prototype, it would be considered as a
1001: global variable, whose value would be permanently changed. It is not
1002: mandatory to declare in this way all the parameters you use, but beware of
1003: side effects!)
1004:
1005: The function \kbd{zero(r)} will then count the number of zeros (we have simply
1006: made the change of variable $z = r*\exp(i*t)$).
1007:
1008: Now type \kbd{\b{p} 9} (otherwise the computation would take too long, and
1009: anyway we know that the result is an integer), then \kbd{zero(1)}, \kbd{zero(1.5)}.
1010: The result tells us that there are no zeros inside the unit disk, but that
1011: there are two (necessarily complex conjugate) whose modulus is between $1$
1012: and $1.5$. For the sake of completeness, let us compute them. Let $z$ be such
1013: a zero, and write $z=x+iy$ with $x$ and $y$ real. Then the equation $e^z-z=0$
1014: implies, after elementary transformations, that $e^{2x}=x^2+y^2$ and that
1015: $e^x\cos(y)=x$. Hence $y=\pm\sqrt{e^{2x}-x^2}$ and hence
1016: $e^x\cos(\sqrt{e^{2x}-x^2})=x$. Therefore, type
1017: \bprog
1018: fun(x) =
1019: {
1020: local(u);
1021: u = exp(x);
1022: u*cos(sqrt(u^2 - x^2)) - x
1023: }
1024: @eprog
1025:
1026: Then \kbd{fun(0)} is positive while \kbd{fun(1)} is negative. Come back to
1027: precision 28 and type
1028:
1029: \centerline{\tt x0 = solve(x=0,1, fun(x))}
1030:
1031: This quickly gives us the value of \kbd{x}, and we can then type
1032:
1033: \centerline{\tt z = x0 + I*sqrt(exp(2*x0) - x0\pow 2)}
1034:
1035: \noindent which (together with its complex conjugate) is the required zero.
1036: As a check, type \kbd{exp(z) - z}, and also \kbd{abs(z)}.
1037:
1038: Of course you can integrate over contours which are more complicated than
1039: circles, but you must perform yourself the changes of variable as we have
1040: done above to reduce the integral to a number of integrals on line segments.
1041: \smallskip
1042: %
1043: The example above also shows the use of the \kbd{solve} function. To use
1044: \kbd{solve} on functions of a complex variable, it is necessary to reduce the
1045: problem to a real one. For example, to find the first complex zero of the
1046: Riemann zeta function as above, we could try typing
1047:
1048: \kbd{solve(t=14,15, real(zeta(1/2 + I*t)))},
1049:
1050: \noindent but this would not work because the real part is positive for
1051: \kbd{t=14} and \kbd{t=15}. As it happens, the imaginary part works. Type
1052:
1053: \kbd{solve(t=14,15, imag(zeta(1/2 + I*t)))},
1054:
1055: \noindent and this now works. We could also narrow the search interval and
1056: type for instance
1057:
1058: \kbd{solve(t=14,14.2, real(zeta(1/2 + I*t)))}
1059:
1060: \noindent which would also work.
1061:
1062: \section{Functions Related to Polynomials and Power Series}
1063:
1064: First a word of warning to the unwary: it is essential to understand the
1065: crucial difference between exact and inexact objects (see Section~1). This
1066: is especially important in the case of polynomials. Let's immediately take
1067: a plunge into these problems. Type
1068:
1069: \centerline{\tt gcd(x\pow 2 - 1, x\pow 2 - 3*x + 2)}
1070:
1071: No problem, the result is \kbd{x - 1} as expected. But now type
1072:
1073: \centerline{\tt gcd(x\pow 2 - 1., x\pow2 - 3.*x + 2.)}
1074:
1075: You are lucky, the result is almost correct except for a bizarre factor of
1076: 3 which comes from the way PARI does the computation. In any case, it is still
1077: essentially a reasonable result. But now type
1078: \kbd{gcd(x - Pi, x\pow 2 - 6*zeta(2))}.
1079: Although this should be equal to \kbd{x - Pi}, PARI finds a
1080: constant as a result. This is because the notion of GCD of non-exact
1081: polynomials doesn't make much sense. However, type
1082: \kbd{polresultant(x - Pi, x\pow 2 - 6*zeta(2))}.
1083: The result is extremely close to zero, showing that indeed the GCD is
1084: non-trivial, without telling us what it is. This being said, we will usually
1085: use polynomials (and power series) with exact coefficients in our
1086: examples.\smallskip
1087:
1088: Set \kbd{pol = polcyclo(15)}. This gives the $15$-th cyclotomic polynomial,
1089: which is of degree $\varphi(15)=8$. Now, type \kbd{r = polroots(pol)}. You
1090: have the 8 complex roots of pol given to 28 significant digits. To see them
1091: better, type \b{b}. As you see, they are given as pairs of complex conjugate
1092: roots, in a random order. The only ordering done by the function
1093: \kbd{polroots} concerns the real roots, which are given first, and in
1094: increasing order.
1095:
1096: The roots of \kbd{pol} are by definition the primitive $15$-th roots of unity.
1097: To check this, simply type \kbd{rc = r\pow 15}. Why, we get an error message!
1098: Well, fair enough, vectors cannot be multiplied (even less raised to a power)
1099: that easily. However, type \kbd{rc = r\pow 15.} with a $.$ at the end. Now it
1100: works, because powering to a non-integer (here real) exponent is a
1101: transcendental function and hence is applied termwise. Note that the fact that
1102: $15.$ is a real number which is representable exactly as an integer has
1103: nothing to do with the problem.
1104:
1105: We see that the components of the result are very close to 1. It is however
1106: tedious to look at all these real and imaginary parts. It would be impossible
1107: if we had many more. Let's do it automatically. Type
1108: \kbd{rr = round(real(rc))}, then \kbd{sqrt(norml2(rc-rr))}. We see that
1109: \kbd{rr} is indeed all 1's, and that the L2-norm of \kbd{rc - rr} is around
1110: $10^{-27}$, reasonable enough when we work with 28 significant digits. Note
1111: that the function \kbd{norml2}, contrary to what it may seem, does not give the
1112: L2 norm but its square, hence we must take the square root (well, this is not
1113: absolutely necessary in this case!).
1114: %
1115: \smallskip
1116: Now type \kbd{pol = x\pow 5 + x\pow 4 + 2*x\pow 3 - 2*x\pow 2 - 4*x - 3},
1117: then \kbd{factor(pol)}. The polynomial \kbd{pol} factors over $\Q$ (or $\Z$)
1118: as a product of two factors. Type \kbd{fun(p)= factorpadic(pol,p,10)}. This
1119: creates a function \kbd{fun(p)} which factors \kbd{pol} over $\Q_p$ to $p$-adic
1120: precision 10. If now we type \kbd{factor(poldisc(pol))}, we learn that the
1121: primes dividing the discriminant are $11$, $23$ and $37$. Type \kbd{fun(5)},
1122: \kbd{fun(11)}, \kbd{fun(23)}, and \kbd{fun(37)} to see different splittings.
1123:
1124: Similarly, we can type \kbd{lf(p)= lift(factormod(pol,p))}, and
1125: \kbd{lf(2)}, \kbd{lf(11)}, \kbd{lf(23)} and \kbd{lf(37)} which show the
1126: different factorizations, this time over $\F_p$. In fact, even better: type
1127: successively
1128: \bprog
1129: pol2 = x^3 + x^2 + x - 1
1130: fq = factorff(pol2, 3, t^3 + t^2 + t - 1)
1131: centerlift(lift(fq))
1132: @eprog
1133: %
1134: This factors the polynomial \kbd{pol2} over the finite field $\F_3(\theta)$
1135: where $\theta$ is a root of the polynomial \kbd{t\pow 3 + t\pow 2 + t - 1}.
1136: This is of course a form of the field $\F_{27}$. We know that
1137: Gal$(\F_{27}/\F_3)$ is cyclic of order 3 generated by the Frobenius
1138: homomorphism $u\mapsto u^3$, and the roots that we have found give the action
1139: of the powers of the Frobenius on \kbd{t} (if you do not know what I am
1140: talking about, please try some more examples, it's not so hard to figure out).
1141:
1142: Similarly, type \kbd{pol3 = x\pow 4 - 4*x\pow 2 + 16} and
1143: \kbd{fn = factornf(pol3,t\pow 2 + 1)}, and we get the factorization of the
1144: polynomial \kbd{pol3} over the number field defined by \kbd{t\pow 2 + 1},
1145: i.e.~over $\Q(i)$. To see the result even better, type \kbd{lift(fn)},
1146: remembering that \kbd{t} stands for the generator of the number field
1147: (here equal to $i=\sqrt{-1}$).
1148:
1149: Note that it is possible, although ill advised, to use the same variable
1150: for the polynomial and the number field. You may for example type
1151: \kbd{fn2 = factornf(pol3, x\pow 2 + 1)}, and the result is correct. However,
1152: the PARI object thus created may give unreasonable results. For example,
1153: if you type \kbd{lift(fn2)} in the example above, you will get a strange
1154: object, with a symbol such as \kbd{x*x} typed. This is because PARI knows
1155: that the dummy variable \kbd{x} is not the same as the explicit variable
1156: \kbd{x}, but since it must print it when you lift, it has to do something.
1157: \smallskip
1158: %
1159: To summarize, in addition to being able to factor integers, you can
1160: factor polynomials over $\C$ and $\R$ (this is the function \kbd{polroots}),
1161: over $\F_p$ (the function \kbd{factormod}, over $\F_{p^k}$ (the function
1162: \kbd{factorff}), over $\Q_p$ (the function \kbd{factorpadic}), over $\Q$ or
1163: $\Z$ (the function \kbd{factor}), and over number fields (the functions
1164: \kbd{factornf} and \kbd{nffactor}). Note however that \kbd{factor} itself
1165: will try to guess intelligently over which ring you want to factor: set
1166: \kbd{pol = x\pow2+1}
1167: and try to \kbd{factor} successively \kbd{pol}, \kbd{pol * 1.},
1168: \kbd{pol * (1+0.*I)}, \kbd{pol * Mod(1,2)},
1169: \kbd{pol * Mod(1,Mod(1,3)*(t\pow2+1))}.
1170:
1171: In the present version \vers{}, it is {\it not} possible to factor over
1172: other rings than the ones mentioned above, for example GP cannot factor
1173: multivariate polynomials. \smallskip
1174:
1175: Other functions related to factoring are \kbd{padicappr}, \kbd{polrootsmod},
1176: \kbd{polrootspadic}, \kbd{polsturm}. Play with them a little.
1177:
1178: Now let's type \kbd{polsym(pol3,20)}, where \kbd{pol3} is the same
1179: polynomial as above. This gives the sum of the $k$-th powers of the roots
1180: of \kbd{pol3} up to $k=20$, of course computed using Newton's formula and
1181: not using \kbd{polroots}. You notice that every odd sum is zero (this is
1182: trivial since the polynomial is even), but also that the signs follow a
1183: regular pattern and that the (non-zero) absolute values are powers of 2.
1184: This is true: prove it, and more precisely find an explicit formula for the
1185: $k$-th symmetric power not involving (non-rational) algebraic numbers.
1186: \medskip
1187:
1188: Now let's play a little with power series. We have already done so a little
1189: at the beginning. Type
1190:
1191: \centerline{\tt 8*x + prod(n=1,39, if(n\%4, 1 - x\pow n, 1),
1192: 1 + O(x\pow 40))\pow 8}
1193:
1194: You obtain a power series which has apparently only even powers of \kbd{x}
1195: appearing. This is surprising, but can be proved using the theory of modular
1196: forms. Note that we have initialized the product to \kbd{1 + O(x\pow 40)} and
1197: not simply to 1 otherwise the whole computation would have been done with
1198: polynomials, and this would first have been slightly slower and also totally
1199: useless since the coefficients of \kbd{x\pow 40} and above are irrelevant
1200: anyhow if we stop the product at \kbd{n=39}.
1201:
1202: While we are on the subject of modular forms (which, together with taylor
1203: series expansions of common functions, are another great source of power
1204: series), type \kbd{\b{ps} 122} (which is a shortcut for
1205: \kbd{default(seriesprecision, 122)}), then \kbd{d = x * eta(x)\pow 24}. This
1206: gives the first 122 terms of the (Fourier) series expansion of the modular
1207: discriminant function $\Delta$ of Ramanujan, its coefficients giving by
1208: definition the Ramanujan $\tau$ function which has a number of marvelous
1209: properties (look at any book on modular forms for explanations). We would like
1210: to see its properties modulo 2. Type \kbd{d\%2}. Hmm, apparently PARI
1211: doesn't like to
1212: reduce coefficients of power series (or polynomials for that matter) directly.
1213: Can we do it without writing a little program? No problem. Type instead
1214: \kbd{lift(Mod(1,2) * d)} and now this works like a charm.
1215:
1216: The pattern in the result is clear. Of course, it now remains to prove it
1217: (again modular forms, see Antwerp III or your resident modular forms guru).
1218: Similarly, type \kbd{centerlift(Mod(1,3) * d)}. This time the pattern is
1219: less clear, but nonetheless there is one. Refer to Anwerp III again.
1220:
1221: \section{Working with Elliptic Curves}
1222:
1223: Now we are getting to more complicated objects. Just as with number fields
1224: which we will meet later on, the first thing to do is to initialize them.
1225: That's because a lot of data will be needed repeatedly, and it's much more
1226: convenient to have it ready once and for all. Here, this is done with the
1227: function \kbd{ellinit} (try to guess what we'll use for number fields\dots).
1228:
1229: So type \kbd{e = ellinit([6,-3,9,-16,-14])}. This computes a number of things
1230: about the elliptic curve defined by the affine equation
1231: %
1232: $$ y^2+6xy+9y = x^3-3x^2-16x-14\enspace. $$
1233: %
1234: It's not that clear what all these funny numbers mean, except that we
1235: recognize the first few of them as the coefficients we just input. To
1236: retrieve meaningful information from such complicated objects (and number
1237: fields will be much worse), you are advised to use the so-called {\it member
1238: functions}. Type \kbd{?.} to get a complete list. Whenever \kbd{ell} appears
1239: in the right hand side, we can apply the corresponding function to an object
1240: output by \kbd{ellinit} (I'm sure you know how the other \kbd{init} functions
1241: will be called now, don't you? Oh, by the way, there is no \kbd{clgpinit}
1242: function).
1243:
1244: Let's try it. We see that the discriminant \kbd{e.disc} is equal to 37,
1245: hence the conductor of the curve is 37. Of course in general it is not so
1246: trivial. In fact, the equation of the curve is clearly not minimal, so type
1247: \kbd{r = ellglobalred(e)}. The first component \kbd{r[1]} tells us that the
1248: conductor is 37 as we already knew. The second component is a 4-component
1249: vector which will allow us to get the minimal equation: simply type
1250: \kbd{e = ellchangecurve(e, r[2])} and the new \kbd{e} is now our minimal
1251: equation with corresponding data. You can for the moment ignore the third
1252: component \kbd{r[3]} (for the impatient reader, this is the product of the
1253: local Tamagawa numbers, $c_p$).
1254:
1255: The new \kbd{e} tells us that the minimal equation is $y^2+y = x^3-x$.
1256: Let us now play a little with points on \kbd{e}. Type \kbd{q = [0,0]}, which is
1257: clearly on the curve (type \kbd{ellisoncurve(e, q)} to check). Well, \kbd{q}
1258: may be a torsion point. Type \kbd{ellheight(e, q)}, which computes the
1259: canonical Neron-Tate height of \kbd{q}. This is non-zero, hence \kbd{q} is
1260: not torsion. To see this even better, type
1261:
1262: \kbd{for(k = 1, 20, print(ellpow(e, q,k)))}
1263:
1264: \noindent and we see the characteristic parabolic explosion of the size of
1265: the points. As a further check, type
1266: \kbd{ellheight(e, ellpow(e, q,20)) / ellheight(e, q)}. We indeed find
1267: $400=20^2$ as it should be. You can also type \kbd{ellorder(e, q)} which
1268: returns 0, telling you that \kbd{q} is non-torsion.
1269:
1270: Notice how all those \kbd{ell}--prefixed functions take our elliptic curve as
1271: a first argument? This will be true with number fields as well: whatever
1272: object was initialized by an $ob$--\kbd{init} function will have to be used as
1273: a first argument of all the $ob$--prefixed functions. Conversely, you won't be
1274: able to use any such high-level function before you correctly initialize the
1275: relevant object. \smallskip
1276:
1277: Ok, let's try another curve. Type \kbd{e = ellinit([0,-1,1,0,0])}. This
1278: corresponds to the equation $y^2+y = x^3-x^2$. Again from the discriminant
1279: we see that the conductor is equal to 11, and if you type \kbd{ellglobalred(e)}
1280: you will see that the equation for \kbd{e} is minimal. Type \kbd{q = [0,0]}
1281: which is clearly a point on \kbd{e}, and \kbd{ellheight(e, q)}. This time we
1282: obtain a value which is very close to zero, hence \kbd{q} must be a torsion
1283: point. Indeed, typing \kbd{for(k=1,5, print(ellpow(e, q,k)))} we see that
1284: \kbd{q} is a point of order 5 (note that the point at infinity is represented
1285: as \kbd{[0]}). More simply, you can type \kbd{ellorder(e, q)}.\smallskip
1286:
1287: Let's try still another curve. Type \kbd{e = ellinit([0,0,1,-7,6])} to get
1288: the curve $y^2+y=x^3-7x+6$. Typing \kbd{ellglobalred(e)} shows that this is a
1289: minimal equation and that the conductor, equal to the discriminant, is 5077.
1290: There are some trivial integral points on this curve, but let's try to be
1291: more systematic.
1292:
1293: First, let's study the torsion points. Typing \kbd{elltors(e)} shows that the
1294: torsion subgroup is trivial, so we won't have to worry about torsion points.
1295: Next, the member \kbd{e.roots} gives us the 3 roots of the minimal
1296: equation over $\C$, i.e.~$Y^2=X^3-7X+25/4$ (set $(X,Y)=(x,y+1/2)$) so if
1297: $(x,y)$ is a real point on the curve, $x$ must be at least equal to the
1298: smallest root, i.e.~$x\ge-3$. Finally, if $(x,y)$ is on the curve, its
1299: opposite is clearly $(x,-y-1)$. So we are going to use the \kbd{ellordinate}
1300: function and type (for instance in \kbd{points.gp} which you can read in with
1301: \kbd{\b{r} points} as we saw before)
1302: \bprog
1303: {
1304: v=[];
1305: for (x = -3, 1000,
1306: s = ellordinate(e,x);
1307: if (length(s), \\ @com we could use \kbd{if (!s,)} instead
1308: v = concat(v, [[x,s[1]]])
1309: )
1310: ); v
1311: }
1312: @eprog
1313:
1314: \noindent By the way, this is how you insert a comment in a script:
1315: everything following a double backslash (up to the first newline character)
1316: is ignored. If you want comments which span many lines, you can brace them
1317: between \kbd{/* ... */} pairs. Everything in between will be ignored as well.
1318: For instance as a header for the file \kbd{points.gp} you could insert the
1319: following:
1320: \bprog
1321: /* Finds rational points on the elliptic curve e, using the naivest
1322: * algorithm I could think of right now (TO BE IMPROVED).
1323: * e should have rational coefficients.
1324: * TODO: Make that into a usable function.
1325: */
1326: @eprog
1327:
1328: (I hope you did not waste your time copying this nonsense, did you?)
1329:
1330: We thus get a large number (18) of integral points. Together with their
1331: opposites and the point at infinity, this makes a total of 37 integral
1332: points, which is large for a curve having such a small conductor. So we
1333: suspect (if we don't know already, since this curve is quite famous!) that
1334: the rank of this curve must be high. Let's try and put some order into this
1335: (note that we work only with the integral points, but in general rational
1336: points should also be considered).
1337:
1338: Type \kbd{hv = ellheight(e, v)}. This gives the vector of canonical heights.
1339: Let us order the points according to their height. For this, type
1340:
1341: \kbd{iv = vecindexsort(hv)}, then \kbd{hv = vecextract(hv,iv)} and
1342: \kbd{v = vecextract(v,iv)}.
1343:
1344: It seems reasonable to take the numbers with smallest height as generators of
1345: the Mordell-Weil group. Let's try the first 4: type
1346:
1347: \kbd{m = ellheightmatrix(e, vecextract(v,[1,2,3,4])); matdet(m)}
1348:
1349: Since the curve has no torsion, the determinant being close to zero implies
1350: that the first four points are dependent. To find the dependency, it is
1351: enough to find the kernel of the matrix \kbd{m}. So type \kbd{matker(m)}:
1352: we indeed get a non-trivial kernel, and the coefficients are (close to)
1353: integers as they should. Typing \kbd{elladd(e, v[1],v[3])} does indeed show
1354: that it is equal to \kbd{v[4]}.
1355:
1356: Taking any other four points, we would in fact always find a dependency.
1357: Let's find them all. Type \kbd{vp = [v[1],v[2],v[3]]\til;
1358: m = ellheightmatrix(e,vp);
1359: matdet(m)}. This is now clearly non-zero so the first 3 points
1360: are linearly independent, showing that the rank of the curve is at least
1361: equal to 3 (it is in fact equal to 3, and \kbd{e} is the curve of smallest
1362: conductor having rank 3). We would like to see whether the other points are
1363: dependent. For this, we use the function \kbd{ellbil}. Indeed, if \kbd{Q} is
1364: some point which is dependent on \kbd{v[1],v[2]} and \kbd{v[3]}, then
1365: \kbd{matsolve(m, ellbil(e, vp,Q))} will by definition give the coefficients
1366: of the dependence relation. If these coefficients are close to integers, then
1367: there is a dependency, otherwise not. This is much safer than using the
1368: \kbd{matker} function. Thus, type
1369:
1370: \centerline{\tt w = vector(18,k, matsolve(m, ellbil(e, vp,v[k])))}
1371:
1372: We ``see'' that the coefficients are all very close to integers, and we can
1373: prove it by typing
1374:
1375: \centerline{\tt wr = round(w); sqrt(norml2(w - wr))}
1376:
1377: \noindent which gives an upper bound on the maximum distance to an integer.
1378: Thus \kbd{wr} is the vector expressing all the components of \kbd{v} on its
1379: first 3. We are thus led to strongly believe that the curve has rank exactly
1380: 3, and this can be proved to be the case.\smallskip
1381:
1382: Let's now explore a few more elliptic curve related functions. Let us keep
1383: our curve \kbd{e} of rank 3, and type
1384: \bprog
1385: v1 = [1,0]; v2 = [2,0];
1386: z1 = ellpointtoz(e, v1)
1387: z2 = ellpointtoz(e, v2)
1388: @eprog
1389:
1390: We thus get the complex parametrization of the curve. To add the points
1391: \kbd{v1} and \kbd{v2}, we should of course type \kbd{elladd(e, v1,v2)},
1392: but we can also type \kbd{ellztopoint(e, z1 + z2)} which of course has the
1393: disadvantage of giving complex numbers, but illustrates how the group law on
1394: \kbd{e} is obtained from the addition law on $\C$.
1395:
1396: Type \kbd{f = x * Ser(ellan(e, 30))}. This gives a power series which
1397: is the Fourier expansion of a modular form of weight 2 for $\Gamma_0(5077)$
1398: (this has been proved directly, but also follows from Wiles' result since
1399: \kbd{e} is semi-stable). In fact, to find the modular parametrization of
1400: the curve, type \kbd{modul = elltaniyama(e)}, then
1401: \kbd{u=modul[1]; v=modul[2];}. Type
1402:
1403: \centerline{\tt (v\pow 2 + v) - (u\pow 3 - 7*u + 6)}
1404:
1405: \noindent to see that this indeed parametrizes the curve.
1406:
1407: Now type \kbd{x * u' / (2*v + 1)}, and we see that this is equal to the
1408: modular form \kbd{f} found above (the quote \kbd{'} tells GP to take the
1409: derivative of the expression with respect to its main variable). The
1410: functions \kbd{u} and \kbd{v}, considered on the upper half plane (with
1411: $x=e^{2i\pi\tau}$), are in fact modular {\it functions} for $\Gamma_0(5077)$.
1412: \smallskip
1413:
1414: Finally, let us come back to the curve defined by typing
1415: \kbd{e = ellinit([0,-1,1,0,0])} where we had seen that the point
1416: \kbd{q = [0,0]} was of order 5. We know that the conductor of this curve is
1417: equal to 11 (type \kbd{ellglobalred(e)}). We want the sign of the functional
1418: equation. Type
1419:
1420: \centerline{\kbd{elllseries(e, 1,-11,1.1)}, \quad then \quad
1421: \kbd{elllseries(e, 1,-11,1)}.}
1422:
1423: Since the values are clearly different, the sign cannot be $-$. In fact
1424: there is an algebraic algorithm which would allow to compute this sign, but
1425: it has not yet been completely written, although in case of conductors prime
1426: to 6 it is very simple.
1427:
1428: Now type \kbd{ls = elllseries(e, 1,11,1)}, and just as a check
1429: \kbd{elllseries(e, 1,11,1.1)}. The values agree (approximately) as they should,
1430: and give the value of the L-function of \kbd{e} at 1. Now according to
1431: the Birch and Swinnerton-Dyer conjecture (which is proved for this curve),
1432: \kbd{ls} is given by the following formula (in this case):
1433: %
1434: \def\sha{\hbox{III}}
1435: $$L(E,1)=\dfrac{\Omega\cdot c\cdot|\sha|}{|E_{\text{tors}}|^2}\enspace,$$
1436: %
1437: where $\Omega$ is the real period of $E$, $c$ is the global Tamagawa number,
1438: product of the local $c_p$ for primes $p$ dividing the conductor, $|\sha|$ is
1439: the order of the Tate-Shafarevich group, and $E_{\text{tors}}$ is the
1440: torsion group of $E$.
1441:
1442: Now we know many of these quantities: $\Omega$ is equal to \kbd{e.omega[1]}
1443: (if there had been 3 real roots instead of 1 in \kbd{e.roots}, $\Omega$ would
1444: be equal to \kbd{2 * e.omega[1]}). The Tamagawa number $c$ is given as the
1445: last component of \kbd{ellglobalred(e)}, and here is equal to 1. We already
1446: know that the torsion subgroup of $E$ contains a point of order 5, and typing
1447: \kbd{torsell(e)} shows that it is of order exactly 5. Hence type
1448: \kbd{ls * 25/e[15]}. This shows that $|\sha|$ must be equal to 1.
1449:
1450: \section{Working in Quadratic Number Fields}
1451:
1452: The simplest of all number fields outside $\Q$ are quadratic fields. Such
1453: fields are characterized by their discriminant, and even better, any non-square
1454: integer $D$ congruent to 0 or 1 modulo 4 is the discriminant of a specific
1455: order in a quadratic field. We can check whether this order is maximal by
1456: using the command \kbd{isfundamental(D)}. Elements of a quadratic field are
1457: of the form $a+b\omega$, where $\omega$ is chosen as $\sqrt{D}/2$ if $D$ is
1458: even and $(1+\sqrt{D})/2$ if $D$ is odd, and are represented in PARI by
1459: quadratic numbers. To initialize working in a quadratic order, one should
1460: start by the command \kbd{w = quadgen($D$)}.
1461:
1462: This sets \kbd{w} equal to $\omega$ as above, and is printed \kbd{w}. Note
1463: however that if several different quadratic orders are used, a printed \kbd{w}
1464: may have several different meanings. For example if you type
1465:
1466: \centerline{\tt w1 = quadgen(-23); w2 = quadgen(-15);}
1467:
1468: Then ask for the value of \kbd{w1} and \kbd{w2}, both will be printed as
1469: \kbd{w}, but of course they are not equal. Hence beware when dealing with
1470: several quadratic orders at once. \smallskip
1471: %
1472: In addition to elements of a quadratic order, we also want to be able to
1473: handle ideals of such orders. In the quadratic case, it is equivalent to
1474: handling binary quadratic forms, and this has been chosen in PARI. For
1475: negative discriminants, quadratic forms are triples $(a,b,c)$ representing
1476: the form $ax^2+bxy+cy^2$. Such a form will be printed as, and can be created
1477: by, \kbd{Qfb($a$,$b$,$c$)}.
1478:
1479: Such forms can be multiplied, divided, powered as many PARI objects using
1480: the usual operations, and they can also be reduced using the function
1481: \kbd{qfbred} (it is not the purpose of this tutorial to explain what all
1482: these things mean). In addition, Shanks's NUCOMP algorithm has been
1483: implemented (functions \kbd{qfbnucomp} and \kbd{qfbnupow}), and this is
1484: usually a little faster.
1485:
1486: Finally, you have at your disposal the functions \kbd{qfbclassno} which
1487: ({\it usually\/}) gives the class number, the function \kbd{qfbhclassno}
1488: which gives the Hurwitz class number, and the much more sophisticated
1489: \kbd{quadclassunit} function which gives the class number and class group
1490: structure.
1491:
1492: Let us see examples of all this at work.
1493:
1494: Type \kbd{qfbclassno(-10007)}. GP tells us that the result is 77. However,
1495: you may have noticed in the explanation above that the result is only
1496: {\it usually\/} correct. This is because the implementers of the algorithm
1497: have been lazy and have not put the complete Shanks algorithm into PARI,
1498: causing it to fail in certain very rare cases. In practice, it is almost
1499: always correct, and the much more powerful \kbd{quadclassunit} program, which
1500: {\it is} complete (at least for fundamental discriminants) can give
1501: confirmation (but now, under the Generalized Riemann Hypothesis!!!).
1502:
1503: So we may be a little suspicious of this class number. Let us check it.
1504: First, we need to find a quadratic form of discriminant $-10007$. Since this
1505: discriminant is congruent to 1 modulo 8, we know that there is an ideal of
1506: norm equal to 2, i.e.~a binary quadratic form $(a,b,c)$ with $a=2$. To
1507: compute it we type \kbd{f = qfbprimeform(-10007, 2)}. OK, now we have a form.
1508: If the class number is correct, the very least is that this form raised to
1509: the power 77 should equal the identity. Let's check this. Type \kbd{f\pow 77}.
1510: We get a form starting with 1, i.e.~the identity, so this test is OK. Raising
1511: \kbd{f} to the powers 11 and 7 does not give the identity, thus we now know
1512: that the order of \kbd{f} is exactly 77, hence the class number is a multiple
1513: of 77. But how can we be sure that it is exactly 77 and not a proper multiple?
1514: Well, type
1515: \bprog
1516: sqrt(10007)/Pi * prodeuler(p=2,500, 1./(1 - kronecker(-10007,p)/p))
1517: @eprog
1518: %
1519: This is nothing else than an approximation to the Dirichlet class number
1520: formula. The function \kbd{kronecker} is the Kronecker symbol, in this case
1521: simply the Legendre symbol. Note also that we have written \kbd{1./(1 - \dots)}
1522: with a dot after the first 1. Otherwise, PARI may want to compute the whole
1523: thing as a rational number, which would be terribly long and useless. In fact
1524: PARI does no such thing in this particular case (\kbd{prodeuler} is always
1525: computed as a real number), but you never know. Better safe than sorry!
1526:
1527: We find 77.77, pretty close to 77, so things seem in order. Explicit bounds
1528: on the prime limit to be used in the Euler product can be given which make
1529: the above reasoning rigorous.
1530:
1531: Let us try the same thing with $D=-3299$. \kbd{qfbclassno} and the Euler
1532: product convince us that the class number must be 27. However, we get stuck
1533: when we try to prove this in the simple-minded way above. Indeed, we type
1534: \kbd{f = qfbprimeform(-3299, 3)} (here 2 is not the norm of a prime ideal but
1535: 3 is), and we see that \kbd{f} raised to the power 9 is equal to the identity.
1536: This is the case for any other quadratic form we choose. So we suspect that the
1537: class group is not cyclic. Indeed, if we list all 9 distinct powers of \kbd{f},
1538: we see that \kbd{qfbprimeform(-3299, 5)} is not on the list (although its cube
1539: is as it must). This implies that the class group is probably equal to a
1540: product of a cyclic group of order 9 by a cyclic group of order 3. The Euler
1541: product plus explicit bounds prove this.
1542:
1543: Another way to check it is to use the \kbd{quadclassunit} function by typing
1544: for example
1545:
1546: \centerline{\tt quadclassunit(-3299)}
1547:
1548: Note that this function cheats a little and could still give a wrong answer,
1549: even assuming GRH (we could get a subgroup and not the whole class group).
1550: If we want to use proven bounds under GRH, we have to type
1551:
1552: \centerline{\tt quadclassunit(-3299,,[1,6])}
1553:
1554: The double comma \kbd{,,} is not a typo, it means we omit an optional second
1555: argument (we would use it to compute the narrow class group, which would be
1556: the same here of course). As we want to use the optional {\it third}
1557: argument, we have to indicate to GP we skipped this one.
1558:
1559: Now, if we believe in GRH, the class group is as we thought (see Chapter 3
1560: for a complete description of this function).
1561:
1562: Note that using the even more general function \kbd{bnfinit} (which handles
1563: general number fields and gives much more complicated results), we could
1564: {\it certify\/} this result (remove the GRH assumption). Let's do it, type
1565: \bprog
1566: bnf = bnfinit(x^2 + 3299); bnfcertify(bnf)
1567: @eprog
1568:
1569: A non-zero result (here 1) means that everything is ok. Good, but what did
1570: we certify after all? Let's have a look at this \kbd{bnf} (just type it!).
1571: Enlightening, isn't it? Recall that the \kbd{init} functions (we've already
1572: seen \kbd{ellinit}) store all kind of technical information which you
1573: certainly don't care about, but which will be put to good use by some higher
1574: level functions. That's why \kbd{bnfcertify} could not be used on the output
1575: of \kbd{quadclassunit}: it needs much more data.
1576:
1577: To extract sensible information from such complicated objects, you must use
1578: one of the many {\it member functions} (remember: \kbd{?.} to get a complete
1579: list). In this case \kbd{bnf.clgp} which extracts the class group structure.
1580: This is much better. Type \kbd{\%.no} to check that this leading 27 is indeed
1581: what we think it is and not some stupid technical parameter. Note that
1582: \kbd{bnf.clgp.no} would work just as well, or even \kbd{bnf.no}!
1583:
1584: As a last check, we can request a relative equation for the Hilbert class
1585: field of $\Q(\sqrt{-3299})$: type \kbd{quadhilbert(-3299)}. It is indeed of
1586: degree 27 so everything fits together.
1587:
1588: \medskip
1589: %
1590: Working in real quadratic fields instead of complex ones, i.e.~with $D>0$, is
1591: not very different.
1592:
1593: The same \kbd{quadgen} function is used to create elements. Ideals are again
1594: represented by binary quadratic forms $(a,b,c)$, this time indefinite. However,
1595: the Archimedean valuations of the number field start to come into play (as
1596: is clear if one considers ideles instead of ideals), hence in fact quadratic
1597: forms with positive discriminant will be represented as a quadruplet
1598: $(a,b,c,d)$ where the quadratic form itself is $ax^2+bxy+cy^2$ with $a$,
1599: $b$ and $c$ integral, and $d$ is the Archimedean component, a real number.
1600: For people familiar with the notion, $d$ represents a ``distance'' as defined
1601: by Shanks and Lenstra.
1602:
1603: To create such forms, one uses the same function as for definite ones, but
1604: you can add a fourth (optional) argument to initialize the distance:
1605:
1606: \centerline{\tt Qfb($a$, $b$, $c$, $d$)}
1607:
1608: If the discriminant of $(a,b,c)$ is negative, $d$ is silently
1609: discarded. If you omit it, this component is set to \kbd{0.} (i.e.~a real zero
1610: to the current precision).
1611:
1612: Again these forms can be multiplied, divided, powered, and they can be
1613: reduced using the function \kbd{qfbred}. This function is in fact a
1614: succession of elementary reduction steps corresponding essentially to a
1615: continued fraction expansion, and a single one of these steps can be achieved
1616: by adding an (optional) flag to the arguments of using this function. Since
1617: handling the fourth component $d$ usually involves computing logarithms, the
1618: same flag may be used to ignore the fourth component. Finally, it is
1619: sometimes useful to operate on forms of positive discriminant without
1620: performing any reduction (this is useless in the negative case), the
1621: functions \kbd{qfbcompraw} and \kbd{qfbpowraw} do exactly that.
1622:
1623: Again, the function \kbd{qfbprimeform} gives a prime form, but the form which
1624: is given corresponds to an ideal of prime norm which is usually not reduced.
1625: If desired, it can be reduced using \kbd{qfbred}.
1626:
1627: Finally, you still have at your disposal the function \kbd{qfbclassno} which
1628: gives the class number (this time {\it guaranteed\/} correct),
1629: \kbd{quadregulator} which gives the regulator, and the much more sophisticated
1630: \kbd{quadclassunit} giving the class group's structure and its generators,
1631: as well as the regulator. The \kbd{qfbclassno} and \kbd{quadregulator}
1632: functions use an algorithm which is $O(\sqrt D)$, hence become very slow for
1633: discriminants of more than 10 digits. \kbd{quadclassunit} can be used on a
1634: much larger range.
1635:
1636: Let us see examples of all this at work and learn some little known number
1637: theory at the same time. First of all, type
1638: \kbd{d = 3 * 3299; qfbclassno(d)}. We see that the class number is 3 (we know
1639: in advance that it must be divisible by 3 from the \kbd{d = -3299} case above
1640: and Scholz's theorem). Let us create a form by typing
1641: \kbd{f = qfbred(qfbprimeform(d,2), 2)} (the last 2 tells \kbd{qfbred} to
1642: ignore the archimedean component). This gives us a prime ideal of norm
1643: equal to 2. Is this ideal principal? Well, one way to check this, which is
1644: not the most efficient but will suffice for now, is to look at the complete
1645: cycle of reduced forms equivalent to \kbd{f}. Type
1646: \bprog
1647: g = f; for(i=1,20, g = qfbred(g, 3); print(g))
1648: @eprog\noindent
1649: (this time the 3 means to do a single reduction step, still not using
1650: Shanks's distance). We see that we come back to the form \kbd{f} without
1651: having the principal form (starting with $\pm1$) in the cycle, so the ideal
1652: corresponding to \kbd{f} is not principal.
1653:
1654: Since the class number is equal to 3, we know however that \kbd{f\pow 3} will
1655: be a principal ideal $\alpha\Z_K$. How do we find $\alpha$? For this, type
1656: \kbd{f3 = qfbpowraw(f, 3)}. This computes the cube of \kbd{f}, without
1657: reducing it. Hence it corresponds to an ideal of norm equal to $8=2^3$, so we
1658: already know that the norm of $\alpha$ is equal to $\pm8$. We need more
1659: information, and this will be given by the fourth component of the form.
1660: Reduce your form until you reach the unit form (you will have to type
1661: \kbd{qfbred(\%,~1)} exactly 6 times).
1662:
1663: Extract the Archimedean component by typing \kbd{c = component(\%, 4)}. By
1664: definition of this distance, we know that
1665: $${\alpha\over{\sigma(\alpha)}}=\pm e^{2c},$$
1666: where $\sigma$ denotes real conjugation in our quadratic field. Thus, if we
1667: type
1668:
1669: \centerline{\tt a = sqrt(8 * exp(2*c))}
1670:
1671: \noindent and then \kbd{sa = 8 / a}, we know that up to sign, \kbd{a} and
1672: \kbd{sa} are numerical approximations of $\alpha$ and $\sigma(\alpha)$. Of
1673: course, $\alpha$ can always be chosen to be positive, and a quick numerical
1674: check shows that the difference of \kbd{a} and \kbd{sa} is close to an
1675: integer, and not the sum, so that in fact the norm of $\alpha$ is equal to
1676: $-8$ and the numerical approximation to $\sigma(\alpha)$ is \kbd{$-$sa}. Thus
1677: we type
1678:
1679: \centerline{\tt p = x\pow 2 - round(a-sa)*x - 8}
1680:
1681: \noindent and this is the characteristic polynomial of $\alpha$. We can check
1682: that the discriminant of this polynomial is a square multiple of \kbd{d}, so
1683: $\alpha$ is indeed in our field. More precisely, solving for $\alpha$ and
1684: using the numerical approximation that we have to resolve the sign ambiguity in
1685: the square root, we get explicitly $\alpha=(15221+153\sqrt d)/2$. Note that
1686: this can also be done automatically using the functions \kbd{polred} and
1687: \kbd{modreverse}, as we will see later in the general number field case, or by
1688: solving a system of 2 linear equations in 2 variables.
1689:
1690: \noindent{\bf Exercise:} now that we have $\alpha$ explicitly, check that it
1691: is indeed a generator of the ideal corresponding to the form \kbd{f3}.
1692:
1693: \medskip Let us now play a little with cycles. Set \kbd{D = 10\pow 7 + 1},
1694: then type
1695:
1696: \centerline{\tt quadclassunit(D,,[1,6])}
1697:
1698: We get as a result a 5-component vector, which tells us that (under GRH) the
1699: class number is equal to 1, and the regulator is approximately
1700: equal to $2641.5$. If you want to certify this, use \kbd{qfbclassno} and
1701: \kbd{quadregulator}, {\it not} \kbd{bnfinit} and \kbd{bnfcertify}, which will
1702: take an absurdly long time (well, about 5 minutes if you are careful and set
1703: the initial precision correctly). Indeed \kbd{bnfcertify} needs the fundamental
1704: unit which is so large that \kbd{bnfinit} will have a (relatively) hard time
1705: computing it: you will need about $R/\log(10)\approx 1147$ digits of precision!
1706: On the other hand, you can try \kbd{quadunit(D)}. Impressive, isn't it? (you
1707: can check that its logarithm is indeed equal to the regulator).
1708:
1709: Now just as an example, let's assume that we want the regulator to 500
1710: decimals, say (without cheating and computing the fundamental unit exactly
1711: first!). I claim that by simply knowing the crude approximation above, this
1712: can be computed with no effort.
1713:
1714: This time, we want to start with the unit form. Since \kbd{D} is odd, we can
1715: type:
1716:
1717: \centerline{\tt u = qfbred(Qfb(1,1,(1 - D)/4), 2)}
1718:
1719: We use the function \kbd{qfbred} with no distance since we want the initial
1720: distance to be equal to~0.
1721:
1722: Now we type \kbd{f = qfbred(u, 1)}. This is the first form encountered along
1723: the principal cycle. For the moment, keep the precision low, for example the
1724: initial default precision. The distance from the identity of \kbd{f} is
1725: around 4.253. Very crudely, since we want a distance of $2641.5$, this should
1726: be encountered approximately at $2641.5/4.253=621$ times the distance of
1727: \kbd{f}. Hence, as a first try, we type \kbd{f\pow 621}. Oops, we overshot,
1728: since the distance is now $3173.02$. Now we can refine our initial estimate and
1729: believe that we should be close to the correct distance if we raise \kbd{f} to
1730: the power $621*2641.5/3173$ which is close to $517$. Now if we compute
1731: \kbd{f\pow 517} we hit the principal form right on the dot. Note that this is
1732: not a lucky accident: we will always land extremely close to the correct target
1733: using this method, and usually at most one reduction correction step is
1734: necessary. Of course, only the distance component can tell us where we are
1735: along the cycle.
1736:
1737: Up to now, we have only worked to low precision. The goal was to obtain this
1738: unknown integer $517$. Note that this number has absolutely no mathematical
1739: significance: indeed the notion of reduction of a form with positive
1740: discriminant is not well defined since there are usually many reduced forms
1741: equivalent to a given form. However, when PARI makes its computations, the
1742: specific order and reductions that it performs are dictated entirely by the
1743: coefficients of the quadratic form itself, and not by the distance component,
1744: hence the precision used has no effect.
1745:
1746: Hence we now start again by setting the precision to (for example) 500,
1747: we retype the definition of \kbd{u} (why is this necessary?), and then
1748: \kbd{f = qfbred(u, 1)} and finally \kbd{f\pow 517}. Of course we know in
1749: advance that we land on the unit form, and the fourth component gives us the
1750: regulator to 500 decimal places with no effort at all.
1751:
1752: In a similar way, we could obtain the so-called {\it compact representation}
1753: of the fundamental unit itself, or $p$-adic regulators. I leave this as
1754: exercises for the interested reader.
1755:
1756: You can try the \kbd{quadhilbert} function on that field but, since the class
1757: number is $1$, the result won't be that exciting. If you try it on our
1758: preceding example ($3*3299$) it should take about five minutes (time for a
1759: coffee break?).
1760:
1761: \section{Working in General Number Fields}
1762:
1763: Note for the present release: this section is a little obsolete since many
1764: new powerful functions are available now. This needs to be rewritten
1765: entirely. \smallskip
1766:
1767: The situation here is of course more difficult. First of all, remembering
1768: what we did with elliptic curves, we need to initialize it (with something
1769: more serious than \kbd{quadgen}). For example assume that we want to work in
1770: the number field $K$ defined by one of the roots of the equation
1771: $x^4+24x^2+585x+1791=0$. This is done by typing
1772: \bprog
1773: T = x^4 + 24*x^2 + 585*x + 1791
1774: nf = nfinit(T)
1775: @eprog
1776:
1777: We get quite a complicated object but, thanks to member functions, we don't
1778: need to know anything about its internal structure (which is dutifully
1779: described in Chapter~3). If you type \kbd{nf.pol}, you will get the
1780: polynomial \kbd{T} which you just input. \kbd{nf.sign} yields the signature
1781: $(r_1,r_2)$ of the field, \kbd{nf.disc} the field discriminant, \kbd{nf.zk}
1782: an integral basis, etc\dots.
1783:
1784: The integral basis is expressed in terms of a generic root \kbd{x} of \kbd{T}
1785: and we notice it's very far from being a power integral basis, which is a
1786: little strange for such a small field. Hum, let's check that: \kbd{poldisc(T)}?
1787: Ooops, small wonder we had such denominators, the index is, well, type
1788: \kbd{sqrt(\% / nf.disc)}. That's $3087$, we don't want to work with such
1789: a badly skewed polynomial!
1790:
1791: So, type \kbd{P = polred(T)}. We see from the third component that the
1792: polynomial $x^4-x^3-21x^2+17x+133$ defines the same field with much smaller
1793: coefficients, so type \kbd{A = P[3]}. The \kbd{polred} function gives a
1794: (usually) simpler polynomial, and also sometimes some information on the
1795: existence of subfields. For example in this case, the second component of
1796: \kbd{polred} tells us that the field defined by $x^2-x+1=0$, i.e.~the field
1797: generated by the cube roots of unity, is a subfield of our number field $K$.
1798: Note this is given as incidental information and that the list of subfields
1799: found in this way is usually far from complete. To get the complete list, you
1800: will have to use the function \kbd{nfsubfields} (we'll do that later on).
1801:
1802: Type \kbd{poldisc(A)}, this is much better, but maybe not optimal yet
1803: (the index is still $7$). Type \kbd{polredabs(A)} (the \kbd{abs} stands for
1804: absolute). Since it seems that we won't get anything better, we'll stick with
1805: \kbd{A} (note however that \kbd{polredabs} finds a smallest generating
1806: polynomial with respect to a bizarre norm which ensures that the index will
1807: be small, but not necessarily minimal). In fact, had you typed
1808: \kbd{nfinit(T, 3)}, \kbd{nfinit} would first have tried to find a good
1809: polynomial defining the same field (i.e.~one with small index) before
1810: proceeding.
1811:
1812: It's not too late, let's redefine our number field: \kbd{NF = nfinit(nf, 3)}.
1813: The output is a two-component vector. The first component is the new
1814: \kbd{nf} (type \kbd{nf = NF[1];}). If you type \kbd{nf.pol}, you notice that GP
1815: indeed replaced your bad polynomial \kbd{T} by a much better one, which
1816: happens to be \kbd{A} (small wonder, \kbd{nfinit} internally called
1817: \kbd{polredabs}!). The second component enables you to switch conveniently to
1818: our new polynomial.
1819:
1820: Namely, call $\theta$ a root of our initial polynomial \kbd{T}, and $\alpha$
1821: a root of the one that \kbd{polred} has found, namely \kbd{A}. These are
1822: algebraic numbers, and as already mentioned are represented as polmods. For
1823: example, in our special case $\theta$ is equal to the polmod
1824:
1825: \centerline{\tt Mod(x, x\pow 4 + 24*x\pow 2 + 585*x + 1791)}
1826:
1827: \noindent while $\alpha$ is equal to the polmod
1828:
1829: \centerline{\tt Mod(x, x\pow 4 - x\pow 3 - 21*x\pow 2 + 17*x + 133)}
1830:
1831: \noindent Here of course we are considering only the algebraic aspect, and
1832: hence ignore completely {\it which} root $\theta$ or $\alpha$ is chosen.
1833:
1834: Now probably you may have a number of elements of your number field which are
1835: expressed as polmods with respect to your old polynomial, i.e.~as
1836: polynomials in $\theta$. Since we are now going to work with $\alpha$
1837: instead, it is necessary to convert these numbers to a representation using
1838: $\alpha$. This is what the second component of \kbd{NF} is for: type
1839: \kbd{NF[2]}, you get
1840:
1841: \centerline{\tt Mod(x\pow 2 + x - 11, x\pow 4 - x\pow 3 - 21*x\pow 2 +
1842: 17*x + 133)}
1843:
1844: \noindent meaning that $\theta = \alpha^2+\alpha-11$, and hence the conversion
1845: from a polynomial in $\theta$ to one in $\alpha$ is easy, using \kbd{subst}
1846: (we could get this polynomial from \kbd{polred} as well, try
1847: \kbd{polred(T, 2)}). If we want to do the reverse, i.e.~go back from a
1848: representation in $\alpha$ to a representation in $\theta$, we use the
1849: function \kbd{modreverse} on this polynomial \kbd{NF[2]}. Try it. The result
1850: has a big denominator (147) essentially because our initial polynomial
1851: \kbd{T} was so bad. By the way to get the 147, you should type
1852: \kbd{denominator(content(NF[2]))}. Trying \kbd{denominator} by itself would not
1853: work since the denominator of a polynomial is defined to be 1 (and its
1854: numerator is itself). The reason for this ``surprising'' behaviour is that we
1855: think of a polynomial as a special case of a rational function. \smallskip
1856:
1857: From now on, we completely forget about \kbd{T}, and use only the polynomial
1858: \kbd{A} defining $\alpha$, and the components of the vector \kbd{nf} which
1859: gives information on our number field $K$. Type
1860:
1861: \centerline{\tt u = Mod(x\pow 3 - 5*x\pow 2 - 8*x + 56, A) / 7}
1862:
1863: This is an element in $K$. There are three essentially equivalent
1864: representations for number field elements: polmod, polynomial, and column
1865: vector giving a decomposition in the integral basis \kbd{nf.zk} ({\it not} on
1866: the power basis $(1,x,x^2,\dots)$). All three are equally valid when the
1867: number field is understood (is given as first argument to the function).
1868: You will be able to use any one of them as long as the function you call
1869: requires an \kbd{nf} argument as well. However, most PARI functions will
1870: return elements as column vectors.
1871:
1872: It's a very important feature of number theoretic functions that, although
1873: they may have a preferred format for input, they will accept a wealth of
1874: other different formats. We already saw this for \kbd{nfinit} which
1875: accepts either a polynomial or an \kbd{nf}. It will be true for ideals,
1876: ideles, congruence subgroups, etc.
1877:
1878: Ok, let's stick with elements for the time being. How does one go from one
1879: representation to the other? Between polynomials and polmods, it's easy:
1880: \kbd{lift} and \kbd{Mod} will do the job. Next, from polmods/polynomials to
1881: column vectors: type \kbd{v = nfalgtobasis(nf, u)}. So $\kbd{u} = \alpha^3-
1882: \alpha^2 - \alpha + 8$, right? Wrong! The coordinates of \kbd{u} are given
1883: with respect to the {\it integral basis}, not the power basis
1884: $(1,\alpha,\alpha^2,\alpha^3)$ (and they don't coincide, type \kbd{nf.zk} if
1885: you forgot what the integral basis looked like). As a polynomial in $\alpha$,
1886: we simply have $\kbd{u} = {1\over7}\alpha^3 -
1887: {5\over7}\alpha^2-{8\over7}\alpha+8$, which is trivially deduced from the
1888: original polmod representation!
1889:
1890: Of course \kbd{v = nfalgtobasis(nf, lift(u))} would work equally well. Indeed
1891: we don't need the polmod information since \kbd{nf} already provides the
1892: defining polynomial. To go back to polmod representation, use
1893: \kbd{nfbasistoalg(nf, v)}. Notice that \kbd{u} is in fact an integer since
1894: \kbd{v} has integer coordinates (try \kbd{denominator(v) == 1}, which is of
1895: course overkill here, but not so in a program).
1896:
1897: Let's try this out. We may for instance compute \kbd{u\pow 3}. Try it. Or, we
1898: can type \kbd{1/u}. Better yet, if we want to know the norm from $K$ to $\Q$
1899: of \kbd{u}, we type \kbd{norm(u)} (what else?). \kbd{trace(u)} works as well.
1900: Notice that none of this would work on polynomials or column vectors since
1901: you don't have the opportunity to supply \kbd{nf}! But we could use
1902: \kbd{nfeltpow(nf,u,3)}, \kbd{nfeltdiv(nf,1,u)} (or \kbd{nfeltpow(nf,u,-1)})
1903: which would work whatever representation was chosen. There is no
1904: \kbd{nfeltnorm} function (\kbd{nfelttrace} does not exist either), but we can
1905: easily write one:
1906: \bprog
1907: nfeltnorm(nf,u) =
1908: {
1909: local(t);
1910: t = type(u);
1911: if (t != "t_POLMOD",
1912: if (t == "t_COL",
1913: u = nfbasistoalg(nf, u)
1914: ,
1915: u = Mod(u, nf.pol)
1916: )
1917: );
1918: norm(u)
1919: }
1920: @eprog
1921:
1922: Notice that this is certainly not foolproof (try it with complex or quadratic
1923: arguments!), but we could refine it if the need arose. In fact there was no
1924: need for this function, since you can consider (\kbd{u}) as a principal
1925: ideal, and just type \kbd{idealnorm(nf,u)} whatever the chosen representation
1926: for \kbd{u}. We'll talk about ideals later on.
1927:
1928: If we want all the symmetric functions of \kbd{u} and not only the norm, we
1929: type \kbd{charpoly(u)} (we could write \kbd{charpoly(u, y)} to tell GP to
1930: use the variable \kbd{y} for the characteristic polynomial). Note that this
1931: gives the characteristic polynomial of \kbd{u}, and not in general the
1932: minimal polynomial. Exercises: how does one (easily) find the minimal
1933: polynomial from this? Find a simpler expression for \kbd{u}.\smallskip
1934:
1935: OK, now let's work on the field itself. The \kbd{nfinit} command already
1936: gave us some information. The field is totally complex (its signature
1937: \kbd{nf.sign} is $[0,2]$), its discriminant \kbd{nf.disc} is $D=18981$ and
1938: $(1,\alpha, \alpha^2, {1\over7}\alpha^3+{2\over7}\alpha^2+{6\over7}\alpha)$
1939: is an integral basis (\kbd{nf.zk}). The Galois group of its Galois closure
1940: can be obtained by typing \kbd{polgalois(A)}. The answer ($[8,-1,1]$) shows
1941: that it is equal to $D_4$, the dihedral group with 8 elements, i.e.~the group
1942: of symmetries of a square. \smallskip
1943:
1944: This implies that the field is ``partially Galois'', i.e.~that there exists
1945: at least one non-trivial field isomorphism which fixes $K$ (exactly one in
1946: this case). To find out which it is, we use the function \kbd{nfgaloisconj}.
1947: This uses the LLL algorithm to find linear relations. So type
1948: \kbd{nfgaloisconj(nf)}. The result tells us that, apart from the trivial
1949: automorphism, the map $\alpha \mapsto (-\alpha^3+5\alpha^2+\alpha-49)/7$ (in
1950: the third component) is a field automorphism. Indeed, if we type
1951: \kbd{s = Mod(\%[3], A); charpoly(s)}, we obtain the polynomial \kbd{A} once
1952: again. \smallskip
1953:
1954: The fixed field of this automorphism is going to be the only non-trivial
1955: subfield of $K$. I seem to recall that \kbd{polred} told us this was the
1956: third cyclotomic field. Let's check this: type \kbd{nfsubfields(nf)}. Indeed,
1957: there's a quadratic subfield, but it's given by \kbd{Q = x\pow 2 + 22*x + 133
1958: } and I don't recognize it. Now \kbd{polred(Q)} proves that this subfield is
1959: indeed the field generated by a cube root of unity. Let's check that \kbd{s}
1960: is of order 2: \kbd{subst(lift(s), x, s)}. Yup, it is. Let's express it as a
1961: matrix:
1962: \bprog
1963: {
1964: v = [;]; b = nf.zk;
1965: for (i=1, 4,
1966: v = concat(v, nfalgtobasis(nf, nfgaloisapply(nf, s, b[i])))
1967: )
1968: }
1969: @eprog
1970:
1971: \kbd{v} gives the action of \kbd{s} on the integral basis. Let's check
1972: \kbd{v\pow2}. That's the identity all right. \kbd{k = matker(v-1)} is indeed
1973: two-dimensional, and \kbd{z = nfbasistoalg(nf, k[,2])} generates the
1974: quadratic subfield. Notice that 1, \kbd{z} and \kbd{u} are $\Q$-linearly
1975: dependent (and in fact $\Z$-linearly as well). Exercise: how would you check
1976: these two assertions in general? (Answer: \kbd{concat}, then respectively
1977: \kbd{matrank} or \kbd{matkerint} (or \kbd{qflll})). \kbd{z = charpoly(z)},
1978: \kbd{z = gcd(z,z')} and \kbd{polred(z)} tell us that we found back the same
1979: subfield again (as we ought to!).
1980:
1981: As a final check, type \kbd{nfrootsof1(nf)}. Again we find that $K$ contains
1982: a cube root of unity, since the torsion subgroup of its unit group
1983: is of order 6. And the given generator happens to be equal to \kbd{u} (so if
1984: you did not do the above exercise as you should have, you now know the answer
1985: anyway).
1986: \smallskip
1987:
1988: \noindent{\bf Additional comment} (you're not supposed to skip this anymore,
1989: but do as you wish):
1990:
1991: Before working with ideals, let us note one more thing. The main part of the
1992: work of \kbd{polred} or \kbd{nfinit} is to compute an integral basis, i.e.~a
1993: $\Z$-basis of the maximal order $\Z_K$ of $K$. For a large polynomial, this
1994: implies factoring the discriminant of the polynomial, which is very often out
1995: of the question. There are two ways in which the situation may be improved:
1996:
1997: 1) First, it is often the case that the polynomial that one considers is of
1998: quite a special type, giving some information on the discriminant. For
1999: example, one may know in advance that the discriminant is a square. Hence we
2000: can ``help'' PARI by giving it that information. More precisely, using the
2001: extra information that we have we may be able to factor the discriminant of
2002: the polynomial. We can then use the function \kbd{addprimes} to inform
2003: PARI of this factorization. Namely, add the primes which are known to
2004: divide the discriminant; this will save PARI some work. Do it only for big
2005: primes (bigger than \kbd{primelimit}, whose value you can get using
2006: \kbd{default})~--- it will be useless otherwise.
2007:
2008: 2) The second way in which the situation may be improved is that often we do
2009: not need the complete information on the maximal order, but only require that
2010: the order be $p$-maximal for a certain number of primes $p$ (but then, we
2011: may not be able to use the functions which require a genuine \kbd{nf}). The
2012: function \kbd{nfbasis} specifically computes the integral basis and is not
2013: much quicker than \kbd{nfinit} so is not very useful in its standard use. But
2014: you can provide a factorization of the discriminant as an optional third
2015: argument. And here we can cheat, and give on purpose an incomplete
2016: factorization involving only the primes we want. For example coming
2017: back to our initial polynomial $T$, the discriminant of the polynomial is
2018: $3^7\cdot7^6\cdot19\cdot37$. If we only want a $7$-maximal order, we simply
2019: type
2020:
2021: \centerline{\tt nfbasis(T, ,[7,6; 1537461,1])}
2022:
2023: \noindent and the factors of 1537461 will not be looked at! (of course in
2024: this example it would be stupid to cheat, but if the discriminant has 2000
2025: digits, this can be a handy trick).\medskip
2026:
2027: We now would like to work with ideals (and even with ideles) and not only
2028: with elements. An ideal can be represented in many different ways. First, an
2029: element of the field (in any of the various guises seen above) will be
2030: considered as a principal ideal. Then the standard representation is a
2031: square matrix giving the Hermite Normal Form of a $\Z$-basis of the ideal
2032: expressed on the integral basis. Standard means that most ideal related
2033: functions will use this representation for their output. Note that, as
2034: mentioned before, we always represent elements on the integral basis and not
2035: on a power basis.
2036:
2037: Prime ideals can be represented in a special form as well (see the
2038: description of \kbd{idealprimedec} in Chapter~3) and all ideal-related
2039: functions will accept them. On the other hand, the function \kbd{idealtwoelt}
2040: can be used to find a two-element $\Z_K$-basis of a given ideal (as $a\Z_K +
2041: b\Z_K$, where $a$ and $b$ belong to $K$), but this is {\it not} a valid
2042: representation for an ideal under GP, and most functions will choke on it (or
2043: worse, take it for something else and output a completely meaningless
2044: result). To be able to use such an ideal, you will first have to convert it
2045: to HNF form.
2046:
2047: Whereas it's very easy to go to HNF form (use \kbd{idealhnf(nf,id)} for valid
2048: ideals, or \kbd{idealhnf(nf,a,b)} for a two-element representation as above),
2049: it's a much more complicated problem to check whether an ideal is principal
2050: and find a generator. In fact an \kbd{nf} does not contain enough data for
2051: this particular task. We'll need a Big Number Field (or \kbd{bnf}) for that
2052: (in particular, we need the class group and fundamental units). More on this
2053: later.\smallskip
2054:
2055: An ``idele'' will be represented as a 2-element vector, the first element
2056: being the corresponding ideal (in any valid form), which summarizes the
2057: non-archimedean information, and the second element a vector of real and
2058: complex numbers with $r_1+r_2$ components, the first $r_1$ being real, the
2059: remaining ones complex. In fact, to avoid certain ambiguities, the first
2060: $r_1$ components are allowed to have an imaginary part which is a multiple of
2061: $\pi$. These $r_1+r_2$ components correspond to the Archimedean places of the
2062: number field $K$. \medskip
2063:
2064: Let us keep our number field $K$ as above, and hence the vector \kbd{nf}. Type
2065:
2066: \centerline{\tt P = idealprimedec(nf,7)}
2067:
2068: This gives the decomposition of the prime number 7 into prime ideals. We have
2069: chosen 7 because it is the index of $\Z[\theta]$ in $\Z_K$, hence is the most
2070: difficult case to treat. The index is given in \kbd{nf[4]} and cannot be
2071: accessed through member functions since it's rarely needed. It is in any case
2072: trivial to compute as \kbd{sqrtint(poldisc(nf.pol) / nf.disc)}.
2073:
2074: The result is a vector with 4 components, showing that 7 is totally split in
2075: the field $K$ into four prime ideals of norm 7 (you can check:
2076: \kbd{idealnorm(nf,P[1])}). Let us take one of these ideals, say the first, so
2077: type \kbd{pr = P[1]}. We would like to have the Hermite Normal Form of this
2078: ideal. No problem: since ideal multiplication always gives the result in HNF,
2079: we will simply multiply our prime ideal by $\Z_K$, which is represented by the
2080: identity matrix. So type
2081:
2082: \centerline{\tt idealmul(nf, matid(4), pr)}
2083:
2084: \noindent or in fact simply
2085:
2086: \centerline{\tt idealmul(nf, 1, pr)}
2087:
2088: \noindent or even simpler yet
2089:
2090: \centerline{\tt idealhnf(nf,pr)}
2091:
2092: \noindent and we have the desired HNF. Let's now perform ideal operations.
2093: For example type
2094:
2095: \centerline{\tt idealmul(nf, pr, idealmul(nf, pr,pr))}
2096:
2097: \noindent or more simply
2098:
2099: \centerline{\tt pr3 = idealpow(nf, pr,3)}
2100:
2101: \noindent to get the cube of the ideal \kbd{pr}. Since the norm of this ideal
2102: is equal to $343=7^3$, to check that it is really the cube of \kbd{pr} and
2103: not of other ideals above 7, we can type
2104:
2105: \centerline{\tt for(i=1,4, print(idealval(nf, pr3,P[i])))}
2106:
2107: \noindent and we see that the valuation at \kbd{pr} is equal to 3, while the
2108: others are equal to zero. We could see this as well from
2109: \kbd{idealfactor(nf, pr3)}.
2110:
2111: Let us now ``idelize'' \kbd{pr3} by typing \kbd{id3 = [pr3, [0,0]]}.
2112: (We need $r_1+r_2=2$ components for the second vector.) Then type
2113: \kbd{r1 = idealred(nf, id3)}. We get a new ideal which is equivalent to the
2114: first (modulo the principal ideals). The Archimedean component is non-trivial
2115: and gives the distance from the reduced ideal to the original one (see H. Cohen
2116: {\it A Course in Computational Algebraic Number Theory}, GTM {\bf 138} for
2117: details, especially Sections 5.8.4 and 6.5). Now, just for fun type
2118:
2119: \centerline{\tt r = r1; for(i=1,3, r = idealred(nf,r,[1,5]); print(r))}
2120:
2121: We see that the third \kbd{r} is equal to the initial \kbd{r1}. This means
2122: that we have found a unit in our field, and it is easy to extract this unit
2123: given the Archimedean information (clearly it would be impossible without).
2124: We first form the difference of the Archimedean contributions of \kbd{r} and
2125: \kbd{r1} by typing
2126:
2127: \centerline{\tt arch = r[2] - r1[2]; l1 = arch[1]; l2 = arch[2];}
2128:
2129: From this, we obtain the logarithmic embedding of the unit by typing
2130:
2131: \centerline{\tt l = real(l1 + l2) / 4;
2132: v1 = [l1,l2,conj(l1),conj(l2)]\til / 2 - [l,l,l,l]\til; }
2133:
2134: This 4-component vector contains by definition the logarithms of the
2135: four complex embeddings of the unit. Since the matrix \kbd{nf[5][1]}
2136: contains the values of the $r_1+r_2$ embeddings of the elements of the
2137: integral basis, we can obtain the representation of the unit on the
2138: integral basis by typing
2139: \bprog
2140: {
2141: m1 = nf[5][1];
2142: m = matrix(4,4,j,k,
2143: if (j<=2,
2144: m1[j,k]
2145: ,
2146: conj(m1[j-2, k])
2147: )
2148: );
2149: v = exp(v1);
2150: au = matsolve(m,v);
2151: vu = round(real(au))
2152: }
2153: @eprog
2154:
2155: Then \kbd{vu} is the representation of the unit on the integral basis.
2156: The closeness of the approximation of \kbd{au} to \kbd{vu} gives us
2157: confidence that we have made no numerical mistake. To be sure that \kbd{vu}
2158: represents a unit, type \kbd{u = nfbasistoalg(nf,vu)}, then typing
2159: \kbd{norm(u)} we see that it is equal to 1 so \kbd{u} is a unit.
2160:
2161: There is of course no reason for \kbd{u} to be a fundamental unit. Let us see
2162: if it is a square. Type \kbd{f1 = factor(subst(charpoly(u,x), x, x\pow 2))}.
2163: We see that the characteristic polynomial of \kbd{u} where \kbd{x} is
2164: replaced by \kbd{x\pow 2} is a product of 2 polynomials of degree 4, hence
2165: \kbd{u} is a square (Exercise: why?).
2166:
2167: We now want to find the square root of \kbd{u}. We can again use the
2168: \kbd{matsolve} function as above. For this we need to take the square
2169: root of each element of the vector \kbd{v}, and hence there are
2170: sign ambiguities. Let's do it anyway. Type \kbd{v = sqrt(v)}. We see that
2171: \kbd{v[1]} and \kbd{v[3]} are conjugates, as well as \kbd{v[2]} and \kbd{v[4]},
2172: so for the moment the signs seem OK. Now try \kbd{au = matsolve(m,v)}. The
2173: numbers obtained are clearly not integers, hence the last remaining sign
2174: change must be performed. Type \kbd{v[1] = -v[1]; v[3] = -v[3]} (they must
2175: stay conjugate) and then again \kbd{au = matsolve(m,v)}. This time the
2176: components are close to integers, so we are done (after typing
2177: \kbd{vu = round(real(au))} as before).
2178:
2179: Anyway, we find that a square root \kbd{u2} of \kbd{-u} is represented by
2180: the vector \kbd{vu=[-4,1,1,-1]\til} on the integral basis, and this is in
2181: fact a fundamental unit.\medskip
2182:
2183: The function \kbd{polred} gives us another method to find \kbd{u2} as
2184: follows: type \kbd{q = polred(f1[1,1], 2)}. We recognize the polynomial
2185: \kbd{A} as the component \kbd{q[3,2]}. To obtain the square root of our unit
2186: we then simply type \kbd{up2 = modreverse(Mod(q[3,1], f1[1,1]))}
2187: (Exercise: why?). We find that \kbd{up2} is represented by the vector
2188: \kbd{[-3,-1,0,0]\til} on the integral basis, which is not the result that we
2189: have found before nor its opposite. Where is the error? (Please think about
2190: this before reading on. There is a mathematical subtlety hidden here.)
2191:
2192: Have you solved the problem? Good! The problem occurs because as mentioned
2193: before (but you may not have noticed since it is not stressed in standard
2194: textbooks) although the number field $K$ is not Galois over $\Q$, there does
2195: exist a non-trivial automorphism, and we have found it above by using the
2196: function \kbd{nfgaloisconj}. Indeed, if we apply this automorphism to
2197: \kbd{up2} (by typing \kbd{nfgaloisapply(nf,s,up2)} where \kbd{s} is the
2198: non-trivial component of \kbd{nfgaloisconj(nf)} computed above), we find the
2199: opposite of \kbd{u2}, which is OK. \smallskip
2200:
2201: Still another method which avoids all sign ambiguities and automorphism
2202: problems is as follows. Type \kbd{r = f1[1,1] \% (x\pow 2 - u)} to find the
2203: remainder of the characteristic polynomial of \kbd{u2} divided by
2204: \kbd{x\pow 2 - u}. This will be a polynomial of degree 1 in \kbd{x} (with
2205: polmod coefficients) and we know that \kbd{u2}, being a root of both
2206: polynomials, will be the root of \kbd{r}, hence can be obtained by typing
2207: \kbd{u2 = -coeff(r,0) / coeff(r,1)}. Indeed, we immediately find the correct
2208: result with no trial and error.\smallskip
2209:
2210: Still another method to find the square root of \kbd{u} is to use
2211: \kbd{nffactor(nf,y\pow 2 + u)}. Except that this won't work as is since the
2212: main variable of the polynomial to be factored must have {\it higher}
2213: priority than the number field variable. This won't be possible here since
2214: \kbd{nf} was defined using the variable \kbd{x} which has the highest possible
2215: priority. So we need to substitute variables around (using \kbd{subst}). I
2216: leave you to work out the details.\smallskip
2217:
2218: Now ideals can be used in a wide variety of formats. We have already seen
2219: HNF-representations, ideles and prime ideals. We can also use algebraic
2220: numbers. For example type \kbd{al = Mod(x\pow 2 - 9, A)}, then
2221: \kbd{ideleprincipal(nf,al)}. We obtain the idele corresponding to \kbd{al}
2222: (see the manual for the exact description). However it is usually not
2223: necessary to compute this explicitly since this is done automatically inside
2224: PARI. For example, you can type
2225:
2226: \centerline{\tt for(i=1,4, print(i ": " idealval(nf,al,P[i])))}
2227:
2228: We see that the valuation is non-zero (equal to 1) at the prime ideals
2229: \kbd{P[2]} and \kbd{P[3]}. In addition, typing \kbd{norm(al)} shows that
2230: \kbd{al} is of norm $49=7^2$ (\kbd{idealnorm(nf,al)} gives the same result of
2231: course, and would be more generic). Let's check this differently. Type
2232: \kbd{P23 = idealmul(nf,P[2],P[3])} and then \kbd{idealhnf(nf,al)}. We see that
2233: the results are the same, hence the product of the two prime ideals \kbd{P[2]}
2234: and \kbd{P[3]} is equal to the principal ideal generated by \kbd{al}. There is
2235: still something to be done with this example as we shall see below after we
2236: introduce Big Number Fields (which will trivialize the examples we have just
2237: seen).
2238:
2239: Essentially all functions that you would want on ideals are available.
2240: We mention here the complete list, referring to Chapter 3 for detailed
2241: explanations:
2242:
2243: \kbd{idealadd}, \kbd{idealaddtoone}, \kbd{idealappr}, \kbd{idealchinese},
2244: \kbd{idealcoprime}, \kbd{idealdiv}, \kbd{idealfactor}, \kbd{idealhnf},
2245: \kbd{idealintersect}, \kbd{idealinv}, \kbd{ideallist}, \kbd{ideallog},
2246: \kbd{idealmin}, \kbd{idealmul}, \kbd{idealnorm}, \kbd{idealpow},
2247: \kbd{idealprimedec}, \kbd{idealprincipal}, \kbd{idealred},
2248: \kbd{idealstar}, \kbd{idealtwoelt}, \kbd{idealval}, \kbd{ideleprincipal},
2249: \kbd{nfisideal}.
2250:
2251: We suggest you play with these functions to get a feel for the algebraic
2252: number theory package. Remember simply that when a matrix (usually in Hermite
2253: normal form) is output, it is always a $\Z$-basis of the result expressed on
2254: the {\it integral basis} \kbd{nf.zk} of the number field, which is usually
2255: {\it not} a power basis. \medskip
2256:
2257: Apart from the above functions you have at your disposal the very powerful
2258: functions \kbd{bnfclassunit} which is of the same type as \kbd{quadclassunit}
2259: seen above, but for general number fields, and hence much slower. See Chapter
2260: 3 for a detailed explanation of its use.
2261:
2262: First type \kbd{setrand(1)}: this resets the random seed (to make sure we get
2263: the exact same results). Now type \kbd{m = bnfclassunit(A)}, where \kbd{A} is
2264: the same polynomial as before. After some work, we get a matrix with one
2265: column, which does not look too horrible (this is because much of the really
2266: useful information has been discarded, this is not an \kbd{init} function).
2267: Let's extract the column with \kbd{v = m[,1]}.
2268:
2269: Then \kbd{v} is a vector with 10 components. We immediately recognize the first
2270: component as being the polynomial \kbd{A} once again. In fact member
2271: functions are still available for \kbd{m}, even though it was not created
2272: by an \kbd{init} function. So, let's try them: \kbd{m.pol} gives \kbd{A},
2273: \kbd{m.sign}, \kbd{m.disc}, \kbd{m.zk}, ok nothing really exciting. But new
2274: ones are available now: \kbd{m.no} tells us the class number is 4,
2275: \kbd{m.cyc} that it's cyclic (of order 4 but that we already knew),
2276: \kbd{m.gen} that it's generated by a prime ideal above seven (in HNF form). If
2277: you play a little bit with the \kbd{idealhnf} function you see that this ideal
2278: is \kbd{P[4]}.
2279:
2280: The regulator \kbd{m.reg} is equal to $3.794\dots$. \kbd{m.tu} tells us that
2281: the roots of unity in $K$ are exactly the sixth roots of 1 and gives a
2282: primitive root. Finally \kbd{m.fu} gives us a fundamental unit (which must be
2283: linked to the unit \kbd{u2} found above since the unit rank is 1).
2284:
2285: To find the relation (without trial and error, because in this case it is
2286: quite easy to see it directly!), let us use the logarithmic embeddings. Type
2287: \kbd{uf = m.fu[1]}, \kbd{uu= m.tu[2]} to get the generators of the unit
2288: group, then
2289: \bprog
2290: cu2 = log(conjvec(u2));
2291: cuf = log(conjvec(Mod(uf,A)));
2292: cuu = log(conjvec(Mod(uu,A)));
2293: @eprog\noindent
2294: to get the (complex) logarithmic embeddings. Then type
2295:
2296: \kbd{lindep(real([cu2[1], cuf[1], cuu[1]]))}
2297:
2298: \noindent to find a linear dependence. Unfortunately, the result $[0,0,1]$ is
2299: to be expected since \kbd{uu} is a root of unity, hence the components of
2300: \kbd{cuu} are pure imaginary. Hence we must not take the real part. However,
2301: in that case the logs are defined only up to a multiple of $2i\pi$. Hence we
2302: type
2303:
2304: \kbd{lindep([cu2[1], cuf[1], cuu[1], 2*I*Pi])}
2305:
2306: The result $[1,-1,-3,0]$ shows that \kbd{u2} is probably equal to
2307: \kbd{uf * uu\pow 3}, which is clear since \kbd{uu} is a sixth root of unity,
2308: and which can be checked directly. This works quite nicely, but we'll see
2309: later on a much better method.
2310:
2311: {\bf Note:} since the fundamental unit obtained depends on the random
2312: seed, you could have obtained a different unit from the one given here, had
2313: you not reset the random seed before the computation. This was the purpose
2314: of the initial \kbd{setrand} instruction (which was otherwise completely
2315: unnecessary of course). \medskip
2316:
2317: If you followed closely what component of \kbd{m} the various member
2318: functions were giving (after possibly some cosmetic change), you must have
2319: noticed that the seventh and tenth components were not given. You may ignore
2320: the last component of \kbd{v} (it gives the accuracy to which the fundamental
2321: unit was computed). The seventh is a technical check number, which must be
2322: close to 1 (between $1/\sqrt2$ and $\sqrt2$ at the very least). This tells
2323: us, under LOTS of assumptions (notably the Generalized Riemann Hypothesis),
2324: that our result is correct (i.e.~neither the class group nor the regulator are
2325: smaller). We'll see shortly how to remove all those assumptions: we need much
2326: more data than was output here. This is a little wasteful since the missing
2327: data {\it were} computed, then discarded to provide a human-readable output.
2328: Thus you will usually not use this function, but its \kbd{init} variant:
2329: \kbd{bnfinit}.
2330: \smallskip
2331:
2332: So type: \kbd{setrand(1); bnf = bnfinit(A);}
2333:
2334: This performs exactly the same computations as \kbd{bnfclassunit}, but keeps
2335: much more information. In particular, you don't {\it want} to see the result,
2336: whence the semicolon (if you really want to, you can have a look, it's only
2337: about three screenful). All the member functions that we used previously
2338: still apply to this \kbd{bnf} object. In fact, \kbd{bnfinit} even recomputed
2339: the information provided by \kbd{nfinit} (we could have typed
2340: \kbd{bnfinit(NF)} to avoid the waste), and it's included in the \kbd{bnf}
2341: structure: \kbd{bnf.nf} should be exactly identical to \kbd{NF}. Thus, all
2342: functions which took an \kbd{nf} as first argument, will equally accept a
2343: \kbd{bnf} (and a \kbd{bnr} as well which contains even more data).
2344:
2345: We are now ready to perform more sophisticated operations in the class group.
2346: First and foremost, we can now certify the result: type
2347: \kbd{bnfcertify(bnf)}. The output (\kbd{1} if all went well) is not that
2348: impressive, but it means that we now know the class group and fundamental
2349: units unconditionally (in particular, the GRH assumption could be removed)!
2350: In this case, the certification process takes a very short time, and you
2351: might wonder why it was not built in as a final check in the \kbd{bnfinit}
2352: function. The answer is that as the regulator gets bigger this process gets
2353: increasingly difficult, and becomes soon impractical, while \kbd{bnfinit}
2354: still happily spits out results. So it makes sense to dissociate the two: you
2355: can always make the check afterwards, if the result is interesting enough
2356: (and looking at the tentative regulator, you know in advance whether the
2357: certification can possibly succeed: if \kbd{bnf.reg} is about 2000, don't
2358: waste your time).
2359:
2360: Ok, now that we feel safe about the \kbd{bnf} output, let's do some real
2361: work. For example, let us take again our prime ideal \kbd{pr} above 7. Since
2362: we know that the class group is of order 4, we deduce that \kbd{pr} raised to
2363: the fourth power must be principal. Type \kbd{pr4 = idealpow(nf, pr, 4)} then
2364: \kbd{vis = bnfisprincipal(bnf,pr4)}. The function \kbd{bnfisprincipal} now
2365: uses all the information contained in \kbd{bnf} and tells us that indeed
2366: \kbd{pr4} is principal. More precisely, the first component, \kbd{[0]}, gives
2367: us the factorization of the ideal in the class group. Here, \kbd{[0]} means
2368: that it is up to equivalence equal to the 0-th power of the generator given
2369: in \kbd{bnf.gen}, in other words that it is a principal ideal.
2370:
2371: The second component gives us the algebraic number $\alpha$ such that
2372: \kbd{pr4$=\alpha\Z_K$}, where $\Z_K$ is the ring of integers of our number
2373: field, $\alpha$ being as usual expressed on the integral basis. To get $\alpha$
2374: as an algebraic number, we type \kbd{alpha = nfbasistoalg(bnf,vis[2])} (note
2375: that we can use a \kbd{bnf} with all the \kbd{nf} functions; but not the other
2376: way round, of course).
2377:
2378: Let us check that the result is correct: first, type \kbd{norm(alpha)}
2379: (\kbd{idealnorm(nf, vis[2])} would have worked directly). It is indeed equal to
2380: $7^4 = 2401$, which is the norm of \kbd{pr4} (it could also have been equal to
2381: $-2401$). This is only a first check. The complete check is obtained by
2382: computing the HNF of the principal ideal generated by \kbd{alpha}. To do this,
2383: type \kbd{idealhnf(nf,alpha)}.
2384:
2385: The result that we obtain is identical to \kbd{pr4}, thus showing that
2386: \kbd{alpha} is correct (not that there was any doubt!). You may ignore the
2387: third component of \kbd{vis} which just tells us that the accuracy to which
2388: \kbd{alpha} was computed was amply sufficient.
2389:
2390: But \kbd{bnfisprincipal} also gives us information for non-principal ideals.
2391: For example, type
2392:
2393: \kbd{vit = bnfisprincipal(bnf, pr)}.
2394:
2395: The component \kbd{vit[1]} is now equal to $[3]$, and tells us that \kbd{pr}
2396: is ideal-equivalent to the cube of the generator \kbd{g} given by
2397: \kbd{bnfinit}. Of course we already knew this since the product of \kbd{P[2]}
2398: and \kbd{P[3]} was principal, as well as the product of all the \kbd{P[$i$]}
2399: (generated by 7), and we noticed that \kbd{P[4]} was of order 4 when we looked
2400: at \kbd{bnf.gen}.
2401:
2402: The second component \kbd{vit[2]} gives us $\alpha$ on the integral basis
2403: such that \kbd{pr$=\alpha \kbd{g}^3$}. Note that if you {\it don't\/} want
2404: this $\alpha$, which may be large and whose computation may take some time,
2405: you can just add a flag (1 in this case, see the online help) to the
2406: arguments of the \kbd{bnfisprincipal} function, so that it only returns the
2407: position of \kbd{pr} in the class group. \smallskip
2408:
2409: Let us now take the example of the principal ideal \kbd{P23} that we have
2410: seen above. We know that \kbd{P23} is principal, but of course we have
2411: forgotten the generator that we had found, so we need another one. For this,
2412: we type \kbd{pp = bnfisprincipal(bnf,P23)}. The first component of the result
2413: is \kbd{[0]}, telling us that the ideal is indeed principal, as we knew
2414: already. But the second component gives us the components of a generator of
2415: \kbd{P23} on the integral basis. Now we remember suddenly that we already had
2416: a generator \kbd{al} for the same ideal. Hence type
2417: \kbd{u3=nfeltdiv(nf,al2,al)}. This must be a unit. It's already obviously an
2418: integer and \kbd{idealnorm(nf,\%)} tells us that \kbd{u3} is indeed a unit.
2419: You can again find out what unit it is as we did above. However, as we
2420: mentioned, this is not really the best method.
2421:
2422: To find the unit, we use explicitly the third component of the vector
2423: \kbd{bnf} given by \kbd{bnfinit}. This contains an $(r+1)\times r$ complex
2424: matrix whose columns represent the complex logarithmic embedding of the
2425: fundamental units. Here $r=r_1+r_2-1$ is the unit rank. We first compute
2426: the component of \kbd{u3} on the torsion-free part of the group of units
2427: by proceeding as follows. Type
2428:
2429: \kbd{me = concat(bnf[3],[2,2]\til)}.
2430:
2431: Indeed, this is a variant of the regulator matrix and is more practical to
2432: use since it is more symmetric and avoids suppressing one row arbitrarily.
2433: Now type
2434:
2435: \kbd{cu3 = ideleprincipal(nf,u3)[2]\til}
2436:
2437: to get the complex logarithmic embedding of \kbd{u3} (as a column vector). We
2438: could of course also have computed this logarithmic embeddings directly using
2439: \kbd{conjvec} as we did above, but then we must take care of the factors of 1
2440: and 2 occurring.
2441:
2442: Then type \kbd{xc = matsolve(real(me), real(cu3))}
2443:
2444: Whatever field we are in, if \kbd{u3} is a unit this {\it must} end with a 0
2445: (approximate of course) because of the ``spurious'' vector $[2,2]$, and the
2446: other components (here only one) give the exponents on the fundamental units.
2447: Here the only other component is the first, with a coefficient of $1$ (we
2448: could type \kbd{round(xc)} to tidy up the result). So we know that \kbd{u} is
2449: equal to \kbd{uf} multiplied by a root of unity.
2450:
2451: To find this root of unity, we type \kbd{xd = cu3 - me*xc} then
2452: \kbd{xu = ideleprincipal(nf,uu)[2]} and finally \kbd{xd[1] / xu[1]}. We find
2453: $3$ as a result, so finally our unit \kbd{u3} must be equal to
2454: \kbd{uu\pow 3 * uf} itself, which is the case.
2455:
2456: Of course, you don't need to do all that: just type \kbd{bnfisunit(bnf,u3)}.
2457: Like the \kbd{bnfisprincipal} function, this gives us the decomposition of
2458: some object (here a unit) on the precomputed generators (here \kbd{bnf.tufu})
2459: of some abelian group of finite type (here the units of $K$). The result
2460: \kbd{[1,Mod(3,6)]} tells us that \kbd{u3} is equal to \kbd{uu\pow 3 * uf} as
2461: before.\smallskip
2462:
2463: Another famous so-called {\it discrete logarithm} problem can be easily
2464: solved with PARI, namely the one associated to the invertible elements modulo
2465: an ideal: $(\Z_K / I)^*$. Just use \kbd{idealstar} (this is an \kbd{init}
2466: function) and \kbd{ideallog}.
2467:
2468: ----- TO BE COMPLETED -----
2469:
2470: \section{GP Programming}
2471:
2472: ----- TO BE WRITTEN -----
2473:
2474: \section{Plotting}
2475:
2476: PARI supports a multitude of high and low-level graphing functions, on a
2477: variety of output devices~: a special purpose window under the \kbd{X
2478: Windows} system, a \kbd{PostScript} file ready for the printer, or a
2479: \kbd{gnuplot} output device (only the first two are available by default).
2480: These functions use a multitude of flags, which are mostly power-of-2. To
2481: simplify understanding we first give these flags symbolic names.
2482: \bprog
2483: /* Generic flags: */
2484: parametric = 1; no_x_axis = 8; points = 64;
2485: recursive = 2; no_y_axis = 16; points_lines = 128;
2486: norescale = 4; no_frame = 32; splines = 256;
2487:
2488: /* Relative positioning of graphic objects: */
2489: nw = 0; se = 4; relative = 1;
2490: sw = 2; ne = 6;
2491:
2492: /* String positioning: */
2493: /* V */ bottom = 0; /* H */ left = 0; /* Fine tuning */ hgap = 16;
2494: vcenter = 4; center = 1; vgap = 32;
2495: top = 8; right = 2;
2496: @eprog
2497: We also decrease drastically the default precision.
2498: \bprog
2499: \p 9
2500: @eprog
2501: This is very important, since plotting involves calculation of functions at
2502: a huge number of points, and a relative precision of 28 significant digits
2503: is an obvious overkill: the output device resolution certainly won't reach
2504: $10^{28} \times 10^{28}$ pixels!
2505:
2506: Start with something really simple:
2507: \bprog
2508: ploth(X = -2, 2, sin(X^7))
2509: @eprog
2510:
2511: You can see the limitations of the ``straightforward'' mode of plotting:
2512: while the first several cycles of \kbd{sin} reach $-1$ and $1$, the cycles
2513: which are closer to the left and right border do not. This is understandable,
2514: since PARI is calculating $\sin(X^7)$ at many (evenly spaced) points, but
2515: these points have no direct relationship to the ``interesting'' points on
2516: the graph of this function. No value close enough to the maxima and minima
2517: are calculated, which leads to wrong turning points of the graph.
2518:
2519: There is a way to fix this: one can ask PARI to use variable step which
2520: smaller at the points where the graph of the function is more curved:
2521: \bprog
2522: ploth(X = -2, 2, sin(X^7),recursive)
2523: @eprog
2524:
2525: \noindent The precision near the edges of the graph is much better now.
2526: However, the recursive plotting (named so since PARI subdivides intervals
2527: until the graph becomes almost straight) has its own pitfalls. Try
2528: \bprog
2529: ploth(X = -2, 2, sin(X*7), recursive)
2530: @eprog
2531:
2532: \noindent Note that the graph looks correct far away, but it has a straight
2533: interval near the origin, and some sharp corners as well. This happens
2534: because the graph is symmetric with respect to the origin, thus the middle 3
2535: points calculated during the initial subdivision of $[-2,2]$ are exactly on
2536: the same line. To PARI this indicates that no further subdivision is needed,
2537: and it plots the graph on this subinterval as a straight line.
2538:
2539: There are many ways to circumvent this. Say, one can make the right limit
2540: 2.1. Or one can ask PARI for an initial subdivision into 16 points instead
2541: of default 15:
2542: \bprog
2543: ploth(X = -2, 2, sin(X*7), recursive, 16)
2544: @eprog
2545:
2546: All these arrangements break the symmetry of the initial subdivision, thus
2547: make the problem go away. Eventually PARI will be able to better detect such
2548: pathological cases, but currently some manual intervention may be required.
2549:
2550: Function \kbd{ploth} has some additional enhancements which allow graphing
2551: in situations when the calculation of the function takes a lot of time. Let
2552: us plot $\zeta({1\over 2} + it)$:
2553: \bprog
2554: ploth(t = 100, 110, real(zeta(0.5+I*t)), /*empty*/, 1000)
2555: @eprog
2556:
2557: \noindent This can take quite some time. (1000 is close to the default for
2558: many plotting devices, we want to specify it explicitely so that the result
2559: do not depend on the output device.) Try the recursive plot:
2560: \bprog
2561: ploth(t = 100, 110, real(zeta(0.5+I*t)), recursive)
2562: @eprog
2563:
2564: It takes approximately the same time. Now try specifying fewer points,
2565: but make PARI approximate the data by a smooth curve:
2566: \bprog
2567: ploth(t = 100, 110, real(zeta(0.5+I*t)), splines, 100)
2568: @eprog
2569:
2570: This takes much less time, and the output is practically the same. How to
2571: compare these two outputs? We will see it shortly. Right now let us
2572: plot both real and complex parts of $\zeta$ on the same graph:
2573: \bprog
2574: f(t) = z=zeta(0.5+I*t); [real(z),imag(z)]
2575: ploth(t = 100, 110, f(t), , 1000)
2576: @eprog
2577:
2578: Note how one half of the roots of the real and imaginary parts coincide.
2579: Why did we define a function \kbd{f(t)}? To avoid calculation of
2580: $\zeta({1\over2} + it)$ twice for the same value of t. Similarly, we can
2581: plot parametric graphs:
2582: \bprog
2583: ploth(t = 100, 110, f(t), parametric, 1000)
2584: @eprog
2585:
2586: Again, one can speed up the calculation with
2587: \bprog
2588: ploth(t = 100, 110, f(t), parametric+splines, 100)
2589: @eprog
2590:
2591: If your plotting device supports it, you may ask PARI to show the points
2592: in which it calculated your function:
2593: \bprog
2594: ploth(t = 100, 110, f(t), parametric+splines+points_lines, 100)
2595: @eprog
2596:
2597: As you can see, the points are very dense on the graph. To see some crude
2598: graph, one can even decrease the number of points to 30. However, if you
2599: decrease the number of points to 20, you can see that the approximation to
2600: the graph now misses zero. Using splines, one can create reasonable graphs
2601: for larger values of t, say with
2602: \bprog
2603: ploth(t = 10000, 10004, f(t), parametric+splines+points_lines, 50)
2604: @eprog
2605:
2606: How can we compare two graphs of the same function plotted by different
2607: methods? Documentation shows that \kbd{ploth} does not provide any direct
2608: method to do so. However, it is possible, and even not very complicated.
2609:
2610: The solution comes from the other direction. PARI has a power mix of high
2611: level plotting function with low level plotting functions, and these functions
2612: can be combined together to obtain many different effects. Return back
2613: to the graph of $\sin(X^7)$. Suppose we want to create an additional
2614: rectangular frame around our graph. No problem!
2615:
2616: First, all low-level graphing work takes place in some virtual drawing
2617: boards (numbered from 0 to 15), called ``rectangles'' (or ``rectwindows'').
2618: So we create an empty ``rectangle'' and name it rectangle 2 (any
2619: number between 0 and 15 would do):
2620: \bprog
2621: plotinit(2)
2622: plotscale(2, 0,1, 0,1)
2623: @eprog
2624: This creates a rectwindow whose size exactly fits the size of the output
2625: device, and makes the coordinate system inside it go from 0 to 1 (for both
2626: $x$ and $y$). Create a rectangular frame along the boundary of this rectangle:
2627: \bprog
2628: plotmove(2, 0,0)
2629: plotbox(2, 1,1)
2630: @eprog
2631: Suppose we want to draw the graph inside a subrectangle of this with upper
2632: and left margins of $0.10$ (so 10\% of the full rectwindow width), and
2633: lower and top margins of $0.02$, just to make it more interesting. That
2634: makes it an $0.88 \times 0.88$ subrectangle; so we create another rectangle
2635: (call it 3) of linear size 0.88 of the size of the initial rectangle and
2636: graph the function in there:
2637: \bprog
2638: plotinit(3, 0.88, 0.88, relative)
2639: plotrecth(3, X = -2, 2, sin(X^7), recursive)
2640: @eprog
2641: (nothing is output yet, these commands only fills the virtual drawing
2642: boards with PARI graphic objects). Finally, output rectangles 2 and 3 on
2643: the same plot, with the required offets (counted from upper-left corner):
2644: \bprog
2645: plotdraw([2, 0,0, 3, 0.1,0.02], relative)
2646: @eprog
2647: \noindent The output misses something comparing to the output of
2648: \kbd{ploth}: there are no coordinates of the corners of the internal
2649: rectangle. If your output device supports mouse operations (only
2650: \kbd{gnuplot} does), you can find coordinates of particular points of the
2651: graph, but it is nice to have something printed on a hardcopy too.
2652:
2653: However, it is easy to put $x$- and $y$-limits on the graph. In the
2654: coordinate system of the rectangle 2 the corners are $(0.1,0.1)$,
2655: $(0.1,0.98)$, $(0.98,0.1)$, $(0.98,0.98)$. We can mark lower $x$-limit by
2656: doing
2657: \bprog
2658: plotmove(2, 0.1,0.1)
2659: plotstring(2, "-2.000", left+top+vgap)
2660: @eprog\noindent
2661: Computing the minimal and maximal $y$-coordinates might be trickier, since
2662: in principle we do not know the range in advance (though for $\sin X^7$ it
2663: is quite easy to guess!). Fortunately, \kbd{plotrecth} returns the $x$- and
2664: $y$-limits.
2665:
2666: Here is the complete program:
2667: \bprog
2668: plotinit(3, 0.88, 0.88, relative)
2669: lims = plotrecth(3, X = -2, 2, sin(X^7), recursive)
2670: \p 3 \\ @com $3$ significant digits for the bounding box are enough
2671: limits = vector(4,i, Str(lims[i]))
2672: plotinit(2); plotscale(2, 0,1, 0,1)
2673: plotmove(2, 0,0); plotbox(2, 1,1)
2674: plotmove(2, 0.1,0.1);
2675: plotstring(2, limits[1], left+top+vgap)
2676: plotstring(2, limits[3], bottom+vgap+right+hgap)
2677: plotmove(2, 0.98,0.1); plotstring(2, limits[2], right+top+vgap)
2678: plotmove(2, 0.1,0.98); plotstring(2, limits[4], right+hgap+top)
2679: plotdraw([2, 0,0, 3, 0.1,0.02], relative)
2680: @eprog
2681:
2682: We started with a trivial requirement: have an additional frame around
2683: the graph, and it took some effort to do so. But at least it was possible,
2684: and PARI did the hardest part: creating the actual graph.
2685: Now do a different thing: plot together the ``exact'' graph of
2686: $\zeta({1/2}+it)$ together with one obtained from splines approximation.
2687: We can emit these graphs into two rectangles, say 0 and 1, then put these
2688: two rectangles together on one plot. Or we can emit these graphs into one
2689: rectangle 0.
2690:
2691: However, a problem arises: note how we
2692: introduced a coordinate system in rectangle 2 of the above example, but we
2693: did not introduce a coordinate system in rectangle 3. Plotting a
2694: graph into rectangle 3 automatically created a coordinate system
2695: inside this rectangle (you could see this coordinate system in action
2696: if your output device supports mouse operations). If we use two different
2697: methods of graphing, the bounding boxes of the graphs will not be exactly
2698: the same, thus outputting the rectangles may be tricky. Thus during
2699: the second plotting we ask \kbd{plotrecth} to use the coordinate system of
2700: the first plotting. Let us add another plotting with fewer
2701: points too:
2702: \bprog
2703: plotinit(0, 0.9,0.9, relative)
2704: plotrecth(0, t=100,110, f(t), parametric, 300)
2705: plotrecth(0, t=100,110, f(t), parametric+splines+points_lines+norescale, 30);
2706: plotrecth(0, t=100,110, f(t), parametric+splines+points_lines+norescale, 20);
2707: plotdraw([0, 0.05,0.05], relative)
2708: @eprog
2709:
2710: This achieves what we wanted: we may compare different ways to plot a graph,
2711: but the picture is confusing: which graph is what, and why there are multiple
2712: boxes around the graph? At least with some output devices one can control
2713: how the output curves look like, so we can use this to distinguish different
2714: graphs. And the mystery of multiple boxes is also not that hard to solve:
2715: they are bounding boxes for calculated points on each graph. We can disable
2716: output of bounding boxes with appropriate options for \kbd{plotrect}.
2717: With these frills the script becomes:
2718: \bprog
2719: plotinit(0, 0.9,0.9, relative)
2720: plotpointtype(-1, 0) \\@com set color of graph points
2721: plotpointsize(0, 0.4) \\@com use tiny markers (if available)
2722: plotrecth(0, t=100,110, f(t), parametric+points, 300)
2723: plotpointsize(0, 1) \\@com normal-size markers
2724: plotlinetype(-1, 1) \\@com set color of graph lines
2725: plotpointtype(-1, 1) \\@com set color of graph points
2726: curve_only = norescale + no_frame + no_x_axis + no_y_axis
2727: plotrecth(0, t=100,110,f(t), parametric+splines+points_lines+curve_only, 30);
2728: plotlinetype(-1, 2) \\@com set color of graph lines
2729: plotpointtype(-1, 2) \\@com set color of graph points
2730: plotrecth(0, t=100,110,f(t), parametric+splines+points_lines+curve_only, 20);
2731: plotdraw([0, 0.05,0.05], relative)
2732: @eprog
2733:
2734: \noindent Plotting axes on the second and third graph would not hurt, but
2735: is not needed either, so we omit them. One can see that the discrepancies
2736: between the exact graph and one based on 30 points exist, but are pretty
2737: small. On the other hand, decreasing the number of points to 20 makes
2738: quite a noticable difference.
2739:
2740: Keep in mind that \kbd{plotlinetype}, \kbd{plotpointtype},
2741: \kbd{plotpointsize} may do nothing on some terminals.
2742:
2743: What if we
2744: want to create a high-resolution hardcopy of the plot? There may be several
2745: possible solutions. First, the display output device may allow a
2746: high-resolution hardcopy itself. Say, PM display (with gnuplot output on
2747: OS/2) pretends that its resolution is $19500\times 12500$, thus the data
2748: PARI sends to it are already high-resolution, and printing is available
2749: through the menubar. Alternatively, with gnuplot output one can switch
2750: the output plotting device to many different hardcopy devices:
2751: \kbd{plotfile("plot.tex")}, \kbd{plotterm("texdraw")}.
2752: After this all the plotting will go into file {\tt plot.tex} with whatever
2753: output conventions gnuplot format {\tt texdraw} provides. To switch output
2754: back to normal, one needs to restore the initial plotting terminal, and
2755: restore the initial output file by doing \kbd{plotfile("-")}.
2756:
2757: One can combine PARI programming capabilities to produce multiple plots:
2758: \bprog
2759: plotfile("manypl1.gif") \\@com Avoid switching STDOUT to binary mode
2760: plotterm("gif=300,200")
2761: wpoints = plothsizes()[1] \\@com $300 \times 200$ is advice only
2762: {
2763: for( k=1,6,
2764: plotfile("manypl" k ".gif");
2765: ploth(x = -1, 3, sin(x^k), , wpoints)
2766: )
2767: }
2768: @eprog
2769:
2770: \noindent This plots 6 graphs of $\sin x^k$, $k=1\dots 6$ into
2771: $300\times 200$ GIF files {\tt manypl1.gif}\dots {\tt manypl6.gif}.
2772:
2773: Additionally, one can ask PARI to output a plot into a PS file: just
2774: use the command \kbd{psdraw} instead of \kbd{plotdraw} in the above examples
2775: (or \kbd{psploth} instead of \kbd{ploth}). See \kbd{psfile} argument
2776: to \kbd{default} for how to change the output file for this operation. Keep
2777: in mind that the precision of PARI PS output will be the same as for the
2778: current output device.
2779:
2780: Now suppose we want to join many different small graphs into one picture.
2781: We cannot use one rectangle for all the output as we did in the example
2782: with $\zeta({1/2}+it)$, since the graphs should go into different places.
2783: Rectangles are a scarce commodity in PARI, since only 16 of them are
2784: user-accessible. Does it mean that we cannot have more than 16 graphs on
2785: one picture? Thanks to an additional operation of PARI plotting engine,
2786: there is no such restrictions. This operation is \kbd{plotcopy}.
2787:
2788: The following script puts 4 different graphs on one plot using 2 rectangles
2789: only, \kbd{A} and \kbd{T}:
2790: \bprog
2791: A = 2; \\@com accumulator
2792: T = 3; \\@com temporary target
2793: plotinit(A); plotscale(A, 0, 1, 0, 1)
2794:
2795: plotinit(T, 0.42, 0.42, relative);
2796: plotrecth(T, x= -5, 5, sin(x), recursive)
2797: plotcopy(T, 2, 0.05, 0.05, relative + nw)
2798:
2799: plotmove(A, 0.05 + 0.42/2, 1 - 0.05/2)
2800: plotstring(A,"Graph", center + vcenter)
2801:
2802: plotinit(T, 0.42, 0.42, relative);
2803: plotrecth(T, x= -5, 5, [sin(x),cos(2*x)], 0)
2804: plotcopy(T, 2, 0.05, 0.05, relative + ne)
2805:
2806: plotmove(A, 1 - 0.05 - 0.42/2, 1 - 0.05/2)
2807: plotstring(A,"Multigraph", center + vcenter)
2808:
2809: plotinit(T, 0.42, 0.42, relative);
2810: plotrecth(T, x= -5, 5, [sin(3*x), cos(2*x)], parametric)
2811: plotcopy(T, 2, 0.05, 0.05, relative + sw)
2812:
2813: plotmove(A, 0.05 + 0.42/2, 0.5)
2814: plotstring(A,"Parametric", center + vcenter)
2815:
2816: plotinit(T, 0.42, 0.42, relative);
2817: plotrecth(T, x= -5, 5, [sin(x), cos(x), sin(3*x),cos(2*x)], parametric)
2818: plotcopy(T, 2, 0.05, 0.05, relative + se)
2819:
2820: plotmove(A, 1 - 0.05 - 0.42/2, 0.5)
2821: plotstring(A,"Multiparametric", center + vcenter)
2822:
2823: plotlinetype(A, 3)
2824: plotmove(A, 0, 0)
2825: plotbox(A, 1, 1)
2826:
2827: plotdraw([A, 0, 0])
2828: \\ psdraw([A, 0, 0], relative) \\ @com if hardcopy is needed
2829: @eprog
2830:
2831: The rectangle \kbd{A} plays the role of accumulator, rectangle \kbd{T} is
2832: used as a target to \kbd{plotrecth} only. Immediately after plotting into
2833: rectangle \kbd{T} the contents is copied to accumulator. Let us explain
2834: numbers which appear in this example: we want to create 4 internal rectangles
2835: with a gap 0.06 between them, and the outside gap 0.05 (of the size of the
2836: plot). This leads to the size 0.42 for each rectangle. We also
2837: put captions above each graph, centered in the middle of each gap. There
2838: is no need to have a special rectangle for captions: they go into the
2839: accumulator too.
2840:
2841: To simplify positioning of the rectangles, the above example uses relative
2842: ``geographic'' notation: the last argument of \kbd{plotcopy} specifies the
2843: corner of the graph (say, northwest) one counts offset from. (Positive
2844: offsets go into the rectangle.)
2845:
2846: To demonstrate yet another useful plotting function, design a program to
2847: plot Taylor polynomials for a $\sin x$ about 0. For simplicity, make the
2848: program good for any function, but assume that a function is odd, so only
2849: odd-numbered Taylor series about 0 should be plotted. Start with defining
2850: some useful shortcuts
2851: \bprog
2852: xlim = 13; ordlim = 25; f(x) = sin(x);
2853: default(seriesprecision,ordlim)
2854: farray(t) = vector((ordlim+1)/2, k, truncate( f(1.*t + O(t^(2*k+1)) )))
2855: FARRAY = farray('t); \\@com\kbd{'t} to make sure \kbd{t} is a free variable
2856: @eprog
2857:
2858: \noindent \kbd{farray(x)} returns a vector of Taylor polynomials for
2859: $f(x)$, which we store in \kbd{FARRAY}. We want to plot $f(x)$ into a
2860: rectangle, then make the rectangle which is 1.2 times higher, and plot
2861: Taylor polynomials into the larger rectangle. Assume that the larger
2862: rectangle takes 0.9 of the final plot.
2863:
2864: First of all, we need to measure the height of the smaller rectangle:
2865: \bprog
2866: plotinit(3, 0.9, 0.9/1.2, 1);
2867: curve_only = no_x_axis+no_y_axis+no_frame;
2868: lims = plotrecth(3, x= -xlim, xlim, f(x), recursive+curve_only,16);
2869: h = lims[4] - lims[3];
2870: @eprog
2871:
2872: \noindent Next step is to create a larger rectangle, and plot the Taylor
2873: polynomials into the larger rectangle:
2874: \bprog
2875: plotinit(4, 0.9,0.9, relative);
2876: plotscale(4, lims[1], lims[2], lims[3] - h/10, lims[4] + h/10)
2877: plotrecth(4, x = -xlim, xlim, subst(FARRAY,t,x), norescale);
2878: @eprog
2879:
2880: Here comes the central command of this example:
2881:
2882: \kbd{plotclip(4)}
2883:
2884: \noindent What does it do? The command \kbd{plotrecth(\dots, norescale)}
2885: scales the graphs according to coordinate system in the
2886: rectangle, but it does not pay any other attension to the size of
2887: the rectangle. Since \kbd{xlim} is 13, the Taylor polynomials take
2888: very large values in the interval \kbd{-xlim...xlim}. In particular,
2889: significant part of the graphs is going to be {\it outside} of the rectangle.
2890: \kbd{plotclip} removes the parts of the plottings which fall off the
2891: rectangle boundary
2892: \bprog
2893: plotinit(2)
2894: plotscale(2, 0.0, 1.0, 0.0, 1.0)
2895: plotmove(2,0.5,0.975)
2896: plotstring(2,"Multiple Taylor Approximations",center+vcenter)
2897: plotdraw([2, 0, 0, 3, 0.05, 0.05 + 0.9/12, 4, 0.05, 0.05], relative)
2898: @eprog
2899:
2900: These commands draw a caption, and combine 3 rectangles (one with the
2901: caption, one with the graph of the function, and one with graph of Taylor
2902: polynomials) together.
2903:
2904: This finishes our short survey of PARI plotting functions, but let us add
2905: some remarks. First of all, for a typical output device the picture is
2906: composed of small colored squares (pixels) (as a very large checkerboard).
2907: Each output rectangle is in fact a union of such squares. Each drop
2908: of paint in the rectangle will color a whole square it is in. Since the
2909: rectangle has a coordinate system, sometimes it is important to know how
2910: this coordinate system is positioned with respect to the boundaries of
2911: these squares.
2912:
2913: The command \kbd{plotscale} describes a range of $x$ and $y$ in the
2914: rectangle. The limit values of $x$ and $y$ in the coordinate system are
2915: coordinates {\it of the centers} of corner squares. In particular,
2916: if ranges of $x$ and $y$ are $[0,1]$, then the segment which connects (0,0)
2917: with (0,1) goes along the {\it middle} of the left column of the rectangle.
2918: In particular, if we made tiny errors in calculation of endpoints of this
2919: segment, this will not change which squares the segment intersects, thus
2920: the resulting picture will be the same. (It is important to know such details
2921: since many calculations in PARI are approximate.)
2922:
2923: Another consideration is that all examples we did in this section were
2924: using relative sizes and positions for the rectangles. This is nice, since
2925: different output devices will have very similar pictures, while we
2926: did not need to care about particular resolution of the output device.
2927: On the other hand,
2928: using relative positions does not guarantie that the pictures will be
2929: similar. Why? Even if two output devices have the same resolution,
2930: the picture may be different. The devices may use fonts of different
2931: size, or may have a different ``unit of length''.
2932:
2933: The information about the device in PARI is encoded in 6 numbers: resolution,
2934: size of a character cell of the font, and unit of length, all separately
2935: for horizontal and vertical direction. These sizes are expressed as
2936: numbers of pixels. To inspect these numbers one may use the function
2937: \kbd{plothsizes}. The ``units of length'' are currently used to calculate
2938: right and top gaps near graph rectangle of \kbd{ploth}, and gaps for
2939: \kbd{plotstring}. Left and bottom gaps near graph rectangle are calculate
2940: using both units of length, and sizes of character boxes (so that there
2941: is enough place to print limits of the graphs).
2942:
2943: What does it show? Using relative sizes during plotting produces
2944: \var{approximately} the same plotting on different devices, but does not
2945: ensure that the plottings ``look the same''. Moreover, ``looking the
2946: same'' is not a desirable target, ``looking tuned for the environment''
2947: will be much better. If you want to produce such fine-tuned plottings,
2948: you need to abandone a relative-size model, and do your plottings in
2949: pixel units. To do this one removes flag \kbd{relative} from the above
2950: examples, which will make size and offset arguments interpreted this way.
2951: After querying sizes with \kbd{plothsizes} one can fine-tune sizes and
2952: locations of subrectangles to the details of an arbitrary plotting
2953: device.
2954:
2955: To check how good your fine-tuning is, you may test your graphs with a
2956: medium-resolution plotting (as many display output devices are), and
2957: with a low-resolution plotting (say, with \kbd{plotterm("dumb")} of gnuplot).
2958: \vfill\eject\bye
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to tutorial.tex CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / doc