Group Abstract Group Abstract

Message Boards Message Boards

Explain SetPrecision and N

Posted 26 days ago

I calculate huge matrices (of the order of 10^4 x 10^4 or more) in C++ (because Mathematica is too slow for this...) in 80-bit (long double) and export it to 128-bit (because Mathematica can only read either 64-bit or 128-bit with BinaryReadList). I import the matrices into Mathematica because it's much easier to make graphics, convergence analysis and so on. So after having loaded them into Mathematica, the numbers have mantissas of about 50 digits (the numbers on the left after the back tick `). But because they were originally calculated only in 80-bits, 50 digits are useless. In order to load them quickly, I need to convert them into compressed text. But before doing this, I would like to reduce the numbers to 20 digits, because this is 80-bit precision. This would reduce space on disk and load them more quickly. However, SetPrecision[z,p] or N[z,p] where z is a number with a mantissa of 50 decimal digits and p=20 is the desired precision does not reduce the digits at all. If I want to reduce it to 16 digits for example, I need to set p=4 or so. Why is this? Normally, when someone talks of a number z having p precision, I understand that it has a mantissa of p decimal digits. Am I missing something?

POSTED BY: Ulrich Utiger
5 Replies
Posted 26 days ago

The problem is that the numbers in the matrices are mostly very small, of the order of 10^-4000, which is also why I calculated them in C++ using long double. So Round[mantissa*10^-4000, 10^-20]=0. I could create a little function with MantissaExponent and separate the numbers in two and round only the mantissa. But this would take much more time than to create the matrices. Mathematica can handle arbitrary precision numbers but this is also slow. Speed is key here as I need to maximize the size of the matrices. I could do the same with C++ and save the file in text format, but then I had to compress the text, otherwise the file would get too voluminous. But if I am well informed, Mathematica would not be able to read the file as Compress and Uncompress use another format.

POSTED BY: Ulrich Utiger

I think guard digit retention will be a problem.

For an example we will start with a number in the ballpark you describe. All variants below seem to be retain an excess of guard digits. They are in fact (mostly) zero bits, but, alas, those still take up memory.

pi50 = N[10^(-4000)*Pi, 50];
pi50 // InputForm
(* Out[353]= 3.1415926535897932384626433832795028841971693993751*10^-4000
Out[354]//InputForm=
3.14159265358979323846264338327950288419716939937510582097494459231`50.*^-4000 *)

I'll cut back to 20 digits. The 17 is because we work here in hex and 17 hex digits exceeds 20 decimal digits.

{m, e} = MantissaExponent[pi50, 16];
pi20 = SetPrecision[Round[m, 16^(-17)], 20]*16^e
InputForm[pi20]
(* Out[358]= 3.1415926535897932384*10^-4000
Out[359]//InputForm=
3.1415926535897932384427420953578946468715017823`20.*^-4000 *)

So we got the desired 20 digits but a slew of guard digits came along for the ride. Here is method 2 and we see the same happening.

pi20b = 2^SetPrecision[Round[Log[2, pi50], 16^(-20)], 24]
InputForm[pi20b]
(* Out[360]= 3.1415926535897932385*10^-4000
Out[361]//InputForm=
3.1415926535897932384626434215219292373261840846`20.0357783006888*^-4000 *)

Method 3 uses an internal function to extract bits. The resulting mantissa remains about the same length as before due to guard digits.

{s, m, m2, e} = NumericalMath`$NumberBits[pi50];
keep = 68;
pi20c = N[FromDigits[m[[1 ;; keep]], 2], 20]*2^(e - keep)
InputForm[pi20c]
(* Out[364]= 3.1415926535897932385*10^-4000
Out[365]//InputForm=
3.1415926535897932384538450330644203682569164055`20.*^-4000 *)

Method 4.

rd = RealDigits[pi50, 16];
keep = 17;
rd2 = {rd[[1, 1 ;; keep]], rd[[2]]};
pi20d = N[FromDigits[rd2[[1]], 16]*16^(-keep + rd2[[2]]), 20]
InputForm[pi20d]
(* Out[374]= 3.1415926535897932384*10^-4000
Out[375]//InputForm=
3.1415926535897932384427420953578946468715017823`20.*^-4000 *)

And method 5.

pi20str = ResourceFunction["RealToHexString"][pi50,{15,20}]
pi20e = ResourceFunction["HexStringToReal"][pi20str, {15, 20}]
InputForm[pi20e]
(* Out[18]= "0c18ead7b2f602e3da88"
Out[19]= 3.1415926535897932384*10^-4000
Out[20]//InputForm=
3.1415926535897932384427420953578946468715017823`20.*^-4000 *)

Some days you just don't get a break.

POSTED BY: Daniel Lichtblau

SetPrecision and N do not reduce the number of digits, they just take note to disregard the digits beyond a position. You may try ResourceFunction["DecimalRound"].

POSTED BY: Gianluca Gorni

Related to Gianluca Gorni's response, you can round to a certain decimal place and then reset precision like so.

pi50 = N[Pi, 50];
InputForm[pi20 = SetPrecision[Round[pi50, 10^(-20)], 20]]
(* Out[72]//InputForm= 3.14159265358979323846`20. *)

This cuts down on the retained guard digits.

POSTED BY: Daniel Lichtblau
Posted 21 days ago

So if the problem is recognized even by a Wolfram engineer, why not introduce an option in SetPrecision or N that removes the superfluous digits instead of keeping them? By the way, I solved the problem by exporting the matrices as gzip compressed text (space separated numbers) and importing them in Mathematica with Import["filename.gz", "Table"]. This way the mantissas keep 20 digits. This takes a bit more time to save the matrices but they are loaded 3 times faster compared to BinaryReadList. In addition, space on disk is halved.

POSTED BY: Ulrich Utiger
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard