Gitlab will migrate to a new storage backend starting 0300 UTC on 2020-04-04. We do not anticipate a maintenance window for this migration. Performance may be impacted over the weekend. Thanks for your patience.

Commit c5ceb603 authored by Phil Jones's avatar Phil Jones

Switch xaxis & x2axis to use Kahan summation.

This drastically reduces rounding errors for certain step sizes.
parent e8207d93
......@@ -319,15 +319,24 @@ void compute_data(void) {
// start running along x-axis -------------------------------------
int x_step_index;
double x_correction = 0;
double x_step = 0;
for (x_step_index = 0; x_step_index <= inter.x1.xsteps; x_step_index++) {
if (x_step_index > 0) {
// Kahan summation of step sizes to prevent accumulation of errors
volatile double tmp = dx - x_correction;
volatile double tmp2 = x_step + tmp;
x_correction = (tmp2 - x_step) - tmp;
x_step = tmp2;
}
if (inter.x1.xtype == FLOG) {
x_axis_point = exp(x_min + x_step_index * dx);
x_axis_point = exp(x_min + x_step);
// in order to avoid rounding error we explictly use the max value as last step
if (x_step_index == inter.x1.xsteps) {
x_axis_point = exp(x_max);
x_axis_point = exp(x_min + x_max);
}
} else {
x_axis_point = x_min + x_step_index*dx;
x_axis_point = x_min + x_step;
// in order to avoid rounding error we explictly use the max value as last step
if (x_step_index == inter.x1.xsteps) {
......@@ -372,18 +381,27 @@ void compute_data(void) {
// start running along second x-axis ----------------------
int y_step_index;
double y_correction = 0;
double y_step = 0;
for (y_step_index = 0; y_step_index <= inter.x2.xsteps; y_step_index++) {
set_progress_action_text("Calculating");
//set_progress_message_text("");
if (y_step_index > 0) {
// Kahan summation of step sizes to prevent accumulation of errors
volatile double tmp = dx - y_correction;
volatile double tmp2 = y_step + tmp;
y_correction = (tmp2 - y_step) - tmp;
y_step = tmp2;
}
if (inter.x2.xtype == FLOG) {
y_axis_point = exp(y_min + y_step_index * dy);
y_axis_point = exp(y_min + y_step);
// in order to avoid rounding error we explictly use the max value as last step
if (y_step_index == inter.x2.xsteps) {
y_axis_point = exp(y_max);
}
} else {
y_axis_point = y_min + y_step_index*dy;
y_axis_point = y_min + y_step;
// in order to avoid rounding error we explictly use the max value as last step
if (y_step_index == inter.x2.xsteps) {
y_axis_point = y_max;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment