Estude fácil! Tem muito documento disponível na Docsity
Ganhe pontos ajudando outros esrudantes ou compre um plano Premium
Prepare-se para as provas
Estude fácil! Tem muito documento disponível na Docsity
Prepare-se para as provas com trabalhos de outros alunos como você, aqui na Docsity
Encontra documentos específicos para os exames da tua universidade
Prepare-se com as videoaulas e exercícios resolvidos criados a partir da grade da sua Universidade
Responda perguntas de provas passadas e avalie sua preparação.
Ganhe pontos para baixar
Ganhe pontos ajudando outros esrudantes ou compre um plano Premium
Arquivo tutorial do PARI
Tipologia: Notas de estudo
1 / 49
Esta página não é visível na pré-visualização
Não perca as partes importantes!
C. Batut, K. Belabas, D. Bernardi, H. Cohen, M. Olivier
Laboratoire A2X, U.M.R. 9936 du C.N.R.S. Universit´e Bordeaux I, 351 Cours de la Lib´eration 33405 TALENCE Cedex, FRANCE e-mail: [email protected]
http://pari.math.u-bordeaux.fr/
Copyright ©c 2000–2006 The PARI Group
Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies.
Permission is granted to copy and distribute modified versions, or translations, of this manual under the conditions for verbatim copying, provided also that the entire resulting derived work is distributed under the terms of a permission notice identical to this one.
PARI/GP is Copyright ©c 2000–2006 The PARI Group
PARI/GP is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY WHATSOEVER.
This booklet is a guided tour and a tutorial to the gp calculator. Many examples will be given, but each time a new function is used, the reader should look at the appropriate section in the User’s Manual to PARI/GP for detailed explanations. This chapter can be read independently, for example to get acquainted with the possibilities of gp without having to read the whole manual. At this point.
So you are sitting in front of your workstation (or terminal, or PC... ), and you type gp to get the program started (or click on the relevant icon, or select some menu item). It says hello in its particular manner, and then waits for you after its prompt, initially? (or something like gp >).
Type 2 + 2. What happens? Maybe not what you expect. First of all, of course, you should tell gp that your input is finished, and this is done by hitting the Return (or Newline, or Enter) key. If you do exactly this, you will get the expected answer. However some of you may be used to other systems like Gap, Macsyma, Magma or Maple. In this case, you will have subconsciously ended the line with a semicolon “;” before hitting Return, since this is how it is done on those systems. In that case, you will simply see gp answering you with a smug expression, i.e. a new prompt and no answer! This is because a semicolon at the end of a line tells gp not to print the result (it is still stored in the result history). You will certainly want to use this feature if the output is several pages long.
Try 27 * 37. Wow! even multiplication works. Actually, maybe those spaces are not necessary after all. Let’s try 27*37. Seems to be ok. We will still insert them in this document since it makes things easier to read, but as gp does not care about them, you don’t have to type them all.
Now this session is getting lengthy, so the second thing one needs to learn is to quit. Each system has its quit signal. In gp, you can use quit or \q (backslash q), the q being of course for quit. Try it.
Now you’ve done it! You’re out of gp, so how do you want to continue studying this tutorial? Get back in please.
Let’s get to more serious stuff. I seem to remember that the decimal expansion of 1/7 has some interesting properties. Let’s see what gp has to say about this. Type
1 / 7
What? This computer is making fun of me, it just spits back to me my own input, that’s not what I want!
Now stop complaining, and think a little. Mathematically, 1/7 is an element of the field Q of rational numbers, so how else but 1/7 can the computer give the answer to you? Well maybe 2/ 14 or 7−^1 , but why complicate matters?. Seriously, the basic point here is that PARI, hence gp, will almost always try to give you a result which is as precise as possible (we will see why “almost” later). Hence since here the result can be represented exactly, that’s what it gives you.
But I still want the decimal expansion of 1/7. No problem. Type one of the following:
1./ 7 1 / 7. 1./ 7. 1 / 7 + 0. etc...
Immediately a number of decimals of this fraction appear, 28 on most systems, 38 on the others, and the repeating pattern is 142857. The reason is that you have included in the operations numbers like 0., 1. or 7. which are imprecise real numbers, hence gp cannot give you an exact result.
Why 28 / 38 decimals by the way? Well, it is the default initial precision. This has been chosen so that the computations are very fast, and gives already 12 decimals more accuracy than conventional double precision floating point operations. The precise value depends on a technical reason: if your machine supports 64-bit integers (the standard C library can handle integers up to 264 ), the default precision is 38 decimals, and 28 otherwise. As the latter is probably the case, we will assume it henceforth. Of course, you can extend the precision (almost) as much as you like as we will see in a moment.
I’m getting bored, why don’t we get on with some more exciting stuff? Well, try exp(1). Presto, comes out the value of e to 28 digits. Try log(exp(1)). Well, we get a floating point number and not an exact 1, but pretty close! That’s what you lose by working numerically.
What could we try now? Hum, pi? The answer is not that enlightening. Pi? Ah. This works better. But let’s remember that gp distinguishes between uppercase and lowercase letters. pi was as meaningless to it as stupid garbage would have been: in both cases gp will just create a variable with that funny unknown name you just used. Try it! Note that it is actually equivalent to type stupidgarbage: all spaces are suppressed from the input. In the 27 * 37 example it was not so conspicuous as we had an operator to separate the two operands. This has important consequences for the writing of gp scripts. More about this later.
By the way, you can ask gp about any identifier you think it might know about: just type it, prepending a question mark “?”. Try ?Pi and ?pi for instance. On most systems, an extended online help should be available: try doubling the question mark to check whether it’s the case on yours: ??Pi. In fact the gp header already gave you that information if it was the case, just before the copyright message. As well, if it says something like “readline enabled” then you should have a look at the readline introduction in the User’s Manual before you go on: it will be much easier to type in examples and correct typos after you’ve done that.
Now try exp(Pi * sqrt(163)). Hmmm, we suspect that the last digit may be wrong, can this really be an integer? This is the time to change precision. Type \p 50, then try exp(Pi * sqrt(163)) again. We were right to suspect that the last decimal was incorrect, since we get quite a few nines in its place, but it is now convincingly clear that this is not an integer. Maybe it’s a bug in PARI, and the result is really an integer? Type
(log(%) / Pi)^
immediately after the preceding computation; % means the result of the last computed expression. More generally, the results are numbered %1, %2,... including the results that you do not want to see printed by putting a semicolon at the end of the line, and you can evidently use all these quantities in any further computations. The result seems to be indistinguishable from 163, hence it does not seem to be a bug.
In fact, it is known that exp(π ∗
n) not only is not an integer or a rational number, but is even a transcendental number when n is a non-zero rational number.
So gp is just a fancy calculator, able to give me more decimals than I will ever need? Not so, gp is incredibly more powerful than an ordinary calculator, independently of its arbitrary precision possibilities.
Note that the default command reacts differently according to the number of input argu- ments. This is not an uncommon behaviour for gp functions. You can see this from the online help, or the complete description in Chapter 3: any argument surrounded by braces {} in the function prototype is optional, which really means that a default argument will be supplied by gp. You can then check out from the text what effect a given value will have, and in particular the default one.
Try the following: starting in precision 28, type first default(format, "e0.100"), then exp(1). Where are my 100 significant digits? Well, default(format,) only changes the output format, but not the default precision. On the other hand, the \p command changes both the precision and the output format.
Another thing you better get used to pretty fast is error messages. Try typing 1/0. Couldn’t be clearer. Taking again our universal example in precision 28, type
floor(exp(75))
floor is the mathematician’s integer part, not to be confused with truncate, which is the computer scientist’s: floor(-3.4) is equal to −4 whereas truncate(-3.4) is equal to −3. You get a more cryptic error message, which you would immediately understand if you had read the additional comments of the preceding section. Since you were told not to read them, here’s the explanation: gp is unable to compute the integer part of exp(75) given only 28 decimals of accuracy, since it has 33 digits.
Some error messages are more cryptic and sometimes not so easy to understand. For instance, try log(x). It simply tells you that gp does not understand what log(x) is, although it does know the log function, as ?log will readily tell us.
Now let’s try sqrt(-1) to see what error message we get now. Haha! gp even knows about complex numbers, so impossible to trick it that way. Similarly, try typing log(-2), exp(I*Pi), I^I... So we have a lot of real and complex analysis at our disposal. There always is a specific branch of multivalued complex transcendental functions which is taken, specified in the manual. Again, beware that I and i are not the same thing. Compare I^2 with i^2 for instance.
Just for fun, let’s try 6*zeta(2) / Pi^2. Pretty close, no?
Now gp didn’t seem to know what log(x) was, although it did know how to compute numerical values of log. This is annoying. Maybe it knows the exponential function? Let’s give it a try. Type exp(x). What’s this? If you had any experience with other computer algebra systems, the answer should have simply been exp(x) again. But here the answer is the Taylor expansion of the function around x = 0, to 16 terms. 16 is the default seriesprecision, which can be changed by typing \ps n or default(seriesprecision, n) where n is the number of terms that you want in your power series. Note the O(x^16) which ends the series, and which is trademark of this type of object in gp. It is the familiar “big–oh” notation of analysis.
You thus automatically get the Taylor expansion of any function that can be expanded around 0, and incidentally this explains why we weren’t able to do anything with log(x) which is not defined at 0. (In fact gp knows about Laurent series, but log(x) is not meromorphic either at 0.) If we try log(1+x), then it works. But what if we wanted the expansion around a point different from 0? Well, you’re able to change x into x−a, aren’t you? So for instance you can type log(x+2) to have the expansion of log around x = 2. As exercises you can try
cos(x) cos(x)^2 + sin(x)^ exp(cos(x)) gamma(1 + x) exp(exp(x) - 1) 1 / tan(x)
for different values of serieslength (change it using \ps newvalue).
Let’s try something else: type (1 + x)^3. No O(x) here, since the result is a polynomial. Haha, but I have learnt that if you do not take exponents which are integers greater or equal to 0, you obtain a power series with an infinite number of non-zero terms. Let’s try. Type (1 + x)^(-3) (the parentheses around -3 are not necessary but make things easier to read). Surprise! Contrary to what we expected, we don’t get a power series but a rational function. Again this is for the same reason that 1 / 7 just gave you 1/7: the result being exact, PARI doesn’t see any reason to make it non-exact.
But I still want that power series. To obtain it, you can do as in the 1/7 example and type
(1 + x)^(-3) + O(x^16) (1 + x)^(-3) * (1 + O(x^16)) (1 + x + O(x^16))^(-3), etc...
(Not on this example, but there is a difference between the first 2 methods. Do you spot it?) Better yet, use the series constructor which transforms any object into a power series, using the current seriesprecision, and simply type
Ser( (1 + x)^(-3) )
Now try (1 + x)^(1/2): we obtain a power series, since the result is an object which PARI does not know how to represent exactly. (We could teach PARI about algebraic functions, but then take (1 + x)^Pi as another example.) This gives us still another solution to our preceding exercise: we can type (1 + x)^(-3.). Since -3. is not an exact quantity, PARI has no means to know that we are dealing with a rational function, and will instead give you the power series, this time with real instead of integer coefficients.
To summarize, in this section we have seen that in addition to integers, real numbers and rational numbers, PARI can handle complex numbers, polynomials, rational functions and power series. A large number of functions exist which handle these types, but in this tutorial we will only look at a few.
Of course PARI knows about vectors (rows and columns are distinguished, even though math- ematically there is no difference) and matrices. Type for example [1,2,3,4]. This gives the row vector whose coordinates are 1, 2, 3 and 4. If you want a column vector, type [1,2,3,4]~, the tilde meaning of course transpose. You don’t see much difference in the output, except for the tilde at the end. However, now type \b: lo and behold, the column vector appears as a proper vertical thingy now. The \b command is used mainly for this purpose. The length of a vector is given by, well length of course. The shorthand “cardinal” notation #v for length(v) is also available, for instance v[#v] is the last element of v.
Type m = [a,b,c; d,e,f]. You have just entered a matrix with 2 rows and 3 columns. Note that the matrix is entered by rows and the rows are separated by semicolons “;”. The matrix is printed naturally in a rectangle shape. If you want it printed horizontally just as you typed it, type \a, or if you want this type of printing to be the permanent default type default(output, 0). Type default(output, 1) if you want to come back to the original output mode.
Now type m[1,2], m[1,], m[,2]. Are explanations necessary? (In an expression such as m[j,k], the j always refers to the row number, and the k to the column number, and the first index is always 1, never 0. This default cannot be changed.)
Even better, type m[1,2] = 5; m. The semicolon also allows us to put several instructions on the same line; the final result is the output of the last statement on the line. Now type m[1,] = [15,-17,8]. No problem. Finally type m[,2] = [j,k]. You have an error message since you have typed a row vector, while m[,2] is a column vector. If you type instead m[,2] = [j,k]~ it works.
Type now h = mathilbert(20). You get the so-called “Hilbert matrix” whose coefficient of row i and column j is equal to (i + j − 1)−^1. Incidentally, the matrix h takes too much room. If you don’t want to see it, simply type a semi-colon “;” at the end of the line (h = mathilbert(20);). This is an example of a “precomputed” matrix, built into PARI. We will see a more general construction later.
What is interesting about Hilbert matrices is that first their inverses and determinants can be computed explicitly (and the inverse has integer coefficients), and second they are numerically very unstable, which make them a severe test for linear algebra packages in numerical analysis. Of course with PARI, no such problem can occur: since the coefficients are given as rational numbers, the computation will be done exactly, so there cannot be any numerical error. Try it. Type d = matdet(h). The result is a rational number (of course) of numerator equal to 1 and denominator having 226 digits. How do I know, by the way? Well, type sizedigit(1/d). Or #Str(1/d). (The length of the character string representing the result.)
Now type hr = 1.* h; (do not forget the semicolon, we don’t want to see the result!), then dr = matdet(hr). You notice two things. First the computation, is much faster than in the rational case. (If your computer is too fast for you to notice, try again with h = mathilbert(40), or some larger value.) The reason for this is that PARI is handling real numbers with 28 digits of accuracy, while in the rational case it is handling integers having up to 226 decimal digits.
The second, more important, fact is that the result is terribly wrong. If you compare with 1.∗d computed earlier, which is the correct answer, you will see that only 2 decimals agree! (None agree if you replaced 20 by 40 as suggested above.) This catastrophic instability is as already mentioned one of the characteristics of Hilbert matrices. In fact, the situation is worse than that. Type norml2(1/h - 1/hr) (the function norml2 gives the square of the L^2 norm, i.e. the sum of the squares of the coefficients). The result is larger than 10^50 , showing that some coefficients of 1/hr are wrong by as much as 10^24 (the largest error is in fact equal to 4. 244 · 1024 for the coefficient
of row 15 and column 15, which is a 28 digit integer). To obtain the correct result after rounding for the inverse, we have to use a default precision of 56 digits (try it).
Although vectors and matrices can be entered manually, by typing explicitly their elements, very often the elements satisfy a simple law and one uses a different syntax. For example, as- sume that you want a vector whose i-th coordinate is equal to i^2. No problem, type for example vector(10,i, i^2) if you want a vector of length 10. Similarly, if you type
matrix(5,5, i,j, 1 / (i+j-1))
you will get the Hilbert matrix of order 5, hence the mathilbert function is in fact redundant. The i and j represent dummy variables which are used to number the rows and columns respectively (in the case of a vector only one is present of course). You must not forget, in addition to the dimensions of the vector or matrix, to indicate explicitly the names of these variables. You may omit the variables and the final expression to get zero entries, as in matrix(10,20).
Warning. The letter I is reserved for the complex number equal to the square root of −1. Hence it is forbidden to use it as a variable. Try typing vector(10,I, I^2), the error message that you get clearly indicates that gp does not consider I as a variable. There are two other reserved variable names: Pi and Euler. All function names are forbidden as well. On the other hand there is nothing special about i, pi or euler.
When creating vectors or matrices, it is often useful to use boolean operators and the if() statement. Indeed, an if expression has a value, which is of course equal to the evaluated part of the if. So for example you can type
matrix(8,8, i,j, if ((i-j)%2, 1, 0))
to get a checkerboard matrix of 1 and 0. Note however that a vector or matrix must be created first before being used. For example, it is forbidden to write
for (i = 1, 5, v[i] = 1/i)
if the vector v has not been created beforehand, for example by a v = vector(5) command. Another useful way to create vectors and matrices is to extract them from larger ones, using vecextract. For instance, if h is the 20 × 20 Hilbert matrix as above,
vecextract(h, "11..20", "11..20")}
is its lower right quadrant.
The last PARI types which we have not yet played with are closely linked to number theory. People not interested in number theory can skip ahead.
The first is the type “integer–modulo”. Let us see an example. Type n = 10^15 + 3. We want to know whether this number is prime or not. Of course we could make use of the built-in facilities of PARI, but let us do otherwise. We first trial divide by the built-in table of primes. We slightly cheat here and use a variant of the function factor which does exactly this. So type factor(n, 200000). The last argument tells factor to trial divide up to the given bound and stop at this point. Set it to 0 to trial divide by the full set of built-in primes, which goes up to 500000 by default.
As for all factoring function, the result is a 2 column matrix: the first column gives the primes and the second their exponents. Here we get a single row, telling us that if primes stopped at 200000 as we made factor believe, n would be prime. (Or is that a contradiction?) More seriously, n is not divisible by any prime up to 200000.
of the type “complex”. To start, we must specify which quadratic field we want to work in. For this, we use the function quadgen applied to the discriminant d (as opposed to the radicand) of the quadratic field. This returns a number (always printed as w) equal to (d + a)/2 where a is equal to 0 or 1 according to whether d is even or odd. The behavior of quadgen is a little special: although its result is always printed as w, the variable w itself is not set to that value. Hence it is necessary to write systematically w = quadgen(d) using the variable name w (or w1 etc. if you have several quadratic fields), otherwise things will get confusing.
So type w = quadgen(-163), then charpoly(w) which asks for the characteristic polynomial of w. The result shows what w will represent. You may ask for 1.*w to see which root of the quadratic has been taken, but this is rarely necessary. We can now play in the field Q(
Type for example w^10, norm(3 + 4*w), 1 / (4+w). More interesting, type a = Mod(1,23) * w then b = a^264. This is a generalization of Fermat’s theorem to quadratic fields. If you do not want to see the modulus 23 all the time, type lift(b).
Another example: type p = x^2 + wx + 5w + 7, then norm(p). We thus obtain the quartic equation over Q corresponding to the relative quadratic extension over Q(w) defined by p.
On the other hand, if you type wr = sqrt(w^2), do not expect to get back w. Instead, you get the numerical value, the function sqrt being considered as a “transcendental” function, even though it is algebraic. Type algdep(wr,2): this looks for algebraic relations involving the powers of w up to degree 2. This is one way to get w back. Similarly, type algdep(sqrt(3*w + 5), 4). See the user’s manual for the function algdep.
The fourth number-theoretic type is the type “polynomial–modulo”, i.e. polynomial modulo another polynomial. This type is used to work in general algebraic extensions, for example elements of number fields (if the base field is Q), or elements of finite fields (if the base field is Z/pZ for a prime p). In a sense it is a generalization of the type quadratic number. The syntax used is the same as for intmods. For example, instead of typing w = quadgen(-163), you can type
w = Mod(x, quadpoly(-163))
Then, exactly as in the quadratic case, you can type w^10, norm(3 + 4w), 1 / (4+w), a = Mod(1,23)w, b = a^264, obtaining of course the same results. (Type lift(... ) if you don’t want to see the polynomial x^2 - x + 41 repeated all the time.) Of course, you can work in any degree, not only quadratic. For the latter, the corresponding elementary operations will be slower than with quadratic numbers. Start the timer, then compare
w = quadgen(-163); W = Mod(x, quadpoly(-163)); a = 2 + w; A = 2 + W; b = 3 + w; B = 3 + W; for (i=1,10^5, a+b) for (i=1,10^5, A+B) for (i=1,10^5, ab) for (i=1,10^5, AB) for (i=1,10^5, a/b) for (i=1,10^5, A/B)
Don’t retype everything, use the arrow keys!
There is however a slight difference in behavior. Keeping our polmod w, type 1.*w. As you can see, the result is not the same. Type sqrt(w). Here, we obtain a vector with 2 components, the two components being the principal branch of the square root of all the possible embeddings
of w in C. More generally, if w was of degree n, we would get an n-component vector, and similarly for all transcendental functions.
We have at our disposal the usual arithmetic functions, plus a few others. Type a = Mod(x, x^3 - x - 1) defining a cubic extension. We can for example ask for b = a^5. Now assume we want to express a as a polynomial in b. This is possible since b is also a generator of the same field. No problem, type modreverse(b). This gives a new defining polynomial for the same field, i.e. the characteristic polynomial of b, and expresses a in terms of this new polmod, i.e. in terms of a. We will see this in more detail in the number field section.
Since PARI is aimed at number theorists, it is not surprising that there exists a large number of arithmetic functions; see the list by typing ?4. We have already seen several, such as factor. Note that factor handles not only integers, but also univariate polynomials. Type for example factor(x^200 - 1). You can also ask to factor a polynomial modulo a prime p (factormod) and even in a finite field which is not a prime field (factorff).
Evidently, you have functions for computing GCD’s (gcd), extended GCD’s (bezout), solving the Chinese remainder theorem (chinese) and so on.
In addition to the factoring facilities, you have a few functions related to primality testing such as isprime, ispseudoprime, precprime, and nextprime. As previously mentioned, only strong pseudoprimes are produced by the latter two (they pass the ispseudoprime test); the more sophisticated primality tests in isprime, being so much slower, are not applied by default.
We also have the usual multiplicative arithmetic functions: the M¨obius μ function (moebius), the Euler φ function (eulerphi), the ω and Ω functions (omega and bigomega), the σk functions (sigma), which compute sums of k-th powers of the positive divisors of a given integer, etc...
You can compute continued fractions. For example, type \p 1000, then contfrac(exp(1)): you obtain the continued fraction of the base of natural logarithms, which as you can see obeys a very simple pattern. Can you prove it?
In many cases, one wants to perform some task only when an arithmetic condition is satisfied. gp gives you the following functions: isprime as mentioned above, Z issquare, isfundamental to test whether an integer is a fundamental discriminant (i.e. 1 or the discriminant of a quadratic field), and the forprime, fordiv and sumdiv loops. Assume for example that we want to compute the product of all the divisors of a positive integer n. The easiest way is to write
p = 1; fordiv(n,d, p *= d); p
(There is a simple formula for this product in terms of n and the number of its divisors: find and prove it!) The notation p *= d is just a shorthand for p = p * d.
If we want to know the list of primes p less than 1000 such that 2 is a primitive root modulo p, one way would be to write:
forprime(p=3,1000, if (znprimroot(p) == 2, print(p)))
Note that this assumes that znprimroot returns the smallest primitive root, and this is indeed the case. Had we not known about this, we could have written
forprime(p=3,1000, if (znorder(Mod(2,p)) == p-1, print(p)))
Can we do better? Perhaps, but then we’d better drop the requirement that the basis be triangular, since the latter is essentially canonical. Type M = H * qflll(H): its column give an LLL-reduced basis for Λ. (Of course, qflll(H) itself gives the base change matrix.) The LLL algorithm outputs a nice basis for a lattice given by an arbitrary basis, where nice means the basis vectors are almost orthogonal and short, with precise guarantees on their relations to the shortest vectors. Not really spectacular on this example, though.
Let us try something else, there should be an integer relation between the components of u = [log(15), log(5), log(3)]. So input u then type
m = matid(3); m[3,] = round(u * 10^25); v = qflll(m)[,1] \ first vector of the LLL-reduced basis u * v
Pretty close. In fact, lindep automates this kind of search for integer relations, using more sophis- ticated algorithms like PSLQ.
Let’s come back to Λ above, and our LLL basis in M. Type
G = M~*M \ Gram matrix m = qfminim(G, norml2(M[,1]), 100)
This enumerates the vectors in Λ which are shorter than the first LLL basis vector, at most 100 of them. There are m[ 1 ] = 6 such vectors, and m[3] gives half of them (-m[3] would complete the lot): they are the first 3 basis vectors! So these are optimally short, at least with respect to the Euclidean length. Let us try
m = qfminim(G, norml2(M[,4]), 100, 2)
(The flag 2 instructs qfminim to use a different enumeration strategy, which is much faster when we expect more short vectors than we want to store. Without the flag, this example requires several hours. This is an exponential time algorithm, after all!) This time, we find a slew of short vectors; matrank(m[3]) says the 100 we have are all included in a 2-dimensional space. Let us try
m = qfminim(G, norml2(M[,4]) - 1, 100000, 2)
this time we find 50840 vectors of the requested length, spanning a 4-dimensional space, which is actually generated by M[,1], M[,2] M[,3] and M[,5].
All the elementary transcendental functions and several higher transcendental functions are provided: Γ function, incomplete Γ function, error function, exponential integral, Bessel functions (H^1 , H^2 , I, J, K, N ), confluent hypergeometric functions, Riemann ζ function, polylogarithms, Weber functions, theta functions. More will be written if the need arises.
In this type of functions, the default precision plays an essential role. In almost all cases transcendental functions work in the following way. If the argument is exact, the result is computed using the current default precision. If the argument is not exact, the precision of the argument is used for the computation. A note of warning however: even in this case the printed value is the current real format, usually the same as the default precision. In the present chapter we assume that your machine works with 32-bit long integers. If it is not the case, we leave it to you as a good exercise to make the necessary modifications.
Let’s assume that we have 28 decimals of default precision (this is what we get automatically at the start of a gp session on 32-bit machines). Type e = exp(1). We get the number e = 2. 718... to 28 decimals. Let us check how many correct decimals we really have. Change the precision to a substantially higher value, for example by typing \p 50. Then type e, then exp(1) once again. This last value is the correct value of the mathematical constant e to 50 decimals, while the variable e shows the value that was computed to 28 decimals. Clearly they coincide to exactly 29 significant digits.
So 28 digits are printed, but how many significant digits are actually contained in the variable e? Type #e which indicates we have exactly 3 mantissa words. Since 3 ln(2^32 )/ ln(10) ≈ 28 .9 we see that we have 28 or 29 significant digits (on 32-bit machines).
Come back to 28 decimals (\p 28). If we type exp(1.) you can check that we also obtain 28 decimals. However, type f = exp(1 + 1E-30). Although the default precision is still 28, you can check using the method above that we have in fact 59 significant digits! The reason is that 1
Now type cos(1E-15). The result is printed as 1. 0000.. ., but is of course not exactly equal to 1. Using #%, we see that the result has 7 mantissa words, giving us the possibility of having 67 correct significant digits. In fact (look in precision 100), only 60 are correct. PARI gives you as much as it can, and since 6 mantissa words would have given you only 57 digits, it uses 7. But why does it give so precise a result? Well, it is the same reason as before. When x is close to 1, cos(x) is close to 1 − x^2 /2, hence the precision is going to be approximately the same as this quantity, which will be 1 − 0. 5 ∗ 10 −^30 where 0. 5 ∗ 10 −^30 is considered with 28 significant digit accuracy, hence the result will have approximately 28 + 30 = 58 significant digits.
Unfortunately, this philosophy cannot go too far. For example, when you type cos(0), gp should give you exactly 1. Since it is reasonable for a program to assume that a transcendental function never gives you an exact result, gp gives you 1. 000... to one more mantissa word than the current precision.
Let’s see some more transcendental functions at work. Type gamma(10). No problem (type 9! to check). Type gamma(100). The number is now written in exponential format because the default accuracy is too small to give the correct result (type 99! to get the complete answer). To get the complete value, there are two solutions. The first and most natural one is to increase the precision. Since gamma(100) has 156 decimal digits, type \p 170 to be on the safe side, then gamma(100) once again. After some work, the printed result is this time perfectly correct.
However, this is probably not the best way to proceed. Come back first to the default precision (type \p 28). As the gamma function increases very rapidly, one usually uses its logarithm. Type lngamma(100). This time the result has a reasonable size, and is exactly equal to log(99!).
Try gamma(1/2 + 10*I). No problem, we have the complex gamma function. Now type
t = 1000; z = gamma(1 + It) * t^(-1/2) * exp(Pi/2t) / sqrt(2*Pi),
then norm(z). The latter is very close to 1, in accordance with the complex Stirling formula.
Let’s play now with the Riemann zeta function. First turn on the timer (type #). Type zeta(2), then Pi^2/6. This seems correct. Type zeta(3). All this takes essentially no time at all. However, type zeta(3.). Although the result is the same, you will notice that the time
Now type 2 * agm(1,I) / (1+I). As you see, the complex AGM also works, although one must be careful with its definition. The result found is almost identical to the previous one. Do you see why?
Finally, type agm(1, 1 + 7 + O(7^10)). So we also have p-adic AGM. Note however that since the square root of a p-adic number is not in general an element of the same p-adic field, only certain p-adic AGMs can be computed. In addition, when p = 2, the congruence restriction is that agm(a,b) can be computed only when a/b is congruent to 1 modulo 16, and not 8 as could be expected.
Now type ?3. This gives you the list of all transcendental functions. Instead of continuing with more examples, we suggest that you experiment yourself with the list of functions. In each case, try integer, real, complex and p-adic arguments. You will notice that some have not been implemented (or do not have a reasonable definition).
Although not written to be a numerical analysis package, PARI can nonetheless perform some numerical computations. Since linear algebra and polynomial computations are treated somewhere else, this section focuses on solving equations and various methods of summation.
You of course know the formula π = 4(1 − 13 + 15 − 17 + · · ·) which is deduced from the power series expansion of atan(x). You also know that π cannot be computed from this formula, since the convergence is so slow. Right? Wrong! Type \p 100 to make it more interesting, then 4 * sumalt(k=0, (-1)^k/(2*k + 1)). In a split second, we get π to 100 significant digits (type Pi to check).
Similarly, try sumpos(k=1, 1 / k^2). Although once again the convergence is slow, the sum- mation is rather fast; compare with the exact result Pi^2/6. This is less impressive because a bit slower than for alternating sums, but still useful.
Even better, sumalt can be used to sum divergent series! Type
zet(s) = sumalt(k=1, (-1)^(k-1) / k^s) / (1 - 2^(1-s))
Then for positive values of s different from 1, zet(s) is equal to zeta(s) and the series converges, albeit slowly; sumalt doesn’t care however. For negative s, the series diverges, but zet(s) still gives the correct result! (Namely, the value of a suitable analytic continuation.) Try zet(-1), zet(-2), zet(-1.5), and compare with the corresponding values of zeta. You should not push the game too far: zet(-100), for example, gives a completely wrong answer.
Try zet(I), and compare with zeta(I). Even (some) complex values work, although the sum is not alternating any more! Similarly, try
sumalt(n=1, (-1)^n / (n+I))
More traditional functions are the numerical integration functions. Come back to \p 28 since these routines are very slow when working with too many significant digits. Try intnum(t=1,2, 1/t) and presto! you get 26 decimals of log(2). Look at Chapter 3 to see the available integration functions.
With PARI, however, you can go further since complex types are allowed. For example, assume that we want to know the location of the zeros of the function h(z) = ez^ − z. We use Cauchy’s
theorem, which tells us that the number of zeros in a disk of radius r centered around the origin is equal to 1 2 iπ
Cr
h′(z) h(z)
dz ,
where Cr is the circle of radius r centered at the origin. The function we want to integrate is
fun(z) = local(u = exp(z)); (u-1) / (u-z)
(Here u is a local variable to the function f: whenever a function is called, gp fills its argument list with the actual arguments given, and initializes the other declared parameters and local variables to 0. It will then restore their former values upon exit. If we had not declared u in the function prototype, it would be considered as a global variable, whose value would be permanently changed. It is not mandatory to declare in this way all parameters, but beware of side effects!)
Type now: zero(r) = r/(2Pi) * intnum(t=0, 2Pi, real( fun(rexp(It)) * exp(I*t) ))
The function zero(r) will count the number of zeros of fun whose modulus is less than r: we simply made the change of variable z = r ∗ exp(i ∗ t), and took the real part to avoid integrating the imaginary part. Actually, there is a built-in function intcirc to integrate over a circle, yielding the much simpler:
zero(r) = intcirc(z=0, r, fun(z)) (This is a little faster than the previous implementation, and no less accurate.)
We may type \p 9 since we know that the result is a small integer (but the computations should be instantaneous even at \p 100 or so), then zero(1), zero(1.5). The result tells us that there are no zeros inside the unit disk, but that there are two (necessarily complex conjugate) in the annulus 1 < |z| < 1 .5. For the sake of completeness, let us compute them. Let z = x + iy be such a zero, with x and y real. Then the equation ez^ − z = 0 implies, after elementary transformations, that e^2 x^ = x^2 +y^2 and that ex^ cos(y) = x. Hence y = ±
e^2 x^ − x^2 and hence ex^ cos(
e^2 x^ − x^2 ) = x. Therefore, type
fun(x) = local(u = exp(x)); u * cos(sqrt(u^2 - x^2)) - x
Then fun(0) is positive while fun(1) is negative. Come back to precision 28 and type
x0 = solve(x=0,1, fun(x)) z = x0 + Isqrt(exp(2x0) - x0^2)
which is the required zero. As a check, type exp(z) - z.
Of course you can integrate over contours which are more complicated than circles, but you must perform yourself the changes of variable as we have done above to reduce the integral to a number of integrals on line segments.
The example above also shows the use of the solve function. To use solve on functions of a complex variable, it is necessary to reduce the problem to a real one. For example, to find the first complex zero of the Riemann zeta function as above, we could try typing
solve(t=14,15, real( zeta(1/2 + I*t) )),
but this does not work because the real part is positive for t = 14 and 15. As it happens, the imaginary part works. Type
