[m-users.] Too Slow

Peter Wang novalazy at gmail.com
Fri May 8 17:54:14 AEST 2015


Hi,

On Thu, 07 May 2015 08:20:57 +0100, matthias.guedemann at googlemail.com (Matthias Güdemann) wrote:
> > If anyone wanted to contribute a binding to libgmp that we can place
> > in the extras directory that would be okay with me *hint hint*.
> 
> I created a pull request for the inclusion of libtommath. Please review.
> 
> https://github.com/Mercury-Language/mercury/pull/24
> 

diff --git a/extras/mp_int/mp_int.m b/extras/mp_int/mp_int.m
new file mode 100644
index 0000000..7b3fccf
--- /dev/null
+++ b/extras/mp_int/mp_int.m
@@ -0,0 +1,1283 @@
+%---------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sw=4 et wm=0 tw=0
+%---------------------------------------------------------------------------%
+%
+% File: mp_int.m.
+% Main author: Matthias Güdemann.
+% Stability: low.
+%
+% This module implements a binding to libtomath.
+%
+% This library provides a portable ISO C implementation of multi precision
+% integers. libtommath is in the public domain and its source code is available
+% from https://github.com/libtom/libtommath.
+%
+% To use the provided binding, one needs the compiled library and the .h
+% include files, see README.txt for the details.
+%
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+:- module mp_int.
+
+:- interface.
+
+:- type mp_int.
+
+    % Addition.
+    %
+:- func '+'(mp_int, mp_int) = mp_int.
+
+    % Subtraction.
+    %
+:- func '-'(mp_int, mp_int) = mp_int.
+
+    % Unary minus.
+    %
+:- func '-'(mp_int) = mp_int.
+
+    % Absolute Value.
+    %
+:- func 'abs'(mp_int) = mp_int.

lowercase value

No quotes required for abs.

Maybe move this so that the basic arithmetic operators appear together.

+
+    % Multiplication.
+    %
+:- func '*'(mp_int, mp_int) = mp_int.
+
+    % Truncating integer division, e.g., (-10) // 3 = (-3).
+    %
+:- func '//'(mp_int, mp_int) = mp_int.
+:- func 'quotient'(mp_int, mp_int) = mp_int.

No quotes required for quotient.  Any reason to export it?

