16
$\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.

Improve this question
$\endgroup$
4
  • 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, 2015 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$
    – kmingulov
    Aug 24, 2015 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, 2015 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$
    – kmingulov
    Aug 24, 2015 at 19:00

2 Answers 2

Reset to default
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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ In this way, I can transfer reals from GMP to Mathematica. But how can I transfer reals from Mathematica to GMP? $\endgroup$
    – kmingulov
    Feb 19, 2016 at 2:37
  • $\begingroup$ See my updated response. $\endgroup$
    – Simon Byrne
    Feb 19, 2016 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.

Improve this answer
$\endgroup$
9
  • $\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, 2016 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, 2016 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, 2016 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, 2016 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, 2016 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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