Commit 379f897a authored by Andreas Enge's avatar Andreas Enge

Initial commit.

parents
Andreas Enge <andreas.enge@inria.fr>
This diff is collapsed.
DIR=/usr/local
PARI=${DIR}/pari-dev
GMP=${DIR}/gmp-5.1.3
MPFR=${DIR}/mpfr-3.1.2
MPC=${DIR}/mpc-1.0.2
# DO NOT CHANGE ANYTHING BEYOND THIS LINE #
VERSION=0.0.1
DISTDIR=pari-gnump-${VERSION}
INCLUDE=-I${PARI}/include -I${GMP}/include -I${MPFR}/include -I${MPC}/include
LIB=-L${PARI}/lib -L${GMP}/lib -L${MPFR}/lib -L${MPC}/lib -lpari -lmpc -lmpfr -lgmp
all : libpari-gnump.so
check : tests
${PARI}/bin/gp < example.gp
./tests
dist : clean
mkdir ${DISTDIR}
cp AUTHOR COPYING example.gp Makefile pari-gnump.c pari-gnump.h pari-gnump-user.c pari-gnump-user.h README tests.c ${DISTDIR}
tar zcvf ${DISTDIR}.tar.gz ${DISTDIR}/
tests : libpari-gnump.so tests.o
gcc -o tests tests.o libpari-gnump.so ${LIB} -Wl,-rpath=${PARI}/lib -Wl,-rpath=${GMP}/lib -Wl,-rpath=${MPFR}/lib -Wl,-rpath=${MPC}/lib -Wl,-rpath=.
tests.o : tests.c pari-gnump.h
gcc -c tests.c ${INCLUDE}
libpari-gnump.so : pari-gnump.o pari-gnump-user.o
gcc -o libpari-gnump.so $^ -shared ${LIB} -Wl,-rpath=${PARI}/lib
%.o : %.c %.h
gcc $< -c -fPIC ${INCLUDE}
clean :
-rm -rf *.o libpari-gnump.so tests ${DISTDIR} ${DISTDIR}.tgz
Pari-gnump allows to easily switch between number types from the GNU
multiprecision ecosystem (mpz, mpq, mpfr, mpc) and corresponding types in
libpari, as well as to use functions from the GNU libraries in GP.
INSTALLATION
------------
Adapt the location and version numbers of Pari/GP, GNU MP (GMP), GNU MPFR
and GNU MPC in the first lines of Makefile.
If you wish to use functions relying on GMP, GNU MPFR or GNU MPC in your GP
session, you need to add them to pari-gnump-user.h and pari-gnump-user.c
(see below).
Execute
make
This creates the library libpari-gnump.so. There is no need to install this
library, it may simply be copied to the location of your Pari/GP project.
To test if everything went well, execute
make check
NUMBER TRANSFORMATIONS
----------------------
The exported functions are given in pari-gnump.h. For each GNU
multiprecision type X from mpz, mpq, mpfr and mpc, there are functions
X_set_GEN
X_get_GEN
The first one takes an argument x of type X and g of type GEN, and assigns
the value of g to x if the types are compatible. The second one takes as
argument x of type X and returns a GEN with the same value, allocated on the
Pari stack.
The functions pari_mpfr_init2, pari_mpc_init2 and pari_mpc_init3 work
exactly as their counterparts without the pari prefix, except that they
allocate the mantissae on the Pari stack instead of using the GMP memory
allocation functions. They should not be freed with calls to mpfr_clear or
mpc_clear, but with the usual Pari stack handling ("avma").
pari_mpfr_init_set_GEN and pari_mpc_init_set_GEN combine such an
initialisation with an assignment; as precision, they use the precision of
the source variable if it is of a floating point type, or the default
precision if it is an integer or a rational number. If you are unsure what
is meant by this paragraph, you probably do not wish to use these functions.
USER FUNCTIONS IN GP
--------------------
If you wish to use functions relying on GMP, GNU MPFR or GNU MPC in a GP
session, you need to write a wrapper function in pari-gnump-user.c and
declare it in pari-gnump-user.h. The function should take GEN arguments,
transform them to mpz, mpq, mpfr or mpc, compute with them in the desired
way, and transform the result back into a GEN. By calling
make
this function will be automatically included into the library
libpari-gnump.so. The functions must then be made available inside your GP
session using the "install" command.
As examples, pari-gnump-user.c implements correctly rounded multiplications
of real and complex numbers using MPFR and MPC, as well as the erf and zeta
functions from MPFR. The file example.gp shows how to install and use these
functions in a GP script.
/*
Copyright © 2014 Andreas Enge <andreas.enge@inria.fr>
This file is part of pari-gnump.
Pari-gnump is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.
Pari-gnump is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with Pari-gnump. If not, see <http://www.gnu.org/licenses/>.
*/
install ("pari_mpc_mul", "GGp", "mpc_mul", "./libpari-gnump.so");
install ("pari_mpfr_mul", "GGp", "mpfr_mul", "./libpari-gnump.so");
install ("pari_mpfr_erf", "Gp", "mpfr_erf", "./libpari-gnump.so");
install ("pari_mpfr_zeta", "Gp", "mpfr_zeta", "./libpari-gnump.so");
mpc_mul (1+2*I, 3-I)
(1+2*I) * (3-I)
mpfr_mul (Pi, Pi)
\p 200
mpc_mul (Pi, Pi)
\p 100000
pi = Pi;
mpfr_mul (pi, pi) - pi*pi
\p 38
mpfr_erf (1)
mpfr_zeta (2)
/*
Copyright © 2014 Andreas Enge <andreas.enge@inria.fr>
This file is part of pari-gnump.
Pari-gnump is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.
Pari-gnump is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with Pari-gnump. If not, see <http://www.gnu.org/licenses/>.
*/
#include "pari-gnump-user.h"
/* The following functions perform floating point computations using the
mpfr and mpc libraries. They take as additional parameter a long prec,
coding a pari floating point precision which is used for the result.
Floating point arguments are first transformed into variables of mpfr
type, using their own precision; the parameter prec is also used to
determine the floating point precision of exact arguments of type t_INT
or t_FRAC, for which there is no intrinsic way of deriving a precision.
If such a function is used inside GP, it should be installed with a
parameter code of "p" for the prec parameter, which is then automagically
replaced by the current default precision.
*/
/****************************************************************************/
GEN pari_mpc_mul (GEN x, GEN y, long prec)
/* Returns the product of x and y.
Install into GP with
install ("pari_mpc_mul", "GGp", "mpc_mul", "./libpari-gnump.so");
*/
{
mpfr_prec_t p = bit_accuracy (prec);
mpc_t z, z1, z2;
pari_mpc_init2 (z, p);
pari_mpc_init_set_GEN (z1, x, p);
pari_mpc_init_set_GEN (z2, y, p);
mpc_mul (z, z1, z2, MPC_RNDNN);
return mpc_get_GEN (z);
}
/****************************************************************************/
GEN pari_mpfr_mul (GEN x, GEN y, long prec)
/* Returns the product of x and y.
Install into GP with
install ("pari_mpfr_mul", "GGp", "mpfr_mul", "./libpari-gnump.so");
*/
{
mpfr_prec_t p = bit_accuracy (prec);
mpfr_t z, z1, z2;
pari_mpfr_init2 (z, p);
pari_mpfr_init_set_GEN (z1, x, p);
pari_mpfr_init_set_GEN (z2, y, p);
mpfr_mul (z, z1, z2, MPFR_RNDN);
return mpfr_get_GEN (z);
}
/****************************************************************************/
GEN pari_mpfr_erf (GEN x, long prec)
/* Returns the error function at x.
Install into GP with
install ("pari_mpfr_erf", "Gp", "mpfr_erf", "./libpari-gnump.so");
*/
{
mpfr_prec_t p = bit_accuracy (prec);
mpfr_t z, z1;
pari_mpfr_init2 (z, p);
pari_mpfr_init_set_GEN (z1, x, p);
mpfr_erf (z, z1, MPFR_RNDN);
return mpfr_get_GEN (z);
}
/****************************************************************************/
GEN pari_mpfr_zeta (GEN x, long prec)
/* Returns the Riemann zeta function at x.
Install into GP with
install ("pari_mpfr_zeta", "Gp", "mpfr_zeta", "./libpari-gnump.so");
*/
{
mpfr_prec_t p = bit_accuracy (prec);
mpfr_t z, z1;
pari_mpfr_init2 (z, p);
pari_mpfr_init_set_GEN (z1, x, p);
mpfr_zeta (z, z1, MPFR_RNDN);
return mpfr_get_GEN (z);
}
/****************************************************************************/
/*
Copyright © 2014 Andreas Enge <andreas.enge@inria.fr>
This file is part of pari-gnump.
Pari-gnump is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.
Pari-gnump is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with Pari-gnump. If not, see <http://www.gnu.org/licenses/>.
*/
#include "pari-gnump.h"
GEN pari_mpc_mul (GEN x, GEN y, long prec);
GEN pari_mpfr_mul (GEN x, GEN y, long prec);
GEN pari_mpfr_erf (GEN x, long prec);
GEN pari_mpfr_zeta (GEN x, long prec);
This diff is collapsed.
/*
Copyright © 2014 Andreas Enge <andreas.enge@inria.fr>
This file is part of pari-gnump.
Pari-gnump is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.
Pari-gnump is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with Pari-gnump. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <limits.h>
#include <assert.h>
#include <mpc.h>
#include <pari/pari.h>
void mpz_set_GEN (mpz_ptr z, GEN x);
GEN mpz_get_GEN (mpz_srcptr z);
void mpq_set_GEN (mpq_ptr q, GEN x);
GEN mpq_get_GEN (mpq_srcptr q);
void pari_mpfr_init2 (mpfr_ptr f, mpfr_prec_t prec);
void pari_mpfr_init_set_GEN (mpfr_ptr f, GEN x, mpfr_prec_t default_prec);
int mpfr_set_GEN (mpfr_ptr f, GEN x, mpfr_rnd_t rnd);
GEN mpfr_get_GEN (mpfr_srcptr f);
void pari_mpc_init2 (mpc_ptr c, mpfr_prec_t prec);
void pari_mpc_init3 (mpc_ptr c, mpfr_prec_t prec_re, mpfr_prec_t prec_im);
void pari_mpc_init_set_GEN (mpc_ptr c, GEN x, mpfr_prec_t default_prec);
int mpc_set_GEN (mpc_ptr c, GEN x, mpc_rnd_t rnd);
GEN mpc_get_GEN (mpc_srcptr c);
/*
Copyright © 2014 Andreas Enge <andreas.enge@inria.fr>
This file is part of pari-gnump.
Pari-gnump is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.
Pari-gnump is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with Pari-gnump. If not, see <http://www.gnu.org/licenses/>.
*/
#include "pari-gnump.h"
/****************************************************************************/
/* */
/* Test functions */
/* */
/****************************************************************************/
int test_z_to_GEN (mpz_srcptr z)
{
int ret;
mpz_t z2;
GEN x;
mpz_init (z2);
x = mpz_get_GEN (z);
mpz_set_GEN (z2, x);
ret = mpz_cmp (z, z2);
mpz_clear (z2);
return (ret);
}
/****************************************************************************/
int test_GEN_to_z (GEN x)
{
int ret;
mpz_t z;
GEN x2;
mpz_init (z);
mpz_set_GEN (z, x);
x2 = mpz_get_GEN (z);
mpz_clear (z);
return (cmpii (x, x2));
}
/****************************************************************************/
int test_q_to_GEN (mpq_srcptr q)
{
int ret;
mpq_t q2;
GEN x;
mpq_init (q2);
x = mpq_get_GEN (q);
mpq_set_GEN (q2, x);
ret = mpq_cmp (q, q2);
mpq_clear (q2);
return (ret);
}
/****************************************************************************/
int test_GEN_to_q (GEN x)
{
int ret;
mpq_t q;
GEN x2;
mpq_init (q);
mpq_set_GEN (q, x);
x2 = mpq_get_GEN (q);
mpq_clear (q);
return (gcmp (x, x2));
}
/****************************************************************************/
int test_f_to_GEN (mpfr_srcptr f)
{
int ret;
mpfr_t f2;
GEN x;
mpfr_init2 (f2, mpfr_get_prec (f));
x = mpfr_get_GEN (f);
mpfr_set_GEN (f2, x, MPFR_RNDN);
ret = mpfr_cmp (f, f2);
mpfr_clear (f2);
return (ret);
}
/****************************************************************************/
int test_f_to_GEN_stack (mpfr_srcptr f)
{
int ret;
mpfr_t f2;
GEN x;
pari_mpfr_init2 (f2, mpfr_get_prec (f));
x = mpfr_get_GEN (f);
mpfr_set_GEN (f2, x, MPFR_RNDN);
ret = mpfr_cmp (f, f2);
return (ret);
}
/****************************************************************************/
int test_GEN_to_f (GEN x)
{
int ret;
mpfr_t f;
GEN x2;
mpfr_init (f);
mpfr_set_GEN (f, x, MPFR_RNDN);
x2 = mpfr_get_GEN (f);
mpfr_clear (f);
return (gcmp (x, x2));
}
/****************************************************************************/
int test_GEN_to_f_stack (GEN x)
{
int ret;
mpfr_t f;
GEN x2;
pari_mpfr_init_set_GEN (f, x, 53);
x2 = mpfr_get_GEN (f);
return (gcmp (x, x2));
}
/****************************************************************************/
int test_c_to_GEN (mpc_srcptr c)
{
int ret;
mpc_t c2;
GEN x;
mpc_init3 (c2, mpfr_get_prec (mpc_realref (c)),
mpfr_get_prec (mpc_imagref (c)));
x = mpc_get_GEN (c);
mpc_set_GEN (c2, x, MPC_RNDNN);
ret = mpc_cmp (c, c2);
mpc_clear (c2);
return (ret);
}
/****************************************************************************/
int test_c_to_GEN_stack (mpc_srcptr c)
{
int ret;
mpc_t c2;
GEN x;
pari_mpc_init3 (c2, mpfr_get_prec (mpc_realref (c)),
mpfr_get_prec (mpc_imagref (c)));
x = mpc_get_GEN (c);
mpc_set_GEN (c2, x, MPC_RNDNN);
ret = mpc_cmp (c, c2);
return (ret);
}
/****************************************************************************/
int test_GEN_to_c (GEN x)
{
int ret;
mpc_t c;
GEN x2;
mpc_init2 (c, 100);
mpc_set_GEN (c, x, MPC_RNDNN);
x2 = mpc_get_GEN (c);
mpc_clear (c);
return (gequal (x, x2) != 1);
}
/****************************************************************************/
int test_GEN_to_c_stack (GEN x)
{
int ret;
mpc_t c;
GEN x2;
pari_mpc_init_set_GEN (c, x, 53);
x2 = mpc_get_GEN (c);
return (gequal (x, x2) != 1);
}
/****************************************************************************/
void test ()
{
mpz_t z, z2;
mpq_t q, q2;
mpfr_t f, f2;
mpc_t c, c2;
gmp_randstate_t rand;
GEN x;
mpfr_prec_t prec;
pari_init (1000000, 0);
gmp_randinit_default (rand);
mpz_init (z);
mpz_set_si (z, 3);
mpz_mul_2exp (z, z, 200);
assert (test_z_to_GEN (z) == 0);
mpz_neg (z, z);
assert (test_z_to_GEN (z) == 0);
mpz_set_si (z, 0);
assert (test_z_to_GEN (z) == 0);
mpz_clear (z);
x = shifti (stoi (3), 200);
assert (test_GEN_to_z (x) == 0);
x = negi (x);
assert (test_GEN_to_z (x) == 0);
assert (test_GEN_to_z (gen_0) == 0);
mpq_init (q);
mpq_set_si (q, 3, 1);
mpq_div_2exp (q, q, 200);
assert (test_q_to_GEN (q) == 0);
mpq_neg (q, q);
assert (test_q_to_GEN (q) == 0);
mpq_set_ui (q, 3, 1);
assert (test_q_to_GEN (q) == 0);
mpq_set_ui (q, 0, 1);
assert (test_q_to_GEN (q) == 0);
mpq_clear (q);
x = gmul2n (stoi (3), -200);
assert (test_GEN_to_q (x) == 0);
x = gneg (x);
assert (test_GEN_to_q (x) == 0);
assert (test_GEN_to_q (stoi (3)) == 0);
assert (test_GEN_to_q (gen_0) == 0);
mpfr_init (f);
mpfr_set_si (f, 3, MPFR_RNDN);
mpfr_div_2ui (f, f, 200, MPFR_RNDN);
assert (test_f_to_GEN (f) == 0);
assert (test_f_to_GEN_stack (f) == 0);
mpfr_neg (f, f, MPFR_RNDN);
assert (test_f_to_GEN (f) == 0);
assert (test_f_to_GEN_stack (f) == 0);
mpfr_set_ui (f, 0, MPFR_RNDN);
assert (test_f_to_GEN (f) == 0);
assert (test_f_to_GEN_stack (f) == 0);
for (prec = 2; prec < 200; prec++) {
mpfr_set_prec (f, prec);
mpfr_urandomb (f, rand);
assert (test_f_to_GEN (f) == 0);
assert (test_f_to_GEN_stack (f) == 0);
}
mpfr_clear (f);
x = gmul2n (stor (3, 6), -200);
assert (test_GEN_to_f (x) == 0);
x = gneg (x);
assert (test_GEN_to_f (x) == 0);
assert (test_GEN_to_f (stor (3, 6)) == 0);
assert (test_GEN_to_f (real_0_bit (-230)) == 0);
assert (test_GEN_to_f (real_0_bit (230)) == 0);
assert (test_GEN_to_f (gen_0) == 0);
x = gmul2n (stor (3, 6), -200);
assert (test_GEN_to_f_stack (x) == 0);
x = gneg (x);
assert (test_GEN_to_f_stack (x) == 0);
assert (test_GEN_to_f_stack (stor (3, 6)) == 0);
assert (test_GEN_to_f_stack (real_0_bit (-230)) == 0);
assert (test_GEN_to_f_stack (real_0_bit (230)) == 0);
mpc_init2 (c, 123);
mpc_set_ui_ui (c, 3, 5, MPC_RNDNN);
mpc_div_2ui (c, c, 100, MPC_RNDNN);
assert (test_c_to_GEN (c) == 0);
assert (test_c_to_GEN_stack (c) == 0);
mpc_set_ui_ui (c, 0, 5, MPC_RNDNN);
assert (test_c_to_GEN (c) == 0);
assert (test_c_to_GEN_stack (c) == 0);
mpc_set_ui_ui (c, 3, 0, MPC_RNDNN);
assert (test_c_to_GEN (c) == 0);
assert (test_c_to_GEN_stack (c) == 0);
for (prec = 2; prec < 200; prec++) {
mpc_set_prec (c, prec);
mpc_urandom (c, rand);
assert (test_c_to_GEN (c) == 0);
assert (test_c_to_GEN_stack (c) == 0);
}
x = gmul2n (mkcomplex (stor (3, 6), stor (5, 6)), -200);
assert (test_GEN_to_c (x) == 0);
assert (test_GEN_to_c_stack (x) == 0);
x = gmul2n (mkcomplex (gen_0, stor (5, 6)), -200);
assert (test_GEN_to_c (x) == 0);
assert (test_GEN_to_c_stack (x) == 0);
x = gmul2n (mkcomplex (stor (3, 6), gen_0), -200);
assert (test_GEN_to_c (x) == 0);
assert (test_GEN_to_c_stack (x) == 0);
mpc_clear (c);
pari_close ();
gmp_randclear (rand);
}
/****************************************************************************/
int main ()
{
test ();
printf ("All library tests passed.\n");
return 0;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment