From 305f9f200c13c74a75c9f5e13f0732e5196f3406 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Antoine=20Min=C3=A9?= Date: Tue, 21 Jan 2020 22:23:01 +0100 Subject: [PATCH] Add support for some more mpz functions adds: divisible, congruent, jacobi, legendre, krobecker, remove, fac (several versions), primorial, bin, fib, lucnum regression tests are still missing --- caml_z.c | 233 ++++++++++++++++++++++++++++++++++++++++++++++++++++--- z.mlip | 72 +++++++++++++++++ z.mlp | 14 ++++ 3 files changed, 307 insertions(+), 12 deletions(-) diff --git a/caml_z.c b/caml_z.c index 20ae46e..11fa47f 100644 --- a/caml_z.c +++ b/caml_z.c @@ -2784,6 +2784,33 @@ CAMLprim value ml_z_root(value a, value b) CAMLreturn(r); } +CAMLprim value ml_z_rootrem(value a, value b) +{ + CAMLparam2(a,b); + CAMLlocal3(r1,r2,r3); + Z_DECL(a); + mpz_t ma, mr1, mr2; + intnat mb = Long_val(b); + if (mb <= 0) + caml_invalid_argument("Z.rootrem: exponent must be positive"); + Z_ARG(a); + if (!(mb & 1) && sign_a) + caml_invalid_argument("Z.rootrem: even root of a negative number"); + ml_z_mpz_init_set_z(ma, a); + mpz_init(mr1); + mpz_init(mr2); + mpz_rootrem(mr1, mr2, ma, mb); + r1 = ml_z_from_mpz(mr1); + r2 = ml_z_from_mpz(mr2); + r3 = caml_alloc_small(2, 0); + Field(r3,0) = r1; + Field(r3,1) = r2; + mpz_clear(ma); + mpz_clear(mr1); + mpz_clear(mr2); + CAMLreturn(r3); +} + CAMLprim value ml_z_perfect_power(value a) { CAMLparam1(a); @@ -2847,19 +2874,200 @@ CAMLprim value ml_z_invert(value base, value mod) CAMLreturn(r); } +CAMLprim value ml_z_divisible(value a, value b) +{ + CAMLparam2(a,b); + mpz_t ma, mb; + int r; + ml_z_mpz_init_set_z(ma, a); + ml_z_mpz_init_set_z(mb, b); + r = mpz_divisible_p(ma, mb); + mpz_clear(ma); + mpz_clear(mb); + CAMLreturn(Val_bool(r)); +} + +CAMLprim value ml_z_congruent(value a, value b, value c) +{ + CAMLparam3(a,b,c); + mpz_t ma, mb, mc; + int r; + ml_z_mpz_init_set_z(ma, a); + ml_z_mpz_init_set_z(mb, b); + ml_z_mpz_init_set_z(mc, c); + r = mpz_congruent_p(ma, mb, mc); + mpz_clear(ma); + mpz_clear(mb); + mpz_clear(mc); + CAMLreturn(Val_bool(r)); +} + +CAMLprim value ml_z_jacobi(value a, value b) +{ + CAMLparam2(a,b); + mpz_t ma, mb; + int r; + ml_z_mpz_init_set_z(ma, a); + ml_z_mpz_init_set_z(mb, b); + r = mpz_jacobi(ma, mb); + mpz_clear(ma); + mpz_clear(mb); + CAMLreturn(Val_int(r)); +} + +CAMLprim value ml_z_legendre(value a, value b) +{ + CAMLparam2(a,b); + mpz_t ma, mb; + int r; + ml_z_mpz_init_set_z(ma, a); + ml_z_mpz_init_set_z(mb, b); + r = mpz_legendre(ma, mb); + mpz_clear(ma); + mpz_clear(mb); + CAMLreturn(Val_int(r)); +} + +CAMLprim value ml_z_kronecker(value a, value b) +{ + CAMLparam2(a,b); + mpz_t ma, mb; + int r; + ml_z_mpz_init_set_z(ma, a); + ml_z_mpz_init_set_z(mb, b); + r = mpz_kronecker(ma, mb); + mpz_clear(ma); + mpz_clear(mb); + CAMLreturn(Val_int(r)); +} + +CAMLprim value ml_z_remove(value a, value b) +{ + CAMLparam2(a,b); + CAMLlocal1(r); + mpz_t ma, mb, mr; + int i; + ml_z_mpz_init_set_z(ma, a); + ml_z_mpz_init_set_z(mb, b); + mpz_init(mr); + i = mpz_remove(mr, ma, mb); + r = caml_alloc_small(2, 0); + Field(r,0) = ml_z_from_mpz(mr); + Field(r,1) = Val_int(i); + mpz_clear(ma); + mpz_clear(mb); + mpz_clear(mr); + CAMLreturn(r); +} + +CAMLprim value ml_z_fac(value a) +{ + CAMLparam1(a); + CAMLlocal1(r); + mpz_t mr; + intnat ma = Long_val(a); + if (ma < 0) + caml_invalid_argument("Z.fac: non-positive argument"); + mpz_init(mr); + mpz_fac_ui(mr, ma); + r = ml_z_from_mpz(mr); + mpz_clear(mr); + CAMLreturn(r); +} + +CAMLprim value ml_z_fac2(value a) +{ + CAMLparam1(a); + CAMLlocal1(r); + mpz_t mr; + intnat ma = Long_val(a); + if (ma < 0) + caml_invalid_argument("Z.fac2: non-positive argument"); + mpz_init(mr); + mpz_2fac_ui(mr, ma); + r = ml_z_from_mpz(mr); + mpz_clear(mr); + CAMLreturn(r); +} + +CAMLprim value ml_z_facM(value a, value b) +{ + CAMLparam2(a,b); + CAMLlocal1(r); + mpz_t mr; + intnat ma = Long_val(a), mb = Long_val(b); + if (ma < 0 || mb < 0) + caml_invalid_argument("Z.facM: non-positive argument"); + mpz_init(mr); + mpz_mfac_uiui(mr, ma, mb); + r = ml_z_from_mpz(mr); + mpz_clear(mr); + CAMLreturn(r); +} + +CAMLprim value ml_z_primorial(value a) +{ + CAMLparam1(a); + CAMLlocal1(r); + mpz_t mr; + intnat ma = Long_val(a); + if (ma < 0) + caml_invalid_argument("Z.primorial: non-positive argument"); + mpz_init(mr); + mpz_primorial_ui(mr, ma); + r = ml_z_from_mpz(mr); + mpz_clear(mr); + CAMLreturn(r); +} + +CAMLprim value ml_z_bin(value a, value b) +{ + CAMLparam2(a,b); + CAMLlocal1(r); + mpz_t ma; + intnat mb = Long_val(b); + if (mb < 0) + caml_invalid_argument("Z.bin: non-positive argument"); + ml_z_mpz_init_set_z(ma, a); + mpz_bin_ui(ma, ma, mb); + r = ml_z_from_mpz(ma); + mpz_clear(ma); + CAMLreturn(r); +} + +CAMLprim value ml_z_fib(value a) +{ + CAMLparam1(a); + CAMLlocal1(r); + mpz_t mr; + intnat ma = Long_val(a); + if (ma < 0) + caml_invalid_argument("Z.fib: non-positive argument"); + mpz_init(mr); + mpz_fib_ui(mr, ma); + r = ml_z_from_mpz(mr); + mpz_clear(mr); + CAMLreturn(r); +} + +CAMLprim value ml_z_lucnum(value a) +{ + CAMLparam1(a); + CAMLlocal1(r); + mpz_t mr; + intnat ma = Long_val(a); + if (ma < 0) + caml_invalid_argument("Z.lucnum: non-positive argument"); + mpz_init(mr); + mpz_lucnum_ui(mr, ma); + r = ml_z_from_mpz(mr); + mpz_clear(mr); + CAMLreturn(r); +} + + + /* XXX should we support the following? - mpz_divisible_p - mpz_congruent_p - mpz_powm_sec - mpz_rootrem - mpz_jacobi - mpz_legendre - mpz_kronecker - mpz_remove - mpz_fac_ui - mpz_bin_ui - mpz_fib_ui - mpz_lucnum_ui mpz_scan0, mpz_scan1 mpz_setbit, mpz_clrbit, mpz_combit, mpz_tstbit mpz_odd_p, mpz_even_p @@ -2867,6 +3075,7 @@ CAMLprim value ml_z_invert(value base, value mod) */ + /*--------------------------------------------------- CUSTOMS BLOCKS ---------------------------------------------------*/ diff --git a/z.mlip b/z.mlip index 4416900..4c70e3f 100644 --- a/z.mlip +++ b/z.mlip @@ -192,6 +192,20 @@ external divexact: t -> t -> t = divexact@ASM Can raise a [Division_by_zero]. *) +external divisible: t -> t -> bool = "ml_z_divisible" +(** [divisible a b] returns [true] if [a] is exactly divisible by [b]. + Unlike the other division functions, [b = 0] is accepted + (only 0 is considered divisible by 0). +*) + +external congruent: t -> t -> t -> bool = "ml_z_congruent" +(** [congruent a b c] returns [true] if [a] is congruent to [b] modulo [c]. + Unlike the other division functions, [c = 0] is accepted + (only equal numbers are considered equal congruent 0). +*) + + + (** {1 Bit-level operations} *) @@ -465,6 +479,57 @@ external nextprime: t -> t = "ml_z_nextprime" The result is only prime with very high probability. *) +external jacobi: t -> t -> int = "ml_z_jacobi" +(** [jacobi a b] returns the Jacobi symbol [(a/b)]. *) + +external legendre: t -> t -> int = "ml_z_legendre" +(** [legendre a b] returns the Legendre symbol [(a/b)]. *) + +external kronecker: t -> t -> int = "ml_z_kronecker" +(** [kronecker a b] returns the Kronecker symbol [(a/b)]. *) + +external remove: t -> t -> t * int = "ml_z_remove" +(** [remove a b] returns [a] after removing all the occurences of the + factor [b]. + Also returns how many occurrences were removed. + *) + +external fac: int -> t = "ml_z_fac" +(** [fac n] returns the factorial of [n] ([n!]). + Raises an [Invaid_argument] if [n] is non-positive. +*) + +external fac2: int -> t = "ml_z_fac2" +(** [fac2 n] returns the double factorial of [n] ([n!!]). + Raises an [Invaid_argument] if [n] is non-positive. +*) + + external facM: int -> int -> t = "ml_z_facM" +(** [facM n m] returns the [m]-th factorial of [n]. + Raises an [Invaid_argument] if [n] or [m] is non-positive. +*) + +external primorial: int -> t = "ml_z_primorial" +(** [primorial n] returns the product of all positive prime numbers less + than or equal to [n]. + Raises an [Invaid_argument] if [n] is non-positive. +*) + +external bin: t -> int -> t = "ml_z_bin" +(** [bin n k] returns the binomial coefficient [n] over [k]. + Raises an [Invaid_argument] if [k] is non-positive. +*) + +external fib: int -> t = "ml_z_fib" +(** [fib n] returns the [n]-th Fibonacci number. + Raises an [Invaid_argument] if [n] is non-positive. +*) + +external lucnum: int -> t = "ml_z_lucnum" +(** [lucnum n] returns the [n]-th Lucas number. + Raises an [Invaid_argument] if [n] is non-positive. +*) + (** {1 Powers} *) @@ -493,6 +558,13 @@ external root: t -> int -> t = "ml_z_root" Otherwise, an [Invalid_argument] is raised. *) +external rootrem: t -> int -> t * t = "ml_z_rootrem" +(** [rootrem x n] computes the [n]-th root of [x] and the remainder + [x-root**n]. + [n] must be positive and, if [n] is even, then [x] must be nonnegative. + Otherwise, an [Invalid_argument] is raised. + *) + external perfect_power: t -> bool = "ml_z_perfect_power" (** True if the argument has the form [a^b], with [b>1] *) diff --git a/z.mlp b/z.mlp index 73ec52e..09377ab 100644 --- a/z.mlp +++ b/z.mlp @@ -80,6 +80,7 @@ external pow: t -> int -> t = "ml_z_pow" external powm_sec: t -> t -> t -> t = "ml_z_powm_sec" external divexact: t -> t -> t = divexact@ASM external root: t -> int -> t = "ml_z_root" +external rootrem: t -> int -> t * t = "ml_z_rootrem" external invert: t -> t -> t = "ml_z_invert" external perfect_power: t -> bool = "ml_z_perfect_power" external perfect_square: t -> bool = "ml_z_perfect_square" @@ -88,6 +89,19 @@ external nextprime: t -> t = "ml_z_nextprime" external hash: t -> int = "ml_z_hash" @NOALLOC external to_bits: t -> string = "ml_z_to_bits" external of_bits: string -> t = "ml_z_of_bits" +external divisible: t -> t -> bool = "ml_z_divisible" +external congruent: t -> t -> t -> bool = "ml_z_congruent" +external jacobi: t -> t -> int = "ml_z_jacobi" +external legendre: t -> t -> int = "ml_z_legendre" +external kronecker: t -> t -> int = "ml_z_kronecker" +external remove: t -> t -> t * int = "ml_z_remove" +external fac: int -> t = "ml_z_fac" +external fac2: int -> t = "ml_z_fac2" +external facM: int -> int -> t = "ml_z_facM" +external primorial: int -> t = "ml_z_primorial" +external bin: t -> int -> t = "ml_z_bin" +external fib: int -> t = "ml_z_fib" +external lucnum: int -> t = "ml_z_lucnum" let zero = of_int 0 let one = of_int 1