# Precision of Float64 in Array

#1

Always in the study of my wrapper for ACME.jl, I need to convert 1d and 2d array of float64 from julia (libjulia/c-api) to armadillo C++ (arma::vec / arma::mat).

In julia console, when I dump the array I need I get :

``````julia> model.c
1×2 Array{Float64,2}:
-9.4e-8  0.0
``````

This is my code to convert the array to armadillo matrix (assuming that type is always double precision aka Float64), the code is not optimized, more theorical :

``````static inline arma::mat make_arma_mat (jl_value_t* array_ptr)
{
arma::mat container;

if (jl_is_array(array_ptr))
{
jl_array_t* array_jl = (jl_array_t*)array_ptr;

// Get array pointer
double* raw = (double*)jl_array_ptr(array_jl);

// Get number of dimensions
int ndims = jl_array_ndims(array_jl);
if (ndims != 2)
throw std::exception("make_arma_mat: container is not a 2d-array");

// Get the size of the i-th dim
size_t nrows = jl_array_dim(array_jl, 0);
size_t ncols = jl_array_dim(array_jl, 1);

container = arma::mat(nrows, ncols);

// Fill array with data
for (size_t c = 0; c<ncols; c++)
{
for (size_t r = 0; r<nrows; r++)
{
double value = raw[r+nrows*c];
container(r,c) = value;
}
}
}

return container;
}
``````

The code work fine for me, I get a correct matrix except with the values precision.

This is the resulted armadillo matrix dump :

``````c
0        0
``````

I don’t understand why ? A simple test show that the libjulia float64 can be box and unbox with good resolution :

``````jl_init (NULL);

double cpp_64 = 1.8e-15;
jl_value_t* arg = jl_box_float64(cpp_64);
double julia_64 = jl_unbox_float64(arg);

std::cout << julia_64 << std::endl;
``````

Output :

``````    1.8e-15
``````

The array ‘`eltype`’ value is 8, I cast the raw array ptr to ‘`double`’, maybe this is not the right way to get the full precision?

#2

This looks like a problem with the precision of your computations, not the Float64 type?

#3

I just get a raw pointer on a double array (see my code on the first post) :

``````double* raw = (double*)jl_array_ptr(array_jl);
``````

And the resulting double value (without any computation) is wrong :

``````double value = raw[r+nrows*c];
``````

I don’t use the Float64 type directly.

Maybe it’s not the best way to get array values?

Note : the problem come with negative exponants values :

``````-9.4e-8
7.2e-14
1.0e-15
``````

Become `0.0000000000` but others values seems to be correct.

#4

Are you sure this isn’t a printing precision issue? See here.

#5

This is not a armadillo precision issue because the problem come inside the MSVC debugger when I breakpoint to the `double value`.

``````double i_need = -9.4e-8;
double i_get = raw[r+nrows*c];
container(r,c) = i_get;
``````

The debugger show me:

``````i_need -9.3999999999999995e-08
i_get   0.00000000000000000
``````

The problem come not only when I get double by casting a Float64 array, but even when I dump some expr (expression) and need to unbox Float64. The following code just dump an expression recursively with Symbol, Int32, Float64 …

``````static inline void dump_expr (jl_expr_t* expr, std::string indent = "")
{
if (jl_is_expr(expr))
{
std::cout << indent << jl_symbol_name(expr->head) << std::endl;
size_t nargs = jl_expr_nargs(expr);
for (int i = 0; i < nargs; ++i)
{
jl_value_t* arg = jl_exprarg(expr, i);
if (jl_is_expr(arg))
{
jl_expr_t* subexpr = (jl_expr_t*)jl_exprarg(expr, i);
if (subexpr != nullptr)
dump_expr(subexpr, indent + "    ");
}
else
if (jl_is_symbol(arg))
{
const char* subexpr = jl_symbol_name((jl_sym_t*)arg);
std::cout << indent << "Symbol: " << subexpr << std::endl;
}
else
if (jl_is_int32(arg))
{
int subexpr = jl_unbox_int32(arg);
std::cout << indent << "Int32: " << subexpr << std::endl;
}
else
if (jl_is_float64(arg))
{
double subexpr = jl_unbox_float64(arg);
std::cout << indent << "Float64: " << subexpr << std::endl;
}
else
{
std::cout << indent << std::string("NEW TYPE = ") << jl_typeof_str(arg) << std::endl;
}
}
}
}
``````

