15
$\begingroup$

I'm developing a Mathematica package, inside of which there is a number-crunching function. It uses a list of arbitrary-precision numbers as input and it takes a lot of time to do its job.

So I've tried to implement this function in C using GMP (GNU Multiple Precision Library) and OpenMP. This implementation works faster and I decided to include it into my package using WSTP. But as far as I know, one cannot send arbitrary-precision number from Mathematica into C via WSTP. So my question is:

How can I efficiently send an arbitrary-precision number from Mathematica into C and vice versa?

My first idea was to transfer a number as a char string (GMP can create reals from strings), but this is apparently not the best way.

$\endgroup$
  • 3
    $\begingroup$ "Apparently not the best way"--why? As far as I know, it is the best way (simply because it is the only direct way) to do it using MathLink. But certainly a more efficient way would be preferable, and in that sense, it is a good question. Perhaps one can try partitioning the number into a list of digits in some large base and transfer a numerical array, but I'm not sure if GMP can easily reconstruct a number represented this way. $\endgroup$ – Oleksandr R. Aug 23 '15 at 13:06
  • $\begingroup$ @OleksandrR. Well, I thought, I can represent my number as an array of bytes, where each byte encodes two digits. In C, to create a char string from this structure and feed it into GMP is not a problem. But is it possible to create such structure in Mathematica? $\endgroup$ – Kirill Mingulov Aug 24 '15 at 5:34
  • 1
    $\begingroup$ Do you mean binary-coded decimal (BCD)? I am not sure if there is any way to produce this in Mathematica. For that matter, a byte array is not so easy either. The only native arrays are of signed 64-bit integers. There is the undocumented RawArray object, but I don't know if it can be passed through MathLink, even if you can somehow represent a floating-point number in this way. $\endgroup$ – Oleksandr R. Aug 24 '15 at 13:48
  • $\begingroup$ @OleksandrR. Yes, I meant BCD. If there is no way to create such array in Mathematica, then I probably have to use strings. $\endgroup$ – Kirill Mingulov Aug 24 '15 at 19:00
4
$\begingroup$

You might be able to pass the GMP internal array of limbs directly as a list, and then call FromDigits with base = 2^(limb size). This won't require any work on behalf of GMP, though I don't know anything about Mathematica's internals to say how efficient it would be on their end.

UPDATE: For the reverse, you could call IntegerDigits.

$\endgroup$
  • $\begingroup$ In this way, I can transfer reals from GMP to Mathematica. But how can I transfer reals from Mathematica to GMP? $\endgroup$ – Kirill Mingulov Feb 19 '16 at 2:37
  • $\begingroup$ See my updated response. $\endgroup$ – Simon Byrne Feb 19 '16 at 16:39
2
$\begingroup$

I'm not sure if this answers your question but it's a bit long for a comment. One efficient approach is to convert the parts before and after radix into a base that is a power of 2 e.g. 2^16, and then process that list so each bigit ("bignum digit") is encoded as a hex string.

Here is an example.

Map[StringDelete[ToString[BaseForm[#, 16]], "\n" ~~ __] &, 
 IntegerDigits[123425342345234523466674, 2^16]]

(* Out[136]= {"1a22", "e6b6", "f80", "2a35", "abb2"} *)

Let's see what this is internally:

In[137]:= InputForm[%]
Out[137]//InputForm=
{"1a22", "e6b6", "f80", "2a35", "abb2"}

A benefit to this approach is that it scales linearly* in the size of the input, whereas use of base 10 would throw in a logarithmic factor.

* I have not tested this, but it should behave as advertised.

$\endgroup$
  • $\begingroup$ "whereas use of base 10 would throw in a logarithmic factor" <- I don't understand this comment. The number of digits required to represent a value scales logarithmically with the value itself, regardless of the base, no? What would be the advantage of this compared to Simon's suggestion, which seems to do the same but transfer the data as an integer array instead of a list of strings? $\endgroup$ – Szabolcs Feb 18 '16 at 15:56
  • $\begingroup$ This question is unrelated, but it came up while I was trying to understand how MathLink transfers various types of data. Is it true that the three MathLink functions MLGetType, MLBytesToGet and MLGetData always need to be used in this sequence and always together? If I leave out either MLGetType or MLBytesToGet then MLGetData doesn't give me data. $\endgroup$ – Szabolcs Feb 18 '16 at 15:59
  • $\begingroup$ (To clarify: I do see in the documentation that other functions can be called after MLGetType. What I am asking about: is it required to call MLGetType just before calling MLBytesToGet and MLGetData? And is it required to call MLBytesToGet before MLGetData?) $\endgroup$ – Szabolcs Feb 18 '16 at 16:03
  • $\begingroup$ What I am really curious about is: It is clear that we can send large integers or large floating point numbers through MathLink. 1. Can can we send such a value to the kernel so it comes out as something that doesn't need post-processing (a plain large integer or large real)? 2. And the reverse: if such a value is sent by the kernel, how can we identify and read it in a C program using the MathLink API? I do see that we can read it as a string, so I guess my question is: is it the only way? $\endgroup$ – Szabolcs Feb 18 '16 at 16:15
  • $\begingroup$ @Szabolcs "Size of input", when working with integers, by convention is log(integer) (use whatever base, as long as it is fixed). So the complexity I refer to for a given number n, would be O(log n) for what I showed, and O(log n log log n) if I first convert to say blocks of 10^5. The extra factor comes from explicit divide-anbd-conquer conversion of binary-based bigits to base power-of-10. $\endgroup$ – Daniel Lichtblau Feb 18 '16 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.