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