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)')
85.6215190762676
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 operations cause all arguments to be converted to PARI:
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 in the PARI interface¶
In the PARI interface, “real precision” refers to the precision of real numbers, so it is the floating-point precision. This is a non-trivial issue, since there are various interfaces for different things.
Internal representation and conversion between Sage and PARI¶
Real numbers in PARI have a precision associated to them, which is always a multiple of the CPU wordsize. So, it is a multiple of 32 of 64 bits. When converting from Sage to PARI, the precision is rounded up to the nearest multiple of the wordsize:
sage: x = 1.0
sage: x.precision()
53
sage: pari(x)
1.00000000000000
sage: pari(x).bitprecision()
64
With a higher precision:
sage: x = RealField(100).pi()
sage: x.precision()
100
sage: pari(x).bitprecision()
128
When converting back to Sage, the precision from PARI is taken:
sage: x = RealField(100).pi()
sage: y = pari(x).sage()
sage: y
3.1415926535897932384626433832793333156
sage: parent(y)
Real Field with 128 bits of precision
So pari(x).sage()
is definitely not equal to x
since it has
28 bogus bits.
Therefore, some care must be taken when juggling reals back and forth
between Sage and PARI. The correct way of avoiding this is to convert
pari(x).sage()
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). 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.sage(); 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
Output precision for printing¶
Even though PARI reals have a precision, not all significant bits are
printed by default. The maximum number of digits when printing a PARI
real can be set using the methods
PariInstance.set_real_precision_bits()
or
PariInstance.set_real_precision()
.
We create a very precise approximation of pi and see how it is printed in PARI:
sage: pi = pari(RealField(1000).pi())
The default precision is 15 digits:
sage: pi
3.14159265358979
With a different precision:
sage: _ = pari.set_real_precision(50)
sage: pi
3.1415926535897932384626433832795028841971693993751
Back to the default:
sage: _ = pari.set_real_precision(15)
sage: pi
3.14159265358979
Input precision for function calls¶
When we talk about precision for PARI functions, we need to distinguish three kinds of calls:
- Using the string interface, for example
pari("sin(1)")
. - Using the library interface with exact inputs, for example
pari(1).sin()
. - Using the library interface with inexact inputs, for example
pari(1.0).sin()
.
In the first case, the relevant precision is the one set by the methods
PariInstance.set_real_precision_bits()
or
PariInstance.set_real_precision()
:
sage: pari.set_real_precision_bits(150)
sage: pari("sin(1)")
0.841470984807896506652502321630298999622563061
sage: pari.set_real_precision_bits(53)
sage: pari("sin(1)")
0.841470984807897
In the second case, the precision can be given as the argument
precision
in the function call, with a default of 53 bits.
The real precision set by
PariInstance.set_real_precision_bits()
or
PariInstance.set_real_precision()
is irrelevant.
In these examples, we convert to Sage to ensure that PARI’s real precision is not used when printing the numbers. As explained before, this artificically increases the precision to a multiple of the wordsize.
sage: s = pari(1).sin(precision=180).sage(); print(s); print(parent(s))
0.841470984807896506652502321630298999622563060798371065673
Real Field with 192 bits of precision
sage: s = pari(1).sin(precision=40).sage(); print(s); print(parent(s))
0.841470984807896507
Real Field with 64 bits of precision
sage: s = pari(1).sin().sage(); print(s); print(parent(s))
0.841470984807896507
Real Field with 64 bits of precision
In the third case, the precision is determined only by the inexact
inputs and the precision
argument is ignored:
sage: pari(1.0).sin(precision=180).sage()
0.841470984807896507
sage: pari(1.0).sin(precision=40).sage()
0.841470984807896507
sage: pari(RealField(100).one()).sin().sage()
0.84147098480789650665250232163029899962
Elliptic curve functions¶
An elliptic curve given with exact \(a\)-invariants is considered an exact object. Therefore, you should set the precision for each method call individually:
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
TESTS:
Check that output from PARI’s print command is actually seen by Sage (trac ticket #9636):
sage: pari('print("test")')
test
Check that default()
works properly:
sage: pari.default("debug")
0
sage: pari.default("debug", 3)
sage: pari(2^67+1).factor()
IFAC: cracking composite
49191317529892137643
IFAC: factor 6713103182899
is prime
IFAC: factor 7327657
is prime
IFAC: prime 7327657
appears with exponent = 1
IFAC: prime 6713103182899
appears with exponent = 1
IFAC: found 2 large prime (power) factors.
[3, 1; 7327657, 1; 6713103182899, 1]
sage: pari.default("debug", 0)
sage: pari(2^67+1).factor()
[3, 1; 7327657, 1; 6713103182899, 1]
-
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)
For more information about how precision works in the PARI interface, see
sage.libs.pari.pari_instance
.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.
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.-
List
(x=None)¶ Create an empty list or convert \(x\) to a list.
EXAMPLES:
sage: pari.List(range(5)) List([0, 1, 2, 3, 4]) sage: L = pari.List() sage: L List([]) sage: L.listput(42, 1) 42 sage: L List([42]) sage: L.listinsert(24, 1) 24 sage: L List([24, 42])
-
PARI_ONE
¶
-
PARI_TWO
¶
-
PARI_ZERO
¶
-
allocatemem
(s=0, sizemax=0, silent=False)¶ Change the PARI stack space to the given size
s
(or double the current size ifs
is \(0\)) and change the maximum stack size tosizemax
.PARI tries to use only its current stack (the size which is set by
s
), but it will increase its stack if needed up to the maximum size which is set bysizemax
.The PARI stack is never automatically shrunk. You can use the command
pari.allocatemem(10^6)
to reset the size to \(10^6\), 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.
INPUT:
s
- an integer (default: 0). A non-zero argument is the size in bytes of the new PARI stack. If \(s\) is zero, double the current stack size.sizemax
- an integer (default: 0). A non-zero argument is the maximum size in bytes of the PARI stack. Ifsizemax
is 0, the maximum of the current maximum ands
is taken.
EXAMPLES:
sage: pari.allocatemem(10^7) PARI stack size set to 10000000 bytes, maximum size set to 67108864 sage: pari.allocatemem() # Double the current size PARI stack size set to 20000000 bytes, maximum size set to 67108864 sage: pari.stacksize() 20000000 sage: pari.allocatemem(10^6) PARI stack size set to 1000000 bytes, maximum size set to 67108864
The following computation will automatically increase the PARI stack size:
sage: a = pari('2^100000000') *** _^s: Warning: increasing stack size to...
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
Setting a small maximum size makes this fail:
sage: pari.allocatemem(10^6, 2^22) PARI stack size set to 1000000 bytes, maximum size set to 4194304 sage: a = pari('2^100000000') Traceback (most recent call last): ... PariError: _^s: the PARI stack overflows (current size: 1000000; maximum size: 4194304) You can use pari.allocatemem() to change the stack size and try again
TESTS:
Do the same without using the string interface and starting from a very small stack size:
sage: pari.allocatemem(1, 2^26) PARI stack size set to 1024 bytes, maximum size set to 67108864 sage: a = pari(2)^100000000 *** Warning: increasing stack size to... sage: pari.stacksize() 16777216
We do not allow
sizemax
less thans
:sage: pari.allocatemem(10^7, 10^6) Traceback (most recent call last): ... ValueError: the maximum size (10000000) should be at least the stack size (1000000)
-
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
-
double_to_gen
(x)¶
-
euler
(precision=0)¶ Euler’s constant \(\gamma = 0.57721...\). Note that
Euler
is one of the few reserved names which cannot be used for user variables.
-
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 \(P\) be a polynomial with integer coefficients. Determines the reduction of the (proper, smooth) genus 2 curve \(C/\QQ\), defined by the hyperelliptic equation \(y^2 = P\). The special syntax
genus2red([P,Q])
is also allowed, where the polynomials \(P\) and \(Q\) have integer coefficients, to represent the model \(y^2 + Q(x)y = P(x)\).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", []]]]]
-
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.
See also
get_real_precision_bits()
to get the precision in bits.EXAMPLES:
sage: pari.get_real_precision() 15
-
get_real_precision_bits
()¶ Return the current PARI default real precision in bits.
This is used both for creation of new objects from strings and for printing. It determines the number of digits in which real numbers 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.
See also
get_real_precision()
to get the precision in decimal digits.EXAMPLES:
sage: pari.get_real_precision_bits() 53
-
get_series_precision
()¶
-
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
(*args, **kwds)¶ Deprecated: Use
prime()
instead. See trac ticket #20216 for details.
-
one
()¶ EXAMPLES:
sage: pari.one() 1
-
pari_version
()¶
-
pi
(precision=0)¶ The constant \(\pi\) (\(3.14159...\)). Note that
Pi
is one of the few reserved names which cannot be used for user variables.
-
polchebyshev
(n, v=None)¶ Chebyshev polynomial of the first kind of degree \(n\), in the 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_eval
(*args, **kwds)¶ Deprecated: Use
polcyclo()
instead. See trac ticket #20217 for details.
-
polsubcyclo
(n, d, v=None)¶ polsubcyclo(n, d, v=x): return the pari list of polynomial(s) defining the sub-abelian extensions of degree \(d\) of the cyclotomic field \(\QQ(\zeta_n)\), where \(d\) divides \(\phi(n)\).
EXAMPLES:
sage: pari.polsubcyclo(8, 4) [x^4 + 1] sage: pari.polsubcyclo(8, 2, 'z') [z^2 + 2, z^2 - 2, z^2 + 1] 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
(*args, **kwds)¶ Deprecated: Use
primes()
instead. See trac ticket #20216 for details.
-
primes
(n=None, end=None)¶ Return a pari vector containing the first \(n\) primes, the primes in the interval \([n, end]\), or the primes up to \(end\).
INPUT:
Either
n
– integer
or
n
– list or tuple \([a, b]\) defining an interval of primes
or
n, end
– start and end point of an interval of primes
or
end
– end point for the list of primes
OUTPUT: a PARI list of prime numbers
EXAMPLES:
sage: pari.primes(3) [2, 3, 5] sage: pari.primes(10) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] sage: pari.primes(20) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71] sage: len(pari.primes(1000)) 1000 sage: pari.primes(11,29) [11, 13, 17, 19, 23, 29] sage: pari.primes((11,29)) [11, 13, 17, 19, 23, 29] sage: pari.primes(end=29) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] sage: pari.primes(10^30, 10^30 + 100) [1000000000000000000000000000057, 1000000000000000000000000000099]
TESTS:
sage: pari.primes(0) [] sage: pari.primes(-1) [] sage: pari.primes(end=1) [] sage: pari.primes(end=-1) [] sage: pari.primes(3,2) []
-
primes_up_to_n
(n)¶
-
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.
See also
set_real_precision_bits()
to set the precision in bits.EXAMPLES:
sage: pari.set_real_precision(60) 15 sage: pari('1.2') 1.20000000000000000000000000000000000000000000000000000000000 sage: pari.set_real_precision(15) 60
-
set_real_precision_bits
(n)¶ Sets the PARI default real precision in bits.
This is used both for creation of new objects from strings and for printing. It determines the number of digits in which real numbers 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.
See also
set_real_precision()
to set the precision in decimal digits.EXAMPLES:
sage: pari.set_real_precision_bits(200) sage: pari('1.2') 1.20000000000000000000000000000000000000000000000000000000000 sage: pari.set_real_precision_bits(53)
-
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
()¶ Return the current size of the PARI stack, which is \(10^6\) by default. However, the stack size is automatically doubled when needed up to some maximum.
See also
stacksizemax()
to get the maximum stack sizeallocatemem()
to change the current or maximum stack size
EXAMPLES:
sage: pari.stacksize() 16000000 sage: pari.allocatemem(2^18, silent=True) sage: pari.stacksize() 262144
-
stacksizemax
()¶ Return the maximum size of the PARI stack, which is determined at startup in terms of available memory. Usually, the PARI stack size is (much) smaller than this maximum but the stack will be increased up to this maximum if needed.
See also
stacksize()
to get the current stack sizeallocatemem()
to change the current or maximum stack size
EXAMPLES:
sage: pari.allocatemem(2^18, 2^26, silent=True) sage: pari.stacksizemax() 67108864
-
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)
-
zero
()¶ EXAMPLES:
sage: pari.zero() 0
-
class
sage.libs.pari.pari_instance.
PariInstance_auto
¶ Bases:
object
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 \(G = \sum_{n >= 0}((-1)^n)/((2n+1)^2) = 0.91596...\). Note that
Catalan
is one of the few reserved names which cannot be used for user variables.
-
Euler
(precision=0)¶ Euler’s constant \(\gamma = 0.57721...\). Note that
Euler
is one of the few reserved names which cannot be used for user variables.
-
I
()¶ The complex number \(\sqrt{-1}\).
-
Pi
(precision=0)¶ The constant \(\pi\) (\(3.14159...\)). 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 \(f\), 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 \(B_x\), where \(B_0 = 1\), \(B_1 = -1/2\), \(B_2 = 1/6\),..., expressed as a rational number. The argument \(x\) should be of type integer.
-
bernpol
(n, v=None)¶ Bernoulli polynomial \(B_n\) in variable \(v\).
? bernpol(1) %1 = x - 1/2 ? bernpol(3) %2 = x^3 - 3/2*x^2 + 1/2*x
-
bernreal
(x, precision=0)¶ Bernoulli number \(B_x\), as
bernfrac
, but \(B_x\) is returned as a real number (with the current precision).
-
bernvec
(x)¶ This routine is obsolete, kept for backward compatibility only.
-
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)¶ Given a prime \(N < 500\), return a vector \([P,t]\) where \(P(x,y)\) is a modular equation of level \(N\), i.e. a bivariate polynomial with integer coefficients; \(t\) indicates the type of this equation: either canonical (\(t = 0\)) or Atkin (\(t = 1\)). This function requires the
seadata
package and its only use is to give access to the package contents. Seepolmodular
for a more general and more flexible function.Let \(j\) be the \(j\)-invariant function. The polynomial \(P\) satisfies the functional equation,
\[P(f,j) = P(f \| W_N, j \| W_N) = 0\]for some modular function \(f = f_N\) (hand-picked for each fixed \(N\) to minimize its size, see below), where \(W_N(\tau) = -1 / (N \tau)\) is the Atkin-Lehner involution. These two equations allow to compute the values of the classical modular polynomial \(\Phi_N\), such that \(\Phi_N(j(\tau), j(N\tau)) = 0\), while being much smaller than the latter. More precisely, we have \(j(W_N(\tau)) = j(N \tau)\); the function \(f\) is invariant under \(\Gamma_0(N)\) and also satisfies
- for Atkin type: \(f \| W_N = f\);
- for canonical type: let \(s = 12/\mathrm{gcd}(12,N-1)\), then \(f \| W_N = N^s / f\). In this case, \(f\) has a simple definition: \(f(\tau) = N^s (\eta(N \tau) / \eta(\tau) )^{2 s}\), where \(\eta\) is Dedekind’s eta function.
The following GP function returns values of the classical modular polynomial by eliminating \(f_N(\tau)\) in the above functional equation, for \(N <= 31\) or \(N\in{41,47,59,71}\).
classicaleqn(N, X='X, Y='Y)= { my([P,t] = ellmodulareqn(N), Q, d); if (poldegree(P,'y) > 2, error("level unavailable in classicaleqn")); if (t == 0, \\ Canonical 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); , \\ Atkin Q = subst(P,'y,Y); d = (X-Y)^(N+1)); polresultant(subst(P,'y,X), Q) / d; }
-
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 \(x\). The expression \(x!\) gives a result which is an integer, while \(factorial(x)\) gives a real number.
-
fibonacci
(x)¶ \(x-th\) 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 \(s = 1\) (default) and totally complex if \(s = 2\). The output is a vector [
pol
,den
] wherepol
is the polynomial of degree \(a\)den
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 \(b\) and \(s\) are omitted, return the number of isomorphism classes of groups of order \(a\).
-
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 \(top-avma\), 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 functionGEN addii(GEN x, GEN y)
adds two PARI integers, and is not directly accessible under
gp
(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
(usuallylong 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, \(-1\) 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 \(n\)-point Gauss-Legendre integration of a smooth function \(f\) lon a compact interval \([a,b]\) at current
realprecision
. If \(n\) is omitted, make a default choice \(n ~ realprecision\), suitable for analytic functions on \([-1,1]\). The error is bounded by\[((b-a)^{2n+1} (n!)^4)/((2n+1)[(2n)!]^3) f^{(2n)} (\xi) , a < \xi < b\]so, if the interval length increases, \(n\) 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 \(n\) too much, see
intnumgauss
.
-
kill
(sym)¶ Restores the symbol
sym
to its “undefined” status, and deletes any help messages attached 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 attached 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 attached 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);
-
localbitprec
(p)¶ Set the real precision to \(p\) bits in the dynamic scope. All computations are performed as if
realbitprecision
was \(p\): transcendental constants (e.g.Pi
) and conversions from exact to floating point inexact data use \(p\) bits, as well as iterative routines implicitly using a floating point accuracy as a termination criterion (e.g.solve
orintnum
). Butrealbitprecision
itself is unaffected and is “unmasked” when we exit the dynamic (not lexical) scope. In effect, this is similar tomy(bit = default(realbitprecision)); default(realbitprecision,p); ... default(realbitprecision, bit);
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,
realbitprecision
will not be restored as intended.Such
localbitprec
statements can be nested, the innermost one taking precedence as expected. Beware thatlocalbitprec
follows the semantic oflocal
, notmy
: a subroutine called fromlocalbitprec
scope uses the local accuracy:? f()=bitprecision(1.0); ? f() %2 = 128 ? localbitprec(1000); f() %3 = 1024
Note that the bit precision of data (
1.0
in the above example) increases by steps of 64 (32 on a 32-bit machine) so we get \(1024\) instead of the expected \(1000\);localbitprec
bounds the relative error exactly as specified in functions that support that granularity (e.g.lfun
), and rounded to the next multiple of 64 (resp. 32) everywhere else.Warning. Changing
realbitprecision
orrealprecision
in programs is deprecated in favor oflocalbitprec
andlocalprec
. Think about therealprecision
andrealbitprecision
defaults as interactive commands 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,1000); Pi %1 = 3.141592653589793239 ? \p realprecision = 1001 significant digits (1000 digits displayed)
Indeed,
realprecision
itself is ignored withinlocalprec
scope, soPi
is computed to a low accuracy. And when we leave thelocalprec
scope,realprecision
only regains precedence, it is not “restored” to the original value.
-
localprec
(p)¶ Set the real precision to \(p\) in the dynamic scope. All computations are performed as if
realprecision
was \(p\): transcendental constants (e.g.Pi
) and conversions from exact to floating point inexact data use \(p\) decimal 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)¶ \(x\) being a
long
, creates the Hilbert matrixof order \(x\), i.e. the matrix whose coefficient (\(i\),:math:\(j\)) is \(1/ (i+j-1)\).
-
matid
(n)¶ Creates the \(n x n\) identity matrix.
-
matpascal
(n, q=None)¶ Creates as a matrix the lower triangular Pascal triangle of order \(x+1\) (i.e. with binomial coefficients up to \(x\)). If \(q\) is given, compute the \(q\)-Pascal triangle (i.e. using \(q\)-binomial coefficients).
-
numtoperm
(n, k)¶ Generates the \(k\)-th permutation (as a row vector of length \(n\)) of the numbers \(1\) to \(n\). The number \(k\) is taken modulo \(n!\), i.e. inverse function of
permtonum
. The numbering used is the standard lexicographic ordering, starting at \(0\).
-
oo
()¶ Returns an object meaning \(+ oo\), for use in functions such as
intnum
. It can be negated (-oo
represents \(- oo\)), and compared to real numbers (t_INT
,t_FRAC
,t_REAL
), with the expected meaning: \(+ oo\) is greater than any real number and \(- oo\) is smaller.
-
partitions
(k, a=None, n=None)¶ Returns the vector of partitions of the integer \(k\) as a sum of positive integers (parts); for \(k < 0\), it returns the empty set
[]
, and for \(k = 0\) the trivial partition (no parts). A partition is given by at_VECSMALL
, where parts are sorted in nondecreasing order:? partitions(3) %1 = [Vecsmall([3]), Vecsmall([1, 2]), Vecsmall([1, 1, 1])]
correspond to \(3\), \(1+2\) and \(1+1+1\). The number of (unrestricted) partitions of \(k\) is given by
numbpart
:? #partitions(50) %1 = 204226 ? numbpart(50) %2 = 204226
Optional parameters \(n\) and \(a\) are as follows:
- \(n = nmax\) (resp. \(n = [nmin,nmax]\)) restricts partitions to length less than \(nmax\) (resp. length between \(nmin\) and \(nmax\)), where the length is the number of nonzero entries.
- \(a = amax\) (resp. \(a = [amin,amax]\)) restricts the parts to integers less than \(amax\) (resp. between \(amin\) and \(amax\)).
? 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 \(amin <= 0\), in which case \(nmin\) is ignored and \(X\) is of constant length \(nmax\):
? 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 \(n-th\) Chebyshev polynomial of the first kind \(T_n\) (\(flag = 1\)) or the second kind \(U_n\) (\(flag = 2\)), evaluated at \(a\) (
'x
by default). Both series of polynomials satisfy the 3-term relation\[P_{n+1} = 2xP_n - P_{n-1},\]and are determined by the initial conditions \(U_0 = T_0 = 1\), \(T_1 = x\), \(U_1 = 2x\). In fact \(T_n' = n U_{n-1}\) and, for all complex numbers \(z\), we have \(T_n(\cos z) = \cos (nz)\) and \(U_{n-1}(\cos z) = \sin(nz)/\sin z\). If \(n >= 0\), then these polynomials have degree \(n\). For \(n < 0\), \(T_n\) is equal to \(T_{-n}\) and \(U_n\) is equal to \(-U_{-2-n}\). In particular, \(U_{-1} = 0\).
-
polcyclo
(n, a=None)¶ \(n\)-th cyclotomic polynomial, evaluated at \(a\) (
'x
by default). The integer \(n\) must be positive.Algorithm used: reduce to the case where \(n\) is squarefree; to compute the cyclotomic polynomial, use \(\Phi_{np}(x) = \Phi_n(x^p)/\Phi(x)\); to compute it evaluated, use \(\Phi_n(x) = \prod_{d \| n} (x^d-1)^{\mu(n/d)}\). In the evaluated case, the algorithm assumes that \(a^d - 1\) is either \(0\) or invertible, for all \(d \| n\). If this is not the case (the base ring has zero divisors), use
subst(polcyclo(n),x,a)
.
-
polhermite
(n, a=None)¶ \(n-th\) Hermite polynomial \(H_n\) evaluated at \(a\) (
'x
by default), i.e.\[H_n(x) = (-1)^n e^{x^2} (d^n)/(dx^n)e^{-x^2}.\]
-
pollegendre
(n, a=None)¶ \(n-th\) Legendre polynomial evaluated at \(a\) (
'x
by default).
-
polmodular
(L, inv=0, x=None, y=None, derivs=0)¶ Return the modular polynomial of prime level \(L\) in variables \(x\) and \(y\) for the modular function specified by
inv
. Ifinv
is 0 (the default), use the modular \(j\) function, ifinv
is 1 use the Weber-\(f\) function, and ifinv
is 5 use \(\gamma_2 = \sqrt[3]{j}\). Seepolclass
for the full list of invariants. If \(x\) is given asMod(j, p)
or an element \(j\) of a finite field (as at_FFELT
), then return the modular polynomial of level \(L\) evaluated at \(j\). If \(j\) is from a finite field andderivs
is non-zero, then return a triple where the last two elements are the first and second derivatives of the modular polynomial evaluated at \(j\).? polmodular(3) %1 = x^4 + (-y^3 + 2232*y^2 - 1069956*y + 36864000)*x^3 + ... ? polmodular(7, 1, , 'J) %2 = x^8 - J^7*x^7 + 7*J^4*x^4 - 8*J*x + J^8 ? polmodular(7, 5, 7*ffgen(19)^0, 'j) %3 = j^8 + 4*j^7 + 4*j^6 + 8*j^5 + j^4 + 12*j^2 + 18*j + 18 ? polmodular(7, 5, Mod(7,19), 'j) %4 = Mod(1, 19)*j^8 + Mod(4, 19)*j^7 + Mod(4, 19)*j^6 + ... ? u = ffgen(5)^0; T = polmodular(3,0,,'j)*u; ? polmodular(3, 0, u,'j,1) %6 = [j^4 + 3*j^2 + 4*j + 1, 3*j^2 + 2*j + 4, 3*j^3 + 4*j^2 + 4*j + 2] ? subst(T,x,u) %7 = j^4 + 3*j^2 + 4*j + 1 ? subst(T',x,u) %8 = 3*j^2 + 2*j + 4 ? subst(T'',x,u) %9 = 3*j^3 + 4*j^2 + 4*j + 2
-
polsubcyclo
(n, d, v=None)¶ Gives polynomials (in variable \(v\)) defining the sub-Abelian extensions of degree \(d\) of the cyclotomic field \(\mathbb{Q}(\zeta_n)\), where \(d \| \phi(n)\).
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 \(flag = 0\) or is omitted: \(m-th\) polylogarithm of \(x\), i.e. analytic continuation of the power series \(Li_m(x) = \sum_{n >= 1}x^n/n^m\) (\(x < 1\)). Uses the functional equation linking the values at \(x\) and \(1/x\) to restrict to the case \(\|x\| <= 1\), then the power series when \(\|x\|^2 <= 1/2\), and the power series expansion in \(\log(x)\) otherwise.
Using \(flag\), computes a modified \(m-th\) polylogarithm of \(x\). We use Zagier’s notations; let \(\Re_m\) denote \(\Re\) or \(\Im\) depending on whether \(m\) is odd or even:
If \(flag = 1\): compute \(~ D_m(x)\), defined for \(\|x\| <= 1\) by
\[\Re_m(\sum_{k = 0}^{m-1} ((-\log\|x\|)^k)/(k!)Li_{m-k}(x) +((-\log\|x\|)^{m-1})/(m!)\log\|1-x\|).\]If \(flag = 2\): compute \(D_m(x)\), defined for \(\|x\| <= 1\) by
\[\Re_m(\sum_{k = 0}^{m-1}((-\log\|x\|)^k)/(k!)Li_{m-k}(x) -(1)/(2)((-\log\|x\|)^m)/(m!)).\]If \(flag = 3\): compute \(P_m(x)\), defined for \(\|x\| <= 1\) by
\[\Re_m(\sum_{k = 0}^{m-1}(2^kB_k)/(k!)(\log\|x\|)^kLi_{m-k}(x) -(2^{m-1}B_m)/(m!)(\log\|x\|)^m).\]These three functions satisfy the functional equation \(f_m(1/x) = (-1)^{m-1}f_m(x)\).
-
polzagier
(n, m)¶ Creates Zagier’s polynomial \(P_n^{(m)}\) used in the functions
sumalt
andsumpos
(with \(flag = 1\)), see “Convergence acceleration of alternating series”, Cohen et al., Experiment. Math., vol. 9, 2000, pp. 3–12.If \(m < 0\) or \(m >= n\), \(P_n^{(m)} = 0\). We have \(P_n := P_n^{(0)}\) is \(T_n(2x-1)\), where \(T_n\) is the Legendre polynomial of the second kind. For \(n > m > 0\), \(P_n^{(m)}\) is the \(m\)-th difference with step \(2\) of the sequence \(n^{m+1}P_n\); in this case, it satisfies
\[2 P_n^{(m)}(sin^2 t) = (d^{m+1})/(dt^{m+1})(\sin(2t)^m \sin(2(n-m)t)).\]
-
prime
(n)¶ The \(n-th\) prime number
? prime(10^9) %1 = 22801763489
Uses checkpointing and a naive \(O(n)\) 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.
-
self
()¶ Return the calling function or closure as a
t_CLOSURE
object. This is useful for defining anonymous recursive functions.? (n->if(n==0,1,n*self()(n-1)))(5) %1 = 120
-
stirling
(n, k, flag=1)¶ Stirling number of the first kind \(s(n,k)\) (\(flag = 1\), default) or of the second kind \(S(n,k)\) (flag = 2), where \(n\), \(k\) are non-negative integers. The former is \((-1)^{n-k}\) times the number of permutations of \(n\) symbols with exactly \(k\) cycles; the latter is the number of ways of partitioning a set of \(n\) elements into \(k\) non-empty subsets. Note that if all \(s(n,k)\) are needed, it is much faster to compute
\[\sum_k s(n,k) x^k = x(x-1)...(x-n+1).\]Similarly, if a large number of \(S(n,k)\) are needed for the same \(k\), one should use
\[\sum_n S(n,k) x^n = (x^k)/((1-x)...(1-kx)).\](Should be implemented using a divide and conquer product.) Here are simple variants for \(n\) 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 \(v\) (of all existing variables if \(v\) 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 \(x\), 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 \(v\) (of all existing variables if \(v\) 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 variable \(t\) has correct priority. While in the second case, two different variables are displayed ast
, one with higher priority than \(z\) (created in the first line) and another one with lower priority (created byvarlower
).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 \(x\), 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.
default_bitprec
()¶ Return the default precision in bits.
EXAMPLES:
sage: from sage.libs.pari.pari_instance import default_bitprec sage: default_bitprec() 64
-
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) 50 sage: [(n, prec_dec_to_bits(n)) for n in range(10, 100, 10)] [(10, 34), (20, 67), (30, 100), (40, 133), (50, 167), (60, 200), (70, 233), (80, 266), (90, 299)]
-
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