PARI C-library interface¶
AUTHORS:
- William Stein (2006-03-01): updated to work with PARI 2.2.12-beta
- William Stein (2006-03-06): added newtonpoly
- Justin Walker: contributed some of the function definitions
- Gonzalo Tornaria: improvements to conversions; much better error handling.
- Robert Bradshaw, Jeroen Demeyer, William Stein (2010-08-15): Upgrade to PARI 2.4.3 (trac ticket #9343)
- Jeroen Demeyer (2011-11-12): rewrite various conversion routines (trac ticket #11611, trac ticket #11854, trac ticket #11952)
- Peter Bruin (2013-11-17): split off this file from gen.pyx (trac ticket #15185)
- Jeroen Demeyer (2014-02-09): upgrade to PARI 2.7 (trac ticket #15767)
- Jeroen Demeyer (2014-09-19): upgrade to PARI 2.8 (trac ticket #16997)
- Jeroen Demeyer (2015-03-17): automatically generate methods from
pari.desc
(trac ticket #17631 and trac ticket #17860)
EXAMPLES:
sage: pari('5! + 10/x')
(120*x + 10)/x
sage: pari('intnum(x=0,13,sin(x)+sin(x^2) + x)')
83.8179442684285 # 32-bit
84.1818153922297 # 64-bit
sage: f = pari('x^3-1')
sage: v = f.factor(); v
[x - 1, 1; x^2 + x + 1, 1]
sage: v[0] # indexing is 0-based unlike in GP.
[x - 1, x^2 + x + 1]~
sage: v[1]
[1, 1]~
Arithmetic obeys the usual coercion rules:
sage: type(pari(1) + 1)
<type 'sage.libs.pari.gen.gen'>
sage: type(1 + pari(1))
<type 'sage.libs.pari.gen.gen'>
GUIDE TO REAL PRECISION AND THE PARI LIBRARY
The default real precision in communicating with the PARI library is the same as the default Sage real precision, which is 53 bits. Inexact Pari objects are therefore printed by default to 15 decimal digits (even if they are actually more precise).
Default precision example (53 bits, 15 significant decimals):
sage: a = pari(1.23); a
1.23000000000000
sage: a.sin()
0.942488801931698
Example with custom precision of 200 bits (60 significant decimals):
sage: R = RealField(200)
sage: a = pari(R(1.23)); a # only 15 significant digits printed
1.23000000000000
sage: R(a) # but the number is known to precision of 200 bits
1.2300000000000000000000000000000000000000000000000000000000
sage: a.sin() # only 15 significant digits printed
0.942488801931698
sage: R(a.sin()) # but the number is known to precision of 200 bits
0.94248880193169751002382356538924454146128740562765030213504
It is possible to change the number of printed decimals:
sage: R = RealField(200) # 200 bits of precision in computations
sage: old_prec = pari.set_real_precision(60) # 60 decimals printed
sage: a = pari(R(1.23)); a
1.23000000000000000000000000000000000000000000000000000000000
sage: a.sin()
0.942488801931697510023823565389244541461287405627650302135038
sage: pari.set_real_precision(old_prec) # restore the default printing behavior
60
Unless otherwise indicated in the docstring, most Pari functions that return inexact objects use the precision of their arguments to decide the precision of the computation. However, if some of these arguments happen to be exact numbers (integers, rationals, etc.), an optional parameter indicates the precision (in bits) to which these arguments should be converted before the computation. If this precision parameter is missing, the default precision of 53 bits is used. The following first converts 2 into a real with 53-bit precision:
sage: R = RealField()
sage: R(pari(2).sin())
0.909297426825682
We can ask for a better precision using the optional parameter:
sage: R = RealField(150)
sage: R(pari(2).sin(precision=150))
0.90929742682568169539601986591174484270225497
Warning regarding conversions Sage - Pari - Sage: Some care must be taken when juggling inexact types back and forth between Sage and Pari. In theory, calling p=pari(s) creates a Pari object p with the same precision as s; in practice, the Pari library’s precision is word-based, so it will go up to the next word. For example, a default 53-bit Sage real s will be bumped up to 64 bits by adding bogus 11 bits. The function p.python() returns a Sage object with exactly the same precision as the Pari object p. So pari(s).python() is definitely not equal to s, since it has 64 bits of precision, including the bogus 11 bits. The correct way of avoiding this is to coerce pari(s).python() back into a domain with the right precision. This has to be done by the user (or by Sage functions that use Pari library functions in gen.pyx). For instance, if we want to use the Pari library to compute sqrt(pi) with a precision of 100 bits:
sage: R = RealField(100)
sage: s = R(pi); s
3.1415926535897932384626433833
sage: p = pari(s).sqrt()
sage: x = p.python(); x # wow, more digits than I expected!
1.7724538509055160272981674833410973484
sage: x.prec() # has precision 'improved' from 100 to 128?
128
sage: x == RealField(128)(pi).sqrt() # sadly, no!
False
sage: R(x) # x should be brought back to precision 100
1.7724538509055160272981674833
sage: R(x) == s.sqrt()
True
Elliptic curves and precision: If you are working with elliptic curves, you should set the precision for each method:
sage: e = pari([0,0,0,-82,0]).ellinit()
sage: eta1 = e.elleta(precision=100)[0]
sage: eta1.sage()
3.6054636014326520859158205642077267748
sage: eta1 = e.elleta(precision=180)[0]
sage: eta1.sage()
3.60546360143265208591582056420772677481026899659802474544
Number fields and precision: TODO
TESTS:
Check that output from PARI’s print command is actually seen by Sage (trac ticket #9636):
sage: pari('print("test")')
test
-
class
sage.libs.pari.pari_instance.
PariInstance
¶ Bases:
sage.libs.pari.pari_instance.PariInstance_auto
Initialize the PARI system.
INPUT:
size
– long, the number of bytes for the initial PARI stack (see note below)maxprime
– unsigned long, upper limit on a precomputed prime number table (default: 500000)
Note
In Sage, the PARI stack is different than in GP or the PARI C library. In Sage, instead of the PARI stack holding the results of all computations, it only holds the results of an individual computation. Each time a new Python/PARI object is computed, it it copied to its own space in the Python heap, and the memory it occupied on the PARI stack is freed. Thus it is not necessary to make the stack very large. Also, unlike in PARI, if the stack does overflow, in most cases the PARI stack is automatically increased and the relevant step of the computation rerun.
This design obviously involves some performance penalties over the way PARI works, but it scales much better and is far more robust for large projects.
Note
If you do not want prime numbers, put
maxprime=2
, but be careful because many PARI functions require this table. If you get the error message “not enough precomputed primes”, increase this parameter.-
allocatemem
(s=0, silent=False)¶ Change the PARI stack space to the given size (or double the current size if
s
is).
If
and insufficient memory is avaible to double, the PARI stack will be enlarged by a smaller amount. In any case, a
MemoryError
will be raised if the requested memory cannot be allocated.The PARI stack is never automatically shrunk. You can use the command
pari.allocatemem(10^6)
to reset the size to, which is the default size at startup. Note that the results of computations using Sage’s PARI interface are copied to the Python heap, so they take up no space in the PARI stack. The PARI stack is cleared after every computation.
It does no real harm to set this to a small value as the PARI stack will be automatically doubled when we run out of memory. However, it could make some calculations slower (since they have to be redone from the start after doubling the stack until the stack is big enough).
INPUT:
s
- an integer (default: 0). A non-zero argument should be the size in bytes of the new PARI stack. Ifis zero, try to double the current stack size.
EXAMPLES:
sage: pari.allocatemem(10^7) PARI stack size set to 10000000 bytes sage: pari.allocatemem() # Double the current size PARI stack size set to 20000000 bytes sage: pari.stacksize() 20000000 sage: pari.allocatemem(10^6) PARI stack size set to 1000000 bytes
The following computation will automatically increase the PARI stack size:
sage: a = pari('2^100000000')
a
is now a Python variable on the Python heap and does not take up any space on the PARI stack. The PARI stack is still large because of the computation ofa
:sage: pari.stacksize() 16000000 sage: pari.allocatemem(10^6) PARI stack size set to 1000000 bytes sage: pari.stacksize() 1000000
TESTS:
Do the same without using the string interface and starting from a very small stack size:
sage: pari.allocatemem(1) PARI stack size set to 1024 bytes sage: a = pari(2)^100000000 sage: pari.stacksize() 16777216
-
complex
(re, im)¶ Create a new complex number, initialized from re and im.
-
debugstack
()¶ Print the internal PARI variables
top
(top of stack),avma
(available memory address, think of this as the stack pointer),bot
(bottom of stack).EXAMPLE:
sage: pari.debugstack() # random top = 0x60b2c60 avma = 0x5875c38 bot = 0x57295e0 size = 1000000
-
default
(variable, value=None)¶
-
double_to_gen
(x)¶
-
euler
(precision=0)¶ Return Euler’s constant to the requested real precision (in bits).
EXAMPLES:
sage: pari.euler() 0.577215664901533 sage: pari.euler(precision=100).python() 0.577215664901532860606512090082...
-
factorial
(n)¶ Return the factorial of the integer n as a PARI gen.
EXAMPLES:
sage: pari.factorial(0) 1 sage: pari.factorial(1) 1 sage: pari.factorial(5) 120 sage: pari.factorial(25) 15511210043330985984000000
-
genus2red
(P, P0=None)¶ Let
be a polynomial with integer coefficients. Determines the reduction of the (proper, smooth) genus 2 curve
, defined by the hyperelliptic equation
. The special syntax
genus2red([P,Q])
is also allowed, where the polynomialsand
have integer coefficients, to represent the model
.
EXAMPLES:
sage: x = polygen(QQ) sage: pari.genus2red([-5*x^5, x^3 - 2*x^2 - 2*x + 1]) [1416875, [2, -1; 5, 4; 2267, 1], x^6 - 240*x^4 - 2550*x^3 - 11400*x^2 - 24100*x - 19855, [[2, [2, [Mod(1, 2)]], []], [5, [1, []], ["[V] page 156", [3]]], [2267, [2, [Mod(432, 2267)]], ["[I{1-0-0}] page 170", []]]]]
This is the old deprecated syntax:
sage: pari.genus2red(x^3 - 2*x^2 - 2*x + 1, -5*x^5) doctest:...: DeprecationWarning: The 2-argument version of genus2red() is deprecated, use genus2red(P) or genus2red([P,Q]) instead See http://trac.sagemath.org/16997 for details. [1416875, [2, -1; 5, 4; 2267, 1], x^6 - 240*x^4 - 2550*x^3 - 11400*x^2 - 24100*x - 19855, [[2, [2, [Mod(1, 2)]], []], [5, [1, []], ["[V] page 156", [3]]], [2267, [2, [Mod(432, 2267)]], ["[I{1-0-0}] page 170", []]]]]
-
get_debug_level
()¶ Set the debug PARI C library variable.
-
get_real_precision
()¶ Returns the current PARI default real precision.
This is used both for creation of new objects from strings and for printing. It is the number of digits IN DECIMAL in which real numbers are printed. It also determines the precision of objects created by parsing strings (e.g. pari(‘1.2’)), which is not the normal way of creating new pari objects in Sage. It has no effect on the precision of computations within the pari library.
EXAMPLES:
sage: pari.get_real_precision() 15
-
get_series_precision
()¶
-
getrand
()¶ Returns PARI’s current random number seed.
OUTPUT:
GEN of type t_VECSMALL
EXAMPLES:
sage: a = pari.getrand() sage: a.type() 't_INT'
-
init_primes
(M)¶ Recompute the primes table including at least all primes up to M (but possibly more).
EXAMPLES:
sage: pari.init_primes(200000)
We make sure that ticket trac ticket #11741 has been fixed:
sage: pari.init_primes(2^30) Traceback (most recent call last): ... ValueError: Cannot compute primes beyond 436273290
-
matrix
(m, n, entries=None)¶ matrix(long m, long n, entries=None): Create and return the m x n PARI matrix with given list of entries.
-
new_with_bits_prec
(s, precision)¶ pari.new_with_bits_prec(self, s, precision) creates s as a PARI gen with (at most) precision bits of precision.
-
nth_prime
(n)¶
-
pari_version
()¶
-
pi
(precision=0)¶ Return the value of the constant pi to the requested real precision (in bits).
EXAMPLES:
sage: pari.pi() 3.14159265358979 sage: pari.pi(precision=100).python() 3.1415926535897932384626433832...
-
polchebyshev
(n, v=-1)¶ polchebyshev(n, v=x): Chebyshev polynomial of the first kind of degree n, in variable v.
EXAMPLES:
sage: pari.polchebyshev(7) 64*x^7 - 112*x^5 + 56*x^3 - 7*x sage: pari.polchebyshev(7, 'z') 64*z^7 - 112*z^5 + 56*z^3 - 7*z sage: pari.polchebyshev(0) 1
-
polcyclo
(n, v=-1)¶ polcyclo(n, v=x): cyclotomic polynomial of degree n, in variable v.
EXAMPLES:
sage: pari.polcyclo(8) x^4 + 1 sage: pari.polcyclo(7, 'z') z^6 + z^5 + z^4 + z^3 + z^2 + z + 1 sage: pari.polcyclo(1) x - 1
-
polcyclo_eval
(n, v)¶ polcyclo_eval(n, v): value of the nth cyclotomic polynomial at value v.
EXAMPLES:
sage: pari.polcyclo_eval(8, 2) 17 sage: cyclotomic_polynomial(8)(2) 17
-
pollegendre
(n, v=-1)¶ pollegendre(n, v=x): Legendre polynomial of degree n (n C-integer), in variable v.
EXAMPLES:
sage: pari.pollegendre(7) 429/16*x^7 - 693/16*x^5 + 315/16*x^3 - 35/16*x sage: pari.pollegendre(7, 'z') 429/16*z^7 - 693/16*z^5 + 315/16*z^3 - 35/16*z sage: pari.pollegendre(0) 1
-
polsubcyclo
(n, d, v=-1)¶ polsubcyclo(n, d, v=x): return the pari list of polynomial(s) defining the sub-abelian extensions of degree
of the cyclotomic field
, where
divides
.
EXAMPLES:
sage: pari.polsubcyclo(8, 4) [x^4 + 1] sage: pari.polsubcyclo(8, 2, 'z') [z^2 - 2, z^2 + 1, z^2 + 2] sage: pari.polsubcyclo(8, 1) [x - 1] sage: pari.polsubcyclo(8, 3) []
-
poltchebi
(*args, **kwds)¶ Deprecated: Use
polchebyshev()
instead. See trac ticket #18203 for details.
-
prime_list
(n)¶ prime_list(n): returns list of the first n primes
To extend the table of primes use pari.init_primes(M).
INPUT:
n
- C long
OUTPUT:
gen
- PARI list of first n primes
EXAMPLES:
sage: pari.prime_list(0) [] sage: pari.prime_list(-1) [] sage: pari.prime_list(3) [2, 3, 5] sage: pari.prime_list(10) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] sage: pari.prime_list(20) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71] sage: len(pari.prime_list(1000)) 1000
-
primes_up_to_n
(n)¶ Return the primes <= n as a pari list.
EXAMPLES:
sage: pari.primes_up_to_n(1) [] sage: pari.primes_up_to_n(20) [2, 3, 5, 7, 11, 13, 17, 19]
-
read
(filename)¶ Read a script from the named filename into the interpreter. The functions defined in the script are then available for use from Sage/PARI. The result of the last command in
filename
is returned.EXAMPLES:
Create a gp file:
sage: import tempfile sage: gpfile = tempfile.NamedTemporaryFile(mode="w") sage: gpfile.file.write("mysquare(n) = {\n") sage: gpfile.file.write(" n^2;\n") sage: gpfile.file.write("}\n") sage: gpfile.file.write("polcyclo(5)\n") sage: gpfile.file.flush()
Read it in Sage, we get the result of the last line:
sage: pari.read(gpfile.name) x^4 + x^3 + x^2 + x + 1
Call the function defined in the gp file:
sage: pari('mysquare(12)') 144
-
set_debug_level
(level)¶ Set the debug PARI C library variable.
-
set_real_precision
(n)¶ Sets the PARI default real precision in decimal digits.
This is used both for creation of new objects from strings and for printing. It is the number of digits IN DECIMAL in which real numbers are printed. It also determines the precision of objects created by parsing strings (e.g. pari(‘1.2’)), which is not the normal way of creating new pari objects in Sage. It has no effect on the precision of computations within the pari library.
Returns the previous PARI real precision.
EXAMPLES:
sage: pari.set_real_precision(60) 15 sage: pari('1.2') 1.20000000000000000000000000000000000000000000000000000000000 sage: pari.set_real_precision(15) 60
-
set_series_precision
(n)¶
-
setrand
(seed)¶ Sets PARI’s current random number seed.
INPUT:
seed
– either a strictly positive integer or a GEN of typet_VECSMALL
as output bygetrand()
This should not be called directly; instead, use Sage’s global random number seed handling in
sage.misc.randstate
and callcurrent_randstate().set_seed_pari()
.EXAMPLES:
sage: pari.setrand(50) sage: a = pari.getrand() sage: pari.setrand(a) sage: a == pari.getrand() True
TESTS:
Check that invalid inputs are handled properly (trac ticket #11825):
sage: pari.setrand("foobar") Traceback (most recent call last): ... PariError: incorrect type in setrand (t_POL)
-
stacksize
()¶ Returns the current size of the PARI stack, which is
by default. However, the stack size is automatically doubled when needed. It can also be set directly using
pari.allocatemem()
.EXAMPLES:
sage: pari.stacksize() 1000000
-
vector
(n, entries=None)¶ vector(long n, entries=None): Create and return the length n PARI vector with given list of entries.
EXAMPLES:
sage: pari.vector(5, [1, 2, 5, 4, 3]) [1, 2, 5, 4, 3] sage: pari.vector(2, [x, 1]) [x, 1] sage: pari.vector(2, [x, 1, 5]) Traceback (most recent call last): ... IndexError: length of entries (=3) must equal n (=2)
-
class
sage.libs.pari.pari_instance.
PariInstance_auto
¶ Bases:
sage.structure.parent_base.ParentWithBase
Part of the
PariInstance
class containing auto-generated functions.You must never use this class directly (in fact, Sage may crash if you do), use the derived class
PariInstance
instead.-
Catalan
(precision=0)¶ Catalan’s constant
. Note that
Catalan
is one of the few reserved names which cannot be used for user variables.
-
Euler
(precision=0)¶ Euler’s constant
. Note that
Euler
is one of the few reserved names which cannot be used for user variables.
-
I
()¶ The complex number
.
-
Pi
(precision=0)¶ The constant
(
). Note that
Pi
is one of the few reserved names which cannot be used for user variables.
-
addhelp
(sym, str)¶ Changes the help message for the symbol
sym
. The string str is expanded on the spot and stored as the online help forsym
. It is recommended to document global variables and user functions in this way, althoughgp
will not protest if you don’t.You can attach a help text to an alias, but it will never be shown: aliases are expanded by the
?
help operator and we get the help of the symbol the alias points to. Nothing prevents you from modifying the help of built-in PARI functions. But if you do, we would like to hear why you needed it!Without
addhelp
, the standard help for user functions consists of its name and definition.gp> f(x) = x^2; gp> ?f f = (x)->x^2
Once addhelp is applied to
, the function code is no longer included. It can still be consulted by typing the function name:
gp> addhelp(f, "Square") gp> ?f Square gp> f %2 = (x)->x^2
-
bernfrac
(x)¶ Bernoulli number
, where
,
,
,..., expressed as a rational number. The argument
should be of type integer.
-
bernpol
(n, v=None)¶ Bernoulli polynomial
in variable
.
? bernpol(1) %1 = x - 1/2 ? bernpol(3) %2 = x^3 - 3/2*x^2 + 1/2*x
-
bernreal
(x, precision=0)¶ Bernoulli number
, as
bernfrac
, butis returned as a real number (with the current precision).
-
bernvec
(x)¶ Creates a vector containing, as rational numbers, the Bernoulli numbers
,
,...,
. This routine is obsolete. Use
bernfrac
instead each time you need a Bernoulli number in exact form.Note. This routine is implemented using repeated independent calls to
bernfrac
, which is faster than the standard recursion in exact arithmetic. It is only kept for backward compatibility: it is not faster than individual calls tobernfrac
, its output uses a lot of memory space, and coping with the index shift is awkward.
-
default
(key=None, val=None)¶ Returns the default corresponding to keyword key. If val is present, sets the default to val first (which is subject to string expansion first). Typing
default()
(or\d
) yields the complete default list as well as their current values. Seedefaults
(in the PARI manual) for an introduction to GP defaults,gp_defaults
(in the PARI manual) for a list of available defaults, andmeta
(in the PARI manual) for some shortcut alternatives. Note that the shortcuts are meant for interactive use and usually display more information thandefault
.
-
ellmodulareqn
(N, x=None, y=None)¶ Return a vector [
eqn
,:math:] where
eqn
is a modular equation of level, i.e. a bivariate polynomial with integer coefficients;
indicates the type of this equation: either canonical (
) or Atkin (
). This function currently requires the package
seadata
to be installed and is limited to,
prime.
Let
be the
-invariant function. The polynomial
eqn
satisfies the following functional equation, which allows to compute the values of the classical modular polynomialof prime level
, such that
, while being much smaller than the latter:
- for canonical type:
, where
;
- for Atkin type:
.
In both cases,
is a suitable modular function (see below).
The following GP function returns values of the classical modular polynomial by eliminating
in the above two equations, for
or
.
classicaleqn(N, X='X, Y='Y)= { my(E=ellmodulareqn(N), P=E[1], t=E[2], Q, d); if(poldegree(P,'y)>2,error("level unavailable in classicaleqn")); if (t == 0, my(s = 12/gcd(12,N-1)); Q = 'x^(N+1) * substvec(P,['x,'y],[N^s/'x,Y]); d = N^(s*(2*N+1)) * (-1)^(N+1); , Q = subst(P,'y,Y); d = (X-Y)^(N+1)); polresultant(subst(P,'y,X), Q) / d; }
More precisely, let
be the Atkin-Lehner involution; we have
and the function
also satisfies:
- for canonical type:
;
- for Atkin type:
.
Furthermore, for an equation of canonical type,
is the standard
-quotient
where
is Dedekind’s eta function, which is invariant under
.
- for canonical type:
-
extern
(str)¶ The string str is the name of an external command (i.e. one you would type from your UNIX shell prompt). This command is immediately run and its output fed into
gp
, just as if read from a file.
-
externstr
(str)¶ The string str is the name of an external command (i.e. one you would type from your UNIX shell prompt). This command is immediately run and its output is returned as a vector of GP strings, one component per output line.
-
factorial
(x, precision=0)¶ Factorial of
. The expression
gives a result which is an integer, while
gives a real number.
-
fibonacci
(x)¶ Fibonacci number.
-
galoisgetpol
(a, b=0, s=1)¶ Query the galpol package for a polynomial with Galois group isomorphic to GAP4(a,b), totally real if
(default) and totally complex if
. The output is a vector [
pol
,den
] wherepol
is the polynomial of degreeden
is the denominator ofnfgaloisconj(pol)
. Pass it as an optional argument togaloisinit
ornfgaloisconj
to speed them up:
? [pol,den] = galoisgetpol(64,4,1); ? G = galoisinit(pol); time = 352ms ? galoisinit(pol, den); \\ passing 'den' speeds up the computation time = 264ms ? % == %` %4 = 1 \\ same answer
If
and
are omitted, return the number of isomorphism classes of groups of order
.
-
getabstime
()¶ Returns the CPU time (in milliseconds) elapsed since
gp
startup. This provides a reentrant version ofgettime
:my (t = getabstime()); ... print("Time: ", getabstime() - t);
For a version giving wall-clock time, see
getwalltime
.
-
getenv
(s)¶ Return the value of the environment variable
s
if it is defined, otherwise return 0.
-
getheap
()¶ Returns a two-component row vector giving the number of objects on the heap and the amount of memory they occupy in long words. Useful mainly for debugging purposes.
-
getrand
()¶ Returns the current value of the seed used by the pseudo-random number generator
random
. Useful mainly for debugging purposes, to reproduce a specific chain of computations. The returned value is technical (reproduces an internal state array), and can only be used as an argument tosetrand
.
-
getstack
()¶ Returns the current value of
, i.e. the number of bytes used up to now on the stack. Useful mainly for debugging purposes.
-
gettime
()¶ Returns the CPU time (in milliseconds) used since either the last call to
gettime
, or to the beginning of the containing GP instruction (if insidegp
), whichever came last.For a reentrant version, see
getabstime
.For a version giving wall-clock time, see
getwalltime
.
-
getwalltime
()¶ Returns the time (in milliseconds) elapsed since the UNIX Epoch (1970-01-01 00:00:00 (UTC)).
my (t = getwalltime()); ... print("Time: ", getwalltime() - t);
-
input
()¶ Reads a string, interpreted as a GP expression, from the input file, usually standard input (i.e. the keyboard). If a sequence of expressions is given, the result is the result of the last expression of the sequence. When using this instruction, it is useful to prompt for the string by using the
print1
function. Note that in the present version 2.19 ofpari.el
, when usinggp
under GNU Emacs (seeemacs
(in the PARI manual)) one must prompt for the string, with a string which ends with the same prompt as any of the previous ones (a"? "
will do for instance).
-
install
(name, code, gpname=None, lib=None)¶ Loads from dynamic library lib the function name. Assigns to it the name gpname in this
gp
session, with prototype code (see below). If gpname is omitted, uses name. If lib is omitted, all symbols known togp
are available: this includes the whole oflibpari.so
and possibly others (such aslibc.so
).Most importantly,
install
gives you access to all non-static functions defined in the PARI library. For instance, the function \kbd{GEN addii(GEN x, GEN y)} adds two PARI integers, and is not directly accessible undergp
(it is eventually called by the+
operator of course):? install("addii", "GG") ? addii(1, 2) %1 = 3
It also allows to add external functions to the
gp
interpreter. For instance, it makes the functionsystem
obsolete:? install(system, vs, sys,/*omitted*/) ? sys("ls gp*") gp.c gp.h gp_rl.c
This works because
system
is part oflibc.so
, which is linked togp
. It is also possible to compile a shared library yourself and provide it to gp in this way: usegp2c
, or do it manually (see themodules_build
variable inpari.cfg
for hints).Re-installing a function will print a warning and update the prototype code if needed. However, it will not reload a symbol from the library, even if the latter has been recompiled.
Prototype. We only give a simplified description here, covering most functions, but there are many more possibilities. The full documentation is available in
libpari.dvi
, see??prototype
- First character
i
,l
,v
: return type int / long / void. (Default:GEN
) - One letter for each mandatory argument, in the same order as they appear
in the argument list:
G
(GEN
),&
(GEN*
),L
(long
),s
(char *
),n
(variable). p
to supplyrealprecision
(usuallylong prec
in the argument list),P
to supplyseriesprecision
(usually kbd{long precdl}).
We also have special constructs for optional arguments and default values:
DG
(optionalGEN
,NULL
if omitted),D&
(optionalGEN*
,NULL
if omitted),Dn
(optional variable,if omitted),
For instance the prototype corresponding to
long issquareall(GEN x, GEN *n = NULL)
is
lGD&
.Caution. This function may not work on all systems, especially when
gp
has been compiled statically. In that case, the first use of an installed function will provoke a Segmentation Fault (this should never happen with a dynamically linked executable). If you intend to use this function, please check first on some harmless example such as the one above that it works properly on your machine.- First character
-
intnumgaussinit
(n=0, precision=0)¶ Initialize tables for
-point Gauss-Legendre integration of a smooth function
lon a compact interval
at current
realprecision
. Ifis omitted, make a default choice
, suitable for analytic functions on
. The error is bounded by
so, if the interval length increases,
should be increased as well.
? T = intnumgaussinit(); ? intnumgauss(t=-1,1,exp(t), T) - exp(1)+exp(-1) %1 = -5.877471754111437540 E-39 ? intnumgauss(t=-10,10,exp(t), T) - exp(10)+exp(-10) %2 = -8.358367809712546836 E-35 ? intnumgauss(t=-1,1,1/(1+t^2), T) - Pi/2 %3 = -9.490148553624725335 E-22 ? T = intnumgaussinit(50); ? intnumgauss(t=-1,1,1/(1+t^2), T) - Pi/2 %5 = -1.1754943508222875080 E-38 ? intnumgauss(t=-5,5,1/(1+t^2), T) - 2*atan(5) %6 = -1.2[...]E-8
On the other hand, we recommend to split the integral and change variables rather than increasing
too much, see
intnumgauss
.
-
kill
(sym)¶ Restores the symbol
sym
to its “undefined” status, and deletes any help messages associated tosym
usingaddhelp
. Variable names remain known to the interpreter and keep their former priority: you cannot make a variable “less important” by killing it!? z = y = 1; y %1 = 1 ? kill(y) ? y \\ restored to ``undefined'' status %2 = y ? variable() %3 = [x, y, z] \\ but the variable name y is still known, with y > z !
For the same reason, killing a user function (which is an ordinary variable holding a
t_CLOSURE
) does not remove its name from the list of variable names.If the symbol is associated to a variable — user functions being an important special case —, one may use the quote operator
a = 'a
to reset variables to their starting values. However, this will not delete a help message associated toa
, and is also slightly slower thankill(a)
.? x = 1; addhelp(x, "foo"); x %1 = 1 ? x = 'x; x \\ same as 'kill', except we don't delete help. %2 = x ? ?x foo
On the other hand,
kill
is the only way to remove aliases and installed functions.? alias(fun, sin); ? kill(fun); ? install(addii, GG); ? kill(addii);
-
localprec
(p)¶ Set the real precision to
in the dynamic scope. All computations are performed as if
realprecision
was: transcendental constants (e.g.
Pi
) and conversions from exact to floating point inexact data usedecimal digits, as well as iterative routines implicitly using a floating point accuracy as a termination criterion (e.g.
solve
orintnum
). Butrealprecision
itself is unaffected and is “unmasked” when we exit the dynamic (not lexical) scope. In effect, this is similar tomy(prec = default(realprecision)); default(realprecision,p); ... default(realprecision, prec);
but is both less cumbersome, cleaner (no need to manipulate a global variable, which in fact never changes and is only temporarily masked) and more robust: if the above computation is interrupted or an exception occurs,
realprecision
will not be restored as intended.Such
localprec
statements can be nested, the innermost one taking precedence as expected. Beware thatlocalprec
follows the semantic oflocal
, notmy
: a subroutine called fromlocalprec
scope uses the local accuracy:? f()=precision(1.); ? f() %2 = 38 ? localprec(19); f() %3 = 19
Warning. Changing
realprecision
itself in programs is now deprecated in favor oflocalprec
. Think about therealprecision
default as an interactive command for thegp
interpreter, best left out of GP programs. Indeed, the above rules imply that mixing both constructs yields surprising results:? \p38 ? localprec(19); default(realprecision,100); Pi %1 = 3.141592653589793239 ? \p realprecision = 115 significant digits (100 digits displayed)
Indeed,
realprecision
itself is ignored withinlocalprec
scope, soPi
is computed to a low accuracy. And when we leavelocalprec
scope,realprecision
only regains precedence, it is not “restored” to the original value.
-
mathilbert
(n)¶ being a
long
, creates the Hilbert matrixof order, i.e. the matrix whose coefficient (
,:math:
) is
.
-
matid
(n)¶ Creates the
identity matrix.
-
matpascal
(n, q=None)¶ Creates as a matrix the lower triangular Pascal triangle of order
(i.e. with binomial coefficients up to
). If
is given, compute the
-Pascal triangle (i.e. using
-binomial coefficients).
-
numtoperm
(n, k)¶ Generates the
-th permutation (as a row vector of length
) of the numbers
to
. The number
is taken modulo
, i.e. inverse function of
permtonum
. The numbering used is the standard lexicographic ordering, starting at.
-
partitions
(k, a=None, n=None)¶ Returns the vector of partitions of the integer
as a sum of positive integers (parts); for
, it returns the empty set
[]
, and forthe trivial partition (no parts). A partition is given by a
t_VECSMALL
, where parts are sorted in nondecreasing order:? partitions(3) %1 = [Vecsmall([3]), Vecsmall([1, 2]), Vecsmall([1, 1, 1])]
correspond to
,
and
. The number of (unrestricted) partitions of
is given by
numbpart
:? #partitions(50) %1 = 204226 ? numbpart(50) %2 = 204226
Optional parameters
and
are as follows:
(resp.
) restricts partitions to length less than
(resp. length between
and
), where the length is the number of nonzero entries.
(resp.
) restricts the parts to integers less than
(resp. between
and
).
? partitions(4, 2) \\ parts bounded by 2 %1 = [Vecsmall([2, 2]), Vecsmall([1, 1, 2]), Vecsmall([1, 1, 1, 1])] ? partitions(4,, 2) \\ at most 2 parts %2 = [Vecsmall([4]), Vecsmall([1, 3]), Vecsmall([2, 2])] ? partitions(4,[0,3], 2) \\ at most 2 parts %3 = [Vecsmall([4]), Vecsmall([1, 3]), Vecsmall([2, 2])]
By default, parts are positive and we remove zero entries unless
, in which case
is ignored and
is of constant length
:
? partitions(4, [0,3]) \\ parts between 0 and 3 %1 = [Vecsmall([0, 0, 1, 3]), Vecsmall([0, 0, 2, 2]),\ Vecsmall([0, 1, 1, 2]), Vecsmall([1, 1, 1, 1])]
-
polchebyshev
(n, flag=1, a=None)¶ Returns the
Chebyshev polynomial of the first kind
(
) or the second kind
(
), evaluated at
(
'x
by default). Both series of polynomials satisfy the 3-term relationand are determined by the initial conditions
,
,
. In fact
and, for all complex numbers
, we have
and
. If
, then these polynomials have degree
. For
,
is equal to
and
is equal to
. In particular,
.
-
polcyclo
(n, a=None)¶ -th cyclotomic polynomial, evaluated at
(
'x
by default). The integermust be positive.
Algorithm used: reduce to the case where
is squarefree; to compute the cyclotomic polynomial, use
; to compute it evaluated, use
. In the evaluated case, the algorithm assumes that
is either
or invertible, for all
. If this is not the case (the base ring has zero divisors), use
subst(polcyclo(n),x,a)
.
-
polhermite
(n, a=None)¶ Hermite polynomial
evaluated at
(
'x
by default), i.e.
-
pollegendre
(n, a=None)¶ Legendre polynomial evaluated at
(
'x
by default).
-
polmodular
(L, inv=0, x=None, y=None, compute_derivs=0)¶ Return the modular polynomial of level
in variables
and
for the
function. If
is given as
Mod(j, p)
or an elementof a prime finite field, then return the modular polynomial of level
evaluated at
modulo
. If
is from a finite field and compute_derivs is non-zero, then return a triple where the last two elements are the first and second derivatives of the modular polynomial evaluated at
.
? polmodular(3) %1 = x^4 + (-y^3 + 2232*y^2 - 1069956*y + 36864000)*x^3 + [...]
-
polsubcyclo
(n, d, v=None)¶ Gives polynomials (in variable
) defining the sub-Abelian extensions of degree
of the cyclotomic field
, where
.
If there is exactly one such extension the output is a polynomial, else it is a vector of polynomials, possibly empty. To get a vector in all cases, use
concat([], polsubcyclo(n,d))
.The function
galoissubcyclo
allows to specify exactly which sub-Abelian extension should be computed.
-
poltchebi
(n, v=None)¶ Deprecated alias for
polchebyshev
-
polylog
(m, x, flag=0, precision=0)¶ One of the different polylogarithms, depending on flag:
If
or is omitted:
polylogarithm of
, i.e. analytic continuation of the power series
(
). Uses the functional equation linking the values at
and
to restrict to the case
, then the power series when
, and the power series expansion in
otherwise.
Using
, computes a modified
polylogarithm of
. We use Zagier’s notations; let
denote
or
depending on whether
is odd or even:
If
: compute
, defined for
by
If
: compute
, defined for
by
If
: compute
, defined for
by
These three functions satisfy the functional equation
.
-
polzagier
(n, m)¶ Creates Zagier’s polynomial
used in the functions
sumalt
andsumpos
(with), see “Convergence acceleration of alternating series”, Cohen et al., Experiment. Math., vol. 9, 2000, pp. 3–12.
If
or
,
. We have
is
, where
is the Legendre polynomial of the second kind. For
,
is the
-th difference with step
of the sequence
; in this case, it satisfies
-
prime
(n)¶ The
prime number
? prime(10^9) %1 = 22801763489
Uses checkpointing and a naive
algorithm.
-
read
(filename=None)¶ Reads in the file filename (subject to string expansion). If filename is omitted, re-reads the last file that was fed into
gp
. The return value is the result of the last expression evaluated.If a GP
binary file
is read using this command (seewritebin
(in the PARI manual)), the file is loaded and the last object in the file is returned.In case the file you read in contains an
allocatemem
statement (to be generally avoided), you should leaveread
instructions by themselves, and not part of larger instruction sequences.
-
readstr
(filename=None)¶ Reads in the file filename and return a vector of GP strings, each component containing one line from the file. If filename is omitted, re-reads the last file that was fed into
gp
.
-
readvec
(filename=None)¶ Reads in the file filename (subject to string expansion). If filename is omitted, re-reads the last file that was fed into
gp
. The return value is a vector whose components are the evaluation of all sequences of instructions contained in the file. For instance, if file contains1 2 3
then we will get:
? \r a %1 = 1 %2 = 2 %3 = 3 ? read(a) %4 = 3 ? readvec(a) %5 = [1, 2, 3]
In general a sequence is just a single line, but as usual braces and
\
may be used to enter multiline sequences.
-
stirling
(n, k, flag=1)¶ Stirling number of the first kind
(
, default) or of the second kind
(flag = 2), where
,
are non-negative integers. The former is
times the number of permutations of
symbols with exactly
cycles; the latter is the number of ways of partitioning a set of
elements into
non-empty subsets. Note that if all
are needed, it is much faster to compute
Similarly, if a large number of
are needed for the same
, one should use
(Should be implemented using a divide and conquer product.) Here are simple variants for
fixed:
/* list of s(n,k), k = 1..n */ vecstirling(n) = Vec( factorback(vector(n-1,i,1-i*'x)) ) /* list of S(n,k), k = 1..n */ vecstirling2(n) = { my(Q = x^(n-1), t); vector(n, i, t = divrem(Q, x-i); Q=t[1]; simplify(t[2])); }
-
system
(str)¶ str is a string representing a system command. This command is executed, its output written to the standard output (this won’t get into your logfile), and control returns to the PARI system. This simply calls the C
system
command.
-
varhigher
(name, v=None)¶ Return a variable name whose priority is higher than the priority of
(of all existing variables if
is omitted). This is a counterpart to
varlower
.? Pol([x,x], t) *** at top-level: Pol([x,x],t) *** ^------------ *** Pol: incorrect priority in gtopoly: variable x <= t ? t = varhigher("t", x); ? Pol([x,x], t) %3 = x*t + x
This routine is useful since new GP variables directly created by the interpreter always have lower priority than existing GP variables. When some basic objects already exist in a variable that is incompatible with some function requirement, you can now create a new variable with a suitable priority instead of changing variables in existing objects:
? K = nfinit(x^2+1); ? rnfequation(K,y^2-2) *** at top-level: rnfequation(K,y^2-2) *** ^-------------------- *** rnfequation: incorrect priority in rnfequation: variable y >= x ? y = varhigher("y", x); ? rnfequation(K, y^2-2) %3 = y^4 - 2*y^2 + 9
Caution 1. The name is an arbitrary character string, only used for display purposes and need not be related to the GP variable holding the result, nor to be a valid variable name. In particular the name can not be used to retrieve the variable, it is not even present in the parser’s hash tables.
? x = varhigher("#"); ? x^2 %2 = #^2
Caution 2. There are a limited number of variables and if no existing variable with the given display name has the requested priority, the call to
varhigher
uses up one such slot. Do not create new variables in this way unless it’s absolutely necessary, reuse existing names instead and choose sensible priority requirements: if you only need a variable with higher priority than, state so rather than creating a new variable with highest priority.
\\ quickly use up all variables ? n = 0; while(1,varhigher("tmp"); n++) *** at top-level: n=0;while(1,varhigher("tmp");n++) *** ^------------------- *** varhigher: no more variables available. *** Break loop: type 'break' to go back to GP prompt break> n 65510 \\ infinite loop: here we reuse the same 'tmp' ? n = 0; while(1,varhigher("tmp", x); n++)
-
varlower
(name, v=None)¶ Return a variable name whose priority is lower than the priority of
(of all existing variables if
is omitted). This is a counterpart to
varhigher
.New GP variables directly created by the interpreter always have lower priority than existing GP variables, but it is not easy to check whether an identifier is currently unused, so that the corresponding variable has the expected priority when it’s created! Thus, depending on the session history, the same command may fail or succeed:
? t; z; \\ now t > z ? rnfequation(t^2+1,z^2-t) *** at top-level: rnfequation(t^2+1,z^ *** ^-------------------- *** rnfequation: incorrect priority in rnfequation: variable t >= t
Restart and retry:
? z; t; \\ now z > t ? rnfequation(t^2+1,z^2-t) %2 = z^4 + 1
It is quite annoying for package authors, when trying to define a base ring, to notice that the package may fail for some users depending on their session history. The safe way to do this is as follows:
? z; t; \\ In new session: now z > t ... ? t = varlower("t", 'z); ? rnfequation(t^2+1,z^2-2) %2 = z^4 - 2*z^2 + 9 ? variable() %3 = [x, y, z, t]
? t; z; \\ In new session: now t > z ... ? t = varlower("t", 'z); \\ create a new variable, still printed "t" ? rnfequation(t^2+1,z^2-2) %2 = z^4 - 2*z^2 + 9 ? variable() %3 = [x, y, t, z, t]
Now both constructions succeed. Note that in the first case,
varlower
is essentially a no-op, the existing variablehas correct priority. While in the second case, two different variables are displayed as
t
, one with higher priority than(created in the first line) and another one with lower priority (created by
varlower
).Caution 1. The name is an arbitrary character string, only used for display purposes and need not be related to the GP variable holding the result, nor to be a valid variable name. In particular the name can not be used to retrieve the variable, it is not even present in the parser’s hash tables.
? x = varlower("#"); ? x^2 %2 = #^2
Caution 2. There are a limited number of variables and if no existing variable with the given display name has the requested priority, the call to
varlower
uses up one such slot. Do not create new variables in this way unless it’s absolutely necessary, reuse existing names instead and choose sensible priority requirements: if you only need a variable with higher priority than, state so rather than creating a new variable with highest priority.
\\ quickly use up all variables ? n = 0; while(1,varlower("x"); n++) *** at top-level: n=0;while(1,varlower("x");n++) *** ^------------------- *** varlower: no more variables available. *** Break loop: type 'break' to go back to GP prompt break> n 65510 \\ infinite loop: here we reuse the same 'tmp' ? n = 0; while(1,varlower("tmp", x); n++)
-
version
()¶ Returns the current version number as a
t_VEC
with three integer components (major version number, minor version number and patchlevel); if your sources were obtained through our version control system, this will be followed by further more precise arguments, including e.g. agit
commit hash.This function is present in all versions of PARI following releases 2.3.4 (stable) and 2.4.3 (testing).
Unless you are working with multiple development versions, you probably only care about the 3 first numeric components. In any case, the
lex
function offers a clever way to check against a particular version number, since it will compare each successive vector entry, numerically or as strings, and will not mind if the vectors it compares have different lengths:if (lex(version(), [2,3,5]) >= 0, \\ code to be executed if we are running 2.3.5 or more recent. , \\ compatibility code );
On a number of different machines,
version()
could return either of%1 = [2, 3, 4] \\ released version, stable branch %1 = [2, 4, 3] \\ released version, testing branch %1 = [2, 6, 1, 15174, ""505ab9b"] \\ development
In particular, if you are only working with released versions, the first line of the gp introductory message can be emulated by
[M,m,p] = version(); printf("GP/PARI CALCULATOR Version %s.%s.%s", M,m,p);
If you are working with many development versions of PARI/GP, the 4th and/or 5th components can be profitably included in the name of your logfiles, for instance.
Technical note. For development versions obtained via
git
, the 4th and 5th components are liable to change eventually, but we document their current meaning for completeness. The 4th component counts the number of reachable commits in the branch (analogous tosvn
‘s revision number), and the 5th is thegit
commit hash. In particular,lex
comparison still orders correctly development versions with respect to each others or to released versions (provided we stay within a given branch, e.g.master
)!
-
-
sage.libs.pari.pari_instance.
prec_bits_to_dec
(prec_in_bits)¶ Convert from precision expressed in bits to precision expressed in decimal.
EXAMPLES:
sage: from sage.libs.pari.pari_instance import prec_bits_to_dec sage: prec_bits_to_dec(53) 15 sage: [(32*n, prec_bits_to_dec(32*n)) for n in range(1, 9)] [(32, 9), (64, 19), (96, 28), (128, 38), (160, 48), (192, 57), (224, 67), (256, 77)]
-
sage.libs.pari.pari_instance.
prec_bits_to_words
(prec_in_bits)¶ Convert from precision expressed in bits to pari real precision expressed in words. Note: this rounds up to the nearest word, adjusts for the two codewords of a pari real, and is architecture-dependent.
EXAMPLES:
sage: from sage.libs.pari.pari_instance import prec_bits_to_words sage: prec_bits_to_words(70) 5 # 32-bit 4 # 64-bit
sage: [(32*n, prec_bits_to_words(32*n)) for n in range(1, 9)] [(32, 3), (64, 4), (96, 5), (128, 6), (160, 7), (192, 8), (224, 9), (256, 10)] # 32-bit [(32, 3), (64, 3), (96, 4), (128, 4), (160, 5), (192, 5), (224, 6), (256, 6)] # 64-bit
-
sage.libs.pari.pari_instance.
prec_dec_to_bits
(prec_in_dec)¶ Convert from precision expressed in decimal to precision expressed in bits.
EXAMPLES:
sage: from sage.libs.pari.pari_instance import prec_dec_to_bits sage: prec_dec_to_bits(15) 49 sage: [(n, prec_dec_to_bits(n)) for n in range(10, 100, 10)] [(10, 33), (20, 66), (30, 99), (40, 132), (50, 166), (60, 199), (70, 232), (80, 265), (90, 298)]
-
sage.libs.pari.pari_instance.
prec_dec_to_words
(prec_in_dec)¶ Convert from precision expressed in decimal to precision expressed in words. Note: this rounds up to the nearest word, adjusts for the two codewords of a pari real, and is architecture-dependent.
EXAMPLES:
sage: from sage.libs.pari.pari_instance import prec_dec_to_words sage: prec_dec_to_words(38) 6 # 32-bit 4 # 64-bit sage: [(n, prec_dec_to_words(n)) for n in range(10, 90, 10)] [(10, 4), (20, 5), (30, 6), (40, 7), (50, 8), (60, 9), (70, 10), (80, 11)] # 32-bit [(10, 3), (20, 4), (30, 4), (40, 5), (50, 5), (60, 6), (70, 6), (80, 7)] # 64-bit
-
sage.libs.pari.pari_instance.
prec_words_to_bits
(prec_in_words)¶ Convert from pari real precision expressed in words to precision expressed in bits. Note: this adjusts for the two codewords of a pari real, and is architecture-dependent.
EXAMPLES:
sage: from sage.libs.pari.pari_instance import prec_words_to_bits sage: prec_words_to_bits(10) 256 # 32-bit 512 # 64-bit sage: [(n, prec_words_to_bits(n)) for n in range(3, 10)] [(3, 32), (4, 64), (5, 96), (6, 128), (7, 160), (8, 192), (9, 224)] # 32-bit [(3, 64), (4, 128), (5, 192), (6, 256), (7, 320), (8, 384), (9, 448)] # 64-bit
-
sage.libs.pari.pari_instance.
prec_words_to_dec
(prec_in_words)¶ Convert from precision expressed in words to precision expressed in decimal. Note: this adjusts for the two codewords of a pari real, and is architecture-dependent.
EXAMPLES:
sage: from sage.libs.pari.pari_instance import prec_words_to_dec sage: prec_words_to_dec(5) 28 # 32-bit 57 # 64-bit sage: [(n, prec_words_to_dec(n)) for n in range(3, 10)] [(3, 9), (4, 19), (5, 28), (6, 38), (7, 48), (8, 57), (9, 67)] # 32-bit [(3, 19), (4, 38), (5, 57), (6, 77), (7, 96), (8, 115), (9, 134)] # 64-bit