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