[To main | To resources page ]
This document is maintained by
Karim Belabas, the current PARI/GP maintainer (feedback welcome). It
strives to be useful for all versions of PARI/GP but, for the time being, it
documents the behaviour of the latest versions in the stable and
testing branches. And
may hint at features introduced in the development version, available via
Subversion only.
*** isprime: Warning: IFAC: untested integer declared prime.Can I really trust the output of isprime()?
PARI/GP is a widely used computer algebra system designed for fast computations in number theory (factorizations, algebraic number theory, elliptic curves...), but also contains a large number of other useful functions to compute with mathematical entities such as matrices, polynomials, power series, algebraic numbers, etc., and a lot of transcendental functions.
PARI is a C library (C++ - compatible). Writing C (or C++)
programs, and linking them to libpari is the only way to access
the full set of PARI routines (you can install most of them
under gp, but this is not as portable). It will produce much faster code for
low-level applications.
gp is an interactive shell giving access to PARI functions, easier to use. Low-level scripts are typically 10 times slower than directly written C code. High-level scripts, whose running time is dominated by a few calls to PARI functions, incur no penalty.
GP is the name of gp's scripting language.
It is the GP-to-C compiler. It compiles GP scripts to the C language, easing the task of writing PARI programs. It can automagically compile them to object code and load the resulting functions in gp. Low-level gp2c-compiled scripts typically run 3 or 4 times faster. gp2c currently only understands a subset of the GP language, and some functionalities are restricted to Unix operating systems.
It is a library, which provide line-editing facilities (cursor movement, command history, completion, etc). All applications built using this library, for instance the bash shell, benefit from a consistent and convenient text-only user interface.
It is strongly advised to make sure readline is present on your system before building gp yourself (readline is included in the Windows binary we distribute). In fact, before running Configure so that the latter gets a chance to detect it. The gp shell will function without line editing, but typing will be tedious.
GNU MP, or GMP for short, is a multiprecision library, which in particular provides asymptotically fast routines for arithmetic operations between integers. PARI multiprecision kernel can use native implementations or make good use of GMP if it is available. We advise you to download and install GMP before Configuring PARI. Make sure you download version 4.2 or more recent, because of this problem.
See also How do I enable the GMP kernel.
Please use PARI/GP's Bug Tracking System. Others may see and comment on your request, and you will be individually notified whenever a solution is committed to the Subversion server.
If you have used PARI/GP in the preparation of a paper, then yes, we would appreciate that.
When an author formally cites the use of a computer algebra system or package in the bibliography of a paper, then system developers receive credit in citation indices. This practice encourages development of systems and packages and helps provide recognition for this work that is similar to that enjoyed by those who publish research articles. It also makes your experimental results reproducible.
See also SIGSAM's statement of policy and bibliography of computer algebra systems.
You may cite PARI in the following form (BibTeX format)
@preamble("\usepackage{url}")
@manual{PARI2,
organization = "{The PARI~Group}",
title = "{PARI/GP, version {\tt 2.3.4}}",
year = 2008,
address = "Bordeaux",
note = "available from \url{http://pari.math.u-bordeaux.fr/}"
}
(please change the version number to match the version you actually used), which would translate to something like
\bibitem{PARI2}
PARI/GP, version {\tt 2.3.4}, Bordeaux, 2008, \url{http://pari.math.u-bordeaux.fr/}.
PARI/GP is free software. It is released under GNU's General Public License (GPL for short) as published by the Free Software Fundation, either version 2 of the License, or (at your option) any later version.
The short answer to the second question is no.
We are quite happy with the GPL. The original license (for all the 1.xx versions) was proprietary, too casually worded, and led to various conflicts. Starting from version 2.0, we wanted a well known license for the whole PARI/GP package, that would leave contributors secure with the future use of their code and, as far as possible, not prevent anyone from helping out. It was soon agreed that PARI/GP should become free software.
The GPL was a natural choice. It is certainly well-known, and it satisfied every developer involved in the project, as well as their respective employers. For the latter, the alternatives would have been more restrictive, certainly not closer to the LGPL, BSD, or the Artistic license. Most free Computer Algebra Systems also use the GPL, presumably for related reasons.
There are two sites devoted to PARI:
The testing distribution includes two User's Manuals
(for the GP calculator and the libpari C library), a Tutorial, a too
short Developer's guide,
and a Reference Card. Online or downloadable versions are available at http://pari.math.u-bordeaux.fr/doc.html.
For more information about the algorithms used, check out Henri Cohen's two books (Graduate Texts in Mathematics 138 and 193, Springer).
No, English only. PARI is written by a small group of volunteers, and maintaining an adequate documentation in a single language has already proven to be a daunting task.
Please do not harass us with translation requests or lectures on linguistic imperialism. Some developers happen to be native French speakers, and are indeed aware of a number of disturbing linguistic issues. But this does not imply we have unlimited time to translate existing documents, not to mention creating and maintaining new ones.
PARI is free software, and you are welcome to undertake or sponsor a translation effort to the language of your choice. Please notify us if you do.
Current stable release is pari-2.3.4.
The last testing snapshot is pari-2.4.2.
Current development version, available via Subversion only is 2.4.3.
We simultaneously maintain two versions (or branches) of PARI:
stable = odd minor version number, e.g 2.3.x:
only modified to fix big problems, which are curable with minor changes.
Neither minor improvements, nor complicated changes get incorporated. The
primary goal is not to upset working scripts or code. Thus you can update
with confidence from one patchlevel to another, within the same version
branch, e.g. 2.3.x to 2.3.y. This branch usually lags far behind
testing: as a rule, it contains many more known bugs and will
lack some useful features.
testing = even minor version number, e.g 2.4.x: this is a development version, and nothing is guaranteed, although everything usually works. Large sections of code are modified, features appear or are cancelled, interfaces change according to the feedback we receive, many bugs are fixed and new ones sneak in. Eventually, the testing version coalesces into a stable one, and a new testing branch is started.
It depends on your needs; we describe a few plausible scenarios below. First, always use the most recent version in a given branch.
If stable can run your scripts without problems, and you do not
need any of the newer features, just use stable!
If you just want to give PARI/GP a try or have already run into problems with
stable, then try out testing. In general, it
contains fewer bugs, more features, and many functions may be orders of
magnitudes faster.
If you are teaching a class and need to program in GP, you probably
want to use testing: it is easier and cleaner to program, in
particular due to the availability of my (and closures, for
more advanced courses), and generally better documented.
If you want to write new programs with libpari, use
testing: the documentation is much more accurate and many
more useful functions are publicly available, usually with more consistent
interfaces.
If you run into problems with testing, please use
our Subversion server: your problem is probably already
solved. If not, please report your problem using our Bug Tracking System, then
update as soon as the problem is solved; usually at most one or two days
later. (People submitting a bug report are individually notified
whenever a solution to their problem is committed to the Subversion server.)
The version number has three components M.m.p, e.g 2.1.4:
stable versions, even ones to testing versions.
The idea is that patchlevel is bumped whenever a series of modifications
has been completed, and minor version number is increased when a
testing
version is deemed robust enough to become the new stable version.
A new testing branch is then created.
This has a number of implications in terms of interface stability.
When it is ready. We have no definite schedule: if you are bothered by a particular problem or want to check out new improvements, please use the development version available via Subversion.
Please do it yourself: send an email to
pari-users-unsubscribe@list.cr.yp.to
(resp. pari-dev, pari-announce). This should be enough to cancel your subscription; just make sure you use the same email account as when you originally subscribed.
Thousands of administrative requests are sent each month to the various mailing lists hosted by the list maintainer. Handling even a small fraction of these by hand is out of the question. The three PARI lists are fully automated, managed by the EZMLM software.
The PARI mailing lists only support the subscribe, unsubscribe and get.nnn commands.
The EZMLM manual describes other commands (like index or thread) and variants (get.aaa_bbb) that are only available to lists managed by EZMLM-IDX. This is not the case of the PARI lists.
Everything that goes to the list must be confirmed by qsecretary. You are not being singled out for mis-treatment.
To prevent unsolicited emails (spam) and virii to be propagated by the PARI lists, the mailing list manager requires an authenticated confirmation before any post is transmitted. You only have to REPLY to the secretary and your post will go through.
This is certainly an annoyance, but the list maintainer had to choose between two evils, and this one is a minor annoyance to posters, whereas spam was a major one for all subscribers.
[ Further details ] [ pymsgauth, a secure answering machine ]
Depending on the list, one of
pari-announce-list-owner@list.cr.yp.to
pari-dev-list-owner@list.cr.yp.to
pari-users-list-owner@list.cr.yp.to
Please report them directly to Math::Pari author and maintainer Ilya Zakharevich.
Please report them directly to the maintainer of the pariemacs package Olivier Ramaré.
Many web browsers silently expand anything the web server tells them is gzip'd. And will nevertheless offer you to save it with the .gz ending, now inappropriate. So you are probably looking at xxx under the name of xxx.gz.
Under Unix, try.
file xxx.gz
It should find out what the file format really is. The file sizes should also
be a tell-tale.
It is gone. We now use Subversion. We still provide a read-only CVS frontend to our Subversion repository for backward compatibility, but new users should use Subversion. If you were using CVS to fetch sources, it should still work transparently. But please consider upgrading to Subversion.
Subversion is an open-source revision control system. For the developers of PARI/GP, Subversion provides network-transparent source control. For ordinary users it provides a very convenient way to obtain patched versions (source code only) in between releases. Subversion clients are available for all major platforms: Unix, MacOS, Windows.
When a Subversion client connects to the central server for the first time, it copies the source code on your hard drive, exactly as if you had fetched and extracted a source release. PARI can then be compiled and installed in the usual way. Later connections upgrade your local copy using latest available sources, fetching only modified files. It is easy to restrict to a specific version branch, to released snapshots, or to revert your local repository to an older copy.
More details on Subversion setup for PARI.
First, there is no guarantee that bleeding edge sources from the Subversion server will indeed compile. You may be unlucky, and fetch a broken version. More probably, new code modules have been added, which were not compiled in. Just try to run Configure after updating, new files will then be taken into account by the newly generated Makefiles.
If you get an error message of the form command not
found when trying to Configure pari, '.' is not in your search path
for commands. You need to type
./Configure
in the toplevel directory.
You need to install GNU readline first, then run Configure. Check out the Configure output to make sure the readline library is detected. You can help Configure with the
--with-readline
--with-readline-include
--with-readlin-lib
flags if you installed your library in an exotic place, or you installed many and Configure picks up the wrong one.
Linux distributions have separate readline and
readline-devel packages. If only readline is
installed, Configure will complain: readline is there so that
precompiled binaries can function properly, you need both of them installed
to compile programs with readline support.
Configure may also complain about a missing
libncurses.so
in which case, install the ncurses-devel package. (Some
distributions let you install readline-devel without
ncurses-devel, which is a bug in their package dependency
handling.)
variant of this FAQ. It has been reported that in
SuSE distributions, one needs to install both xorg-x11-libs and
xorg-x11.
make clean
./Configure (with proper options)
make gp
Sometimes, a previous build (e.g. interrupted or after an update from
Subversion) interferes with subsequent compilation. make clean
deletes all cached information from the build directory.
This was specifically reported on Mac OS X, while trying to compile GNU readline. The exact message was:
make: *** No rule to make target `/Volumes/Macintosh', needed by `readline.o'.
Stop.
This was caused by a directory named "Macintosh HD" up the directory chain between the user's home directory and the root. In particular the environment variable $HOME contained spaces, something that readline-5 configuration cannot handle. Some versions of PARI could not handle them either, but this should no longer be the case.
If renaming the offending directories is not an option, a simple solution was to do the following
ln -s "$HOME/readline-5.0" /tmp/tmpreadlinethen cd to /tmp/tmpreadline and restart the standard procedure. The above creates a derivation (symbolic link) from /tmp/tmpreadline (doesn't contain spaces) to the directory we intended to work in. Adapt to your specific needs.
This yacc file must be processed with a recent
GNU bison, e.g. version 2.3 is fine. The distribution already
contains the compiled files parse.c and parse.h,
hence this problem is specific to Subversion: people downloading the
distribution do not need bison.
On the sparcv9 architecture, a 32bit binary is produced by default. To override this, you need to type in something like (GNU cc)
env CC="gcc -m64" ./Configure --kernel=none-none
or (Forte cc)
env CC="cc -xarch=v9" CFLAGS="-xarch=v9" ./Configure --kernel=none-none
You may want to specify kernel=none-gmp instead, but it is crucial to use the portable native kernel, and not the sparcv9 one, because the latter is 32bit-specific. (Since sparc64 does not support multiplication of 64bit operands in hardware, there was no point in porting the sparc32 PARI assembler kernel.)
You may run into further trouble if the linker picks up some 32-bit libraries instead of their 64-bit equivalent. This might happen for user-installed libraries like libgcc or libgmp. In this case, try to unset LD_LIBRARY_PATH (or reset it to a suitable value), and possibly set LD_RUN_PATH to the 64-bit directory.
This came up under Fedora Core 6 as Bug #621. Andrew van Herick reported the problem and worked around it by adjusting SELinux security policy.
Using Fedora's graphical interfaces, this is done as follows:
System->Administration->Security Level and Firewall.
Modify SELinux Policy->Memory Protection
enable the box labeled
Allow all unconfined executables to use libraries requiring text
relocation that are not labeled textrel_shlib_t
We do not try to support broken compilers. You may try to compile the faulty code module using a lower optimization level (e.g. -O instead of -O3). Or switching to a stable version of a reliable compiler, e.g. gcc-2.95.3.
gcc-2.96 is an unofficial version of gcc which was unfortunately shipped with many Linux distributions. These initial releases are badly broken on a number of architectures, ia64 for instance. Upgrade!
As in Bug #565. First, try to upgrade your version of XCode. In particular, check the build number for gcc (e.g. in gp headers, Apple Computer Inc. build xxx). Builds < 5300 are very buggy, build 5363 works fine.
Check the relevant 'diff' file mentioned at the end of make
bench output. If the only difference comes from a failed installation
of addii, it only means that install() is not
supported on your system; so you will not be able to use this feature.
Besides that, you can safely ignore the Warning.
On most platforms, make install builds and installs a
shared PARI library libpari.so (meaning that needed
functions are loaded from the library at runtime, not copied into the
executable at link time). By default this library goes to /usr/local/lib/.
Unfortunately, this is often not enough to make your system aware of the new
library.
Since the link command in your Makefile probably
specifies something like
$(CC) ... -L$(LIBDIR) -lpari -lm
with LIBDIR pointing at the directory containing
libpari.so, the link phase succeeds. Unfortunately, this does
not register in the binary where the shared library can be found, and
running the binary fails with the above error. You may check with ldd
my_prog what are the shared libraries needed by your program,
and which ones are found by the dynamic lib loader: this should fail
with a message similar to the above.
ldconfig to make your system aware
of the new library. (Requires root powers.)-Wl,-rpath "my exotic libpath" in the link command
of your Makefile. See
examples/Makefile for how to do it on your system. (This is the solution chosen for the gp binary.)LD_LIBRARY_PATH. This can be as easy as adding
export LD_LIBRARY_PATH=/usr/local/libin a
.bash_profile. Another possibility is to run
env LD_LIBRARY_PATH=/usr/local/lib my_proginstead of
my_prog.
This is possible only if the readline library is compiled in (the gp header gives this information on startup). In this case, you need to create a .inputrc file. You can create this in your home directory under Unix; or any place you want provided you then set the environment variable INPUTRC to the full path of the .inputrc file.
In this file, you can associate actions to sequences of keypresses, e.g.
"\e[5C": forward-word
"\e[5D": backward-word
associates Esc + [ + 5 + C, the sequence produced by hitting Control-RightArrow on my system, to the action forward-word: move the cursor to the right by one word. Complete details can be found in readline's manual.
To determine what sequence is associated to a given combination of key presses, you can hit Control-V, then press the relevant keys. gp will print the corresponding sequence, with Esc coded as ^[.
The function simplify is applied to gp's result, before storing
it in the history, but after the command is fully evaluated, in particular
after any variable assignment takes place. Here's an example:
? a = x+y-x
%1 = y
? variable(a)
%2 = x
? variable(%1)
%3 = y
a is still a polynomial of degree 0 in x, but %1 has been simplified. You can use simplify yourself whenever you want to ensure an object has the expected "simplest possible" type.
Use return(), without arguments.
For instance, the following happens to work
init() = f(x) = x
but this one does not
init() = f(x) = x; g(x) = 2 * x
In fact, calling the latter defines the function f(x) whose
definition is
x; g(x) = 2 * x
which is not what was intended. In fact, the GP grammar states that the end of a function definition is the end of the longest expression sequence following the prototype. And the above is unfortunately a valid expression sequence. The solution is simple:
init(a) = ( f(x) = x ); ( g(x) = 2 * x )
When a call to init is evaluated,
the code within parentheses is evaluated first, defining separately
f then
g.
Use indirect sorting. More generally, to sort a vector according to the
values of some hashing function f, you can use
f(x) = abs(x) SORT(v) = vecextract(v, vecsort(vector(#v, i, f(v[i])),, 1))
(In the case at hand, one can use abs(v) instead of building
explicitly a vector of values.)
It does. Read the following from a file (or copy-paste it)
fun() =
{
print("enter an integer a: ");
a = input();
print("a^2 is "a^2)
}
Then typing fun() from the gp prompt waits for some user input.
The reason why the simpler
print("enter an integer a: ");
a = input();
print("a^2 is "a^2)
does not work, is that input() reads expression from
the program standard input (stdin). As you are reading the
script into gp, the script itself is fed into stdin. So input
gets some data, specifically the return value of
print("a^2 is "a^2)
which is 0, then returns immediately. In the meantime, the expanded
"a^2 is"a^2 has been printed, as "a^2 is a^2"
since at this point the value of a has not been modified.
If you re-read the same commands into gp, it will now print
a^2 is 0.
Indeed, read does not work since it reads in the whole file, and
all intermediate values are discarded. The simplest solution is to write
all data into a vector or list and save the latter.
Another solution is to use the readvec function. It returns a
vector whose entries contain the results of all expression sequences from
the file. Let file contain
1+1
2+2
3+3
Then
? readvec(a)
%1 = [2, 4, 6]
A few amusing hacks: if you read in a file using \r
(not with read), all history entries will be filled,
and you can later access them individually. The following construction could
be helpful:
for (i = first, last,
x = eval( Str("%", i) ); \\ the i-th history entry
...
)
On a Unix system, you can use the following trick (courtesy of Markus Endres) to input the n-th line of a given file:
getline(n, file) = extern("sed -n " n "p" file);
On a Unix system, you can use the following beautiful hack (courtesy of
Bill Allombert) using named fifo files: create a fifo and start
reading from it. In another shell, write your data to it. For example:
sh1$ mkfifo fifo
sh1$ gp
? for(i=1, 10, print(i, read("fifo"))
sh2$ perl -ne 'system("echo '\''$_'\'' > fifo")' < data
The possibility of handling file objects in GP without resorting to such extremities is in the TODO list.
Use readvec.
You may use the following function, courtesy of Martin Larsen on
pari-user:
hextodec(s) =
{ local(v=Vec(s), a=10,b=11,c=12,d=13,e=14,f=15,
A=10,B=11,C=12,D=13,E=14,F=15, h);
for(i=1,#v,h = (h<<4) + eval(v[i])); h
}
? hextodec("ff01")
%1 = 65281
Here is a worked out example for Cremona's mwrank.
n = 4;
for (i = 1, n!, print( numtoperm(n, i) ))
for (i = 0, 2^#S - 1, print( vecextract(S, i) ))
forvec( v=vector(k, i,[1,#S]), print(vecextract(S,v)), 2)
aux(n, m, v) =
{
v = concat(v,m);
if(n,
for(i=1, min(m,n), aux(n-i, i, v))
,
print(v)
);
}
partitions(n) = for(i=1, n, aux(n-i, i, []))
Use fordiv, it is no longer specific to integers.
Use pari-2.2.6 or newer: the limit is gone. Lists can be about as large as vectors (roughly 224 entries).
If you are using an older version, this is not straightforward, since the limit is hardcoded into the object representation.
Not straightforward, this is hardcoded into the object representation.
x[1],
x[2] etc. instead of x1, x2).The maximal recursion depth is governed by the maximal stack space allocated to the gp process (not to be confused with gp's stack). Use the following commands from the shell, before starting gp, to increase the process stack size to 20000 kbytes, say:
If your shell is sh-based (bash, ksh, sh, zsh):
ulimit -s 20000
If your shell is csh-based (csh, tcsh):
limit stacksize 20000
Indeed
? system("sleep 5")
time = 0 ms
? for(i=1,10^4, write("testtimer.txt","12345678901"))
time = 164 ms
where the second command actually requires about 10 seconds on my machine.
This is intentional. PARI supports four timing functions clock, times, ftime and getrusage. The last one is chosen by default, since it is the most reliable, and it reports user time, not system or wallclock time. In particular, it should not depend on the system load.
If you really want wallclock time, use
./Configure --time=ftime
make clean && make gp
The time elapsed in commands called via extern is not taken into account by gp's timer. You may want to use
extern("time cmd")
instead of
extern("cmd")
A patch was proposed to take into account the children of the main gp process, but the current behaviour was deemed preferable.
There is no builtin command for this. On a Unix system, the following should work:
timeout(T,E) = extern( Str("ulimit -t ", T, "; echo \"", E, "\" | gp -q -f") )
timeout(T,"E") runs command E for T seconds,
then aborts it. Returns the output of E if it had enough time to run
to completion, and an empty output otherwise. The quotes are essential
around "E", otherwise the function arguments are fully evaluated
before the timer can be set.
? timeout(10, "1 + 1") \\ allow 10 seconds to compute 1 + 1 time = 0 ms. %1 = 2 ? timeout(10, "bnfinit(x^30 - 2)") ? \\ empty output
The GMP kernel is faster than the native one. Here are detailed benchmarks for integer and floating point arithmetic, comparing the two kernels.
The only drawback of using the GMP kernel is that it is sometimes slightly slower for tiny inputs, and that the GMP-enabled PARI library is not compatible with the native one. In particular, all programs built using a native shared library must be rebuilt if it is replaced by a GMP-enabled one.
It depends on what you mean by "handle".
On a 32 bit machines, PARI cannot even represent integers larger than 2^(2^29), or floating point numbers whose absolute binary exponent is larger than 2^30. But the main problem is that PARI's native kernel does not include FFT-based algorithms that would enable it to process huge numbers in a decent time. I would say mantissas of 10000 decimal digits are a sensible practical limit for intensive number crunching with PARI's native kernel.
Using the GMP kernel, almost linear algorithms are available, and only the physical limitations of PARI's representation of integers and floats remain.
Karatsuba multiplication was not enabled by default for type t_REAL until version 2.2.4 included and a quadratic algorithm was used instead. Starting from version 2.2.5, this is no longer the case. Note that all division algorithms in the native kernel are quadratic, though.
Almost linear algorithms are available through the GMP kernel.
All division algorithms in the current PARI kernel are quadratic.
Almost linear algorithms are available through the GMP kernel.
Type
Configure --with-gmp
or possibly
Configure --with-gmp --builddir=some/exotic/directory
if you do not want to mess up with the binaries built using the native kernel.
For instance
(1 << (1<<25)) % (1<<(1<<24) + 1)
This is a general problem in GMP before version 4.2. It may allocate a huge
amount of memory in the process stack space (using alloca).
The best solution is to upgrade to GMP version 4.2 or higher.
Otherwise, one can increase the maximum size of the process stack segment before starting gp from the shell. Alternatively, you can configure an old GMP with
configure --enable-alloca=malloc-noreentrant
but this will slow down GMP.
No.
For gp, Pi is a rational approximation to pi, not exactly pi. With default precision it is within 10-28 of pi. In particular, its sine is not 0.
In fact cos(Pi) is not -1 either, but rather a real number close to -1:
? cos(Pi) + 1
%1 = 1.2744735289059618216 E-57
You cannot. E.g
? simplify((a+b)*(c+d))
%1 = (c + d)*a + (c + d)*b
The reason is that PARI only knows about univariate polynomials (whose coefficient may themselves be polynomials, thus emulating multivariate algebra). The result above is really a simplified form of
( c^1 + (d^1 + 0*d^0)*c^0 )*a^1 + ( (c^1 + (d^1 + 0*d^0)*c^0)*b^1 + 0*b^0 )*a^0
There is no flat or sparse representation of polynomials, and in particular no way to obtain parentheses-free expanded expressions like
a*c + a*d + b*c + b*d
On the other hand, the following routine (courtesy of Michael Somos) comes close to what was requested:
polmonomials(v)=
{ local(w);
if (type(v) != "t_POL",
if(v, [v], [])
,
w = [];
for(n=0, poldegree(v),
w = concat(w, variable(v)^n * polmonomials(polcoeff(v,n)))
); w
)
}
For instance
? polmonomials((a+b)*(c+d))
%1 = [d*b, c*b, d*a, c*a]
Assuming p has exact coefficients, then the output of polroots, is exact to the current accuracy. This means that for each (non-zero) x in the output vector, there exists a complex root z of p such that the relative error |x - z| / |x| is small: its base-10 logarithm is less than -realprecision. This is guaranteed by the algorithm, due to Shönhage and refined by Gourdon, and any other behaviour is a bug worth reporting. On the other hand, the value p(x) itself may be large.
For instance, with 28 digits of accuracy:
? \p28
? p = x^3 - 10^30;
? a = polroots(p)[2]; subst(p, x, a)
%2 = -41.31322575 + 16.13098018*I \\ does not look like 0 !
? \p 100
? b = polroots(p)[2]; subst(p, x, b)
%3 = 0.E-76 + 0.E-76*I \\ b fares better
? abs(a - b) / abs(a)
%4 = 0.E-28 \\ but a is indeed close to b
No.
The integer factorizer implements a variety of algorithms: Trial Division
(all primes up to certain bound, plus user-defined private primetable, via
addprimes), Pure powers, Shanks's SQUFOF, Pollard-Brent rho,
Lenstra's ECM (with Brent and Montgomery improvements), MPQS
(self-initializing, single large prime variation). Output "primes" are BPSW-pseudoprimes.
It cannot compare with special-purpose high-performance implementations, especially non portable ones.
PARI's ECM is slower than GMP-ECM (portable), written by Paul Zimmermann on top of GMP. [comparative benches needed here].
PARI's Quadratic Sieve routinely (less than 1 minute) factors general numbers in the 60 digits range. A few minutes are needed for 70 digits, a few hours for 80 and it has not really been tested in the 90+ range. Besides PARI, we are aware of two other free MPQS implementations:
On my P4 1.6GHz, PARI's MPQS is 2 to 4 times faster than PPMPQS-2.8 in the
60 digits range, 2 times faster for 70 digits, and gets increasingly
slower afterwards. E.g is already 3 times slower for 80 digits. We have not
tested mpqs4linux.
The factors output are not proven primes if they are larger than 1015. They are Baillie-Pomerance-Selfridge-Wagstaff pseudoprimes, among which no composite has been found yet, although it is expected infinitely many exist. If you need a proven factorization, use
fa = factor(N);
isprime( fa[,1] )
Another possibility is to set
default(factor_proven, 1);
This ensures that a true primality proof is computed whenever we do
not explicitly ask for partial factorization. This will slow down
PARI if largish primes turn up often.
*** isprime: Warning: IFAC: untested integer declared prime.Can I really trust the output of isprime()?
Yes! During tests of Pocklington-Lehmer type, isprime() calls factor() which may produce this message, but it does go on to prove recursively that the prime factors encountered really are prime.
ispseudoprime(N) uses Baillie-Pomerance-Selfridge-Wagstaff
deterministic test, which is a strong Rabin-Miller test for base 2, followed
by a strong Lucas test for the sequence (P, -1), where P is the smallest
positive integer such that P^2 - 4 is not a square modulo N. It has been
checked that no composite smaller than 1015 passes this test.
ispseudoprime(N, k) with k > 0, performs
a strong Rabin-Miller test for k randomly chosen bases (with
end-matching to catch roots of -1). This variant is about as fast as
ispseudoprime(N) when k = 3, and slower
afterwards.
isprime(N) first performs an ispseudoprime test.
If N is larger than 1015, Selfridge-Pocklington-Lehmer p-1 test is
performed if N-1 is smooth (factored part is larger than N1/3 --
Konyagin-Pomerance or more involved algorithms based on lattice reduction
are not implemented). Otherwise we use APR-CL (Adleman-Pomerance-Rumely,
Cohen-Lenstra).
PARI's APR-CL can prove primality of numbers with 1000 digits in a few hours. Some examples on a 2.8GHz Pentium IV:
| N | # digits | time |
| 3^2833 - 2^2833 | 1352 | 4 hours |
| (772! + 1)/773 | 1893 | 17 hours |
In its present form (January 2006), the Agrawal-Kayal-Saxena primality test is not practical for large numbers. Efficiency not being a primary concern, it is a simple exercise in GP programming to implement the algorithm.
You get messages of the form
IFAC: trying Lenstra-Montgomery ECM
ECM: working on 64 curves at a time; initializing for up to 50 rounds...
ECM: time = 0 ms
ECM: dsn = 12, B1 = 1800, B2 = 198000, gss = 128*420
ECM: time = 68920 ms, B1 phase done, p = 1801, setting up for B2
ECM: time = 1210 ms, entering B2 phase, p = 2017
ECM: time = 42990 ms
[...]
dsn = 72, B1 = 1073500, B2 = 118085000, gss = 1024*420 :
and the last line is repeated seemingly forever.
The short answer is: ECM is still making progress, although the diagnostics do not show it (detailed explanation).
This problem is specific to version 2.2.10 and occurs while factoring smallish integers with MPQS. For instance
factor(590772364016942232621279991);
*** factor: Warning: MPQS: sizing marginal, index1 too large or too many primes in A.
This diagnostic was left in place to detect suboptimal behaviour with respect to the chosen tunings. This is not a bug and the result is correct. The associated tuning problem is fixed in 2.2.11.
J is not an ideal. It is a matrix in HNF form, but the Z-module generated by its columns is not an OK-module (exercise).
? K = nfinit(x^2+1); J = [5, 4; 0, 1]; nfisideal(K, J)
%1 = 0
In fact, the two ideals above 5 are
? P = idealprimedec(K, 5);
? idealhnf(K, P[1])
%3 =
[5 3]
[0 1]
? idealhnf(K, P[2])
%4 =
[5 2]
[0 1]
PARI routines transform suitable objects (element in K, prime ideal) to
an ideal (matrix in HNF), but do not check that matrices of the right
dimension do indeed correspond to ideals, which would be quite costly. If you
are unsure about how to input ideals, use
idealhnf systematically, or ideals produced by other PARI
routines.
If you insist on finding the OK-module generated by the columns of your matrix, you may use something like
IDEAL(nf, J) =
{ local(M,i,j);
J = mathnf(J); M = [;];
for (i = 1, poldegree(nf.pol),
for (j = 1, #J,
M = concat(M, nfeltmul(nf, J[,j], nf.zk[i]))
)); mathnf(M)
}
For better efficiency, J is put in HNF form first. You can also use directly
idealhnf if your matrix does not have maximal rank, yielding the
faster variant:
IDEAL(nf, J) =
{
J = mathnf(J);
mathnfmodid( idealhnf(nf, vecextract(J, "2..")), J[1,1] )
}
As written, this assumes J is contained in OK and does not work if J has Z-rank 0 or 1, but you get the idea.
The simplest solution is to use something like
nfinitpartial(P) = nfinit( [P, nfbasis(P,1)] )
which computes nfinit from P and a tentative integer
basis (obtained by considering small primes individually, and possibly lazy
localization with respect to a large composite cofactor). This integer basis
may be wrong, namely generate an order which is not maximal at a large prime
number.
You may also try to replace P by either
polredabs(P, 16)
or
polred(P, 1)
Both commands produce an integral monic polynomial generating the same number
field as P, hopefully simpler. polredabs may be very
expensive when that number field has many subfields, so use
polred in that case; it produces worse polynomials, but is much
faster.
This should succeed whenever the index of Z[X]/(P) in the maximal order is mostly divisible by small primes, or if the field discriminant is smooth. If the index has large prime divisors, on the other hand, you will make little progress. Try to factor the discriminant using
factorint(poldisc(P))
at \g4 debugging level, and record the large prime numbers which
are discovered by the factorization engine. Then use
addprimes(p)
for each such p, or addprimes(L) for a list L of
such primes. Of course, if you know in advance non trivial factors of the
discriminant, you should tell the factoring engine about them using
addprimes.
It is allowed to fuel the addprimes engine with composite
numbers, but you then run the risk of triggering the following message
*** impossible inverse modulo: Mod(...)
This stops the computation, but yields a non-trivial factor of the modulus (take the gcd of the number and the modulus).
The routine bnfinit assumes the truth of the Generalized Riemann
Hypothesis, in order to have a small set of generators for the class group.
From 2.4.0 on, bnfinit assumes no more than the GRH.
It is marginally slower than before, but no longer outputs wrong results
(that we know of !).
generate the class group of K. This is not known even under the GRH, which implies the above with 12 in place of 0.3 (a famous result of Eric Bach); the constant can be further optimized depending on the signature and assuming the discriminant is large enough, but not down to 0.3 in any case. Also, there are counter-examples even for relatively large discriminants : unlucky random seeds actually give a wrong result (see TODO file).
So, technically, any result computed using a bnf from
bnfinit is conditional, and depends on a heuristic strengthening
of a result already depending on the GRH. One can make sure that no more
than the GRH is used with
bnfinit(P,, [0.3, 12]) /* notice the double comma after P */
instead of
bnfinit(P)
This is slower, of course. At worst about 40 times slower when the original
result would have been right. The routine quadclassunit has the
same problem as bnfinit (0.3 being replaced by 0.2).
If you want to use no more than the GRH, use
quadclassunit(D,, [0.2, 6])
You can also use
bnf = bnfinit(P);
bnfcertify(bnf)
to certify the bnf unconditionally (if the discriminant is
large, this will be slow, and may overflow the possibilities of the
implementation). Results obtained using a certified bnf do
not depend on the GRH.
Besides all functions using a bnf argument, the routines
bnfclassunit, bnfclgp, bnfreg are
trivial wrappers around bnfinit and make the same assumptions.
All these routines predate the introduction of member functions and are now
obsolete: bnfinit gives much more information in about the same
time, is about as easy to use (using member functions), and has the further
advantage of being an init function.
Note that output from quadclassunit cannot be certified.
Use bnfinit instead if certification is important. Finally,
quadunit and quadregulator do not
assume the GRH (and are much slower than quadclassunit).
If E/F_p has complex multiplication by a principal imaginary quadratic order, we use a very fast explicit formula (quadratic in log p). Otherwise we rely on Shanks-Mestre's baby-step giant step method, whose run-time is unfortunately exponential. Hence this naive algorithm becomes unreasonable when p has about 30 decimal digits.
To go beyond that, you are advised to try PARI version 2.4.3 or higher. Then download and install the optional package SEAdata, implementing the Schoof-Elkies-Atkin algorithm. It should work at least up to 400 decimal digits, way beyond cryptographic sizes. This data package replaces the old SEA script mentioned below.
You may want to try the external ellSEA package (1.5MB), which implements Schoof's polynomial-time algorithm (with improvements by Elkies and Atkin). It should work at least up to about 150 digits, way beyond cryptographic sizes.
You may use SAGE
(online demo calculator).
You may also call mwrank directly from gp, using extern.
Many routines expect operands of a specified type and shape, for instance nf,bnf,bnr,ell,modpr structures. For reasons of efficiency, only minimal checking is done, and it is possible to fool the routine with an object of correct type and length, but definitely unexpected content. (It is easy to do it on purpose, but difficult to generate an example randomly.)
When a subtle undetected violation is performed, one of the following may occur: a regular exception is raised, the PARI stack overflows, a SIGSEGV or SIGBUS signal is generated, or we enter an infinite loop. Since signals are trapped, a Control-C will get you out of any infinite loop, and the garbage collector is rather robust, it is exceedingly rare to lose a session.
Bugs in this class (Bad Input Bugs, or BIBs for short) are never fixed whenever a noticeable code complication or performance penalty would occur. All PARI objects and internals are exposed and at your disposal, which is both convenient and unsafe. There are plans to add more semantic to PARI types so that e.g., an ell structure would record the fact that it is actually defined over Z, and a minimal integral Weierstrass model. This would automatically cure problems such as
? a6; ellap( ellinit([a1,a2,a3,a4,a6]), 2 );
*** ellap: bug in GP (Segmentation Fault), please report
? intnum(x = 0, 10^6, x^2 * exp(-x^2))
%1 = 0e-55
? intnum(x=0, Pi, sin(2^11 * x)^2)
%2 = 3.788253725 E-46
The first integral should be indistinguishable from the integral to infinity, which is sqrt(Pi)/4, and the second is exactly Pi/2.
Numeric integration evaluates the function at regularly spaced points in the integration interval, then uses a standard quadrature rule. This is tried for a few decreasing values of the subdivision diameter until the value appears to stabilize (in fact, was accurately predicted by interpolation from the four last values). This gives incorrect results for functions of period a large power of 2 as in the second example. Or if most of the weight is concentrated in a very small interval as in the first example, where the function is essentially 0 in most of the interval.
For the time being, the user must break the integral into harmless chunks, or perform relevant change of variables herself. E.g make sure the function varies slowly, move singularities to the boundary or (better) remove them, restrict the interval to a region where the function is not negligible.
Starting from version 2.2.9, we use the "double exponential method" instead of standard Romberg, and get better results:
? intnum(x = 0, 10^6, x^2 * exp(-x^2))
%1 = 0.4431134622969907111829434676
? intnum(x = 0, [[1], 2], x^2 * exp(-x^2)) / (sqrt(Pi) / 4)
%2 = 1.000000000000000000000000000
? intnum(x=0, Pi, sin(2^11 * x)^2)
%3 = 1.482067574055964668166042567 \\ not so good
? intnum(x=0, Pi, sin(2^11 * x)^2, 11) / (Pi/2)
%4 = 1.000000000000000000000000000 \\ slower but perfect
It still helps a lot to split the integral into pieces, and indicate the growth rate and oscillation pattern of the function at their boundary, and possibly to use more sampling points for oscillating functions. Otherwise, loss of accuracy may be expected.
For instance
? Mod(1,2) * Mod(1,3)
%1 = Mod(0, 1)
? Mod(1,2)*x + Mod(1,3)
%2 = Mod(1,2)*x
Elementary operations involving INTMODs and POLMODs with different moduli p and q first project both operands to the quotient modulo the ideal (p,q). In the above, the moduli 2 and 3 are coprime and we obtain the 0 ring. This explains the first example. For the second, notice that the result of Mod(1,2)*x is a polynomial in (Z/2Z)[x] whose constant coefficient is Mod(0,2).
Another example is
? Mod(0,2)^2
%3 = Mod(0, 2)
whereas one might expect Mod(0,4) instead. Modular object operate with respect to a fixed modulus, which is never increased (it may be replaced by a strict divisor as explained above).
To keep track of accuracy, you may use PADICs, but operations will be much slower than with INTMODs. If efficiency is a primary concern, use INTs, keeping track of denominators separately, and reduce with relevant moduli to avoid coefficient explosion.
It depends on what you mean by that. PARI was not written to handle huge objects (e.g. millions of digits, dimension in the thousands). It can do it, but different algorithms should be used for these ranges, which are not currently implemented. PARI is quite good for heavy-duty computations in algebraic number theory.
A self-installing binary is available from the download page. The single pending Windows-specific bug concerns the online help which does not work on some configurations (we lack details at this point). See the README.
I am assuming that you started the GP.EXE program by
(double) clicking on the PARI link which was created on the windows
Desktop. Check the Properties of the link in the right-click menu. If you
did not tamper with it at installation time, it states that the program
starts in directory
C:\Program Files\PARI
This is where GP expects to find its input files. That is if you type in
\r foo
this will work if and only if we have a file foo or
foo.gp in
C:\Program Files\PARI
Change this if you want to store your GP files in some other directory.
Another possibility: some Windows editors (e.g. NotePad) automatically
append a .txt suffix to the files they create. Make sure
your file is really called foo and not foo.txt.
Right click on the GP icon and change Properties->Misc:
select 'Quick Edit' and deselect 'Fast Pasting'.
See also the GAP help pages on the subject.
You are running a keyboard driver for some non-English languages. These
drivers catch the ^ character to produce the circumflex accent
and do not pass it properly to PARI. You may try to type a SPACE
after the ^ dead key, to produce a true ^;
it may or may not fix the problem.
This depends on your Windows installation and the compilers used when
building the binary. It may work properly, or it may emit a beep and
prints a tilde character '~'. In the latter case, no fix is
known: use Backspace instead, or readline commands.
Current Cygwin versions do not exhibit the problem anymore (at least from april 2006 on). The next release of the self-extractor will fix the problem.
Cygwin is a UNIX environment for Windows, developed by Red Hat Software. It consists of two parts:
The binary gp.exe we provide only needs the DLL to function. If
you wish to compile gp yourself (as described in the manual for Unix
systems), you need the Cygwin building tools.
If the error looks like
process_begin: CreateProcess((null), /cygdrive/c/mingw/bin/gcc -c ...) failed.
make (e=3): The system cannot find the path specified.
c:\mingw\bin\make.exe: *** [kernel.o] Error 3
then you did not install the full Cygwin distribution
(specifically, you are missing Cygwin's gcc here). In the
example above, there are already versions of gcc and
make available from the mingw compiler,
adding to the confusion.
Using recent Cygwin development tools (downloaded March 20 2003) and Subversion version of PARI, a DLL is produced which at least works under Cygwin. I do not know whether it still works in other build environments. Please test and report!
This version of gp is not maintained anymore. You can get a 6 years-old legacy binary gp-2.0.14.sit.bin from the FTP archives. This is a sit archive, which may be extracted using StuffIt expander.
You can get gp from the Fink distribution. But Fink only provides the stable version of pari, and you may want to try new features.
Alternatively you can install Apple's Developer Tools, which should be on one of your installation's CDs. The current name of the tool suite is Xcode. Then you can compile pari yourself:
You probably want to read this FAQ, which solves a very common problem.
Towards the end, Configure prints the message
###
### Editline wrapper detected, building without readline support
###
You can still compile gp, but line editing, history and command completion will not work.
This one is tough, since we must solve three problems. First the immediate one (readline missing), which is standard and not OS X-specific, then two problems with the normal fix. (The fault lies with Apple in both cases, and cannot be fixed within PARI.) The good news is that the problem will then be fixed once and for all on your system.
First things first: although a libreadline.dylib library is part of OS 10.4, it is a fake: a simple link to a replacement library, editline, which is only partially compatible with readline. So we must install the true readline library. But the linker shipped with OS 10.4 will insist on using the default libreadline.dylib instead of the correct one (ours), unless we provide it as a shared library (this is a bug). Unfortunately, in the current readline-5.0, building a shared library on OS 10.4 fails (flags changed their meanings in Apple's gcc in between versions). So we must fix that also.
This requires a recent version of PARI, for instance the current stable version pari-2.3.2. Your can also try the current svn snapshot (bleeding edge release).
Download http://ftp.gnu.org/pub/gnu/readline/readline-5.0.tar.gz.
Extract the archive (tar zxf readline-5.0.tar.gz), then copy-paste the following in a Terminal:
cd readline-5.0
./configure
sed -e 's/-dynamic/-dynamiclib/' shlib/Makefile > shlib/Makefile.ok
mv shlib/Makefile.ok shlib/Makefile
make && sudo make install
This installs a proper shared library for readline, under /usr/local (thus it does not disturb your vendor's libraries and won't cause problems in other parts of the system). Now rejoice: you have solved the readline problem once and for all.
Now, you can compile and install PARI, making sure that the proper readline is used:
./Configure --with-readline=/usr/local
make && sudo make install
Under MacOS 9 and below, you can't. With MacOS X, you can install an XFree server and X11 developer's libraries, then compile gp as explained above. It should pick up the X11 installation. Another possibility is to download the FLTK library and try to
Configure --with-fltk=path to your installed fltk library
| Karim Belabas |
|
| 2009-06-15 0:49:37 |