For example, I need to get 2 expression lines with Float64, this is a julia console dump :

``````res[1 + 1] = 1.8e-15 * (ex - 1) - i # C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl, line 239:
J[1 + 1,2 + 1] = 7.2e-14ex # C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl, line 240:
``````

My libjulia/c-api dump is :

``````block
line
Int32: 238
Symbol: C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl
=
ref
Symbol: res
call
Symbol: +
Int32: 1
Int32: 1
call
Symbol: -
call
Symbol: *
Float64: 0
call
Symbol: -
Symbol: ex
Int32: 1
Symbol: i
line
Int32: 239
Symbol: C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl
=
ref
Symbol: J
call
Symbol: +
Int32: 1
Int32: 1
call
Symbol: +
Int32: 2
Int32: 1
call
Symbol: *
Float64: 0
Symbol: ex
``````

The two Float64 are 0 value where the values must be, respectively:

``````Float64 1.8e-15
Float64 7.2e-14
``````

So, this is not really a problem of array or cast of array, but (in my project ofcourse) a global precision problem when the `Float64` use negative exponants because in this dump I use correctly the `jl_unbox_float64` method.

And more, when the Float64 is a positive non decimal value, I get the correct value:

``````Symbol: ex
call
Symbol: exp
call
Symbol: *
Symbol: v
Float64: 40
``````

Only the `ex = exp(v * 40.0)` part of the line dumped with the julia console :

``````let v = q[2 + 1], i = q[2 + 2], ex = exp(v * 40.0) # C:\Users\maxprod\.julia\v0.5\ACME\src\elements.jl, line 238:
``````

Really, i’m dubious!

#6

It is difficult to figure out what is going on because your examples are not self-contained and don’t show how you are calling the code from Julia.

That said, here is a minimal example which works as expected. Hopefully it will help you figure out the problem.

test.c

``````#include "julia.h"

void test(jl_value_t* array_jl) {
jl_array_t* ar = (jl_array_t*)array_jl;
size_t nrows = jl_array_dim(array_jl, 0);
size_t ncols = jl_array_dim(array_jl, 1);

double* raw = (double*)jl_array_data(array_jl);
for (size_t c = 0; c<ncols; c++)
{
for (size_t r = 0; r<nrows; r++)
{
double value = raw[r+nrows*c];
printf("%.*f\n", 20, value);
}
}
}
``````

compile with: `cc -g -I\$JL/src -I\$JL/src/support -I\$JL/usr/include test.c -otest -shared -ljulia -L\$JL/usr/lib`

test:

``````julia> a = rand(2,2); a[1] = 1.8e-15; a[2] = -9.4e-8; a[3] = 7.2e-14;

julia> ccall((:test, :test), Void, (Any,), a)
0.00000000000000180000
-0.00000009400000000000
0.00000000000007200000
0.08850464391479762050
``````

#7

Hi ihnorton. Thank you very much for your answer. You’re right and you’re code too. But when I use extended precision printf like your example, I get :

``````i_need : -0.00000009400000000000
i_get  :  0.00000000000000000000
``````

I’m sure that’s not a bug of julia ofcourse. I’ve tested some precision tips by the libjulia, no problem. But i’m sure at 99% that I use the same call in the libjulia (c-api) and the monitoring console (julia.exe). That’s why i’m dubious about having slight difference only in float64 with this kind of precision.

I will search in my code a little more, and will be back with some details of my julia call to the ACME.jl package if I don’t find solution. The only point that can generate problem is the sample period (1/sample rate) parameter. All other calls are copy/paste of the ACME example code.

One more time, I thank you very much for you help.