Outputting values from CAMPARY

Question

I'm trying to use the CAMPARY library (CudA Multiple Precision ARithmetic librarY). I've downloaded the code and included it in my project. Since it supports both cpu and gpu, I'm starting with cpu to understand how it works and make sure it does what I need. But the intent is to use this with CUDA.

I'm able to instantiate an instance and assign a value, but I can't figure out how to get things back out. Consider:

#include <time.h>
#include "c:\\vss\\CAMPARY\\Doubles\\src_cpu\\multi_prec.h"

int main()
{
    const char *value = "123456789012345678901234567";

    multi_prec<2> a(value);

    a.prettyPrint();
    a.prettyPrintBin();
    a.prettyPrintBin_UnevalSum();
    char *cc = a.prettyPrintBF();
    printf("\n%s\n", cc);
    free(cc);
}

Compiles, links, runs (VS 2017). But the output is pretty unhelpful:

Prec = 2
   Data[0] = 1.234568e+26
   Data[1] = 7.486371e+08

Prec = 2
   Data[0] = 0x1.987bf7c563caap+86;
   Data[1] = 0x1.64fa5c3800000p+29;

0x1.987bf7c563caap+86 + 0x1.64fa5c3800000p+29;

1.234568e+26 7.486371e+08

Printing each of the doubles like this might be easy to do, but it doesn't tell you much about the value of the 128 number being stored. Performing highly accurate computations is of limited value if there's no way to output the results.

In addition to just printing out the value, eventually I also need to convert these numbers to ints (I'm willing to try it all in floats if there's a way to print, but I fear that both accuracy and speed will suffer). Unlike MPIR (which doesn't support CUDA), CAMPARY doesn't have any associated multi-precision int type, just floats. I can probably cobble together what I need (mostly just add/subtract/compare), but only if I can get the integer portion of CAMPARY's values back out, which I don't see a way to do.

CAMPARY doesn't seem to have any docs, so it's conceivable these capabilities are there, and I've simply overlooked them. And I'd rather ask on the CAMPARY discussion forum/mail list, but there doesn't seem to be one. That's why I'm asking here.

To sum up:

1. Is there any way to output the 128bit ( multi_prec<2> ) values from CAMPARY?
2. Is there any way to extract the integer portion from a CAMPARY multi_prec? Perhaps one of the (many) math functions in the library that I don't understand computes this?

Answer 1

There are really only 2 possible answers to this question:

1. There's another (better) multi-precision library that works on CUDA that does what you need.
2. Here's how to modify this library to do what you need.

The only people who could give the first answer are CUDA programmers. Unfortunately, if there were such a library, I feel confident talonmies would have known about it and mentioned it.

As for #2, why would anyone update this library if they weren't a CUDA programmer? There are other, much better multi-precision libraries out there. The ONLY benefit CAMPARY offers is that it supports CUDA. Which means the only people with any real motivation to work with or modify the library are CUDA programmers.

And, as the CUDA programmer with the most vested interest in solving this, I did figure out a solution (albeit an ugly one). I'm posting it here in the hopes that the information will be of value to future CAMPARY programmers. There's not much information out there for this library, so this is a start.

---

The first thing you need to understand is how CAMPARY stores its data. And, while not complex, it isn't what I expected. Coming from MPIR, I assumed that CAMPARY stored its data pretty much the same way: a fixed size exponent followed by an arbitrary number of bits for the mantissa.

But nope, CAMPARY went a different way. Looking at the code, we see:

private:
    double data[prec];

Now, I assumed that this was just an arbitrary way of reserving the number of bits they needed. But no, they really do use prec doubles. Like so:

multi_prec<8> a("2633716138033644471646729489243748530829179225072491799768019505671233074369063908765111461703117249");

    // Looking at a in the VS debugger:

    [0] 2.6337161380336443e+99  const double
    [1] 1.8496577979210756e+83  const double
    [2] 1.2618399223120249e+67  const double
    [3] -3.5978270144026257e+48 const double
    [4] -1.1764513205926450e+32 const double
    [5] -2479038053160511.0 const double
    [6] 0.00000000000000000 const double
    [7] 0.00000000000000000 const double

So, what they are doing is storing the max amount of precision possible in the first double, then the remainder is used to compute the next double and so on until they encompass the entire value, or run out of precision (dropping the least significant bits). Note that some of these are negative, which means the sum of the preceding values is a bit bigger than the actual value and they are correcting it downward.

With this in mind, we return to the question of how to print it.

In theory, you could just add all these together to get the right answer. But kinda by definition, we already know that C doesn't have a datatype to hold a value this size. But other libraries do (say MPIR). Now, MPIR doesn't work on CUDA, but it doesn't need to. You don't want to have your CUDA code printing out data. That's something you should be doing from the host anyway. So do the computations with the full power of CUDA, cudaMemcpy the results back, then use MPIR to print them out:

#define MPREC 8
void ShowP(const multi_prec<MPREC> value)
{
    multi_prec<MPREC> temp(value), temp2;

    // from mpir at mpir.org
    mpf_t mp, mp2;

    mpf_init2(mp, value.getPrec() * 64); // Make sure we reserve enough room
    mpf_init(mp2); // Only needs to hold one double.

    const double *ptr = value.getData();

    mpf_set_d(mp, ptr[0]);

    for (int x = 1; x < value.getPrec(); x++)
    {
        // MPIR doesn't have a mpf_add_d, so we need to load the value into
        // an mpf_t.
        mpf_set_d(mp2, ptr[x]);
        mpf_add(mp, mp, mp2);
    }

    // Using base 10, write the full precision (0) of mp, to stdout.
    mpf_out_str(stdout, 10, 0, mp); 

    mpf_clears(mp, mp2, NULL);
}

Used with the number stored in the multi_prec above, this outputs the exact same value. Yay.

It's not a particularly elegant solution. Having to add a second library just to print a value from the first is clearly sub-optimal. And this conversion can't be all that speedy either. But printing is typically done (much) less frequently than computing. If you do an hour's worth of computing and a handful of prints, the performance doesn't much matter. And it beats the heck out of not being able to print at all.

CAMPARY has a lot of shortcomings (undoced, unsupported, unmaintained). But for people who need mp numbers on CUDA (especially if you need sqrt), it's the best option I've found.