Commit 1012e9cf authored by Mike Furr's avatar Mike Furr

Imported Upstream version 20021123

parents
2002-11-23 <david@picsou.chatons>
* Replaced -ggdb by -g - MacOS X.2 gcc does not generate proper
assembly code if -ggdb is used.
* Bugfix wrt Achim Blumensath: do not store parameter's mpz_t
address into local variables.
2002-11-21 <david@picsou.chatons>
* Fixed Z.legendre and Z.jacobi
* More testing for Z
* Changed Makefile not to build creal by default.
* Changed Makefile to stop if autotest fails.
* Added tests for Q.
* Changed comparison operators in Q.Infixes to / suffix.
2002-07-29 David Monniaux <monniaux@gamma.csl.sri.com>
* Upgraded to GMP 4.1
2002-01-24 David Monniaux <monniaux@basilic.ens.fr>
* Fixed compilation for MacOS X.
* Added hash functions etc...
* Fixed lots of bugs.
2001-12-31 David Monniaux <monniaux@genievre.ens.fr>
* First public release (3.1-3.04-1).
Some usual problems:
*** COMPILATION
**** On my Ultra Sparc, ML GMP won't link with GNU MP.
On the Ultra Sparc under Solaris, GNU MP tries to use Sun's
proprietary compiler cc with option -xtarget=v9 to generate code
optimised for the Ultra Sparc. This unfortunately generates 64-bit
code, which cannot be linked with 32-bit code.
As of version 3.04, Objective Caml cannot be compiled for 64-bit Ultra
Sparc, at least because of some forbidden constructs in some assembler
support files.
The solution is to configure GNU MP with ABI=32.
INSTALLATION INSTRUCTIONS:
Modify Makefile and the beginning of config.h as necessary.
Use GNU Make. BSD Make won't work. On BSD systems, GNU Make is often
available as "gmake".
This diff is collapsed.
# Use GNU Make !
RANLIB= ranlib
OCAML_LIBDIR:= $(shell ocamlc -where)
GMP_INCLUDES= -I/opt/gmp/include -I/users/absint2/local/include -I$(HOME)/packages/gmp/include
GMP_LIBDIR=/opt/gmp/lib
DESTDIR= $(OCAML_LIBDIR)/gmp
RLIBFLAGS= -cclib "-Wl,-rpath $(GMP_LIBDIR)" # Linux, FreeBSD
#RLIBFLAGS= -cclib "-Wl,-R $(GMP_LIBDIR)" # Solaris
# RLIBFLAGS= # MacOS X
LIBFLAGS= -cclib -L. -cclib -L$(GMP_LIBDIR) $(RLIBFLAGS) \
-cclib -lmpfr -cclib -lgmp -cclib -L$(DESTDIR)
CC= gcc
CFLAGS_MISC= -Wall -Wno-unused -g -O3
#CFLAGS_MISC=
CFLAGS_INCLUDE= -I $(OCAML_LIBDIR) $(GMP_INCLUDES)
CFLAGS= $(CFLAGS_MISC) $(CFLAGS_INCLUDE)
OCAMLC= ocamlc -g
OCAMLOPT= ocamlopt
OCAMLFLAGS=
CMODULES= mlgmp_z.c mlgmp_q.c mlgmp_f.c mlgmp_fr.c mlgmp_random.c mlgmp_misc.c
CMODULES_O= $(CMODULES:%.c=%.o)
LIBS= libmlgmp.a gmp.a gmp.cma gmp.cmxa gmp.cmi
PROGRAMS= test_creal test_creal.opt essai essai.opt toplevel\
test_suite test_suite.opt
TESTS= test_suite test_suite.opt
all: $(LIBS) tests
install: all
-mkdir $(DESTDIR)
cp $(LIBS) gmp.mli $(DESTDIR)
tests: $(LIBS) $(TESTS)
./test_suite
./test_suite.opt
%.i: %.c
$(CC) $(CFLAGS) -E $*.c > $*.i
%.cmo: %.ml %.cmi
$(OCAMLC) $(OCAMLFLAGS) -c $*.ml
%.cmx: %.ml %.cmi
$(OCAMLOPT) $(OCAMLFLAGS) -c $*.ml
%.cmo: %.ml
$(OCAMLC) $(OCAMLFLAGS) -c $*.ml
%.cmx: %.ml
$(OCAMLOPT) $(OCAMLFLAGS) -c $*.ml
%.cmi: %.mli
$(OCAMLC) $(OCAMLFLAGS) -c $*.mli
$(CMODULES_O): conversions.c config.h
libmlgmp.a: $(CMODULES_O)
$(AR) -rc $@ $+
$(RANLIB) $@
gmp.cma: gmp.cmo libmlgmp.a
$(OCAMLC) $(OCAMLFLAGS) -a gmp.cmo -cclib -lmlgmp $(LIBFLAGS) -o $@
gmp.a gmp.cmxa: gmp.cmx libmlgmp.a
$(OCAMLOPT) $(OCAMLFLAGS) -a gmp.cmx -cclib -lmlgmp $(LIBFLAGS) -o $@
pretty_gmp.cmo: pretty_gmp.cmi gmp.cmo
toplevel: gmp.cma creal.cmo pretty_gmp.cmo install_pp.cmo creal_pp.cmo install_creal_pp.cmo
ocamlmktop -custom $+ -o $@
essai: gmp.cma essai.cmo
$(OCAMLC) -custom $+ -o $@
essai.opt: gmp.cmxa essai.cmx
$(OCAMLOPT) $+ -o $@
test_creal: gmp.cma creal.cmo test_creal.cmo
$(OCAMLC) -custom $+ -o $@
test_creal.opt: gmp.cmxa creal.cmx test_creal.cmx
$(OCAMLOPT) $+ -o $@
test_suite: gmp.cma test_suite.cmo
$(OCAMLC) -custom $+ -o $@
test_suite.opt: gmp.cmxa test_suite.cmx
$(OCAMLOPT) $+ -o $@
clean:
rm -f *.o *.cm* $(PROGRAMS) *.a
depend:
ocamldep *.ml *.mli > depend
.PHONY: clean
include depend
This directory contains two separate programs:
*** Creal v0.1
Copyright (C) 2000 Jean-Christophe Filliâtre.
This module consists in the files containing the name "creal" and
carrying J.C. Filliâtre's copyright.
Most algorithms are from Valérie Ménissier-Morain Ph.D. thesis
(http://www-calfor.lip6.fr/~vmm/)
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Library General Public License version 2, as
published by the Free Software Foundation.
This program 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 Library General Public License version 2 for more details
This module was adapted for correct linking with the current version
of ML GMP by D. Monniaux.
*** ML GMP 2002/07/29
Copyright (c) 2001,2002 David Monniaux
This package provides an interface between Objective Caml
(http://www.inria.fr) and
- the GNU MP (http://www.swox.com/gmp/) library
- the MPFR (http://www.mpfr.org) library
The current version is meant for
- Objective Caml 3.04
- GNU MP 4.1
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Library General Public License version 2, as
published by the Free Software Foundation, or any more recent version
published by the Free Software Foundation, at your choice.
This program 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 Library General Public License version 2 for more details
(enclosed in the file LGPL.txt).
As a special exception to the GNU Library General Public License, you
may link, statically or dynamically, a "work that uses the Library"
with a publicly distributed version of the Library to produce an
executable file containing portions of the Library, and distribute
that executable file under terms of your choice, without any of the
additional requirements listed in clause 6 of the GNU Library General
Public License. By "a publicly distributed version of the Library",
we mean either the unmodified Library as distributed by the author, or a
modified version of the Library that is distributed under the
conditions defined in clause 3 of the GNU Library General Public
License. This exception does not however invalidate any other reasons
why the executable file might be covered by the GNU Library General
Public License.
Running the "essai" and "essai.opt" programs.
INTEL ARCHITECTURE
* Intel Pentium II, 400 MHz, PC-133, Linux [orion]
essai_static: 2.55 s
* Intel Pentium II w/ 512k cache (dual processor), 450 MHz
essai_static: 2.30 s
* Intel Pentium III, 450 MHz, FreeBSD 4.2 [chaland]
essai.opt: 3.10 s
* Intel Pentium III, 600 MHz, Linux [pleiades]
essai_static: 1.63s
* Intel Pentium III, 700 MHz, OpenBSD [espie]
essai_static: 1.41s
* Intel Pentium III, 1 GHz, FreeBSD 4.4 [airelle]
essai: 1.47 s
essai.opt: 1.41 s
* Intel Celeron, 300 MHz, Linux 2.4.9/glibc 2.2.4 [quatramaran]
essai: 3.86 s
essai.opt: 4.10 s
essai_static: 3.40 s
* AMD K7, 500 MHz w/ 512k cache, Linux [efge]
essai_static: 1.92s
* AMD K7, 600 MHz, Linux [rineau]
essai_static: 1.49 s
* AMD Athlon, 700 MHz, Linux [apo]
essai_static: 1.22 s
* AMD Duron, 700 MHz, Linux [lhabert]
essai_static: 1.30 s
* AMD Duron, 750 MHz (7.5x100), Linux [said]
essai_static: 1.20 s
* AMD Athlon, 750 MHz, Linux [glouglou, Frisch]
essai_static: 1.12 s
* AMD Athlon, 850 MHz, Linux 2.4.16/glibc 2.2.4 [picsou]
essai: 1.23 s
essai.opt: 1.29 s
essai_static: 1.12 s
* AMD Athlon, 850 MHz, Linux 2.4.16/glibc 2.2.4 [picsou w/ BIOS settings]
essai_static: 1.00 s
* AMD Athlon, 900 MHz, Linux [ssecem]
essai_static: 0.95 s
SPARC ARCHITECTURE
* UltraSparc IIi, 400 MHz, Solaris 2.8 [basilic]
essai: 3.98 s
essai.opt: 4.27 s
* UltraSparc IIi, 360 MHz, Solaris 2.8 [mezcal]
essai: 4.86 s
essai.opt: 5.27 s
POWERPC ARCHITECTURE
* ppc750 "G3", 400 Mhz, MacOS X (Darwin 5.2) [mastorna]
essai: 4 s
#define SERIALIZE
#define USE_MPFR
#define NDEBUG
#undef TRACE
#include <gmp.h>
#ifdef USE_MPFR
#include <mpfr.h>
/* If you use version 20011026 of MPFR, use
#define mpfr_get_z_exp mpz_get_fr */
#endif
/* This is the largest prime less than 2^32 */
#define HASH_MODULUS 4294967291UL
#ifdef TRACE
#define trace(x) do { fprintf(stderr, "mlgmp: %s%s\n", MODULE, #x);\
fflush(stderr); } while(0)
#else
#define trace(x)
#endif
#ifdef __GNUC__
#define noreturn __attribute__((noreturn))
#else
#define noreturn
#endif
/* In C99 or recent versions of gcc,
- you can specify which field you want to initialize
- you have "inline". */
#if defined(__GNUC__) || (defined(__STDC__) && __STDC_VERSION__ >= 199901L)
#define field(x) .x =
#else
#define field(x)
#define inline
#endif
#ifdef SERIALIZE
/* Sizes of types on arch 32/ arch 64 */
/* THOSE SIZES ARE A HACK. */
/* __mpz_struct = 2*int + ptr */
#define MPZ_SIZE_ARCH32 12
#define MPZ_SIZE_ARCH64 16
/* __mpq_struct = 2 * __mpz_struct */
#define MPQ_SIZE_ARCH32 (2 * MPZ_SIZE_ARCH32)
#define MPQ_SIZE_ARCH64 (2 * MPZ_SIZE_ARCH64)
/* __mpf_struct = 3 * int + ptr */
#define MPF_SIZE_ARCH32 16
#define MPF_SIZE_ARCH64 24
/* __mpfr_struct = 3 * int + ptr */
#define MPFR_SIZE_ARCH32 16
#define MPFR_SIZE_ARCH64 24
extern void serialize_int_4(int32 i);
extern void serialize_block_1(void * data, long len);
extern uint32 deserialize_uint_4(void);
extern int32 deserialize_sint_4(void);
extern void deserialize_block_1(void * data, long len);
#endif /* SERIALIZE */
/*
* ML GMP - Interface between Objective Caml and GNU MP
* Copyright (C) 2001 David MONNIAUX
*
* This software is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License version 2 published by the Free Software Foundation,
* or any more recent version published by the Free Software
* Foundation, at your choice.
*
* This software 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 Library General Public License version 2 for more details
* (enclosed in the file LGPL).
*
* As a special exception to the GNU Library General Public License, you
* may link, statically or dynamically, a "work that uses the Library"
* with a publicly distributed version of the Library to produce an
* executable file containing portions of the Library, and distribute
* that executable file under terms of your choice, without any of the
* additional requirements listed in clause 6 of the GNU Library General
* Public License. By "a publicly distributed version of the Library",
* we mean either the unmodified Library as distributed by INRIA, or a
* modified version of the Library that is distributed under the
* conditions defined in clause 3 of the GNU Library General Public
* License. This exception does not however invalidate any other reasons
* why the executable file might be covered by the GNU Library General
* Public License.
*/
#include <assert.h>
struct custom_operations _mlgmp_custom_z;
static inline gmp_randstate_t *randstate_val(value val)
{
return ((gmp_randstate_t *) (Data_custom_val(val)));
}
static inline int Int_option_val(value val, int default_val)
{
if (val == Val_int(0)) return default_val;
return Int_val(Field(val, 0));
}
static inline mpz_t * mpz_val (value val)
{
return ((mpz_t *) (Data_custom_val(val)));
}
static inline value alloc_mpz (void)
{
return alloc_custom(&_mlgmp_custom_z,
sizeof(mpz_t),
0,
1);
}
static inline value alloc_init_mpz (void)
{
value r= alloc_mpz();
mpz_init(*mpz_val(r));
return r;
}
#pragma inline(Int_option_val, mpz_val, alloc_mpz, alloc_init_mpz)
struct custom_operations _mlgmp_custom_q;
static inline mpq_t * mpq_val (value val)
{
return ((mpq_t *) (Data_custom_val(val)));
}
static inline value alloc_mpq (void)
{
return alloc_custom(&_mlgmp_custom_q,
sizeof(mpq_t),
0,
1);
}
static inline value alloc_init_mpq (void)
{
value r= alloc_mpq();
mpq_init(*mpq_val(r));
return r;
}
#pragma inline(mpq_val, alloc_mpq, alloc_init_mpq)
struct custom_operations _mlgmp_custom_f;
static inline mpf_t * mpf_val (value val)
{
return ((mpf_t *) (Data_custom_val(val)));
}
static inline value alloc_mpf (void)
{
return alloc_custom(&_mlgmp_custom_f,
sizeof(mpf_t),
0,
1);
}
static inline value alloc_init_mpf (value prec)
{
value r= alloc_mpf();
mpf_init2(*mpf_val(r), Int_val(prec));
return r;
}
struct custom_operations _mlgmp_custom_fr;
#ifdef USE_MPFR
static inline mpfr_t * mpfr_val (value val)
{
return ((mpfr_t *) (Data_custom_val(val)));
}
static inline mp_rnd_t Mode_val (value val)
{
return (mp_rnd_t) (Int_val(val));
}
static inline value alloc_mpfr (void)
{
return alloc_custom(&_mlgmp_custom_fr,
sizeof(mpfr_t),
0,
1);
}
static inline value alloc_init_mpfr (value prec)
{
value r= alloc_mpfr();
mpfr_init2(*mpfr_val(r), Int_val(prec));
return r;
}
#endif
This diff is collapsed.
(*
* Exact real arithmetic (Constructive reals).
* Copyright (C) 2000 Jean-Christophe FILLIATRE
*
* This software is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License version 2, as published by the Free Software Foundation.
*
* This software 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 Library General Public License version 2 for more details
* (enclosed in the file LGPL).
*)
(*i $Id: creal.mli,v 1.1 2001/12/17 08:20:29 monniaux Exp $ i*)
(*s {\bf Constructive reals} are implemented by the following abstract
datatype [t]. If [x] is a constructive real, then the function call
[approx x n] returns an approximation of [x] up to $4^{-n}$, as
an arbitrary precision integer $x_n$ such that $|4^n\cdot x - x_n| < 1$. *)
open Gmp
type t
val approx : t -> int -> Z.t
(*s Basic operations. *)
val add : t -> t -> t
val neg : t -> t
val sub : t -> t -> t
val abs : t -> t
val mul : t -> t -> t
val inv : t -> t
val div : t -> t -> t
val pow_int : t -> int -> t
val sqrt : t -> t
(*s Transcendental functions. [log x y] is $\log_x(y)$. *)
val ln : t -> t
val log : t -> t -> t
val exp : t -> t
val pow : t -> t -> t
(*s Trigonometric functions. *)
val sin : t -> t
val cos : t -> t
val tan : t -> t
val arcsin : t -> t
val arccos : t -> t
val arctan : t -> t
(*s [arctan_reciproqual n] is $\arctan(1/n)$, but is more efficient than
using [arctan]. *)
val arctan_reciproqual : int -> t
(*s Hyperbolic functions. *)
val sinh : t -> t
val cosh : t -> t
val tanh : t -> t
val arcsinh : t -> t
val arccosh : t -> t
val arctanh : t -> t
(*s Some constants. *)
val zero : t
val one : t
val two : t
val pi : t
val pi_over_2 : t
val e : t
(*s Comparisons. [cmp] is absolute comparison: it may not terminate and only
returns [-1] or [+1]. [rel_cmp] is relative comparison, up to $4^{-k}$,
and it returns [-1], [0] or [+1]. *)
val cmp : t -> t -> int
val rel_cmp : int -> t -> t -> int
(*s Coercions. [to_q] and [to_float] expect a precision. [to_float x
n] returns the best floating point representation of the rational
$\ap{x}{n} / 4^n$. [of_string] expects a base as second argument. *)
val of_int : int -> t
val of_z : Z.t -> t
val of_q : Q.t -> t
val of_float : float -> t
val of_string : string -> int -> t
val to_float : t -> int -> float
val to_q : t -> int -> Q.t
(*s Coercion to type [string]. Given a decimal precision [p],
[to_string x p] returns a decimal approximation [d] of [x] with
[p] digits such that $|d - x| < 10^{-p}$. [to_beautiful_string]
returns the same decimal number but with digits packed 5 by 5. *)
val to_string : t -> int -> string
val to_beautiful_string : t -> int -> string
(*s Infix notations. *)
val ( +! ) : t -> t -> t
val ( -! ) : t -> t -> t
val ( *! ) : t -> t -> t
val ( /! ) : t -> t -> t
open Format
let precision = ref 20;;
let pp x = print_string (Creal.to_string x !precision);;
val precision : int ref
val pp : Creal.t -> unit
open Gmp;;
for a = -10 to 10
do
for b = -10 to 10
do
(if (compare
(Z.neg (Z.add (Z.from_int a) (Z.from_int b)))
(Z.from_int (- (a + b)))) <> 0
then Printf.fprintf stderr "A: %d + %d\n" a b)
done
done;;
for a = -10 to 10
do
for b = 0 to 10
do
(if (compare
(Z.neg (Z.add_ui (Z.from_int a) b))
(Z.from_int (- (a + b)))) <> 0
then Printf.fprintf stderr "B: %d + %d\n" a b)
done
done;;
for a = -100 to 100
do
for b = -100 to 100
do
if (b <> 0) then
(let (q, r) = Z.cdiv_qr (Z.from_int a) (Z.from_int b) in
(if (compare
(Z.add (Z.mul q (Z.from_int b)) r)
(Z.from_int a)) <> 0
then Printf.fprintf stderr "C: cdiv_qr %d %d\n" a b);
(let q' = Z.cdiv_q (Z.from_int a) (Z.from_int b) in
if (Z.compare q q') <> 0
then Printf.fprintf stderr "C: cdiv_q %d %d\n" a b);
(let r' = Z.cdiv_r (Z.from_int a) (Z.from_int b) in
if (Z.compare r r') <> 0
then Printf.fprintf stderr "C: cdiv_r %d %d\n" a b))
done
done;;
(if not (Z.equal (Z.mul_2exp (Z.from_int 5) 3) (Z.from_string "40"))
then Printf.fprintf stderr "E: mul_2exp\n");;
(if not (Z.equal_int (Z.fdiv_q_2exp (Z.from_string "-41") 3) ( -6))
then Printf.fprintf stderr "F: fdiv_q_2exp\n");;
(if (Z.sqrtrem (Z.pow_ui (Z.from_string_base ~base: 8 "23") 2)) <>
((Z.from_int 19), (Z.from_int 0))
then Printf.fprintf stderr "G: pow_ui, sqrtrem\n");;
(if (Z.root (Z.pow_ui (Z.from_float 17.) 3) 3) <>
(Z.from_int 17)
then Printf.fprintf stderr "H: pow_ui, root\n");;
(if (Z.perfect_power_p (Z.from_int 179199))
|| not (Z.perfect_square_p (Z.from_int 4225))
then Printf.fprintf stderr "I: perfect powers\n");;
(let a = (Z.from_int 4935) and b = (Z.from_int 2563) in
let (g, s, t) = (Z.gcdext a b)
in if g <> (Z.add (Z.mul a s) (Z.mul b t))
then Printf.fprintf stderr "J: gcdext\n");;
(let modulus = Z.from_string "4098969870980986751"
and x = Z.from_string "1657867867854181" in
let (Some y) = Z.inverse x modulus in
if (Z.modulo (Z.mul x y) modulus) <> Z.one
then Printf.fprintf stderr "K: inverse\n");;
(let prime = Z.nextprime (Z.from_string "109897328754895897328732816617973")
in if not (Z.is_prime prime)
then Printf.fprintf stderr "L: primes\n");;
(if (Z.remove (Z.from_int 16132319) (Z.from_int 7))
<> ((Z.from_int 6719), 4)
then Printf.fprintf stderr "M: remove factor\n");;
let fibo =
let rec fib a b n =
if n <= 1
then a
else fib b (a + b) (pred n)
in fib 1 1;;
(let n=14 in
if (Z.to_int (Z.fib_ui n)) <> (fibo n)
then Printf.fprintf stderr "N: Fibonacci\n");;
let fact =
let rec fac x n =
if n <= 1
then x
else fac (x * n) (pred n)
in fac 1;;
(let n=8 in
if (Z.to_int (Z.fac_ui n)) <> (fact n)
then Printf.fprintf stderr "N: Factorial\n");;
(let n=13 and k=4 in