Annotation of OpenXM_contrib/pari/doc/tutorial.tex, Revision 1.1.1.1
1.1 maekawa 1: % This should be compiled with plain TeX
2: \input parimacro.tex
3: \def\maketitle#1{\currentlabel.\ #1}
4: \begintitle
5: \vskip 2.5truecm
6: \centerline{\mine A Tutorial}
7: \vskip 1.truecm
8: \centerline{\mine for}
9: \vskip 1.truecm
10: \centerline{\mine PARI-GP}
11: \vskip 1.truecm
12: \authors
13: \centerline{last updated March 5, 1999}
14: \centerline{(this document distributed with version \vers)}
15: \endtitle
16:
17: This booklet is intended to be a guided tour and a tutorial to the GP
18: calculator. Many examples will be given, but each time a new function is
19: used, the reader should look at the appropriate section in the User's Manual
20: for detailed explanations. Hence although this chapter can be read
21: independently (for example to get rapidly acquainted with the possibilities
22: of GP without having to read the whole manual), the reader will profit most
23: from it by reading it in conjunction with the reference manual.
24:
25: \section{Greetings!}
26:
27: So you are sitting in front of your workstation (or terminal, or PC,\dots),
28: and you type \kbd{gp} to get the program started (remember to always hit the
29: \kbd{Enter} key and not the \kbd{Return} key on a Macintosh computer).
30:
31: It says hello in its particular manner, and then waits for you after its
32: \kbd{prompt}, initially \kbd{?} (or something like {\bf gp}~\kbd{>}).
33:
34: Type \kbd{2 + 2}. What happens? Maybe not what you expect. First of all, of
35: course, you should tell GP that your input is finished, and this is done by
36: hitting the \kbd{Return} (or \kbd{Newline}) key, or the \kbd{Enter} key on
37: the Mac. If you do exactly this, you will get the expected answer. However
38: some of you may be used to other systems like Macsyma or Maple. In this case,
39: you will have subconsciously ended the line with a semicolon ``\kbd{;}''
40: before hitting \kbd{Return}, since this is how it is done on those systems.
41: In that case, you will simply see GP answering you with a smug expression,
42: i.e.~a new prompt and no answer! This is because a semicolon at the end of a
43: line in GP tells it to keep the result, but not to print it (you will
44: certainly want to use this feature if the output is several pages long).
45:
46: Try \kbd{27 * 37}. Wow! even multiplication works. Actually, maybe those
47: spaces are not necessary after all. Let's try \kbd{27*37}. Yup, seems to be
48: ok. We'll still insert them in this document since it makes things easier to
49: read, but as GP does not care about them, you don't have to type them all!
50:
51: Now this session is getting lengthy, so the second thing one needs to learn
52: is to quit. Each system has its quit signal (to name a few: \kbd{quit},
53: \kbd{exit}, \kbd{system},\dots). In GP, you can use \kbd{quit} or \b{q}
54: (backslash q), the \kbd{q} being of course for quit. Try it.
55:
56: Now you've done it! You're out of GP, so how do you want to continue studying
57: this tutorial? Get back in please (see above).
58:
59: Let's get to more serious stuff. Let's see, I seem to remember that the
60: decimal expansion of $1/7$ has some interesting properties. Let's see what GP
61: has to say about this. Type \kbd{1 / 7}. What? This computer is making fun of
62: me, it just spits back to me my own input, that's not what I want!
63:
64: Now stop complaining, and think a little. This system has been written mainly
65: for pure mathematicians, and not for physicists (although they are welcome to
66: use it). And mathematically, $1/7$ is an element of the field $\Q$ of
67: rational numbers, so how else but $1/7$ can the computer give the answer to
68: you? (well maybe $2/14$, but why complicate matters?). Seriously, the basic
69: point here is that PARI (hence GP) will almost always try to give you a
70: result which is as precise as possible (we will see why ``almost'' later),
71: hence since here you gave an operation whose result can be represented
72: exactly, that's what it gives you.
73:
74: OK, but I still want the decimal expansion of $1/7$. No problem. Type one of
75: the following:
76: %
77: \bprog
78: 1./ 7
79: 1 / 7.
80: 1./ 7.
81: 1 / 7 + 0. \dots
82: \eprog
83: %
84: Immediately a number of decimals of this fraction appear (28 on most systems,
85: 38 on the others), and the repeating pattern is $142857$. The reason is that
86: you have included in the operations numbers like \kbd{0.}, \kbd{1.} or \kbd{7.}
87: which are {\it imprecise\/} real numbers, hence GP cannot give you an exact
88: result.
89:
90: Why 28 / 38 decimals by the way? Well, it is the default initial precision,
91: as indicated when you start GP. This has been chosen so that the
92: computations are very fast, and gives already 12 decimals more accuracy than
93: conventional double precision floating point operations. The precise value
94: depends on a technical reason: if your machine supports 64-bit integers (the
95: standard library can handle integers up to $2^{64}$), the default precision
96: will be 38 decimals, and 28 otherwise. As the latter is most probably the
97: case (if your machine is not a \kbd{DEC alpha}, that is), we'll assume it
98: henceforth.
99:
100: Only large mainframes or supercomputers have 28 digits of precision in their
101: standard libraries, and that is their absolute limit. Not here of course. You
102: can extend the precision (almost) as much as you like as we will see in a
103: moment.
104:
105: I'm getting bored, why don't we get on with some more exciting stuff? Well,
106: try \kbd{exp(1)}. Presto, comes out the value of $e$ to 28 digits. Try
107: \kbd{log(exp(1))}. Well, it's not exactly equal to 1, but pretty close! That's
108: what you lose by working numerically.
109:
110: Let's see, what could we try now. Hum, \kbd{pi}? The answer is not that
111: enlightening. \kbd{Pi}? Ah. This works better. But let's remember that GP
112: distinguishes between uppercase and lowercase letters. \kbd{pi} was as
113: meaningless to it as \kbd{stupid garbage} would have been: in both cases GP
114: will just create a variable with that funny unknown name you just used. Try
115: it! Note that it is actually equivalent to type \kbd{stupidgarbage}: all
116: spaces are suppressed from the input. In the \kbd{27~*~37} example it was
117: not so conspicuous as we had an operator to separate the two operands. This
118: has important consequences for the writing of GP scripts (more about this
119: later).
120:
121: By the way, you can ask GP about any identifier you think it might know
122: about: just type it, prepending a question mark ``\kbd{?}''. Try \kbd{?Pi}
123: and \kbd{?pi} for instance. On some systems, an extended online help might
124: be available: try doubling the question mark to check whether it's the case on
125: yours: \kbd{??Pi}. In fact the GP header already gave you that information if
126: it was the case (just before the copyright message). As well, if it says
127: something like ``\kbd{readline enabled}'' then you should have a look at
128: \secref{se:readline} in the User's Manual before you go on: it will be much
129: easier to type in examples and correct typos after you've done that.
130:
131: Now try \kbd{exp(Pi * sqrt(163))}. Hmmm, since from our last example we
132: suspect that the last digit may be wrong, can this really be an integer?
133: This is the time to change precision. Type \kbd{\b{p} 50}, then try
134: \kbd{exp(Pi * sqrt(163))} again. We were right to suspect that the last
135: decimal was incorrect, since we get even more nines than before, but it is
136: now convincingly clear that this is not an integer. Maybe it's a bug in PARI,
137: and the result is really an integer? Type \kbd{sqr(log(\%) / Pi)} immediately
138: after the preceding computation (\kbd{\%} means the result of the last
139: computed expression). More generally, the results are numbered
140: \kbd{\%1, \%2, \dots} {\it including} the results that you do not want to see
141: printed by putting a semicolon at the end of the line, and you can evidently
142: use all these quantities in any further computations. \kbd{sqr} is the square
143: function (\kbd{sqr(x) = x * x}), not to be confused with \kbd{sqrt} which is
144: the square root function). The result seems to be indistinguishable from $163$,
145: hence it does not seem to be a bug.
146:
147: In fact, it is known that $\exp(\pi*\sqrt{n})$ not only is not an integer
148: or a rational number, but is even a transcendental number when $n$ is a
149: non-zero rational number.
150:
151: So GP is just a fancy calculator, able to give me more decimals than I will
152: ever need? Not so, GP is incredibly more powerful than an ordinary
153: calculator, even independently of its arbitrary precision possibilities.
154:
155: \noindent {\bf Additional comments} (you are supposed to skip this at first,
156: and come back later)
157:
158: 1) If you are a PARI old timer, you will notice pretty soon (or maybe you
159: have already?) that many many things changed between the older 1.39.xx
160: versions and this one. And conspicuously, most function names have been
161: changed. We sincerely think it's for the best since they are much more
162: logical now and the systematic prefixing is very convenient coupled with the
163: automatic completion mechanism: it's now very easy to know what functions are
164: available to deal with, say, elliptic curves since they all share the prefix
165: \kbd{ell}.
166:
167: Of course, this is going to break all your nice old scripts. Well, you can
168: either change the compatibility level (typing \kbd{default(compatible, 3)}
169: will send you back to the stone-age behaviour of good ol' version 1.39.15),
170: or rewrite the scripts. We really advise you to do the latter if they are not
171: too long, since they can now be written much more cleanly than before
172: (especially with the new control statements: \kbd{break}, \kbd{next},
173: \kbd{return}), and besides it'll be as good a way as any to get used to the new
174: names. We {\it might} provide an automatic transcriptor with future versions.
175:
176: To know how a specific function was changed, just type
177: \kbd{whatnow({\rm function})}.
178:
179: 2) It seems that the text implicitly says that as soon as an imprecise number
180: is entered, the result will be imprecise. Is this always true? There is a
181: unique exception: when you multiply an imprecise number by the exact number 0,
182: you will get the exact 0. Compare \kbd{0 * 1.4} and \kbd{0.~*~1.4}. \smallskip
183: %
184: 3) Not only can the number of decimal places of real numbers be large, but the
185: number of digits of integers also. Try \kbd{100!}. It is never necessary to
186: tell GP in advance the size of the integers that it will encounter, since this
187: is adjusted automatically. On the other hand, for many computations with real
188: numbers, it is necessary to specify a default precision (initially 28 digits).
189: \smallskip
190: %
191: 4) Come back to 28 digits of precision (\kbd{\b{p} 28}), and type
192: \kbd{exp(24 * Pi)}. As you can see the result is printed in exponential format.
193: This is because GP never wants you to believe that a result is correct when it
194: is not.
195:
196: We are working with 28 digits of precision, but the integer part of
197: $\exp(24*\pi)$ has 33 decimal digits. Hence if GP had dutifully printed out 33
198: digits, the last few digits would have been wrong. Hence GP wants to print only
199: 28 significant digits, but to do so it has to print in exponential format.
200: \smallskip
201: %
202: 5) There are two ways to avoid this. One is of course to increase the precision
203: to more than 33 decimals. Let's try it. To give it a wide margin, we set the
204: precision to 40 decimals. Then we recall our last result (\kbd{\%}
205: or \kbd{\%n} where \kbd{n} is the number of the result). What? We still have
206: an exponential format! Do you understand why?
207:
208: Again let's try to see what's happening. The number you recalled had been
209: computed only to 28 decimals, and even if you set the precision to 1000
210: decimals, GP knows that your number has only 28 digits of accuracy but an
211: integral part with 33 digits. So you haven't improved things by increasing
212: the precision. Or have you? What if we retype \kbd{exp(24 * Pi)} now that we
213: have 40 digits? Try it. Now we do not have an exponential format, and we see
214: that at 28 decimals the last 6 digits would have been wrong if they had been
215: printed in fixed-point format.
216: \smallskip
217: %
218: 6) {\it Warning\/}. Try the following: starting in precision 28, type first
219: \kbd{default(format, "e0.50")}, then \kbd{exp(24 * Pi)}. Do you understand why
220: the result is so bad, and why there are lots of zeros at the end? Convince
221: yourself by typing \kbd{log(exp(1))}. The moral is that the
222: \kbd{default(format,)} command changes only the output format, but {\it not\/}
223: the default precision. On the other hand, the \b{p} command changes both the
224: precision and the output format.
225:
226: 7) What if I forget what the current precision is and I don't feel like
227: counting all the decimals? Well, you can learn about GP internal variables
228: (and change them!) using \kbd{default}. Type \kbd{default(realprecision)},
229: then \kbd{default(realprecision, 38)}. Huh? In fact this last command is
230: strictly equivalent to \kbd{\b{p} 38}! (admittedly more cumbersome to type).
231: There are more ``defaults'' than just \kbd{format} and \kbd{realprecision}:
232: type \kbd{default} by itself now, they are all there.
233:
234: 8) Note that the \kbd{default} command reacts differently according to the
235: number of input arguments. This is not an uncommon behaviour for GP functions.
236: You can see this from the online help (or the complete description in
237: Chapter~3): any argument surrounded by braces \kbd{\obr\cbr} in the function
238: prototype is optional (which really means that a {\it default} argument will be
239: supplied by GP). You can then check out from the text what effect a given value
240: will have, and in particular the default one.
241:
242: \section{Warming up}
243:
244: Another thing you better get used to pretty fast is error messages. Try
245: typing \kbd{1/0}. Couldn't be clearer. Taking again our universal example in
246: precision 28, type \kbd{floor(exp(24 * Pi))} (\kbd{floor} is the
247: mathematician's integer part, not to be confused with \kbd{truncate}, which is
248: the computer scientist's: \kbd{floor(-3.4)} is equal to $-4$ whereas
249: \kbd{truncate(-3.4)} is equal to $-3$). You get a more cryptic error message,
250: which you would immediately understand if you had read the additional
251: comments of the preceding section. Since I told you not to read them, the
252: explanation is simply that GP is unable to compute the integer part of
253: \kbd{exp(24 * Pi)} given only 28 decimals of accuracy, since it has 33 digits.
254:
255: Some error messages are even much more cryptic than that and are sometimes
256: not so easy to understand (well, it's nothing compared to \TeX's error
257: messages!).
258:
259: For instance, try \kbd{log(x)}. Not really clear, is it? But no matter, it
260: simply tells you that GP simply does not understand what \kbd{log(x)} is
261: (although it does know the \kbd{log} function, as \kbd{?log} will readily
262: tell us).
263:
264: Now let's try \kbd{sqrt(-1)} to see what error message we get now. Haha! GP
265: even knows about complex numbers, so impossible to trick it that way.
266: Similarly, try typing \kbd{log(-2)}, \kbd{exp(I*Pi)}, \kbd{I\pow I},\dots
267: So we have a lot of real and complex analysis at our disposal (note that
268: there is always a specific branch of multivalued complex transcendental
269: functions which is taken, specified in the manual). Again, beware that
270: \kbd{I} and \kbd{i} are not the same thing. Compare \kbd{I\pow2} with
271: \kbd{i\pow2} for instance.
272:
273: Just for fun, let's try \kbd{6*zeta(2) / Pi\pow2}. Pretty close, no?
274:
275: \medskip
276: Now GP didn't seem to know what \kbd{log(x)} was, although it did know how to
277: compute numerical values of \kbd{log}. This is annoying. Maybe it knows the
278: exponential function? Let's give it a try. Type \kbd{exp(x)}. What's this? If
279: you had had any experience with other systems, the answer should have simply
280: been \kbd{exp(x)} again. But here the answer is the Taylor expansion of the
281: function around $\kbd{x}=0$, to 16 terms. 16 is the default
282: \kbd{seriesprecision}, which can be changed by typing \kbd{\b{ps} $n$} or
283: \kbd{default(seriesprecision, $n$)} where $n$ is the number of terms that you
284: want in your power series (note the \kbd{O(x\pow16)} which ends the series,
285: and which is trademark of this type of object in GP. It is the familiar
286: ``big--oh'' notation of analysis).
287:
288: You will thus automatically get the Taylor expansion of any function that can
289: be expanded around $0$, and incidentally this explains why we weren't
290: able to do anything with \kbd{log(x)} which is not defined at 0. But if we
291: try \kbd{log(1+x)}, then it works. But what if we wanted the expansion
292: around a point different from 0? Well, you're able to change $x$ into
293: $x-a$, aren't you? So for instance you can type \kbd{log(x+2)} to
294: have the expansion of \kbd{log} around $\kbd{x}=2$.
295:
296: \begingroup\def\h{\hskip 0em plus 1ex}
297: As exercises, try \kbd{cos(x)}, \h\kbd{cos(x)\pow 2 + sin(x)\pow 2}, \h
298: \kbd{exp(cos(x))}, \h\kbd{exp(exp(x) - 1)},\break
299: \kbd{gamma(1 + x)}, with different values of \kbd{serieslength} if you like.
300: \endgroup
301:
302: Let's try something else: type \kbd{(1 + x)\pow 3}. No \kbd{O(x)} here, since
303: the result is a polynomial. Haha, but I have learnt that if you do not take
304: exponents which are integers greater or equal to 0, you obtain a power series
305: with an infinite number of non-zero terms. Let's try. Type
306: \kbd{(1 + x)\pow (-3)} (the parentheses around \kbd{-3} are not necessary but
307: make things easier to read). Surprise! Contrary to what we expected, we don't
308: get a power series but a rational function. Again this is for the same reason
309: that \kbd{1 / 7} just gave you $1/7$: the result being exact, PARI doesn't see
310: any reason to make it non-exact.
311:
312: But I still want that power series. To obtain it, just do as in the $1/7$
313: example: type
314: %
315: \bprog
316: (1 + x)\pow (-3) + O(x\pow{16})
317: (1 + O(x\pow {16})) * (1 + x)\pow (-3)
318: (1 + x + O(x\pow {16}))\pow (-3) {\rm, etc\dots}
319: \eprog
320: Now try \kbd{(1 + x)\pow (1/2)}: we obtain a power series, since the
321: result is an object which PARI does not know how to represent exactly (we
322: could teach PARI about algebraic functions, but then take \kbd{(1 + x)\pow Pi}
323: as another example). This gives us still another solution to our preceding
324: exercise: we can type \kbd{(1 + x)\pow (-3.)}. Since \kbd{-3.} is not an exact
325: quantity, PARI has no means to know that we are dealing with a rational
326: function, and will instead give you the power series, this time with real
327: instead of integer coefficients.
328:
329: Finally, if you want to be really fancy, you can type
330: \kbd{taylor((1 + x)\pow (-3), x)} (look at the entry for \kbd{taylor} for the
331: description of the syntax), but this is in fact almost never used.
332: \smallskip
333:
334: To summarize, in this section we have seen that in addition to integers, real
335: numbers and rational numbers, PARI can handle complex numbers, polynomials,
336: power series, rational functions. A large number of functions exist which
337: handle these types, but in this tutorial we will only look at a few.
338:
339: \noindent {\bf Additional comments} (as before, you are supposed to skip this
340: at first reading)
341:
342: 1) To be able to duplicate the following example, first type \b{y} to
343: suppress automatic simplification.
344:
345: A complex number has a real part and an imaginary part (who would have
346: guessed?). However, beware that when the imaginary part is the exact integer
347: zero, it is not printed, but the complex number is not converted to a real
348: number. Hence it may {\it look\/} like a real number without being one, and
349: this may cause some confusion in programs which expect real numbers. For
350: example, type \kbd{n = 3 + I - I}. The answer reads \kbd{3}, as expected. But
351: it is still a complex number. Check it with \kbd{type(n)}. Hence if you then
352: type \kbd{(1+x)\pow n}, instead of getting the polynomial
353: \kbd{1 + 3*x + 3*x\pow 2 + x\pow 3} as expected, you obtain a power series.
354: Worse, when you try to apply an arithmetic function, say the Euler totient
355: function (known as \kbd{eulerphi} to GP), you get a sententious error message
356: recalling you that ``arithmetic functions want integer arguments''. You would
357: have guessed yourself, but the message is difficult to understand since 3 looks
358: like a genuine integer!!! (Please read again if this is not clear. It is a
359: common source of incomprehension).
360:
361: Similarly, \kbd{3 + x - x} is not the integer 3 but a constant polynomial
362: (in the variable \kbd{x}), equal to $3 = 3x^0$.
363:
364: If you want the final expression to be in the simplest form possible (for
365: example before applying an arithmetic function, or simply because things will
366: go faster afterwards), apply the function \kbd{simplify} to the result. In
367: fact, by default GP does this automatically at the end of a GP command. Note
368: that it does {\it not\/} simplify in intermediate expressions. This default
369: can be switched off and on by the command \b{y}. This is why I asked you to
370: type this before starting.
371:
372: 2) As already stated, power series expansions are always implicitly around
373: $\kbd{x}=0$. When we wanted them around $\kbd{x}=\kbd{a}$, we replaced
374: \kbd{x} by \kbd{z + a} in the function we wanted to expand. For complicated
375: functions, it may be simpler to use the substitution function \kbd{subst}.
376: For example, if \kbd{p~= 1 / (x\pow 4 + 3*x\pow 3 + 5*x\pow 2 - 6*x + 7)},
377: you may not want to retype this, replacing \kbd{x} by \kbd{z~+ a}, so you can
378: write \kbd{subst(p, x, z+a)} (look up the exact description of the
379: \kbd{subst} function).
380:
381: Now try typing \kbd{p = 1 + x + x\pow 2 + O(x\pow 10)}, then
382: \kbd{subst(p, x, z+1)}. Do you understand why you get an error message?
383:
384: \section{The Remaining PARI Types}
385: Let's talk some more about the basic PARI types.
386:
387: Type \kbd{p = x * exp(-x)}. As expected, you get the power series expansion
388: to 16 terms (if you have not changed the default). Now type
389: \kbd{pr = serreverse(p)}. You are asking here for the {\it reversion\/} of the
390: power series \kbd{p}, in other words the inverse function. This is possible
391: only for power series whose first non-zero coefficient is that of $x^1$. To
392: check the correctness of the result, you can type \kbd{subst(p, x, pr)} or
393: \kbd{ subst(pr, x, p)} and you should get back \kbd{x + O(x\pow 17)}.
394:
395: Now the coefficients of \kbd{pr} obey a very simple formula. First, we would
396: like to multiply the coefficient of \kbd{x\pow n} by \kbd{n!} (in the case of
397: the exponential function, this would simplify things considerably!). The PARI
398: function \kbd{serlaplace} does just that. So type \kbd{ps = serlaplace(pr)}.
399: The coefficients now become integers, which can be immediately recognized by
400: inspection. The coefficient of $x^n$ is now equal to
401: $n^{n-1}$. In other words, we have
402: %
403: $$\kbd{pr} = \sum_{n\ge1}\dfrac{n^{n-1}}{n!} X^{n}.$$
404: %
405: Do you know how to prove this? (If you have never seen this, the proof is
406: difficult.)
407: \smallskip
408: %
409: Of course PARI knows about vectors (rows and columns are distinguished, even
410: though mathematically there is no difference) and matrices. Type for example
411: \kbd{[1,2,3,4]}. This gives the row vector whose coordinates are 1, 2, 3 and
412: 4. If you want a column vector, type \kbd{[1,2,3,4]\~{}}, the tilde meaning
413: of course transpose. You don't see much difference in the output, except for
414: the tilde at the end. However, now type \b{b}: lo and behold, the vector does
415: become a column. The \b{b} command is used mainly for this purpose.
416:
417: Type \kbd{m = [a,b,c; d,e,f]}. You have just entered a matrix with 2 rows and
418: 3 columns. Note that the matrix is entered by {\it rows\/} and the rows are
419: separated by semicolons ``\kbd{;}''. The matrix is printed naturally in a
420: rectangle shape. If you want it printed horizontally just as you typed it,
421: type \b{a}, or if you want this type of printing to be the permanent default
422: type \kbd{default(output, 1)}. Type \kbd{default(output, 0)} if you want to
423: come back to the original default.
424:
425: Now type \kbd{m[1,2]}, \kbd{m[1,]}, \kbd{m[,2]}. Are explanations necessary?
426: (In an expression such as \kbd{m[j,k]}, the \kbd{j} always refers to the
427: row number, and the \kbd{k} to the column number, and the first index is
428: always 1, never 0. This default cannot be changed.)
429:
430: Even better, type \kbd{m[1,2] = 5; m} (the semicolon also allows us to put
431: several instructions on the same line. The final result will be the output of
432: the last statement on the line). Now type \kbd{m[1,] = [15,-17,8]}. No
433: problem. Finally type \kbd{m[,2] = [j,k]}. You have an error message since you
434: have typed a row vector, while \kbd{m[,2]} is a column vector. If you type
435: instead \kbd{m[,2] = [j,k]\~{}} it works. \smallskip
436: %
437: Type now \kbd{h = mathilbert(20)}. You get the so-called ``Hilbert matrix''
438: whose coefficient of row $i$ and column $j$ is equal to $(i+j-1)^{-1}$.
439: Incidentally, the matrix \kbd{h} takes too much room. If you don't want to
440: see it, simply type a semi-colon ``\kbd{;}'' at the end of the line
441: (\kbd{h = mathilbert(20);}). This is an example of a ``precomputed'' matrix,
442: built into PARI. There are only a few. We will see later an example of a much
443: more general construction.
444:
445: What is interesting about Hilbert matrices is that first their inverses and
446: determinants can be computed explicitly (and the inverse has integer
447: coefficients), and second they are numerically very unstable, which make them
448: a severe test for linear algebra packages in numerical analysis. Of course
449: with PARI, no such problem can occur: since the coefficients are given as
450: rational numbers, the computation will be done exactly, so there cannot be
451: any numerical error. Try it. Type \kbd{d~=~matdet(h)} (you have to be a
452: little patient, this is quite a complicated computation). The result is a
453: rational number (of course) of numerator equal to 1 and denominator having
454: 226 decimal digits. How do I know, by the way? I did not count! Instead,
455: simply type \kbd{1.* d}. The result is now in exponential format, and the
456: exponent gives us the answer. Another, more natural, way would be to simply
457: type \kbd{sizedigit(1/d)}.
458:
459: Now type \kbd{hr = 1.* h;} (do not forget the semicolon, we don't want to see
460: all the junk!), then \kbd{dr = matdet(hr)}. You notice two things. First the
461: computation, although not instantaneous, is much faster than in the rational
462: case. The reason for this is that PARI is handling real numbers with 28
463: digits of accuracy, while in the rational case it is handling integers having
464: up to 226 decimal digits.
465:
466: The second more important fact is that the result is terribly wrong. If you
467: compare with \kbd{1.$*$d} computed earlier, which is correct, you will see
468: that only 2 decimals agree! This catastrophic instability is as already
469: mentioned one of the characteristics of Hilbert matrices. In fact, the
470: situation is much worse than that. Type \kbd{norml2(1/h - 1/hr)} (the
471: function \kbd{norml2} gives the square of the $L^2$ norm, i.e.~the sum of the
472: squares of the coefficients). Again be patient since computing \kbd{1/h} will
473: take even more time (not much) than computing \kbd{matdet(h)}. The result is
474: larger than $10^{50}$, showing that some coefficients of \kbd{1/hr} are
475: wrong by as much as $10^{24}$ (the largest error is in fact equal to $7.644
476: \cdot 10^{24}$ for the coefficient of row 15 and column 14, which is a 28
477: digit integer).
478:
479: To obtain the correct result after rounding for the inverse, we have to use a
480: default precision of 56 digits (try it).
481: \smallskip
482:
483: Although vectors and matrices can be entered manually, by typing explicitly
484: their elements, very often the elements satisfy a simple law and one uses a
485: different syntax. For example, assume that you want a vector whose $i$-th
486: coordinate is equal to $i^2$. No problem, type for example
487: \kbd{vector(10,i, i\pow 2)} if you want a vector of length 10. Similarly, if
488: you type
489:
490: \centerline{\kbd{matrix(5,5,i,j, 1/(i+j-1))}}
491:
492: \noindent you will get the Hilbert matrix of order 5 (hence the
493: \kbd{mathilbert} function is redundant). The \kbd{i} and \kbd{j} represent
494: dummy variables which are used to number the rows and columns respectively
495: (in the case of a vector only one is present of course). You must not forget,
496: in addition to the dimensions of the vector or matrix, to indicate explicitly
497: the names of these variables.
498:
499: {\bf Warning.} The letter \kbd{I} is reserved for the complex number equal to
500: the square root of $-1$. Hence it is absolutely forbidden to use it as a
501: variable. Try typing \kbd{vector(10,I, I\pow 2)}, the error message that you
502: get clearly indicates that GP does not consider \kbd{I} as a variable. There
503: are two other reserved variable names: \kbd{Pi} and \kbd{Euler}. All function
504: names are forbidden as well. On the other hand there is nothing special about
505: \kbd{i}, \kbd{pi} or \kbd{euler}.
506:
507: When creating vectors or matrices, it is often useful to use boolean
508: operators and the \kbd{if()} statement (see the section on programming for
509: details on using this statement). Indeed, an \kbd{if} expression has a value,
510: which is of course equal to the evaluated part of the \kbd{if}. So for
511: example you can type
512:
513: \centerline{\kbd{matrix(8,8,i,j, if((i-j)\%2, x, 0))} }
514: \noindent to get a checkerboard matrix of \kbd{x} and \kbd{0}. Note however
515: that a vector or matrix must be {\it created} first before being used. For
516: example, it is forbidden to write
517:
518: \centerline{\kbd{for(i=1,5, v[i]=1/i)}}
519: \noindent if the vector \kbd{v} has not been created beforehand (for example
520: by a \kbd{v = vector(5,j,0)} command).
521:
522: \medskip The last PARI types which we have not yet played with are types
523: which are closely linked to number theory (hence people not interested in
524: number theory can skip ahead).
525:
526: The first is the type ``integer--modulo''. Let us see an example. Type
527: \kbd{n = 10\pow 15 + 3}. We want to know whether this number is prime or not. Of
528: course we could make use of the built-in facilities of PARI, but let us do
529: otherwise. (Note that PARI does not actually prove that a number is prime: any
530: strong pseudoprime for a number of bases is declared to be prime.) We first
531: trial divide by the built-in table of primes. We slightly cheat here and use a
532: variant of the function \kbd{factor} which does exactly this. So type
533: \kbd{factor(n, 200000)} (the last argument tells \kbd{factor} to trial divide
534: up to the given bound and stop at this point. You can set it to 0 to trial
535: divide by the full set of built-in primes, which goes up to $500000$ by
536: default).
537:
538: The result is a 2 column matrix (as for all factoring functions), the first
539: column giving the primes and the second their exponents. Here we get a single
540: row, telling us that \kbd{n} is not divisible by any prime up to $200000$. We
541: could now trial divide further, or even cheat completely and call the PARI
542: function \kbd{factor} without the optional second argument, but before we do
543: this let us see how to get an answer ourselves.
544:
545: By Fermat's little theorem, if $n$ is prime we must have $a^{n-1}\equiv 1
546: \pmod{n}$ for all $a$ not divisible by $n$. Hence we could try this with $a=2$
547: for example. But $2^{n-1}$ is a number with approximately $3\cdot10^{14}$
548: digits, hence impossible to write down, let alone to compute. But instead type
549: \kbd{a = Mod(2,n)}. This creates the number $2$ considered now as an element
550: of the ring $R = \Z/\kbd{n}\Z$. The elements of $R$, called integermods, can
551: always be represented by numbers smaller than \kbd{n}, hence very small.
552: Fermat's theorem can be rewritten
553: %
554: $\kbd{a}^{n-1} = \kbd{Mod(1,n)}$
555: %
556: in the ring $R$, and this can be computed very efficiently. Type
557: \kbd{a\pow (n-1)}. The result is definitely {\it not\/} equal to
558: \kbd{Mod(1,n)}, thus {\it proving\/} that \kbd{n} is not a prime (if we had
559: obtained \kbd{Mod(1,n)} on the other hand, it would have given us a hint that
560: \kbd{n} is maybe prime, but never a proof).
561:
562: To find the factors is another story. One must use less naive techniques than
563: trial division (or be very patient). To finish this example, type
564: \kbd{factor(n)} to see the factors. Since the smallest factor is 14902357,
565: you would have had to be very patient with trial division! Note that none of
566: the factors in the decomposition are {\it proven} primes.
567: \smallskip
568:
569: The second specifically number-theoretic type is the $p$-adic numbers. I have
570: no room for definitions, so please skip ahead if you have no use for such
571: beasts. A $p$-adic number is entered as a rational or integer valued
572: expression to which is added \kbd{O(p\pow n)} (or simply \kbd{O(p)} if
573: $\kbd{n}=1$) where \kbd{p} is the prime and \kbd{n} the $p$-adic precision.
574: Apart from the usual arithmetic operations, you can apply a number of
575: transcendental functions. For example, type \kbd{n = 569 + O(7\pow 8)}, then
576: \kbd{s~=~sqrt(n)}, you obtain one of the square roots of \kbd{n} (if you want
577: to check, type \kbd{s\pow 2 - n}). Type now \kbd{l~=~log(n)}, then \kbd{e =
578: exp(l)}. If you know about $p$-adic logarithms, you will not be surprised that
579: \kbd{e} is not equal to \kbd{n}. Type \kbd{(n/e)\pow 6}: \kbd{e} is in fact
580: equal to \kbd{n} times a $(p-1)$-st root of unity.
581:
582: Incidentally, if you want to get back the integer 569 from the $p$-adic
583: number \kbd{n}, type \kbd{truncate(n)}.
584: \smallskip
585:
586: The third number-theoretic type is the type ``quadratic number''. This type
587: is specially tailored so that we can easily work in a quadratic extension of
588: a base field (usually $\Q$). It is a generalization of the type
589: ``complex''. To start, we must specify which quadratic field we want to work
590: in. For this, we use the function \kbd{quadgen} applied to the
591: {\it discriminant\/} \kbd{d} (as opposed to the radicand) of the quadratic
592: field. This returns a number (always printed as \kbd{w}) equal to
593: $(\kbd{d}+a) / 2$ where $a$ is equal to 0 or 1 according to whether \kbd{d} is
594: even or odd. The behavior of \kbd{quadgen} is a little special: although its
595: result is always printed as \kbd{w}, the variable \kbd{w} itself is not set
596: to that value. Hence it is necessary to write systematically
597: \kbd{w = quadgen(d)} using the variable name \kbd{w} (or \kbd{w1} etc. if you
598: have several quadratic fields), otherwise things will get confusing.
599:
600: So type \kbd{w = quadgen(-163)}, then \kbd{charpoly(w)} which asks for the
601: characteristic polynomial of \kbd{w} (in the variable \kbd{x};
602: you can type \kbd{charpoly(w, y)} to directly express it in terms of some
603: other variable without resorting to the \kbd{subst} function). The result
604: shows what \kbd{w} will represent. You can also ask for \kbd{1.*w} to see
605: which root of the quadratic has been taken, but this is rarely necessary. We
606: can now play in the field $\Q(\sqrt{-163})$. Type for example
607: \kbd{w\pow 10}, \kbd{norm(3 + 4*w)}, \kbd{1 / (4+w)}. More interesting, type
608: \kbd{a = Mod(1,23) * w} then \kbd{b = a\pow 264}. This is a generalization of
609: Fermat's theorem to quadratic fields. If you do not want to see the modulus 23
610: all the time, type \kbd{lift(b)}.
611:
612: Another example: type \kbd{p = x\pow 2 + w*x + 5*w + 7}, then \kbd{norm(p)}. We
613: thus obtain the quartic equation over $\Q$ corresponding to the relative
614: quadratic extension over $\Q(\kbd{w})$ defined by \kbd{p}.
615:
616: On the other hand, if you type \kbd{wr = sqrt(w\pow 2)}, do not expect to get
617: back \kbd{w}. Instead, you get the numerical value, the function \kbd{sqrt}
618: being considered as a ``transcendental'' function, even though it is
619: algebraic. Type \kbd{algdep(wr,2)} (this looks for algebraic relations
620: involving the powers of \kbd{w} up to degree 2). This is one way to get
621: \kbd{w} back. Similarly, type \kbd{algdep(sqrt(3*w + 5), 4)}. See the user's
622: manual for the function \kbd{algdep}.\smallskip
623:
624: The fourth number-theoretic type is the type ``polynomial--modulo'', i.e.
625: polynomial modulo another polynomial. This type is used to work in general
626: algebraic extensions, for example elements of number fields (if the base
627: field is $\Q$), or elements of finite fields (if the base field is
628: $\Z/p\Z$ for a prime $p$, defined by an integermod). In a sense it is a
629: generalization of the type quadratic number. The syntax used is the same as
630: for integermods. For example, instead of typing \kbd{w = quadgen(-163)}, you
631: can type \kbd{w = Mod(x, x\pow 2 - x + 41)}. Then, exactly as in the
632: quadratic case, you can type \kbd{w\pow 10}, \kbd{norm(3 + 4*w)},
633: \kbd{1 / (4+w)}, \kbd{a = Mod(1,23)*w}, \kbd{b = a\pow 264}, obtaining of
634: course the same results (type \kbd{lift(\dots)} if you don't want to see the
635: polynomial \kbd{x\pow 2 - x + 41} repeated all the time). Of course, the basic
636: interest is that you can work in any degree, not only quadratic (of course,
637: even for quadratic moduli, the corresponding elementary operations will be
638: slower than with quadratic numbers).
639:
640: There is however a slight difference in behavior. Keeping our polmod \kbd{w},
641: type \kbd{1.*w}. As you can see, the result is not the same. Type
642: \kbd{sqrt(w)}. Here, we obtain a vector with 2 components, the two components
643: being the principal branch of the square root of all the possible embeddings
644: of \kbd{w} in $\C$ ({\it not\/} the two square roots). More generally, if
645: \kbd{w} was of degree $n$, we would get an $n$-component vector, and similarly
646: for transcendental functions.
647:
648: We have at our disposal the usual arithmetic functions, plus a few others.
649: Type \kbd{a = Mod(x, x\pow 3 - x - 1)} defining a cubic extension. We can for
650: example ask for \kbd{b = a\pow 5}. Now assume we want to express \kbd{a}
651: as a polynomial in \kbd{b}. This is possible since \kbd{b} is also a
652: generator of the same field. No problem, type \kbd{modreverse(b)}. This gives
653: a new defining polynomial for the same field (i.e.~the characteristic
654: polynomial of \kbd{b}), and expresses \kbd{a} in terms of this new polmod,
655: i.e.~in terms of \kbd{a}. We will see this in much more detail in the number
656: field section.
657:
658: \section{Elementary Arithmetic Functions}
659:
660: Since PARI is aimed at number theorists, it is not surprising that there
661: exists a large number of arithmetic functions (see the list in the
662: corresponding section of the users manual). We have already seen several,
663: such as \kbd{factor}. Note that \kbd{factor} handles not only integers, but
664: also (univariate) polynomials. Type for example \kbd{factor(x\pow 15 - 1)}.
665: You can also ask to factor a polynomial modulo a prime $p$ (\kbd{factormod})
666: and even in a finite field which is not a prime field (\kbd{factorff}).
667:
668: Evidently, you have functions for computing GCD's (\kbd{gcd}), extended GCD's
669: (\kbd{bezout}), solving the Chinese remainder theorem (\kbd{chinese}) and so
670: on.
671:
672: In addition to the factoring facilities, you have a few functions related to
673: primality testing such as \kbd{isprime}, \kbd{ispseudoprime},
674: \kbd{precprime}, and \kbd{nextprime}. As previously mentioned, only strong
675: pseudoprimes are produced; there is no sophisticated primality test.
676:
677: We also have the usual multiplicative arithmetic functions: the M\"obius $\mu$
678: function (\kbd{moebius}), the Euler $\phi$ function (\kbd{eulerphi}), the
679: $\omega$ and $\Omega$ functions (\kbd{omega} and \kbd{bigomega}), the
680: $\sigma_k$ functions (\kbd{sigma}), which compute sums of $k$-th powers of the
681: positive divisors of a given integer, etc\dots
682:
683: You can compute continued fractions. For example, type \kbd{\b{p} 1000}, then
684: \kbd{contfrac(exp(1))}: you obtain the continued fraction of the base of
685: natural logarithms, which as you can see obeys a very simple pattern (can you
686: prove it?).
687:
688: In many cases, one wants to perform some task only when an arithmetic
689: condition is satisfied. GP gives you the following functions: \kbd{isprime}
690: as mentioned above, \kbd{issquare}, \kbd{isfundamental} to test whether an
691: integer is a fundamental discriminant (i.e.~$1$ or the discriminant of the
692: ring of integers of a quadratic field), and the \kbd{forprime}, \kbd{fordiv}
693: and \kbd{sumdiv} loops. Assume for example that we want to compute the
694: product of all the divisors of a positive integer \kbd{n}. The easiest way is
695: to write
696:
697: \centerline{\tt p=1; fordiv(n,d, p *= d); p }
698:
699: \noindent
700: (there is a very simple formula for this product: find and prove it). The
701: notation \kbd{p *= d} (inherited from the C programming language) is just a
702: shorthand for \kbd{p = p * d}.
703:
704: If we want to know the list of primes $p$ less than 1000 such that 2 is a
705: primitive root modulo $p$, one way would be to write:
706: %
707: \bprog
708: forprime(p=3,1000, if (znprimroot(p) == 2, print(p)))
709: \eprog\noindent
710: %
711: Note that this assumes that \kbd{znprimroot} returns the smallest primitive
712: root, and this is indeed the case. Had we not known about this, we could
713: have written
714: %
715: \bprog
716: forprime(p=3,1000, if (znorder(Mod(2,p)) == p-1, print(p)))
717: \eprog\noindent
718: %
719: (which is actually faster since we only compute the order of $2$ in $\Z/p\Z$,
720: instead of looking for a generator by trying successive elements whose orders
721: have to be computed as well.)
722:
723: Functions related to quadratic fields, binary quadratic forms and general
724: number fields will be seen in the next sections.
725:
726: \section{Performing Linear Algebra}
727: All the standard linear algebra programs are available of course, and many
728: more. In addition, linear algebra over $\Z$, i.e.~work on lattices, can also
729: be performed. Let us see how this works. First of all, a vector space (more
730: generally a module) will be given by a generating set of vectors (often a
731: basis) which will be representend as {\it column} vectors. This set of vectors
732: will in turn be represented as a matrix: in PARI, we have chosen to consider
733: matrices as row vectors of column vectors. The base field (or ring) can be any
734: ring type PARI supports (except $p$-adics which are currently not correctly
735: handled by the linear algebra package). However, certain operations are
736: specifically written for a real or complex base field, while others are
737: written for $\Z$ as the base ring.
738:
739: ----- TO BE COMPLETED -----
740:
741:
742:
743: \section{Using Transcendental Functions}
744:
745: All the elementary transcendental functions and several higher transcendental
746: functions (gam\-ma function, incomplete gamma function, error function,
747: exponential integral, $K$-bessel functions, confluent hypergeometric functions,
748: Riemann $\zeta$ function, polylogarithms, Weber functions, theta functions)
749: are provided. More may be written if the need arises.
750:
751: In this type of functions, the default precision plays an essential role.
752: In almost all cases transcendental functions work in the following way.
753: If the argument is exact, the result will be computed using the current
754: default precision. If the argument is not exact, the precision of the
755: argument is used for the computation. A note of warning however: even in this
756: case the {\it printed\/} value will be the current real format (usually the
757: same as the default precision). In the present chapter we assume that your
758: machine works with 32-bit long integers. If it is not the case, we leave it
759: to you as a very good exercise to make the necessary modifications.
760:
761: Let's assume that we have 28 decimals of default precision (this is what we
762: get automatically at the start of a GP session on 32-bit machines). Type
763: \kbd{e = exp(1)}. We get the number $e=2.718\dots$ to 28 decimals. Let us check
764: how many correct decimals we really have. The hard (but reasonable) way is to
765: proceed as follows. Change the precision to a substantially higher value, for
766: example by typing \kbd{\b{p} 50}. Then type \kbd{e}, then \kbd{exp(1)} once
767: again. This last value is the correct value of the mathematical constant $e$ to
768: 50 decimals, while the variable \kbd{e} shows the value that was computed to 28
769: decimals. Clearly they coincide to exactly 29 significant digits.
770:
771: A simpler way to see how many decimals we have, is to ask for \kbd{length(e)}
772: which indicates we have exactly $3$ mantissa words. Since
773: $3\ln(2^{32}) / \ln(10)\approx28.9$ we see that we have 28 or 29 significant
774: digits (on 32-bit machines).
775:
776: \smallskip
777: Come back to 28 decimals (\kbd{\b{p} 28}). If we type \kbd{exp(1.)}
778: you can check that we also obtain 28 decimals. However, type
779: \kbd{f = exp(1 + 10.\pow(-30))}. Although the default precision is still 28,
780: you can check using one of the two methods above that we have in fact 59
781: significant digits! The reason is that \kbd{1 + 10.\pow(-30)} is computed
782: according to the PARI philosophy, i.e.~to the best possible precision. Since
783: \kbd{10.} has 29 significant digits and 1 has ``infinite'' precision, the
784: number \kbd{1 + 10.\pow(-30)} will have $59=29+30$ significant digits,
785: hence \kbd{f} also.
786:
787: Now type \kbd{cos(10.\pow(-15))}. The result is printed as $1.0000\dots$, but
788: is of course not exactly equal to 1. Using \kbd{length(\%)}, we see that the
789: result has 7 mantissa words, giving us the possibility of having 67
790: correct significant digits. In fact (look in precision 100), only 60 are
791: correct. PARI gives you as much as it can, and since 6 mantissa words
792: would have given you only 57 digits, it uses 7. But why does it give so
793: precise a result? Well, it is the same reason as before. When $x$ is close
794: to 1, $\cos(x)$ is close to $1-x^2/2$, hence the precision is going to be
795: approximately the same as this quantity, which will be $1-0.5*10^{-30}$ where
796: $0.5*10^{-30}$ is considered with 28 significant digit accuracy, hence the
797: result will have approximately $28+30=58$ significant digits.
798:
799: Unfortunately, this philosophy cannot go too far. For example, when you
800: type \kbd{cos(0)}, GP should give you exactly 1. Since it is reasonable for
801: a program to assume that a transcendental function never gives you an exact
802: result, GP gives you $1.000\dots$ to one more mantissa word than the current
803: precision.
804: \medskip
805: OK, now let's see some more transcendental functions at work. Type
806: \kbd{gamma(10)}. No problem (type \kbd{9!} to check). Type \kbd{gamma(100)}.
807: The number is now written in exponential format because the default
808: accuracy is too small to give the correct result (type \kbd{99!} to get the
809: complete answer).
810: To get the complete value, there are two solutions. The first and most natural
811: one is to increase the precision. Since \kbd{gamma(100)} has 156 decimal
812: digits, type \kbd{\b{p} 170} (to be on the safe side), then \kbd{gamma(100)}
813: once again. After some work, the printed result is this time perfectly
814: correct.
815:
816: However, this is probably not the best way to proceed. Come back first to the
817: default precision (type \kbd{\b{p} 28}). As the gamma function increases
818: very rapidly, one usually uses its logarithm. Type \kbd{lngamma(100)}. This
819: time the result has a reasonable size, and is exactly equal to \kbd{log(99!)}.
820:
821: Try \kbd{gamma(1/2 + 10*I)}. No problem, we have the complex gamma function.
822: Now type
823:
824: \kbd{t = 1000; z = gamma(1 + I*t) * t\pow(-1/2) * exp(Pi/2*t)/sqrt(2*Pi)},
825:
826: \noindent then \kbd{norm(z)}. We see that \kbd{norm(z)} is very close to 1,
827: in accordance with the complex Stirling formula. \smallskip
828: %
829: Let's play now with the Riemann zeta function. First turn on the timer (type
830: \kbd{\#}). Type \kbd{zeta(2)}, then \kbd{Pi\pow 2/6}. This seems correct. Type
831: \kbd{zeta(3)}. All this takes essentially no time at all. However, type
832: \kbd{zeta(3.)}. Although the result is the same, you will notice that the
833: time is substantially larger (if your machine is too fast to see the
834: difference, increase the precision!). This is because PARI uses special
835: formulas to compute \kbd{zeta(n)} when \kbd{n} is a (positive or negative)
836: integer.
837:
838: Type \kbd{zeta(1 + I)}. This also works. Now for fun, let us compute in a
839: very naive way the first complex zero of \kbd{zeta}. We know that it is
840: of the form $1/2 + i*t$ with $t$ between 14 and 15. Thus, we can use the
841: following series of instructions. But instead of typing them directly, write
842: them into a file, say \kbd{zeta.gp}, then type \kbd{\b{r} zeta.gp} under GP to
843: read it in:
844:
845: \bprog\obr
846: \q t1 = 1/2 + 14*I;
847: \q t2 = 1/2 + 15*I; eps = 10\pow-50;
848: \q z1 = zeta(t1);
849: \q until (norm(z2) < eps,
850: \q\q z2 = zeta(t2);
851: \q\q if (norm(z2) < norm(z1),
852: \q\q\q t3 = t1; t1 = t2; t2 = t3; z1 = z2
853: \q\q );
854: \q\q t2 = (t1+t2) / 2.;
855: \q\q print(t1 ": " z1)
856: \q )
857: \cbr\eprog
858:
859: Don't forget the braces: they tell GP that a sequence of instructions is going
860: to span many lines (another, less convenient, way would be to type \b{} at the
861: end of each line to tell the parser that the input was not yet finished).
862: By the way, you don't need to type in the suffix~\kbd{.gp} it will be
863: supplied by GP, if you forget it (the suffix is not mandatory either, it is
864: just more convenient to have all your GP scripts labeled in the same
865: distinctive way).
866:
867: We thus obtain the first zero to 25 significant digits.
868: \medskip
869: %
870: As mentioned at the beginning of this tutorial, some transcendental functions
871: can also be applied to $p$-adic numbers. This is as good a time as any to
872: familiarize yourself with them. Type \kbd{a~=~exp(7 + O(7\pow 10))}, then
873: \kbd{log(a)}. All seems in order. Now type \kbd{b = log(5 + O(7\pow 10))},
874: then \kbd{exp(b)}. Is something wrong, since we don't recover the number we
875: started with? Absolutely not. Type
876: \kbd{exp(b) * teichmuller(5 + O(7\pow 10))},
877: and we indeed recover our initial number. The function \kbd{teichmuller(x)}
878: is the Teichm\"uller character, and is characterized by the fact that it is
879: the unique \hbox{$(p-1)$-st} root of unity (here with $p=7$) which is
880: congruent to \kbd{x} modulo $p$, assuming that \kbd{x} is a $p$-adic
881: unit.\smallskip
882: %
883: Let us come back to real numbers for the moment. Type \kbd{agm(1,sqrt(2))}.
884: This gives the arithmetic-geometric mean of 1 and $\sqrt2$, and is the basic
885: method for computing (complete) elliptic integrals. In fact, type
886:
887: \kbd{Pi/2 / intnum(t=0,Pi/2, 1 / sqrt(1 + sin(t)\pow 2))},
888:
889: \noindent and the result is the same. The elementary transformation
890: \kbd{x = sin(t)} gives the mathematical equality
891: $$\int_0^1 \dfrac{dx}{\sqrt{1-x^4}} = \dfrac{\pi}{2\text{agm}(1,\sqrt2)}
892: \enspace,$$
893: which was one of Gauss's remarkable discoveries in his youth.
894:
895: Now type \kbd{2 * agm(1,I) / (1+I)}. As you see, the complex AGM also works,
896: although one must be careful with its definition. The result found is
897: almost identical to the previous one. Exercise: do you see why?
898:
899: Finally, type \kbd{agm(1, 1 + 7 + O(7\pow 10))}. So we also have $p$-adic
900: AGM. Note however that since the square root of a $p$-adic number is not
901: in general an element of the same $p$-adic field,
902: only certain $p$-adic AGMs can be computed. In addition,
903: when $p=2$, the congruence restriction is that \kbd{agm(a,b)} can be computed
904: only when \kbd{a/b} is congruent to 1 modulo $16$ (not 8 as could be
905: expected).\smallskip
906: %
907: Now type \kbd{?3}. This gives you the list of all transcendental functions.
908: Instead of continuing with more examples, we suggest that you experiment
909: yourself with the list of functions. In each case, try integer, real, complex
910: and $p$-adic arguments. You will notice that some have not been implemented
911: (or do not have a reasonable definition).
912:
913: \section{Using Numerical Tools}
914:
915: Although not written to be a numerical analysis package, PARI can
916: nonetheless perform some numerical computations. We leave for a subsequent
917: section linear algebra and polynomial computations.
918:
919: You of course know the formula $\pi = 4(1-\dfrac13+\dfrac15-\dfrac17+\cdots)$
920: which is deduced from the power series expansion of \kbd{atan(x)}. You also
921: know that $\pi$ cannot be computed from this formula, since the convergence
922: is so slow. Right? Wrong! Type \kbd{\b{p} 100} (just to make it more
923: interesting), then \kbd{4~*~sumalt(k=0, (-1)\pow k/(2*k + 1))}. In a split
924: second (admittedly more than simply typing \kbd{Pi}), we get $\pi$ to 100
925: significant digits (type \kbd{Pi} to check). In version 1.38, the method used
926: was a combination of a method due to Euler for accelerating alternating sums,
927: and a programming trick due to the Dutch mathematician van Wijngaarden (see
928: one of the Numerical Recipes books for an explanation). The method which we
929: presently use is considerably better, and is based on a combination of ideas of
930: F.~Villegas, D.~Zagier and H.~Cohen.
931:
932: Similarly, type \kbd{\b{p} 28} (otherwise the time will be too long) and
933: \kbd{sumpos(k=1, 1 / k\pow 2)}. Although once again the convergence is slow, a
934: similar trick allows to compute the sum when the terms are positive (compare
935: with the exact result \kbd{Pi\pow 2/6}). This is much less impressive because
936: it is much slower, but is still useful.
937:
938: Even better, \kbd{sumalt} can be used to sum divergent series! Type
939:
940: \centerline{\tt zet(s)= sumalt(k=1, (-1)\pow(k-1) / k\pow s) / (1 - 2\pow(1-s))}
941:
942: Then for positive values of \kbd{s} different from 1, \kbd{zet(s)} is equal
943: to \kbd{zeta(s)} and the series cvonverges, albeit slowly (sumalt doesn't
944: care however). For negative \kbd{s}, the series diverges, but \kbd{zet(s)}
945: still gives the correct result! Try \kbd{zet(-1)}, \kbd{zet(-2)},
946: \kbd{zet(-1.5)}, and compare with the corresponding values of \kbd{zeta}.
947: You should not push the game too far: \kbd{zet(-14.5)}, for example,
948: gives a completely wrong answer.
949:
950: Try \kbd{zet(I)}, and compare with \kbd{zeta(I)}. Even (some) complex values
951: work, although the sum is not alternating any more!
952:
953: Similarly, try \kbd{sumalt(n=1, (-1)\pow n/(n+I))}.
954: \medskip
955: %
956: More traditional functions are the numerical integration functions.
957: For example, type \kbd{intnum(t=1,2, 1/t)} and presto! you get 26 decimals
958: of $\log(2)$. Look at Chapter 3 to see the different integration functions
959: which are available.
960:
961: With PARI, however, you can go further since complex types are allowed.
962: For example, assume that we want to know the location of the zeros of the
963: function $h(z)=e^z-z$. We use Cauchy's theorem, which tells us that the
964: number of zeros in a disk of radius $r$ centered around the origin is
965: equal to
966: $$\dfrac{1}{2i\pi}\int_{C_r}\dfrac{h'(z)}{h(z)}\,dz\enspace,$$
967: where $C_r$ is the circle of radius $r$ centered at the origin.
968: Hence type
969:
970: \bprog fun(z) =
971: \obr
972: \q local(u);
973: \q u = exp(z);
974: \q (u-1) / (u-z)
975: \cbr
976: zero(r) = r/(2*Pi) * intnum(t=0, 2*Pi, fun(r*exp(I*t)) * exp(I*t))
977: \eprog
978: \noindent (Here \kbd{u} is a local variable to the function \kbd{f}: whenever
979: a function is called, GP fills its argument list with the actual arguments
980: given, and initializes the other declared parameters and local variables to
981: 0. It will then restore their former values upon exit. If, however, we had
982: not declared \kbd{u} in the function prototype, it would be considered as a
983: global variable, whose value would be permanently changed. It is not
984: mandatory to declare in this way all the parameters you use, but beware of
985: side effects!)
986:
987: The function \kbd{zero(r)} will then count the number of zeros (we have simply
988: made the change of variable $z = r*\exp(i*t)$).
989:
990: Now type \kbd{\b{p} 9} (otherwise the computation would take too long, and
991: anyway we know that the result is an integer), then \kbd{zero(1)}, \kbd{zero(1.5)}.
992: The result tells us that there are no zeros inside the unit disk, but that
993: there are two (necessarily complex conjugate) whose modulus is between $1$
994: and $1.5$. For the sake of completeness, let us compute them. Let $z$ be such
995: a zero, and write $z=x+iy$ with $x$ and $y$ real. Then the equation $e^z-z=0$
996: implies, after elementary transformations, that $e^{2x}=x^2+y^2$ and that
997: $e^x\cos(y)=x$. Hence $y=\pm\sqrt{e^{2x}-x^2}$ and hence
998: $e^x\cos(\sqrt{e^{2x}-x^2})=x$. Therefore, type
999:
1000: \bprog fun(x) =
1001: \obr
1002: \q local(u);
1003: \q u = exp(x);
1004: \q u*cos(sqrt(u\pow2 - x\pow2)) - x
1005: \cbr\eprog
1006:
1007: Then \kbd{fun(0)} is positive while \kbd{fun(1)} is negative. Come back to
1008: precision 28 and type
1009:
1010: \centerline{\tt x0 = solve(x=0,1, fun(x))}
1011:
1012: This quickly gives us the value of \kbd{x}, and we can then type
1013:
1014: \centerline{\tt z = x0 + I*sqrt(exp(2*x0) - x0\pow 2)}
1015:
1016: \noindent which (together with its complex conjugate) is the required zero.
1017: As a check, type \kbd{exp(z) - z}, and also \kbd{abs(z)}.
1018:
1019: Of course you can integrate over contours which are more complicated than
1020: circles, but you must perform yourself the changes of variable as we have
1021: done above to reduce the integral to a number of integrals on line segments.
1022: \smallskip
1023: %
1024: The example above also shows the use of the \kbd{solve} function. To use
1025: \kbd{solve} on functions of a complex variable, it is necessary to reduce the
1026: problem to a real one. For example, to find the first complex zero of the
1027: Riemann zeta function as above, we could try typing
1028:
1029: \kbd{solve(t=14,15, real(zeta(1/2 + I*t)))},
1030:
1031: \noindent but this would not work because the real part is positive for
1032: \kbd{t=14} and \kbd{t=15}. As it happens, the imaginary part works. Type
1033:
1034: \kbd{solve(t=14,15, imag(zeta(1/2 + I*t)))},
1035:
1036: \noindent and this now works. We could also narrow the search interval and
1037: type for instance
1038:
1039: \kbd{solve(t=14,14.2, real(zeta(1/2 + I*t)))}
1040:
1041: \noindent which would also work.
1042:
1043: \section{Functions Related to Polynomials and Power Series}
1044:
1045: First a word of warning to the unwary: it is essential to understand the
1046: crucial difference between exact and inexact objects (see Section~1). This
1047: is especially important in the case of polynomials. Let's immediately take
1048: a plunge into these problems. Type
1049:
1050: \centerline{\tt gcd(x\pow 2 - 1, x\pow 2 - 3*x + 2)}
1051:
1052: No problem, the result is \kbd{x - 1} as expected. But now type
1053:
1054: \centerline{\tt gcd(x\pow 2 - 1., x\pow2 - 3.*x + 2.)}
1055:
1056: You are lucky, the result is almost correct except for a bizarre factor of
1057: 3 which comes from the way PARI does the computation. In any case, it is still
1058: essentially a reasonable result. But now type
1059: \kbd{gcd(x - Pi, x\pow 2 - 6*zeta(2))}.
1060: Although this should be equal to \kbd{x - Pi}, PARI finds a
1061: constant as a result. This is because the notion of GCD of non-exact
1062: polynomials doesn't make much sense. However, type
1063: \kbd{polresultant(x - Pi, x\pow 2 - 6*zeta(2))}.
1064: The result is extremely close to zero, showing that indeed the GCD is
1065: non-trivial, without telling us what it is. This being said, we will usually
1066: use polynomials (and power series) with exact coefficients in our
1067: examples.\smallskip
1068:
1069: Set \kbd{pol = polcyclo(15)}. This gives the $15$-th cyclotomic polynomial,
1070: which is of degree $\varphi(15)=8$. Now, type \kbd{r = polroots(pol)}. You
1071: have the 8 complex roots of pol given to 28 significant digits. To see them
1072: better, type \b{b}. As you see, they are given as pairs of complex conjugate
1073: roots, in a random order. The only ordering done by the function
1074: \kbd{polroots} concerns the real roots, which are given first, and in
1075: increasing order.
1076:
1077: The roots of \kbd{pol} are by definition the primitive $15$-th roots of unity.
1078: To check this, simply type \kbd{rc = r\pow 15}. Why, we get an error message!
1079: Well, fair enough, vectors cannot be multiplied (even less raised to a power)
1080: that easily. However, type \kbd{rc = r\pow 15.} with a $.$ at the end. Now it
1081: works, because powering to a non-integer (here real) exponent is a
1082: transcendental function and hence is applied termwise. Note that the fact that
1083: $15.$ is a real number which is representable exactly as an integer has
1084: nothing to do with the problem.
1085:
1086: We see that the components of the result are very close to 1. It is however
1087: tedious to look at all these real and imaginary parts. It would be impossible
1088: if we had many more. Let's do it automatically. Type
1089: \kbd{rr = round(real(rc))}, then \kbd{sqrt(norml2(rc-rr))}. We see that
1090: \kbd{rr} is indeed all 1's, and that the L2-norm of \kbd{rc - rr} is around
1091: $10^{-27}$, reasonable enough when we work with 28 significant digits. Note
1092: that the function \kbd{norml2}, contrary to what it may seem, does not give the
1093: L2 norm but its square, hence we must take the square root (well, this is not
1094: absolutely necessary in this case!).
1095: %
1096: \smallskip
1097: Now type \kbd{pol = x\pow 5 + x\pow 4 + 2*x\pow 3 - 2*x\pow 2 - 4*x - 3},
1098: then \kbd{factor(pol)}. The polynomial \kbd{pol} factors over $\Q$ (or $\Z$)
1099: as a product of two factors. Type \kbd{fun(p)= factorpadic(pol,p,10)}. This
1100: creates a function \kbd{fun(p)} which factors \kbd{pol} over $\Q_p$ to $p$-adic
1101: precision 10. If now we type \kbd{factor(poldisc(pol))}, we learn that the
1102: primes dividing the discriminant are $11$, $23$ and $37$. Type \kbd{fun(5)},
1103: \kbd{fun(11)}, \kbd{fun(23)}, and \kbd{fun(37)} to see different splittings.
1104:
1105: Similarly, we can type \kbd{lf(p)= lift(factormod(pol,p))}, and
1106: \kbd{lf(2)}, \kbd{lf(11)}, \kbd{lf(23)} and \kbd{lf(37)} which show the
1107: different factorizations, this time over $\F_p$. In fact, even better: type
1108: successively
1109: \bprog
1110: pol2 = x\pow 3 + x\pow 2 + x - 1
1111: fq = factorff(pol2, 3, t\pow 3 + t\pow 2 + t - 1)
1112: centerlift(lift(fq))
1113: \eprog
1114: %
1115: This factors the polynomial \kbd{pol2} over the finite field $\F_3(\theta)$
1116: where $\theta$ is a root of the polynomial \kbd{t\pow 3 + t\pow 2 + t - 1}.
1117: This is of course a form of the field $\F_{27}$. We know that
1118: Gal$(\F_{27}/\F_3)$ is cyclic of order 3 generated by the Frobenius
1119: homomorphism $u\mapsto u^3$, and the roots that we have found give the action
1120: of the powers of the Frobenius on \kbd{t} (if you do not know what I am
1121: talking about, please try some more examples, it's not so hard to figure out).
1122:
1123: Similarly, type \kbd{pol3 = x\pow 4 - 4*x\pow 2 + 16} and
1124: \kbd{fn = factornf(pol3,t\pow 2 + 1)}, and we get the factorization of the
1125: polynomial \kbd{pol3} over the number field defined by \kbd{t\pow 2 + 1},
1126: i.e.~over $\Q(i)$. To see the result even better, type \kbd{lift(fn)},
1127: remembering that \kbd{t} stands for the generator of the number field
1128: (here equal to $i=\sqrt{-1}$).
1129:
1130: Note that it is possible, although ill advised, to use the same variable
1131: for the polynomial and the number field. You may for example type
1132: \kbd{fn2 = factornf(pol3, x\pow 2 + 1)}, and the result is correct. However,
1133: the PARI object thus created may give unreasonable results. For example,
1134: if you type \kbd{lift(fn2)} in the example above, you will get a strange
1135: object, with a symbol such as \kbd{x*x} typed. This is because PARI knows
1136: that the dummy variable \kbd{x} is not the same as the explicit variable
1137: \kbd{x}, but since it must print it when you lift, it has to do something.
1138: \smallskip
1139: %
1140: To summarize, in addition to being able to factor integers, you can
1141: factor polynomials over $\C$ and $\R$ (this is the function \kbd{polroots}),
1142: over $\F_p$ (the function \kbd{factormod}, over $\F_{p^k}$ (the function
1143: \kbd{factorff}), over $\Q_p$ (the function \kbd{factorpadic}), over $\Q$ or
1144: $\Z$ (the function \kbd{factor}), and over number fields (the functions
1145: \kbd{factornf} and \kbd{nffactor}). Note however that \kbd{factor} itself
1146: will try to guess intelligently over which ring you want to factor: set
1147: \kbd{pol = x\pow2+1}
1148: and try to \kbd{factor} successively \kbd{pol}, \kbd{pol * 1.},
1149: \kbd{pol * (1+0.*I)}, \kbd{pol * Mod(1,2)},
1150: \kbd{pol * Mod(1,Mod(1,3)*(t\pow2+1))}.
1151:
1152: In the present version \vers{}, it is {\it not} possible to factor over
1153: other rings than the ones mentioned above, for example GP cannot factor
1154: multivariate polynomials. \smallskip
1155:
1156: Other functions related to factoring are \kbd{padicappr}, \kbd{polrootsmod},
1157: \kbd{polrootspadic}, \kbd{polsturm}. Play with them a little.
1158:
1159: Now let's type \kbd{polsym(pol3,20)}, where \kbd{pol3} is the same polynomial
1160: as above. This gives the sum of the $k$-th powers of the roots of \kbd{pol3}
1161: up to $k=20$, of course computed using Newton's formula and not using
1162: \kbd{polroots}. You notice that every odd sum is zero (this is trivial since the
1163: polynomial is even), but also that the signs follow a regular pattern and
1164: that the (non-zero) absolute values are powers of 2. This is true: prove it,
1165: and more precisely find an explicit formula for the $k$-th symmetric power
1166: not involving (non-rational) algebraic numbers. \smallskip
1167:
1168: You can also try the function \kbd{polkaramul} (for Karatsuba
1169: multiplication). Put the timer on (\kbd{\#}; well, if you've been very
1170: dedicated and followed this tutorial all the way from the beginning, the
1171: timer should already be on), and then set \kbd{pol4 = (x + 1)\pow 200;} (do
1172: not forget the semicolon, you don't want to see the result, just the timing).
1173: Similarly for the next few commands. Then type \kbd{pol5 = pol4 * pol4;}, and
1174: look at the timing. Now type instead \kbd{pol6 = polkaramul(pol4,pol4,5);}.
1175: This should be about twice as fast (it is platform dependent so the gain may
1176: vary). Try with other values than 5. Also type \kbd{ pol5 - pol6} just to be
1177: sure that the result is correct. So if you have polynomials which are that
1178: large, it may be worthwhile to use this function.
1179:
1180: By the way, try \kbd{pol4\pow 2;}. Surprised? Well, if you write down the
1181: multiplication rules for polynomials, you'll notice that it's easier to
1182: compute a square than a general product (you can save about half of the
1183: computations). This is of course taken into account in the powering operation,
1184: but not when doing a multiplication, since GP does not check whether the two
1185: operands are equal.\medskip
1186:
1187: Now let's play a little with power series. We have already done so a little
1188: at the beginning. Type
1189:
1190: \centerline{\tt 8*x + prod(n=1,39, if(n\%4, 1 - x\pow n, 1),
1191: 1 + O(x\pow 40))\pow 8}
1192:
1193: You obtain a power series which has apparently only even powers of \kbd{x}
1194: appearing. This is surprising, but can be proved using the theory of modular
1195: forms. Note that we have initialized the product to \kbd{1 + O(x\pow 40)} and
1196: not simply to 1 otherwise the whole computation would have been done with
1197: polynomials, and this would first have been slightly slower and also totally
1198: useless since the coefficients of \kbd{x\pow 40} and above are irrelevant
1199: anyhow if we stop the product at \kbd{n=39}.
1200:
1201: While we are on the subject of modular forms (which, together with taylor
1202: series expansions of common functions, are another great source of power
1203: series), type \kbd{\b{ps} 122} (which is a shortcut for
1204: \kbd{default(seriesprecision, 122)}), then \kbd{d = x * eta(x)\pow 24}. This
1205: gives the first 122 terms of the (Fourier) series expansion of the modular
1206: discriminant function $\Delta$ of Ramanujan, its coefficients giving by
1207: definition the Ramanujan $\tau$ function which has a number of marvelous
1208: properties (look at any book on modular forms for explanations). We would like
1209: to see its properties modulo 2. Type \kbd{d\%2}. Hmm, apparently PARI
1210: doesn't like to
1211: reduce coefficients of power series (or polynomials for that matter) directly.
1212: Can we do it without writing a little program? No problem. Type instead
1213: \kbd{lift(Mod(1,2) * d)} and now this works like a charm.
1214:
1215: The pattern in the result is clear. Of course, it now remains to prove it
1216: (again modular forms, see Antwerp III or your resident modular forms guru).
1217: Similarly, type \kbd{centerlift(Mod(1,3) * d)}. This time the pattern is
1218: less clear, but nonetheless there is one. Refer to Anwerp III again.
1219:
1220: \section{Working with Elliptic Curves}
1221:
1222: Now we are getting to more complicated objects. Just as with number fields
1223: which we will meet later on, the first thing to do is to initialize them.
1224: That's because a lot of data will be needed repeatedly, and it's much more
1225: convenient to have it ready once and for all. Here, this is done with the
1226: function \kbd{ellinit} (try to guess what we'll use for number fields\dots).
1227:
1228: So type \kbd{e = ellinit([6,-3,9,-16,-14])}. This computes a number of things
1229: about the elliptic curve defined by the affine equation
1230: %
1231: $$ y^2+6xy+9y = x^3-3x^2-16x-14\enspace. $$
1232: %
1233: It's not that clear what all these funny numbers mean, except that we
1234: recognize the first few of them as the coefficients we just input. To
1235: retrieve meaningful information from such complicated objects (and number
1236: fields will be much worse), you are advised to use the so-called {\it member
1237: functions}. Type \kbd{?.} to get a complete list. Whenever \kbd{ell} appears
1238: in the right hand side, we can apply the corresponding function to an object
1239: output by \kbd{ellinit} (I'm sure you know how the other \kbd{init} functions
1240: will be called now, don't you? Oh, by the way, there is no \kbd{clgpinit}
1241: function).
1242:
1243: Let's try it. We see that the discriminant \kbd{e.disc} is equal to 37,
1244: hence the conductor of the curve is 37. Of course in general it is not so
1245: trivial. In fact, the equation of the curve is clearly not minimal, so type
1246: \kbd{r = ellglobalred(e)}. The first component \kbd{r[1]} tells us that the
1247: conductor is 37 as we already knew. The second component is a 4-component
1248: vector which will allow us to get the minimal equation: simply type
1249: \kbd{e = ellchangecurve(e, r[2])} and the new \kbd{e} is now our minimal
1250: equation with corresponding data. You can for the moment ignore the third
1251: component \kbd{r[3]} (for the impatient reader, this is the product of the
1252: local Tamagawa numbers, $c_p$).
1253:
1254: The new \kbd{e} tells us that the minimal equation is $y^2+y = x^3-x$.
1255: Let us now play a little with points on \kbd{e}. Type \kbd{q = [0,0]}, which is
1256: clearly on the curve (type \kbd{ellisoncurve(e, q)} to check). Well, \kbd{q}
1257: may be a torsion point. Type \kbd{ellheight(e, q)}, which computes the
1258: canonical Neron-Tate height of \kbd{q}. This is non-zero, hence \kbd{q} is
1259: not torsion. To see this even better, type
1260:
1261: \kbd{for(k=1,20, print(ellpow(e, q,k)))}
1262:
1263: \noindent and we see the characteristic parabolic explosion of the size of
1264: the points. As a further check, type
1265: \kbd{ellheight(e, ellpow(e, q,20)) / ellheight(e, q)}. We indeed find
1266: $400=20^2$ as it should be. You can also type \kbd{ellorder(e, q)} which
1267: returns 0, telling you that \kbd{q} is non-torsion.
1268:
1269: Notice how all those \kbd{ell}--prefixed functions take our elliptic curve as
1270: a first argument? This will be true with number fields as well: whatever
1271: object was initialized by an $ob$--\kbd{init} function will have to be used as
1272: a first argument of all the $ob$--prefixed functions. Conversely, you won't be
1273: able to use any such high-level function before you correctly initialize the
1274: relevant object. \smallskip
1275:
1276: Ok, let's try another curve. Type \kbd{e = ellinit([0,-1,1,0,0])}. This
1277: corresponds to the equation $y^2+y = x^3-x^2$. Again from the discriminant
1278: we see that the conductor is equal to 11, and if you type \kbd{ellglobalred(e)}
1279: you will see that the equation for \kbd{e} is minimal. Type \kbd{q = [0,0]}
1280: which is clearly a point on \kbd{e}, and \kbd{ellheight(e, q)}. This time we
1281: obtain a value which is very close to zero, hence \kbd{q} must be a torsion
1282: point. Indeed, typing \kbd{for(k=1,5, print(ellpow(e, q,k)))} we see that
1283: \kbd{q} is a point of order 5 (note that the point at infinity is represented
1284: as \kbd{[0]}). More simply, you can type \kbd{ellorder(e, q)}.\smallskip
1285:
1286: Let's try still another curve. Type \kbd{e = ellinit([0,0,1,-7,6])} to get
1287: the curve $y^2+y=x^3-7x+6$. Typing \kbd{ellglobalred(e)} shows that this is a
1288: minimal equation and that the conductor, equal to the discriminant, is 5077.
1289: There are some trivial integral points on this curve, but let's try to be
1290: more systematic.
1291:
1292: First, let's study the torsion points. Typing \kbd{elltors(e)} shows that the
1293: torsion subgroup is trivial, so we won't have to worry about torsion points.
1294: Next, the member \kbd{e.roots} gives us the 3 roots of the minimal
1295: equation over $\C$, i.e.~$Y^2=X^3-7X+25/4$ (set $(X,Y)=(x,y+1/2)$) so if
1296: $(x,y)$ is a real point on the curve, $x$ must be at least equal to the
1297: smallest root, i.e.~$x\ge-3$. Finally, if $(x,y)$ is on the curve, its
1298: opposite is clearly $(x,-y-1)$. So we are going to use the \kbd{ellordinate}
1299: function and type (for instance in \kbd{points.gp} which you can read in with
1300: \kbd{\b{r} points} as we saw before)
1301:
1302: \bprog\obr
1303: \q v=[];
1304: \q for (x = -3, 1000,
1305: \q\q s = ellordinate(e,x);
1306: \q\q if (length(s), \q\q{\rm \b{}\b{} we could use \kbd{if (!s,} instead}
1307: \q\q\q v = concat(v, [[x,s[1]]])
1308: \q\q )
1309: \q ); v
1310: \cbr\eprog
1311:
1312: \noindent By the way, this is how you insert a comment in a script:
1313: everything following a double backslash (up to the first newline character)
1314: is ignored. If you want comments which span many lines, you can brace them
1315: between \kbd{/* ... */} pairs. Everything in between will be ignored as well.
1316: For instance as a header for the file \kbd{points.gp} you could insert the
1317: following:
1318: %
1319: \bprog
1320: /*\q Finds rational points on the elliptic curve e, using the naivest
1321: \q\q algorithm I could think of right now (TO BE IMPROVED).
1322: \q\q e should have rational coefficients.
1323: \q\q TODO: Make that into a usable function.
1324: \ */
1325: \eprog
1326:
1327: (I hope you did not waste your time copying this nonsense, did you?)
1328:
1329: We thus get a large number (18) of integral points. Together with their
1330: opposites and the point at infinity, this makes a total of 37 integral
1331: points, which is large for a curve having such a small conductor. So we
1332: suspect (if we don't know already, since this curve is quite famous!) that
1333: the rank of this curve must be high. Let's try and put some order into this
1334: (note that we work only with the integral points, but in general rational
1335: points should also be considered).
1336:
1337: Type \kbd{hv = ellheight(e, v)}. This gives the vector of canonical heights.
1338: Let us order the points according to their height. For this, type
1339:
1340: \kbd{iv = vecindexsort(hv)}, then \kbd{hv = vecextract(hv,iv)} and
1341: \kbd{v = vecextract(v,iv)}.
1342:
1343: It seems reasonable to take the numbers with smallest height as generators of
1344: the Mordell-Weil group. Let's try the first 4: type
1345:
1346: \kbd{m = ellheightmatrix(e, vecextract(v,[1,2,3,4])); matdet(m)}
1347:
1348: Since the curve has no torsion, the determinant being close to zero implies
1349: that the first four points are dependent. To find the dependency, it is
1350: enough to find the kernel of the matrix \kbd{m}. So type \kbd{matker(m)}:
1351: we indeed get a non-trivial kernel, and the coefficients are (close to)
1352: integers as they should. Typing \kbd{elladd(e, v[1],v[3])} does indeed show
1353: that it is equal to \kbd{v[4]}.
1354:
1355: Taking any other four points, we would in fact always find a dependency.
1356: Let's find them all. Type \kbd{vp = [v[1],v[2],v[3]]\til;
1357: m = ellheightmatrix(e,vp);
1358: matdet(m)}. This is now clearly non-zero so the first 3 points
1359: are linearly independent, showing that the rank of the curve is at least
1360: equal to 3 (it is in fact equal to 3, and \kbd{e} is the curve of smallest
1361: conductor having rank 3). We would like to see whether the other points are
1362: dependent. For this, we use the function \kbd{ellbil}. Indeed, if \kbd{Q} is
1363: some point which is dependent on \kbd{v[1],v[2]} and \kbd{v[3]}, then
1364: \kbd{matsolve(m, ellbil(e, vp,Q))} will by definition give the coefficients
1365: of the dependence relation. If these coefficients are close to integers, then
1366: there is a dependency, otherwise not. This is much safer than using the
1367: \kbd{matker} function. Thus, type
1368:
1369: \centerline{\tt w = vector(18,k, matsolve(m, ellbil(e, vp,v[k])))}
1370:
1371: We ``see'' that the coefficients are all very close to integers, and we can
1372: prove it by typing
1373:
1374: \centerline{\tt wr = round(w); sqrt(norml2(w - wr))}
1375:
1376: \noindent which gives an upper bound on the maximum distance to an integer.
1377: Thus \kbd{wr} is the vector expressing all the components of \kbd{v} on its
1378: first 3. We are thus led to strongly believe that the curve has rank exactly
1379: 3, and this can be proved to be the case.\smallskip
1380:
1381: Let's now explore a few more elliptic curve related functions. Let us keep
1382: our curve \kbd{e} of rank 3, and type
1383: \bprog
1384: v1 = [1,0]; v2 = [2,0];
1385: z1 = ellpointtoz(e, v1)
1386: z2 = ellpointtoz(e, v2)
1387: \eprog
1388:
1389: We thus get the complex parametrization of the curve. To add the points
1390: \kbd{v1} and \kbd{v2}, we should of course type \kbd{elladd(e, v1,v2)},
1391: but we can also type \kbd{ellztopoint(e, z1 + z2)} which of course has the
1392: disadvantage of giving complex numbers, but illustrates how the group law on
1393: \kbd{e} is obtained from the addition law on $\C$.
1394:
1395: Type \kbd{f = x * Ser(ellan(e, 30))}. This gives a power series which
1396: is the Fourier expansion of a modular form of weight 2 for $\Gamma_0(5077)$
1397: (this has been proved directly, but also follows from Wiles' result since
1398: \kbd{e} is semi-stable). In fact, to find the modular parametrization of
1399: the curve, type \kbd{modul = elltaniyama(e)}, then
1400: \kbd{u=modul[1]; v=modul[2];}. Type
1401:
1402: \centerline{\tt (v\pow 2 + v) - (u\pow 3 - 7*u + 6)}
1403:
1404: \noindent to see that this indeed parametrizes the curve.
1405:
1406: Now type \kbd{x * u' / (2*v + 1)}, and we see that this is equal to the
1407: modular form \kbd{f} found above (the quote \kbd{'} tells GP to take the
1408: derivative of the expression with respect to its main variable). The
1409: functions \kbd{u} and \kbd{v}, considered on the upper half plane (with
1410: $x=e^{2i\pi\tau}$), are in fact modular {\it functions} for $\Gamma_0(5077)$.
1411: \smallskip
1412:
1413: Finally, let us come back to the curve defined by typing
1414: \kbd{e = ellinit([0,-1,1,0,0])} where we had seen that the point
1415: \kbd{q = [0,0]} was of order 5. We know that the conductor of this curve is
1416: equal to 11 (type \kbd{ellglobalred(e)}). We want the sign of the functional
1417: equation. Type
1418:
1419: \centerline{\kbd{elllseries(e, 1,-11,1.1)}, \quad then \quad
1420: \kbd{elllseries(e, 1,-11,1)}.}
1421:
1422: Since the values are clearly different, the sign cannot be $-$. In fact
1423: there is an algebraic algorithm which would allow to compute this sign, but
1424: it has not yet been completely written, although in case of conductors prime
1425: to 6 it is very simple.
1426:
1427: Now type \kbd{ls = elllseries(e, 1,11,1)}, and just as a check
1428: \kbd{elllseries(e, 1,11,1.1)}. The values agree (approximately) as they should,
1429: and give the value of the L-function of \kbd{e} at 1. Now according to
1430: the Birch and Swinnerton-Dyer conjecture (which is proved for this curve),
1431: \kbd{ls} is given by the following formula (in this case):
1432: %
1433: \def\sha{\hbox{III}}
1434: $$L(E,1)=\dfrac{\Omega\cdot c\cdot|\sha|}{|E_{\text{tors}}|^2}\enspace,$$
1435: %
1436: where $\Omega$ is the real period of $E$, $c$ is the global Tamagawa number,
1437: product of the local $c_p$ for primes $p$ dividing the conductor, $|\sha|$ is
1438: the order of the Tate-Shafarevich group, and $E_{\text{tors}}$ is the
1439: torsion group of $E$.
1440:
1441: Now we know many of these quantities: $\Omega$ is equal to \kbd{e.omega[1]}
1442: (if there had been 3 real roots instead of 1 in \kbd{e.roots}, $\Omega$ would
1443: be equal to \kbd{2 * e.omega[1]}). The Tamagawa number $c$ is given as the
1444: last component of \kbd{ellglobalred(e)}, and here is equal to 1. We already
1445: know that the torsion subgroup of $E$ contains a point of order 5, and typing
1446: \kbd{torsell(e)} shows that it is of order exactly 5. Hence type
1447: \kbd{ls * 25/e[15]}. This shows that $|\sha|$ must be equal to 1.
1448:
1449: \section{Working in Quadratic Number Fields}
1450:
1451: The simplest of all number fields outside $\Q$ are quadratic fields. Such
1452: fields are characterized by their discriminant, and even better, any non-square
1453: integer $D$ congruent to 0 or 1 modulo 4 is the discriminant of a specific
1454: order in a quadratic field. We can check whether this order is maximal by
1455: using the command \kbd{isfundamental(D)}. Elements of a quadratic field are
1456: of the form $a+b\omega$, where $\omega$ is chosen as $\sqrt{D}/2$ if $D$ is
1457: even and $(1+\sqrt{D})/2$ if $D$ is odd, and are represented in PARI by
1458: quadratic numbers. To initialize working in a quadratic order, one should
1459: start by the command \kbd{w = quadgen($D$)}.
1460:
1461: This sets \kbd{w} equal to $\omega$ as above, and is printed \kbd{w}. Note
1462: however that if several different quadratic orders are used, a printed \kbd{w}
1463: may have several different meanings. For example if you type
1464:
1465: \centerline{\tt w1 = quadgen(-23); w2 = quadgen(-15);}
1466:
1467: Then ask for the value of \kbd{w1} and \kbd{w2}, both will be printed as
1468: \kbd{w}, but of course they are not equal. Hence beware when dealing with
1469: several quadratic orders at once. \smallskip
1470: %
1471: In addition to elements of a quadratic order, we also want to be able to
1472: handle ideals of such orders. In the quadratic case, it is equivalent to
1473: handling binary quadratic forms, and this has been chosen in PARI. For
1474: negative discriminants, quadratic forms are triples $(a,b,c)$ representing
1475: the form $ax^2+bxy+cy^2$. Such a form will be printed as (and can be created
1476: by)
1477:
1478: \centerline{\tt Qfb($a$,$b$,$c$)}
1479:
1480: Such forms can be multiplied, divided, powered as many PARI objects using
1481: the usual operations, and they can also be reduced using the function
1482: \kbd{qfbred} (it is not the purpose of this tutorial to explain what all
1483: these things mean). In addition, Shanks's NUCOMP algorithm has been
1484: implemented (functions \kbd{qfbnucomp} and \kbd{qfbnupow}), and this is
1485: usually a little faster.
1486:
1487: Finally, you have at your disposal the functions \kbd{qfbclassno} which
1488: ({\it usually\/}) gives the class number, the function \kbd{qfbhclassno}
1489: which gives the Hurwitz class number, and the much more sophisticated
1490: \kbd{quadclassunit} function which gives the class number and class group
1491: structure.
1492:
1493: Let us see examples of all this at work.
1494:
1495: Type \kbd{qfbclassno(-10007)}. GP tells us that the result is 77. However,
1496: you may have noticed in the explanation above that the result is only
1497: {\it usually\/} correct. This is because the implementers of the algorithm
1498: have been lazy and have not put the complete Shanks algorithm into PARI,
1499: causing it to fail in certain very rare cases. In practice, it is almost
1500: always correct, and the much more powerful \kbd{quadclassunit} program, which
1501: {\it is} complete (at least for fundamental discriminants) can give
1502: confirmation (but now, under the Generalized Riemann Hypothesis!!!).
1503:
1504: So we may be a little suspicious of this class number. Let us check it.
1505: First, we need to find a quadratic form of discriminant $-10007$. Since this
1506: discriminant is congruent to 1 modulo 8, we know that there is an ideal of
1507: norm equal to 2, i.e.~a binary quadratic form $(a,b,c)$ with $a=2$. To
1508: compute it we type \kbd{f = qfbprimeform(-10007, 2)}. OK, now we have a form.
1509: If the class number is correct, the very least is that this form raised to
1510: the power 77 should equal the identity. Let's check this. Type \kbd{f\pow 77}.
1511: We get a form starting with 1, i.e.~the identity, so this test is OK. Raising
1512: \kbd{f} to the powers 11 and 7 does not give the identity, thus we now know
1513: that the order of \kbd{f} is exactly 77, hence the class number is a multiple
1514: of 77. But how can we be sure that it is exactly 77 and not a proper multiple?
1515: Well, type
1516: %
1517: \bprog
1518: sqrt(10007)/Pi * prodeuler(p=2,500, 1./(1 - kronecker(-10007,p)/p))
1519: \eprog
1520: %
1521: This is nothing else than an approximation to the Dirichlet class number
1522: formula. The function \kbd{kronecker} is the Kronecker symbol, in this case
1523: simply the Legendre symbol. Note also that we have written \kbd{1./(1 - \dots)}
1524: with a dot after the first 1. Otherwise, PARI may want to compute the whole
1525: thing as a rational number, which would be terribly long and useless. In fact
1526: PARI does no such thing in this particular case (\kbd{prodeuler} is always
1527: computed as a real number), but you never know. Better safe than sorry!
1528:
1529: We find 77.77, pretty close to 77, so things seem in order. Explicit bounds
1530: on the prime limit to be used in the Euler product can be given which make
1531: the above reasoning rigorous.
1532:
1533: Let us try the same thing with $D=-3299$. \kbd{qfbclassno} and the Euler
1534: product convince us that the class number must be 27. However, we get stuck
1535: when we try to prove this in the simple-minded way above. Indeed, we type
1536: \kbd{f = qfbprimeform(-3299, 3)} (here 2 is not the norm of a prime ideal but
1537: 3 is), and we see that \kbd{f} raised to the power 9 is equal to the identity.
1538: This is the case for any other quadratic form we choose. So we suspect that the
1539: class group is not cyclic. Indeed, if we list all 9 distinct powers of \kbd{f},
1540: we see that \kbd{qfbprimeform(-3299, 5)} is not on the list (although its cube
1541: is as it must). This implies that the class group is probably equal to a
1542: product of a cyclic group of order 9 by a cyclic group of order 3. The Euler
1543: product plus explicit bounds prove this.
1544:
1545: Another way to check it is to use the \kbd{quadclassunit} function by typing
1546: for example
1547:
1548: \centerline{\tt quadclassunit(-3299)}
1549:
1550: Note that this function cheats a little and could still give a wrong answer,
1551: even assuming GRH (we could get a subgroup and not the whole class group).
1552: If we want to use proven bounds under GRH, we have to type
1553:
1554: \centerline{\tt quadclassunit(-3299,,[1,6])}
1555:
1556: The double comma \kbd{,,} is not a typo, it means we omit an optional second
1557: argument (we would use it to compute the narrow class group, which would be
1558: the same here of course). As we want to use the optional {\it third}
1559: argument, we have to indicate to GP we skipped this one.
1560:
1561: Now, if we believe in GRH, the class group is as we thought (see Chapter 3
1562: for a complete description of this function).
1563:
1564: Note that using the even more general function \kbd{bnfinit} (which handles
1565: general number fields and gives much more complicated results), we could
1566: {\it certify\/} this result (remove the GRH assumption). Let's do it, type
1567: \bprog
1568: bnf = bnfinit(x\pow 2 + 3299); bnfcertify(bnf)
1569: \eprog
1570:
1571: A non-zero result (here 1) means that everything is ok. Good, but what did
1572: we certify after all? Let's have a look at this \kbd{bnf} (just type it!).
1573: Enlightening, isn't it? Recall that the \kbd{init} functions (we've already
1574: seen \kbd{ellinit}) store all kind of technical information which you
1575: certainly don't care about, but which will be put to good use by some higher
1576: level functions. That's why \kbd{bnfcertify} could not be used on the output
1577: of \kbd{quadclassunit}: it needs much more data.
1578:
1579: To extract sensible information from such complicated objects, you must use
1580: one of the many {\it member functions} (remember: \kbd{?.} to get a complete
1581: list). In this case \kbd{bnf.clgp} which extracts the class group structure.
1582: This is much better. Type \kbd{\%.no} to check that this leading 27 is indeed
1583: what we think it is and not some stupid technical parameter. Note that
1584: \kbd{bnf.clgp.no} would work just as well, or even \kbd{bnf.no}!
1585:
1586: As a last check, we can request a relative equation for the Hilbert class
1587: field of $\Q(\sqrt{-3299})$: type \kbd{quadhilbert(-3299)}. It is indeed of
1588: degree 27 so everything fits together.
1589:
1590: \medskip
1591: %
1592: Working in real quadratic fields instead of complex ones, i.e.~with $D>0$, is
1593: not very different.
1594:
1595: The same \kbd{quadgen} function is used to create elements. Ideals are again
1596: represented by binary quadratic forms $(a,b,c)$, this time indefinite. However,
1597: the Archimedean valuations of the number field start to come into play (as
1598: is clear if one considers ideles instead of ideals), hence in fact quadratic
1599: forms with positive discriminant will be represented as a quadruplet
1600: $(a,b,c,d)$ where the quadratic form itself is $ax^2+bxy+cy^2$ with $a$,
1601: $b$ and $c$ integral, and $d$ is the Archimedean component, a real number.
1602: For people familiar with the notion, $d$ represents a ``distance'' as defined
1603: by Shanks and Lenstra.
1604:
1605: To create such forms, one uses the same function as for definite ones, but
1606: you can add a fourth (optional) argument to initialize the distance:
1607:
1608: \centerline{\tt Qfb($a$, $b$, $c$, $d$)}
1609:
1610: If the discriminant of $(a,b,c)$ is negative, $d$ is silently
1611: discarded. If you omit it, this component is set to \kbd{0.} (i.e.~a real zero
1612: to the current precision).
1613:
1614: Again these forms can be multiplied, divided, powered, and they can be
1615: reduced using the function \kbd{qfbred}. This function is in fact a
1616: succession of elementary reduction steps corresponding essentially to a
1617: continued fraction expansion, and a single one of these steps can be achieved
1618: by adding an (optional) flag to the arguments of using this function. Since
1619: handling the fourth component $d$ usually involves computing logarithms, the
1620: same flag may be used to ignore the fourth component. Finally, it is
1621: sometimes useful to operate on forms of positive discriminant without
1622: performing any reduction (this is useless in the negative case), the
1623: functions \kbd{qfbcompraw} and \kbd{qfbpowraw} do exactly that.
1624:
1625: Again, the function \kbd{qfbprimeform} gives a prime form, but the form which
1626: is given corresponds to an ideal of prime norm which is usually not reduced.
1627: If desired, it can be reduced using \kbd{qfbred}.
1628:
1629: Finally, you still have at your disposal the function \kbd{qfbclassno} which
1630: gives the class number (this time {\it guaranteed\/} correct),
1631: \kbd{quadregulator} which gives the regulator, and the much more sophisticated
1632: \kbd{quadclassunit} giving the class group's structure and its generators,
1633: as well as the regulator. The \kbd{qfbclassno} and \kbd{quadregulator}
1634: functions use an algorithm which is $O(\sqrt D)$, hence become very slow for
1635: discriminants of more than 10 digits. \kbd{quadclassunit} can be used on a
1636: much larger range.
1637:
1638: Let us see examples of all this at work and learn some little known number
1639: theory at the same time. First of all, type
1640: \kbd{d = 3 * 3299; qfbclassno(d)}. We see that the class number is 3 (we know
1641: in advance that it must be divisible by 3 from the \kbd{d = -3299} case above
1642: and Scholz's theorem). Let us create a form by typing
1643: \kbd{f = qfbred(qfbprimeform(d,2), 2)} (the last 2 tells \kbd{qfbred} to
1644: ignore the archimedean component). This gives us a prime ideal of norm
1645: equal to 2. Is this ideal principal? Well, one way to check this, which is
1646: not the most efficient but will suffice for now, is to look at the complete
1647: cycle of reduced forms equivalent to \kbd{f}. Type
1648: \bprog
1649: g = f; for(i=1,20, g = qfbred(g, 3); print(g))
1650: \eprog\noindent
1651: (this time the 3 means to do a single reduction step, still not using
1652: Shanks's distance). We see that we come back to the form \kbd{f} without
1653: having the principal form (starting with $\pm1$) in the cycle, so the ideal
1654: corresponding to \kbd{f} is not principal.
1655:
1656: Since the class number is equal to 3, we know however that \kbd{f\pow 3} will
1657: be a principal ideal $\alpha\Z_K$. How do we find $\alpha$? For this, type
1658: \kbd{f3 = qfbpowraw(f, 3)}. This computes the cube of \kbd{f}, without
1659: reducing it. Hence it corresponds to an ideal of norm equal to $8=2^3$, so we
1660: already know that the norm of $\alpha$ is equal to $\pm8$. We need more
1661: information, and this will be given by the fourth component of the form.
1662: Reduce your form until you reach the unit form (you will have to type
1663: \kbd{qfbred(\%,~1)} exactly 6 times).
1664:
1665: Extract the Archimedean component by typing \kbd{c = component(\%, 4)}. By
1666: definition of this distance, we know that
1667: $${\alpha\over{\sigma(\alpha)}}=\pm e^{2c},$$
1668: where $\sigma$ denotes real conjugation in our quadratic field. Thus, if we
1669: type
1670:
1671: \centerline{\tt a = sqrt(8 * exp(2*c))}
1672:
1673: \noindent and then \kbd{sa = 8 / a}, we know that up to sign, \kbd{a} and
1674: \kbd{sa} are numerical approximations of $\alpha$ and $\sigma(\alpha)$. Of
1675: course, $\alpha$ can always be chosen to be positive, and a quick numerical
1676: check shows that the difference of \kbd{a} and \kbd{sa} is close to an
1677: integer, and not the sum, so that in fact the norm of $\alpha$ is equal to
1678: $-8$ and the numerical approximation to $\sigma(\alpha)$ is \kbd{$-$sa}. Thus
1679: we type
1680:
1681: \centerline{\tt p = x\pow 2 - round(a-sa)*x - 8}
1682:
1683: \noindent and this is the characteristic polynomial of $\alpha$. We can check
1684: that the discriminant of this polynomial is a square multiple of \kbd{d}, so
1685: $\alpha$ is indeed in our field. More precisely, solving for $\alpha$ and
1686: using the numerical approximation that we have to resolve the sign ambiguity in
1687: the square root, we get explicitly $\alpha=(15221+153\sqrt d)/2$. Note that
1688: this can also be done automatically using the functions \kbd{polred} and
1689: \kbd{modreverse}, as we will see later in the general number field case, or by
1690: solving a system of 2 linear equations in 2 variables.
1691:
1692: \noindent{\bf Exercise:} now that we have $\alpha$ explicitly, check that it
1693: is indeed a generator of the ideal corresponding to the form \kbd{f3}.
1694:
1695: \medskip Let us now play a little with cycles. Set \kbd{D = 10\pow 7 + 1},
1696: then type
1697:
1698: \centerline{\tt quadclassunit(D,,[1,6])}
1699:
1700: We get as a result a 5-component vector, which tells us that (under GRH) the
1701: class number is equal to 1, and the regulator is approximately
1702: equal to $2641.5$. If you want to certify this, use \kbd{qfbclassno} and
1703: \kbd{quadregulator}, {\it not} \kbd{bnfinit} and \kbd{bnfcertify}, which will
1704: take an absurdly long time (well, about 5 minutes if you are careful and set
1705: the initial precision correctly). Indeed \kbd{bnfcertify} needs the fundamental
1706: unit which is so large that \kbd{bnfinit} will have a (relatively) hard time
1707: computing it: you will need about $R/\log(10)\approx 1147$ digits of precision!
1708: On the other hand, you can try \kbd{quadunit(D)}. Impressive, isn't it? (you
1709: can check that its logarithm is indeed equal to the regulator).
1710:
1711: Now just as an example, let's assume that we want the regulator to 500
1712: decimals, say (without cheating and computing the fundamental unit exactly
1713: first!). I claim that by simply knowing the crude approximation above, this
1714: can be computed with no effort.
1715:
1716: This time, we want to start with the unit form. Since \kbd{D} is odd, we can
1717: type:
1718:
1719: \centerline{\tt u = qfbred(Qfb(1,1,(1 - D)/4), 2)}
1720:
1721: We use the function \kbd{qfbred} with no distance since we want the initial
1722: distance to be equal to~0.
1723:
1724: Now we type \kbd{f = qfbred(u, 1)}. This is the first form encountered along
1725: the principal cycle. For the moment, keep the precision low, for example the
1726: initial default precision. The distance from the identity of \kbd{f} is
1727: around 4.253. Very crudely, since we want a distance of $2641.5$, this should
1728: be encountered approximately at $2641.5/4.253=621$ times the distance of
1729: \kbd{f}. Hence, as a first try, we type \kbd{f\pow 621}. Oops, we overshot,
1730: since the distance is now $3173.02$. Now we can refine our initial estimate and
1731: believe that we should be close to the correct distance if we raise \kbd{f} to
1732: the power $621*2641.5/3173$ which is close to $517$. Now if we compute
1733: \kbd{f\pow 517} we hit the principal form right on the dot. Note that this is
1734: not a lucky accident: we will always land extremely close to the correct target
1735: using this method, and usually at most one reduction correction step is
1736: necessary. Of course, only the distance component can tell us where we are
1737: along the cycle.
1738:
1739: Up to now, we have only worked to low precision. The goal was to obtain this
1740: unknown integer $517$. Note that this number has absolutely no mathematical
1741: significance: indeed the notion of reduction of a form with positive
1742: discriminant is not well defined since there are usually many reduced forms
1743: equivalent to a given form. However, when PARI makes its computations, the
1744: specific order and reductions that it performs are dictated entirely by the
1745: coefficients of the quadratic form itself, and not by the distance component,
1746: hence the precision used has no effect.
1747:
1748: Hence we now start again by setting the precision to (for example) 500,
1749: we retype the definition of \kbd{u} (why is this necessary?), and then
1750: \kbd{f = qfbred(u, 1)} and finally \kbd{f\pow 517}. Of course we know in
1751: advance that we land on the unit form, and the fourth component gives us the
1752: regulator to 500 decimal places with no effort at all.
1753:
1754: In a similar way, we could obtain the so-called {\it compact representation}
1755: of the fundamental unit itself, or $p$-adic regulators. I leave this as
1756: exercises for the interested reader.
1757:
1758: You can try the \kbd{quadhilbert} function on that field but, since the class
1759: number is $1$, the result won't be that exciting. If you try it on our
1760: preceding example ($3*3299$) it should take about five minutes (time for a
1761: coffee break?).
1762:
1763: \section{Working in General Number Fields}
1764:
1765: Note for the present release: this section is a little obsolete since many
1766: new powerful functions are available now. This needs to be rewritten
1767: entirely. \smallskip
1768:
1769: The situation here is of course more difficult. First of all, remembering
1770: what we did with elliptic curves, we need to initialize it (with something
1771: more serious than \kbd{quadgen}). For example assume that we want to work in
1772: the number field $K$ defined by one of the roots of the equation
1773: $x^4+24x^2+585x+1791=0$. This is done by typing
1774: \bprog
1775: T = x\pow 4 + 24*x\pow 2 + 585*x + 1791
1776: nf = nfinit(T)
1777: \eprog
1778:
1779: We get quite a complicated object but, thanks to member functions, we don't
1780: need to know anything about its internal structure (which is dutifully
1781: described in Chapter~3). If you type \kbd{nf.pol}, you will get the
1782: polynomial \kbd{T} which you just input. \kbd{nf.sign} yields the signature
1783: $(r_1,r_2)$ of the field, \kbd{nf.disc} the field discriminant, \kbd{nf.zk}
1784: an integral basis, etc\dots.
1785:
1786: The integral basis is expressed in terms of a generic root \kbd{x} of \kbd{T}
1787: and we notice it's very far from being a power integral basis, which is a
1788: little strange for such a small field. Hum, let's check that: \kbd{poldisc(T)}?
1789: Ooops, small wonder we had such denominators, the index is, well, type
1790: \kbd{sqrt(\% / nf.disc)}. That's $3087$, we don't want to work with such
1791: a badly skewed polynomial!
1792:
1793: So, type \kbd{P = polred(T)}. We see from the third component that the
1794: polynomial $x^4-x^3-21x^2+17x+133$ defines the same field with much smaller
1795: coefficients, so type \kbd{A = P[3]}. The \kbd{polred} function gives a
1796: (usually) simpler polynomial, and also sometimes some information on the
1797: existence of subfields. For example in this case, the second component of
1798: \kbd{polred} tells us that the field defined by $x^2-x+1=0$, i.e.~the field
1799: generated by the cube roots of unity, is a subfield of our number field $K$.
1800: Note this is given as incidental information and that the list of subfields
1801: found in this way is usually far from complete. To get the complete list, you
1802: will have to use the function \kbd{nfsubfields} (we'll do that later on).
1803:
1804: Type \kbd{poldisc(A)}, this is much better, but maybe not optimal yet
1805: (the index is still $7$). Type \kbd{polredabs(A)} (the \kbd{abs} stands for
1806: absolute). Since it seems that we won't get anything better, we'll stick with
1807: \kbd{A} (note however that \kbd{polredabs} finds a smallest generating
1808: polynomial with respect to a bizarre norm which ensures that the index will
1809: be small, but not necessarily minimal). In fact, had you typed
1810: \kbd{nfinit(T, 3)}, \kbd{nfinit} would first have tried to find a good
1811: polynomial defining the same field (i.e.~one with small index) before
1812: proceeding.
1813:
1814: It's not too late, let's redefine our number field: \kbd{NF = nfinit(nf, 3)}.
1815: The output is a two-component vector. The first component is the new
1816: \kbd{nf} (type \kbd{nf = NF[1];}). If you type \kbd{nf.pol}, you notice that GP
1817: indeed replaced your bad polynomial \kbd{T} by a much better one, which
1818: happens to be \kbd{A} (small wonder, \kbd{nfinit} internally called
1819: \kbd{polredabs}!). The second component enables you to switch conveniently to
1820: our new polynomial.
1821:
1822: Namely, call $\theta$ a root of our initial polynomial \kbd{T}, and $\alpha$
1823: a root of the one that \kbd{polred} has found, namely \kbd{A}. These are
1824: algebraic numbers, and as already mentioned are represented as polmods. For
1825: example, in our special case $\theta$ is equal to the polmod
1826:
1827: \centerline{\tt Mod(x, x\pow 4 + 24*x\pow 2 + 585*x + 1791)}
1828:
1829: \noindent while $\alpha$ is equal to the polmod
1830:
1831: \centerline{\tt Mod(x, x\pow 4 - x\pow 3 - 21*x\pow 2 + 17*x + 133)}
1832:
1833: \noindent Here of course we are considering only the algebraic aspect, and
1834: hence ignore completely {\it which} root $\theta$ or $\alpha$ is chosen.
1835:
1836: Now probably you may have a number of elements of your number field which are
1837: expressed as polmods with respect to your old polynomial, i.e.~as
1838: polynomials in $\theta$. Since we are now going to work with $\alpha$
1839: instead, it is necessary to convert these numbers to a representation using
1840: $\alpha$. This is what the second component of \kbd{NF} is for: type
1841: \kbd{NF[2]}, you get
1842:
1843: \centerline{\tt Mod(x\pow 2 + x - 11, x\pow 4 - x\pow 3 - 21*x\pow 2 +
1844: 17*x + 133)}
1845:
1846: \noindent meaning that $\theta = \alpha^2+\alpha-11$, and hence the conversion
1847: from a polynomial in $\theta$ to one in $\alpha$ is easy, using \kbd{subst}
1848: (we could get this polynomial from \kbd{polred} as well, try
1849: \kbd{polred(T, 2)}). If we want to do the reverse, i.e.~go back from a
1850: representation in $\alpha$ to a representation in $\theta$, we use the
1851: function \kbd{modreverse} on this polynomial \kbd{NF[2]}. Try it. The result
1852: has a big denominator (147) essentially because our initial polynomial
1853: \kbd{T} was so bad. By the way to get the 147, you should type
1854: \kbd{denominator(content(NF[2]))}. Trying \kbd{denominator} by itself would not
1855: work since the denominator of a polynomial is defined to be 1 (and its
1856: numerator is itself). The reason for this ``surprising'' behaviour is that we
1857: think of a polynomial as a special case of a rational function. \smallskip
1858:
1859: From now on, we completely forget about \kbd{T}, and use only the polynomial
1860: \kbd{A} defining $\alpha$, and the components of the vector \kbd{nf} which
1861: gives information on our number field $K$. Type
1862:
1863: \centerline{\tt u = Mod(x\pow 3 - 5*x\pow 2 - 8*x + 56, A) / 7}
1864:
1865: This is an element in $K$. There are three essentially equivalent
1866: representations for number field elements: polmod, polynomial, and column
1867: vector giving a decomposition in the integral basis \kbd{nf.zk} ({\it not} on
1868: the power basis $(1,x,x^2,\dots)$). All three are equally valid when the
1869: number field is understood (is given as first argument to the function).
1870: You will be able to use any one of them as long as the function you call
1871: requires an \kbd{nf} argument as well. However, most PARI functions will
1872: return elements as column vectors.
1873:
1874: It's a very important feature of number theoretic functions that, although
1875: they may have a preferred format for input, they will accept a wealth of
1876: other different formats. We already saw this for \kbd{nfinit} which
1877: accepts either a polynomial or an \kbd{nf}. It will be true for ideals,
1878: ideles, congruence subgroups, etc.
1879:
1880: Ok, let's stick with elements for the time being. How does one go from one
1881: representation to the other? Between polynomials and polmods, it's easy:
1882: \kbd{lift} and \kbd{Mod} will do the job. Next, from polmods/polynomials to
1883: column vectors: type \kbd{v = nfalgtobasis(nf, u)}. So $\kbd{u} = \alpha^3-
1884: \alpha^2 - \alpha + 8$, right? Wrong! The coordinates of \kbd{u} are given
1885: with respect to the {\it integral basis}, not the power basis
1886: $(1,\alpha,\alpha^2,\alpha^3)$ (and they don't coincide, type \kbd{nf.zk} if
1887: you forgot what the integral basis looked like). As a polynomial in $\alpha$,
1888: we simply have $\kbd{u} = {1\over7}\alpha^3 -
1889: {5\over7}\alpha^2-{8\over7}\alpha+8$, which is trivially deduced from the
1890: original polmod representation!
1891:
1892: Of course \kbd{v = nfalgtobasis(nf, lift(u))} would work equally well. Indeed
1893: we don't need the polmod information since \kbd{nf} already provides the
1894: defining polynomial. To go back to polmod representation, use
1895: \kbd{nfbasistoalg(nf, v)}. Notice that \kbd{u} is in fact an integer since
1896: \kbd{v} has integer coordinates (try \kbd{denominator(v) == 1}, which is of
1897: course overkill here, but not so in a program).
1898:
1899: Let's try this out. We may for instance compute \kbd{u\pow 3}. Try it. Or, we
1900: can type \kbd{1/u}. Better yet, if we want to know the norm from $K$ to $\Q$
1901: of \kbd{u}, we type \kbd{norm(u)} (what else?). \kbd{trace(u)} works as well.
1902: Notice that none of this would work on polynomials or column vectors since
1903: you don't have the opportunity to supply \kbd{nf}! But we could use
1904: \kbd{nfeltpow(nf,u,3)}, \kbd{nfeltdiv(nf,1,u)} (or \kbd{nfeltpow(nf,u,-1)})
1905: which would work whatever representation was chosen. There is no
1906: \kbd{nfeltnorm} function (\kbd{nfelttrace} does not exist either), but we can
1907: easily write one:
1908:
1909: \bprog nfeltnorm(nf,u)=
1910: \obr
1911: \q local(t);
1912: \q t = type(u);
1913: \q if (t != "t\_POLMOD",
1914: \q\q if (t == "t\_COL",
1915: \q\q\q u = nfbasistoalg(nf, u)
1916: \q\q ,
1917: \q\q\q u = Mod(u, nf.pol)
1918: \q\q )
1919: \q );
1920: \q norm(u)
1921: \cbr\eprog
1922:
1923: Notice that this is certainly not foolproof (try it with complex or quadratic
1924: arguments!), but we could refine it if the need arose. In fact there was no
1925: need for this function, since you can consider (\kbd{u}) as a principal
1926: ideal, and just type \kbd{idealnorm(nf,u)} whatever the chosen representation
1927: for \kbd{u}. We'll talk about ideals later on.
1928:
1929: If we want all the symmetric functions of \kbd{u} and not only the norm, we
1930: type \kbd{charpoly(u)} (we could write \kbd{charpoly(u, y)} to tell GP to
1931: use the variable \kbd{y} for the characteristic polynomial). Note that this
1932: gives the characteristic polynomial of \kbd{u}, and not in general the
1933: minimal polynomial. Exercises: how does one (easily) find the minimal
1934: polynomial from this? Find a simpler expression for \kbd{u}.\smallskip
1935:
1936: OK, now let's work on the field itself. The \kbd{nfinit} command already
1937: gave us some information. The field is totally complex (its signature
1938: \kbd{nf.sign} is $[0,2]$), its discriminant \kbd{nf.disc} is $D=18981$ and
1939: $(1,\alpha, \alpha^2, {1\over7}\alpha^3+{2\over7}\alpha^2+{6\over7}\alpha)$
1940: is an integral basis (\kbd{nf.zk}). The Galois group of its Galois closure
1941: can be obtained by typing \kbd{polgalois(A)}. The answer ($[8,-1,1]$) shows
1942: that it is equal to $D_4$, the dihedral group with 8 elements, i.e.~the group
1943: of symmetries of a square. \smallskip
1944:
1945: This implies that the field is ``partially Galois'', i.e.~that there exists
1946: at least one non-trivial field isomorphism which fixes $K$ (exactly one in
1947: this case). To find out which it is, we use the function \kbd{nfgaloisconj}.
1948: This uses the LLL algorithm to find linear relations. So type
1949: \kbd{nfgaloisconj(nf)}. The result tells us that, apart from the trivial
1950: automorphism, the map $\alpha \mapsto (-\alpha^3+5\alpha^2+\alpha-49)/7$ (in
1951: the third component) is a field automorphism. Indeed, if we type
1952: \kbd{s = Mod(\%[3], A); charpoly(s)}, we obtain the polynomial \kbd{A} once
1953: again. \smallskip
1954:
1955: The fixed field of this automorphism is going to be the only non-trivial
1956: subfield of $K$. I seem to recall that \kbd{polred} told us this was the
1957: third cyclotomic field. Let's check this: type \kbd{nfsubfields(nf)}. Indeed,
1958: there's a quadratic subfield, but it's given by \kbd{Q = x\pow 2 + 22*x + 133
1959: } and I don't recognize it. Now \kbd{polred(Q)} proves that this subfield is
1960: indeed the field generated by a cube root of unity. Let's check that \kbd{s}
1961: is of order 2: \kbd{subst(lift(s), x, s)}. Yup, it is. Let's express it as a
1962: matrix:
1963:
1964: \bprog\obr
1965: \q v = [;]; b = nf.zk;
1966: \q for (i=1, 4,
1967: \q\q v = concat(v, nfalgtobasis(nf, nfgaloisapply(nf, s, b[i])))
1968: \q )
1969: \cbr\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:
2140: \bprog\obr
2141: \q m1 = nf[5][1];
2142: \q m = matrix(4,4,j,k,
2143: \q\q if (j<=2,
2144: \q\q\q m1[j,k]
2145: \q\q ,
2146: \q\q\q conj(m1[j-2, k])
2147: \q\q)
2148: \q );
2149: \q v = exp(v1);
2150: \q au = matsolve(m,v);
2151: \q vu = round(real(au))
2152: \cbr\eprog
2153:
2154: Then \kbd{vu} is the representation of the unit on the integral basis.
2155: The closeness of the approximation of \kbd{au} to \kbd{vu} gives us
2156: confidence that we have made no numerical mistake. To be sure that \kbd{vu}
2157: represents a unit, type \kbd{u = nfbasistoalg(nf,vu)}, then typing
2158: \kbd{norm(u)} we see that it is equal to 1 so \kbd{u} is a unit.
2159:
2160: There is of course no reason for \kbd{u} to be a fundamental unit. Let us see
2161: if it is a square. Type \kbd{f1 = factor(subst(charpoly(u,x), x, x\pow 2))}.
2162: We see that the characteristic polynomial of \kbd{u} where \kbd{x} is
2163: replaced by \kbd{x\pow 2} is a product of 2 polynomials of degree 4, hence
2164: \kbd{u} is a square (Exercise: why?).
2165:
2166: We now want to find the square root of \kbd{u}. We can again use the
2167: \kbd{matsolve} function as above. For this we need to take the square
2168: root of each element of the vector \kbd{v}, and hence there are
2169: sign ambiguities. Let's do it anyway. Type \kbd{v = sqrt(v)}. We see that
2170: \kbd{v[1]} and \kbd{v[3]} are conjugates, as well as \kbd{v[2]} and \kbd{v[4]},
2171: so for the moment the signs seem OK. Now try \kbd{au = matsolve(m,v)}. The
2172: numbers obtained are clearly not integers, hence the last remaining sign
2173: change must be performed. Type \kbd{v[1] = -v[1]; v[3] = -v[3]} (they must
2174: stay conjugate) and then again \kbd{au = matsolve(m,v)}. This time the
2175: components are close to integers, so we are done (after typing
2176: \kbd{vu = round(real(au))} as before).
2177:
2178: Anyway, we find that a square root \kbd{u2} of \kbd{-u} is represented by
2179: the vector \kbd{vu=[-4,1,1,-1]\til} on the integral basis, and this is in
2180: fact a fundamental unit.\medskip
2181:
2182: The function \kbd{polred} gives us another method to find \kbd{u2} as
2183: follows: type \kbd{q = polred(f1[1,1], 2)}. We recognize the polynomial
2184: \kbd{A} as the component \kbd{q[3,2]}. To obtain the square root of our unit
2185: we then simply type \kbd{up2 = modreverse(Mod(q[3,1], f1[1,1]))}
2186: (Exercise: why?). We find that \kbd{up2} is represented by the vector
2187: \kbd{[-3,-1,0,0]\til} on the integral basis, which is not the result that we
2188: have found before nor its opposite. Where is the error? (Please think about
2189: this before reading on. There is a mathematical subtlety hidden here.)
2190:
2191: Have you solved the problem? Good! The problem occurs because as mentioned
2192: before (but you may not have noticed since it is not stressed in standard
2193: textbooks) although the number field $K$ is not Galois over $\Q$, there does
2194: exist a non-trivial automorphism, and we have found it above by using the
2195: function \kbd{nfgaloisconj}. Indeed, if we apply this automorphism to
2196: \kbd{up2} (by typing \kbd{nfgaloisapply(nf,s,up2)} where \kbd{s} is the
2197: non-trivial component of \kbd{nfgaloisconj(nf)} computed above), we find the
2198: opposite of \kbd{u2}, which is OK. \smallskip
2199:
2200: Still another method which avoids all sign ambiguities and automorphism
2201: problems is as follows. Type \kbd{r = f1[1,1] \% (x\pow 2 - u)} to find the
2202: remainder of the characteristic polynomial of \kbd{u2} divided by
2203: \kbd{x\pow 2 - u}. This will be a polynomial of degree 1 in \kbd{x} (with
2204: polmod coefficients) and we know that \kbd{u2}, being a root of both
2205: polynomials, will be the root of \kbd{r}, hence can be obtained by typing
2206: \kbd{u2 = -coeff(r,0) / coeff(r,1)}. Indeed, we immediately find the correct
2207: result with no trial and error.\smallskip
2208:
2209: Still another method to find the square root of \kbd{u} is to use
2210: \kbd{nffactor(nf,y\pow 2 + u)}. Except that this won't work as is since the
2211: main variable of the polynomial to be factored must have {\it higher}
2212: priority than the number field variable. This won't be possible here since
2213: \kbd{nf} was defined using the variable \kbd{x} which has the highest possible
2214: priority. So we need to substitute variables around (using \kbd{subst}). I
2215: leave you to work out the details.\smallskip
2216:
2217: Now ideals can be used in a wide variety of formats. We have already seen
2218: HNF-representations, ideles and prime ideals. We can also use algebraic
2219: numbers. For example type \kbd{al = Mod(x\pow 2 - 9, A)}, then
2220: \kbd{ideleprincipal(nf,al)}. We obtain the idele corresponding to \kbd{al}
2221: (see the manual for the exact description). However it is usually not
2222: necessary to compute this explicitly since this is done automatically inside
2223: PARI. For example, you can type
2224:
2225: \centerline{\tt for(i=1,4, print(i ": " idealval(nf,al,P[i])))}
2226:
2227: We see that the valuation is non-zero (equal to 1) at the prime ideals
2228: \kbd{P[2]} and \kbd{P[3]}. In addition, typing \kbd{norm(al)} shows that
2229: \kbd{al} is of norm $49=7^2$ (\kbd{idealnorm(nf,al)} gives the same result of
2230: course, and would be more generic). Let's check this differently. Type
2231: \kbd{P23 = idealmul(nf,P[2],P[3])} and then \kbd{idealhnf(nf,al)}. We see that
2232: the results are the same, hence the product of the two prime ideals \kbd{P[2]}
2233: and \kbd{P[3]} is equal to the principal ideal generated by \kbd{al}. There is
2234: still something to be done with this example as we shall see below after we
2235: introduce Big Number Fields (which will trivialize the examples we have just
2236: seen).
2237:
2238: Essentially all functions that you would want on ideals are available.
2239: We mention here the complete list, referring to Chapter 3 for detailed
2240: explanations:
2241:
2242: \kbd{idealadd}, \kbd{idealaddtoone}, \kbd{idealappr}, \kbd{idealchinese},
2243: \kbd{idealcoprime}, \kbd{idealdiv}, \kbd{idealfactor}, \kbd{idealhnf},
2244: \kbd{idealintersect}, \kbd{idealinv}, \kbd{ideallist}, \kbd{ideallog},
2245: \kbd{idealmin}, \kbd{idealmul}, \kbd{idealnorm}, \kbd{idealpow},
2246: \kbd{idealprimedec}, \kbd{idealprincipal}, \kbd{idealred},
2247: \kbd{idealstar}, \kbd{idealtwoelt}, \kbd{idealval}, \kbd{ideleprincipal},
2248: \kbd{nfisideal}.
2249:
2250: We suggest you play with these functions to get a feel for the algebraic
2251: number theory package. Remember simply that when a matrix (usually in Hermite
2252: normal form) is output, it is always a $\Z$-basis of the result expressed on
2253: the {\it integral basis} \kbd{nf.zk} of the number field, which is usually
2254: {\it not} a power basis. \medskip
2255:
2256: Apart from the above functions you have at your disposal the very powerful
2257: functions \kbd{bnfclassunit} which is of the same type as \kbd{quadclassunit}
2258: seen above, but for general number fields, and hence much slower. See Chapter
2259: 3 for a detailed explanation of its use.
2260:
2261: First type \kbd{setrand(1)}: this resets the random seed (to make sure we get
2262: the exact same results). Now type \kbd{m = bnfclassunit(A)}, where \kbd{A} is
2263: the same polynomial as before. After some work, we get a matrix with one
2264: column, which does not look too horrible (this is because much of the really
2265: useful information has been discarded, this is not an \kbd{init} function).
2266: Let's extract the column with \kbd{v = m[,1]}.
2267:
2268: Then \kbd{v} is a vector with 10 components. We immediately recognize the first
2269: component as being the polynomial \kbd{A} once again. In fact member
2270: functions are still available for \kbd{m}, even though it was not created
2271: by an \kbd{init} function. So, let's try them: \kbd{m.pol} gives \kbd{A},
2272: \kbd{m.sign}, \kbd{m.disc}, \kbd{m.zk}, ok nothing really exciting. But new
2273: ones are available now: \kbd{m.no} tells us the class number is 4,
2274: \kbd{m.cyc} that it's cyclic (of order 4 but that we already knew),
2275: \kbd{m.gen} that it's generated by a prime ideal above seven (in HNF form). If
2276: you play a little bit with the \kbd{idealhnf} function you see that this ideal
2277: is \kbd{P[4]}.
2278:
2279: The regulator \kbd{m.reg} is equal to $3.794\dots$. \kbd{m.tu} tells us that
2280: the roots of unity in $K$ are exactly the sixth roots of 1 and gives a
2281: primitive root. Finally \kbd{m.fu} gives us a fundamental unit (which must be
2282: linked to the unit \kbd{u2} found above since the unit rank is 1).
2283:
2284: To find the relation (without trial and error, because in this case it is
2285: quite easy to see it directly!), let us use the logarithmic embeddings. Type
2286: \kbd{uf = m.fu[1]}, \kbd{uu= m.tu[2]} to get the generators of the unit
2287: group, then
2288:
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 COMPLETED -----
2473:
2474: \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 / doc