+
+    % Remainder.
+    % X rem Y = X - (X // Y) * Y
+    %
+:- func 'rem'(mp_int, mp_int) = mp_int.

+
+    % Squaring.
+    %
+:- func square(mp_int) = mp_int.
+
+    % Truncating integer division with remainder.
+    %
+:- pred quot_with_rem(mp_int::in, mp_int::in, mp_int::out, mp_int::out)
+    is det.

integer.m has divide_with_rem

+
+    % Multiplication by 2.
+    %
+:- pred mp_mul_2(mp_int::in, mp_int::out) is det.

mul[tiply]_by_2 ?

+
+    % Division by 2.
+    %
+:- pred mp_div_2(mp_int::in, mp_int::out) is det.

div[ide]_by_2 ?

+
+    % Shift Left.
+    %
+:- func shift_left(mp_int, int) = mp_int.
+:- func '<<'(mp_int, int) = mp_int.

left_shift would match int[eger].unchecked_left_shift

Any reason to export shift_left and not just <<?

+
+    % Shift Right.
+    %
+:- func shift_right(mp_int, int) =  mp_int.
+:- func '>>'(mp_int, int) = mp_int.
+

As for shift_left.

+    % zero(X) if X is 0.
+    %
+:- pred zero(mp_int::in) is semidet.

is_zero to match integer.is_zero

+
+    % even(X) if X is even.
+    %
+:- pred even(mp_int::in) is semidet.
+
+    % odd(X) if X is odd.
+    %
+:- pred odd(mp_int::in) is semidet.
+
+    % negative(X) if X is negative.
+    %
+:- pred negative(mp_int::in) is semidet.

is_negative to match is_zero?

+
+    % Greater than.
+    %
+:- pred '>'(mp_int::in, mp_int::in) is semidet.
+
+    % Less than.
+    %
+:- pred '<'(mp_int::in, mp_int::in) is semidet.
+
+    % Greater or equal.
+    %
+:- pred '>='(mp_int::in, mp_int::in) is semidet.
+
+    % Less or equal.
+    %
+:- pred '=<'(mp_int::in, mp_int::in) is semidet.
+
+    % Equal.
+    %
+:- pred mp_eq(mp_int::in, mp_int::in) is semidet.

equal

+
+    % Exponentiation.
+    % Throws exception `math.domain_error` is Y is negative.
+    %
+:- func pow(mp_int, mp_int) = mp_int.
+
+    % Convert mp_int to int.
+    % Fails if not inside [min_int, max_int] interval.
+    %
+:- pred to_int(mp_int::in, int::out) is semidet.
+
+    % As above, but throws exception if value is outside
+    % [min_int, max_int] interval.
+    %
+:- func det_to_int(mp_int) = int.
+
+    % Convert mp_int to a string in given base.
+    %
+:- pred to_base_string(mp_int::in, int::in, string::out) is det.

Show the argument names.  Describe when it throws an exception.

+
+    % Convert mp_int to a string in base 10.
+    %
+:- func to_string(mp_int) = string.
+
+    % Convert string in given base to mp_int. Fails if unsuccesful.
+    %
+:- pred from_string(string::in, int::in, mp_int::out) is semidet.

Show the argument names.

+
+    % Convert string in base 10 to mp_int. Fails if unsuccesful.
+    %
+:- func from_string(string) = mp_int is semidet.

unsuccessful

+
+    % As function above, exits on error/1 instead of failing.
+    %
+:- func det_from_string(string) = mp_int.

As above, but throws an exception instead of failing if unsuccessful.

+
+    % As function above, exits on error/1 instead of failing.
+    %
+:- func det_from_base_string(string, int) = mp_int.

Ditto.

+
+    % Construct mp_int from int.
+    %
+:- func mp_int(int) = mp_int.
+

Convert an int to an mp_int.

+    % Greatest common divisor.
+    %
+:- func gcd(mp_int, mp_int) = mp_int.
+
+    % Least common multiple.
+    %
+:- func lcm(mp_int, mp_int) = mp_int.
+
+    % Jacobi symbol.
+    %
+:- func jacobi(mp_int, mp_int) = int.

Describe when it throws an exception?  (I have no idea what it does.)

+
+    % Modular inverse C = A^(-1) mod B
+    %
+:- func invmod(mp_int, mp_int) = mp_int.

+
+    % Modular exponentiation A = B^C mod D.
+    %
+:- func exptmod(mp_int, mp_int, mp_int) = mp_int.
+

Is this a common abbreviation?

+    % Probabilitic primality test.
+    %
+:- pred is_prime(mp_int::in) is semidet.
+
+    % Probabilistic primality test with given number of rounds. Probability of
+    % reporting composite number for a prime is ca. (1/4)^(-#Rounds)
+    %
+:- pred is_prime(mp_int::in, int::in) is semidet.

+
+    % Square root of mp_int.
+    %
+:- pred sqrt(mp_int::in, mp_int::out) is semidet.

% sqrt(X, Sqrt) is true if Sqrt is the positive square root of X.
% Fails if X is negative.

+
+    % As above, but throws error in case of negative value.
+    %
+:- func det_sqrt(mp_int) = mp_int.
+
+    % Bitwise or.
+    %
+:- func mp_int \/ mp_int = mp_int.
+
+    % Bitwise and.
+    %
+:- func mp_int /\ mp_int = mp_int.
+
+    % Bitwise xor.
+    %
+:- func mp_int `xor` mp_int = mp_int.

I think the less common functions (gcd, lcm, jacobi, primality tests)
should appear later in the interface, perhaps here.

+
+    % Constant 0.
+    %
+:- func zero = mp_int.
+
+    % Constant 1.
+    %
+:- func one = mp_int.
+
+    % Constant 2.
+    %
+:- func two = mp_int.
+
+    % Constant -1.
+    %
+:- func minusone = mp_int.

negative_one to match integer.negative_one.  Place it before zero.

+
+    % Constant 10.
+    %
+:- func ten = mp_int.

+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+:- implementation.
+
+:- import_module bool.
+:- import_module exception.
+:- import_module int.
+:- import_module math.
+:- import_module require.
+
+%---------------------------------------------------------------------------%
+% foreign declarations
+%---------------------------------------------------------------------------%
+
+   % Type declaration for foreign type mp_int*.
+   %
+:- pragma foreign_type("C", mp_int, "mp_int*")
+    where equality is mp_eq, comparison is mp_cmp.
+:- pragma foreign_decl("C",
+                      "#include \"tommath.h\"").
+
+    % Result type to signal success or failure of external functions.
+    %
+:- type mp_result_type --->
+      mp_result_okay
+    ; mp_result_out_of_mem
+    ; mp_result_invalid_input.
+
+    % mapping of libtommath results to Mercury enum.
+:- pragma foreign_enum("C", mp_result_type/0,
+                      [
+                       mp_result_okay          - "MP_OKAY",
+                       mp_result_out_of_mem    - "MP_MEM",
+                       mp_result_invalid_input - "MP_VAL"
+                      ]).
+
+:- type mp_bool_type --->
+      mp_bool_yes
+    ; mp_bool_no.
+
+:- pragma foreign_enum("C", mp_bool_type/0,
+                      [
+                       mp_bool_yes - "MP_YES",
+                       mp_bool_no  - "MP_NO"
+                      ]).
+
+%---------------------------------------------------------------------------%
+% module internal predicate declarations
+%---------------------------------------------------------------------------%
+
+:- pred mp_add(mp_int::in, mp_int::in, mp_int::out) is det.
+:- pred mp_sub(mp_int::in, mp_int::in, mp_int::out) is det.
+:- pred mp_neg(mp_int::in, mp_int::out) is det.
+:- pred mp_abs(mp_int::in, mp_int::out) is det.
+:- pred mp_mul(mp_int::in, mp_int::in, mp_int::out) is det.
+:- pred mp_square(mp_int::in, mp_int::out) is det.

Place the declarations above the clauses.

+
+:- pred mp_init(int::in, mp_int::out) is det.
+
+mp_init(N, Res) :-
+    mp_init(N, Result, Res0),
+    ( Result = mp_result_okay ->
+        Res = Res0
+    ;
+        error("could not initialize mp_int")
+    ).
+
+:- pred mp_init(int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_init(Value::in, Result::out, Mp_Int::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  Mp_Int     = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(Mp_Int);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_set_int(Mp_Int, Value);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").

Most foreign procs could be more succinctly written, e.g.

:- pragma foreign_proc("C",
    mp_init(Value::in, Result::out, Mp_Int::out),
    [will_not_call_mercury, promise_pure, thread_safe],
"
    Mp_Int = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
    Result = mp_init(Mp_Int);
    if (Result == MP_OKAY) {
	Result = mp_set_int(Mp_Int, Value);
    }
").

+
+%---------------------------------------------------------------------------%
+% basic arithmetic
+%---------------------------------------------------------------------------%
+
+mp_add(A, B, C) :-
+    mp_add(A, B, Result, C0),
+    ( Result = mp_result_okay  ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.mp_add: could not add"))
+        )
+    ).

We prefer switches whenever possible.

mp_add(A, B, C) :-
    mp_add(A, B, Result, C0),
    (
	Result = mp_result_okay,
        C = C0
    ;
        Result = mp_result_out_of_mem,
        error("could not initialize mp_int")
    ;
	Result = mp_result_invalid_input,
        throw(math.domain_error("mp_int.mp_add: could not add"))
    ).

BTW, if-then-else chains are usually written like this without another
level of indentation:

    ( Result = mp_result_okay  ->
        C = C0
    ; Result = mp_result_out_of_mem ->
        error("could not initialize mp_int")
    ;
	throw(math.domain_error("mp_int.mp_add: could not add"))
    ).

+:- pred mp_add(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_add(A::in, B::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_add(A, B, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").

+mp_sub(A, B, C) :-
+    mp_sub(A, B, Result, C0),
+    ( Result = mp_result_okay  ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.mp_sub: could not subtract"))
+        )
+    ).
+
+:- pred mp_sub(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_sub(A::in, B::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_sub(A, B, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+mp_neg(A, C) :-
+    mp_neg(A, Result, C0),
+    ( Result = mp_result_okay ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.mp_neg: could not negate value"))
+        )
+    ).
+
+:- pred mp_neg(mp_int::in, mp_result_type::out, mp_int::out)is det.
+:- pragma foreign_proc("C",
+                      mp_neg(A::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_neg(A, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+mp_abs(A, C) :-
+    mp_abs(A, Result, C0),
+    ( Result = mp_result_okay ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.mp_abs: could not compute \
+absolute value"))
+        )
+    ).

Move the string to the next line.  Otherwise, we prefer to join shorter
strings with ++.  The compiler will join them at compile time, at -O3
and above.

+
+:- pred mp_abs(mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_abs(A::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_abs(A, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+abs(A) = Res :- mp_abs(A, Res).
+
+mp_mul(A, B, C) :-
+    mp_mul(A, B, Result, C0),
+    ( Result = mp_result_okay  ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.mp_mul: could not multiply"))
+        )
+    ).
+
+:- pred mp_mul(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_mul(A::in, B::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_mul(A, B, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+mp_mul_2(A, C) :-
+    mp_mul_2(A, Result, C0),
+    ( Result = mp_result_okay ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.mp_mul_2: could not double value"))
+        )
+    ).
+
+:- pred mp_mul_2(mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_mul_2(A::in, Result::out, B::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  B          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(B);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_mul_2(A, B);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+mp_div_2(A, C) :-
+    mp_div_2(A, Result, C0),
+    ( Result = mp_result_okay ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.mp_div_2: could not halve value"))
+        )
+    ).
+
+:- pred mp_div_2(mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_div_2(A::in, Result::out, B::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  B          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(B);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_div_2(A, B);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+quot_with_rem(A, B, Quot, Rem) :-
+    ( zero(B) ->
+        throw(math.domain_error("mp_int.quot_with_rem: division by zero"))
+    ;
+        mp_quot_rem(A, B, Result, Quot0, Rem0),
+        ( Result = mp_result_okay ->
+            Quot = Quot0,
+            Rem = Rem0
+        ;
+            ( Result = mp_result_out_of_mem ->
+                error("could not initialize mp_int")
+            ;
+                throw(math.domain_error("mp_int.quot_with_rem: could not \
+compute quotient and remainder"))
+            )
+        )
+    ).
+
+:- pred mp_quot_rem(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out,
+                   mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_quot_rem(A::in, B::in, Result::out,
+                                  Quot::out, Rem::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult1, initResult2, opResult;
+  Quot        = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  Rem         = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult1 = mp_init(Quot);
+  initResult2 = mp_init(Rem);
+  if (initResult1 == MP_OKAY && initResult2 == MP_OKAY)
+    {
+      opResult = mp_div(A, B, Quot, Rem);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult1 != MP_OKAY ? initResult1 : initResult2;
+").


+rem(A, B) = Res :-
+    ( zero(B) ->
+        throw(math.domain_error("mp_int.rem: division by zero"))
+    ;
+        mp_rem(A, B, Result, Rem0),
+        ( Result = mp_result_okay ->
+            Res = Rem0
+        ;
+            ( Result = mp_result_out_of_mem ->
+                error("could not initialize mp_int")
+            ;
+                throw(math.domain_error("mp_int.rem: could not \
+compute remainder"))
+            )
+        )
+    ).
+
+:- pred mp_rem(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_rem(A::in, B::in, Result::out, Rem::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  Rem        = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(Rem);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_div(A, B, NULL, Rem);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+quotient(A, B) = Res :-
+    ( zero(B) ->
+        throw(math.domain_error("mp_int.quotient: division by zero"))
+    ;
+        mp_quot(A, B, Result, Quot0),
+        ( Result = mp_result_okay ->
+            Res = Quot0
+        ;
+            ( Result = mp_result_out_of_mem ->
+                error("could not initialize mp_int")
+            ;
+                throw(math.domain_error("mp_int.quotient: could not \
+compute quotient"))
+            )
+        )
+    ).
+
+:- pred mp_quot(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out)
+    is det.
+:- pragma foreign_proc("C",
+                      mp_quot(A::in, B::in, Result::out, Quot::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  Quot       = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(Quot);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_div(A, B, Quot, NULL);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+mp_square(A, C) :-
+    mp_square(A, Result, C0),
+    ( Result = mp_result_okay ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.mp_square: could not square"))
+        )
+    ).
+
+:- pred mp_square(mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_square(A::in, Result::out, A_SQ::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  A_SQ       = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(A_SQ);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_sqr(A, A_SQ);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+%---------------------------------------------------------------------------%
+% conversion predicates
+%---------------------------------------------------------------------------%
+
+:- func min_int = mp_int.
+min_int = mp_int(int.min_int).
+
+:- func max_int = mp_int.
+max_int = mp_int(int.max_int).
+
+to_int(A, N) :-
+    ( ( A >= min_int, A =< max_int) ->
+        mp_to_long(A, N)
+    ;
+        fail
+    ).
+
+:- pred mp_to_long(mp_int::in, int::out)
+    is det.
+:- pragma foreign_proc("C",
+                      mp_to_long(A::in, N::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  unsigned long n;
+  n = mp_get_long(A);
+  N = n;
+").

Hmm.

+
+det_to_int(A) = Res :-
+    ( to_int(A, Res0) ->
+        Res0 = Res
+    ;
+        throw(math.domain_error("mp_int.det_to_int: not in int range"))
+    ).
+
+to_base_string(A, Radix, S) :-
+    mp_to_string(A, Radix, Result, S0),
+    ( Result = mp_result_okay ->
+        S = S0
+    ;
+        error("could not convert mp_int to string")
+    ).
+
+:- pred mp_to_string(mp_int::in, int::in, mp_result_type::out, string::out)
+    is det.
+:- pragma foreign_proc("C",
+                      mp_to_string(A::in, Radix::in, Result::out, S::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int length, opResult;
+  opResult = mp_radix_size(A, Radix, &length);
+  if (opResult == MP_OKAY)
+    {
+      MR_allocate_aligned_string_msg(S, length, MR_ALLOC_ID);
+      opResult = mp_toradix(A, S, Radix);
+      Result   = opResult;
+    }
+  else
+    Result     = opResult;
+").
+
+to_string(A) = Res :- to_base_string(A, 10, Res).
+
+from_string(S, Radix, A) :-
+    mp_from_string(S, Radix, Result, A0),
+    ( Result = mp_result_okay ->
+        A = A0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            fail
+        )
+    ).
+
+:- pred mp_from_string(string::in, int::in, mp_result_type::out, mp_int::out)
+    is det.
+:- pragma foreign_proc("C",
+                      mp_from_string(S::in, Radix::in, Result::out, A::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  A          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(A);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_read_radix(A, S, Radix);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+from_string(S) = Res :- from_string(S, 10, Res).
+
+det_from_string(S) = Res :-
+    ( from_string(S, 10, Res0) ->
+        Res = Res0
+    ;
+        error("could not create mp_int from string")
+    ).
+
+det_from_base_string(S, Base) = Res :-
+    ( from_string(S, Base, Res0) ->
+        Res = Res0
+    ;
+        error("could not create mp_int from string")
+    ).
+
+mp_int(N) = Res :-
+    mp_init(abs(N), Res0),
+    ( N < 0 ->
+        Res = -Res0
+    ;
+        Res = Res0
+    ).
+
+%---------------------------------------------------------------------------%
+% bit shifting
+%---------------------------------------------------------------------------%
+
+shift_left(A, N) = Res :-
+    mp_shift_left(A, N, Result, A0),
+    ( Result = mp_result_okay ->
+        Res = A0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.shift_left: could not shift"))
+        )
+    ).
+
+:- pred mp_shift_left(mp_int::in, int::in, mp_result_type::out, mp_int::out)
+    is det.
+:- pragma foreign_proc("C",
+                      mp_shift_left(A::in, N::in, Result::out, B::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  B          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init_copy(B, A);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_lshd(B, N);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+shift_right(A, N) = Res :-
+    mp_shift_right(A, N, Result, A0),
+    ( Result = mp_result_okay ->
+        Res = A0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.shift_right: could not shift"))
+        )
+    ).
+
+:- pred mp_shift_right(mp_int::in, int::in, mp_result_type::out, mp_int::out)
+    is det.
+:- pragma foreign_proc("C",
+                      mp_shift_right(A::in, N::in, Result::out, B::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  B          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init_copy(B, A);
+  if (initResult == MP_OKAY)
+    {
+      mp_rshd(B, N);
+      Result   = MP_OKAY;
+    }
+  else
+    Result     = initResult;
+").
+
+X << N = shift_left(X, N).
+
+X >> N = shift_right(X, N).
+
+%---------------------------------------------------------------------------%
+% test predicates
+%---------------------------------------------------------------------------%
+
+zero(A) :-
+    mp_iszero(A, Result),
+    Result = mp_bool_yes.
+
+:- pred mp_iszero(mp_int::in, mp_bool_type::out) is det.
+:- pragma foreign_proc("C",
+                      mp_iszero(A::in, Result::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int opResult;
+  opResult = mp_iszero(A);
+  Result   = opResult;
+").

You could define the test predicates directly:

:- pragma foreign_proc("C",
    zero(A::in),
    [will_not_call_mercury, promise_pure, thread_safe],
"
    SUCCESS_INDICATOR = mp_iszero(A) ? MR_TRUE : MR_FALSE;
").

+
+even(A) :-
+    mp_iseven(A, Result),
+    Result = mp_bool_yes.
+
+:- pred mp_iseven(mp_int::in, mp_bool_type::out) is det.
+:- pragma foreign_proc("C",
+                      mp_iseven(A::in, Result::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int opResult;
+  opResult = mp_iseven(A);
+  Result   = opResult;
+").
+
+odd(A) :-
+    mp_isodd(A, Result),
+    Result = mp_bool_yes.
+
+:- pred mp_isodd(mp_int::in, mp_bool_type::out) is det.
+:- pragma foreign_proc("C",
+                      mp_isodd(A::in, Result::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int opResult;
+  opResult = mp_isodd(A);
+  Result   = opResult;
+").
+
+negative(A) :-
+    mp_isneg(A, Result),
+    Result = mp_bool_yes.
+
+:- pred mp_isneg(mp_int::in, mp_bool_type::out) is det.
+:- pragma foreign_proc("C",
+                      mp_isneg(A::in, Result::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int opResult;
+  opResult = mp_isneg(A);
+  Result   = opResult;
+").
+
+%---------------------------------------------------------------------------%
+% comparison predicates
+%---------------------------------------------------------------------------%
+
+:- pred mp_cmp(comparison_result::uo, mp_int::in, mp_int::in) is det.
+
+:- pragma foreign_proc("C",
+                      mp_cmp(C::uo, A::in, B::in),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int result;
+  result = mp_cmp(A, B);
+  if (result == MP_LT)
+    C = MR_COMPARE_LESS;
+  else
+    {
+      if (result == MP_GT)
+        C = MR_COMPARE_GREATER;
+      else
+        C = MR_COMPARE_EQUAL;
+    }
+").
+
+A > B :- mp_cmp((>), A, B).
+
+A < B :- mp_cmp((<), A, B).
+
+A >= B :-
+    mp_cmp(C, A, B),
+    ( C = (>); C = (=)).
+
+A =< B :-
+    mp_cmp(C, A, B),
+    ( C = (<); C = (=)).
+
+mp_eq(A, B) :- mp_cmp((=), A, B).
+
+%---------------------------------------------------------------------------%
+
+A + B = C       :- mp_add(A, B, C).
+A - B = C       :- mp_sub(A, B, C).
+-A = C          :- mp_neg(A, C).
+A * B = C       :- mp_mul(A, B, C).
+A // B          = quotient(A, B).
+square(X) = Res :- mp_square(X, Res).
+
+%---------------------------------------------------------------------------%
+% higher level functions
+%---------------------------------------------------------------------------%
+
+pow(A, N) = Res :-
+    ( zero(N) ->
+        Res = one
+    ;
+        ( N rem two = one ->
+            Res = A * pow(A, N - one)
+        ;
+            mp_div_2(N, N0),
+            SQ = pow(A, N0),
+            mp_square(SQ, Res)
+        )
+    ).
+
+%---------------------------------------------------------------------------%
+% number theoretic functions
+%---------------------------------------------------------------------------%
+
+gcd(A, B) = Res :-
+    mp_gcd(A, B, Result, C0),
+    ( Result = mp_result_okay  ->
+        Res = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.gcd: could not compute gcd"))
+        )
+    ).
+
+:- pred mp_gcd(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_gcd(A::in, B::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_gcd(A, B, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+lcm(A, B) = Res :-
+    mp_lcm(A, B, Result, C0),
+    ( Result = mp_result_okay  ->
+        Res = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.lcm: could not compute lcm"))
+        )
+    ).
+
+:- pred mp_lcm(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_lcm(A::in, B::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_lcm(A, B, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+jacobi(A, P) = Res :-
+    mp_jacobi(A, P, Result, C0),
+    ( Result = mp_result_okay  ->
+        Res = C0
+    ;
+        throw(math.domain_error("mp_int.jacobi: could not compute Jacobi \
+symbol of mp_int"))
+    ).
+
+:- pred mp_jacobi(mp_int::in, mp_int::in, mp_result_type::out, int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_jacobi(A::in, P::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int opResult, res;
+  opResult = mp_jacobi(A, P, &res);
+  Result   = opResult;
+  C        = res;
+").
+
+invmod(A, B) = Res :-
+    mp_invmod(A, B, Result, C0),
+    ( Result = mp_result_okay  ->
+        Res = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.invmod: could not compute\
+modular inverse"))
+        )
+    ).
+
+:- pred mp_invmod(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out)
+    is det.
+:- pragma foreign_proc("C",
+                      mp_invmod(A::in, B::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_invmod(A, B, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+exptmod(A, B, C) = Res :-
+    mp_exptmod(A, B, C, Result, D0),
+    ( Result = mp_result_okay  ->
+        Res = D0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.exptmod: could not compute modular \
+exponentiation"))
+        )
+    ).
+
+:- pred mp_exptmod(mp_int::in, mp_int::in, mp_int::in, mp_result_type::out,
+                  mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_exptmod(A::in, B::in, C::in, Result::out, D::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  D          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(D);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_exptmod(A, B, C, D);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+    % Default number of rounds for Miller-Rabin primality test.
+    %
+:- func miller_rabin_rounds = int.
+miller_rabin_rounds = 40.
+
+is_prime(A) :- is_prime(A, miller_rabin_rounds).
+
+is_prime(A, Rounds) :-
+    mp_is_prime(A, Rounds, Result, PResult),
+    ( Result = mp_result_okay ->
+        PResult = 1
+    ;
+        error("could not conduct Miller-Rabin primality test on mp_int")
+    ).
+
+:- pred mp_is_prime(mp_int::in, int::in, mp_result_type::out, int::out)
+    is semidet.
+:- pragma foreign_proc("C",
+                      mp_is_prime(A::in, Rounds::in, Result::out, Value::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int opResult, result;
+  opResult = mp_prime_is_prime(A, Rounds, &result);
+  Value    = result;
+  Result   = opResult;
+").
+
+sqrt(A, Res) :-
+    ( A >= zero ->
+        mp_sqrt(A, Result, C0),
+        ( Result = mp_result_okay  ->
+            Res = C0
+        ;
+            error("could not initialize mp_int")
+        )
+    ;
+        fail
+    ).
+
+:- pred mp_sqrt(mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_sqrt(A::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_sqrt(A, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+det_sqrt(A) = Res :-
+    ( sqrt(A, Res0) ->
+        Res = Res0
+    ;
+        throw(math.domain_error("mp_int.det_sqrt: argument negative"))
+    ).
+
+%---------------------------------------------------------------------------%
+% bitwise operations
+%---------------------------------------------------------------------------%
+
+A /\ B = C :-
+    mp_and(A, B, Result, C0),
+    ( Result = mp_result_okay  ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int./\\: could not compute bitwise \
+AND"))
+        )
+    ).
+
+:- pred mp_and(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_and(A::in, B::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_and(A, B, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+A \/ B = C :-
+    mp_or(A, B, Result, C0),
+    ( Result = mp_result_okay  ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.\\/: could not compute bitwise \
+OR"))
+        )
+    ).
+
+:- pred mp_or(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_or(A::in, B::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_or(A, B, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+A `xor` B = C :-
+    mp_xor(A, B, Result, C0),
+    ( Result = mp_result_okay  ->
+        C = C0
+    ;
+        ( Result = mp_result_out_of_mem ->
+            error("could not initialize mp_int")
+        ;
+            throw(math.domain_error("mp_int.xor: could not compute bitwise \
+XOR"))
+        )
+    ).
+
+:- pred mp_xor(mp_int::in, mp_int::in, mp_result_type::out, mp_int::out) is det.
+:- pragma foreign_proc("C",
+                      mp_xor(A::in, B::in, Result::out, C::out),
+                      [will_not_call_mercury, promise_pure, thread_safe],
+"
+  int initResult, opResult;
+  C          = MR_GC_NEW_ATTRIB(mp_int, MR_ALLOC_ID);
+  initResult = mp_init(C);
+  if (initResult == MP_OKAY)
+    {
+      opResult = mp_xor(A, B, C);
+      Result   = opResult;
+    }
+  else
+    Result     = initResult;
+").
+
+%---------------------------------------------------------------------------%
+% often used constants
+%---------------------------------------------------------------------------%
+
+minusone = mp_int(-1).
+zero     = mp_int(0).
+one      = mp_int(1).
+two      = mp_int(2).
+ten      = mp_int(10).

It would be preferable to statically allocate the mp_ints for these
values.  Something like:

    :- pragma foreign_decl("C", local, "
	static mp_int constant_one;
	...
    ").

    :- initialise init/2.

    :- pred init(io::di, io::uo) is det.

    init(!IO) :-
	...

    :- pragma foreign_proc("C",
	one = (X::out),
	[ ... ],
    "
	X = &constant_one;
    ").

But it can be left for later.

Peter



More information about the users mailing list