@LuisSanz @agriggio
I have spotted a couple of numerical instabilities in the RCD implementation, or at least in the current RT version, which might lead to inf
and nan
values.
- in the lines like
float V_Stat = epssq - 18.0f * cfa[indx] * cfa[indx - w1] - …
the epssq
factor is mis-placed, and should be put either at the end of the expression or in a separate statement, something like
float V_Stat = 0.f - 18.0f * cfa[indx] * cfa[indx - w1] - …
V_Stat += epssq;
The reason is simple: if you add a very small number to a relatively large one (in this order), most likely the small number will not make any difference, because there are no bits available to represent it. Consider this minimal example:
#include <stdio.h>
int main()
{
static const float eps = 1e-5, epssq = 1e-10;
float sum1 = epssq + 0.5f - 0.5f;
float sum2 = 0.5f - 0.5f + epssq;
printf("sum1=%e sum2=%e\n", sum1, sum2);
return 0;
}
At least on my OSX system, this code generates the following output:
sum1=0.000000e+00 sum2=1.000000e-10
You can see that the first sum is identically zero…
- when computing the direction like this
VH_Dir[indx] = V_Stat / (V_Stat + H_Stat);
I encountered cases in which V_Stat = -H_Stat, and the above expression leads to a division by zero. In my code I mitigated the problem by writing
VH_Dir[indx] = V_Stat / (fabs(V_Stat) + fabs(H_Stat));
but I am not sure if this is correct, therefore I would really like to have the option of the experts…