diff --git a/src/backend/nsl/nsl_diff.c b/src/backend/nsl/nsl_diff.c
index 83858da25..08626ae8a 100644
--- a/src/backend/nsl/nsl_diff.c
+++ b/src/backend/nsl/nsl_diff.c
@@ -1,487 +1,487 @@
 /***************************************************************************
     File                 : nsl_diff.c
     Project              : LabPlot
     Description          : NSL numerical differentiation functions
     --------------------------------------------------------------------
     Copyright            : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn)
  ***************************************************************************/
 
 /***************************************************************************
  *                                                                         *
  *  This program is free software; you can redistribute it and/or modify   *
  *  it under the terms of the GNU General Public License as published by   *
  *  the Free Software Foundation; either version 2 of the License, or      *
  *  (at your option) any later version.                                    *
  *                                                                         *
  *  This program is distributed in the hope that it will be useful,        *
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
  *  GNU General Public License for more details.                           *
  *                                                                         *
  *   You should have received a copy of the GNU General Public License     *
  *   along with this program; if not, write to the Free Software           *
  *   Foundation, Inc., 51 Franklin Street, Fifth Floor,                    *
  *   Boston, MA  02110-1301  USA                                           *
  *                                                                         *
  ***************************************************************************/
 
 /* TODO:
 	* add more orders
 */
 
 #include "nsl_diff.h"
 #include "nsl_common.h"
 #include "nsl_sf_poly.h"
 
 const char* nsl_diff_deriv_order_name[] = {i18n("first"), i18n("second"), i18n("third"), i18n("fourth"), i18n("fifth"), i18n("sixth")};
 
 double nsl_diff_first_central(double xm, double fm, double xp, double fp) {
 	return (fp - fm)/(xp - xm);
 }
 
 int nsl_diff_deriv_first_equal(const double *x, double *y, const size_t n) {
 	if (n < 3)
 		return -1;
 
 	double dy=0, oldy=0, oldy2=0;
 	size_t i;
 	for (i=0; i < n; i++) {
 		if (i == 0)	/* forward */
 			dy = (-y[2] + 4.*y[1] - 3.*y[0])/(x[2]-x[0]);
 		else if (i == n-1) {	/* backward */
 			y[i] = (3.*y[i] - 4.*y[i-1] + y[i-2])/(x[i]-x[i-2]);
 			y[i-1] = oldy;
 		}
 		else
 			dy = (y[i+1]-y[i-1])/(x[i+1]-x[i-1]);
 
 		if (i > 1)
 			y[i-2] = oldy2;
 		if (i > 0 && i < n-1)
 			oldy2 = oldy;
 
 		oldy = dy;
 	}
 
 	return 0;
 }
 
 int nsl_diff_first_deriv(const double *x, double *y, const size_t n, int order) {
 	switch (order) {
 	case 2:
 		return nsl_diff_first_deriv_second_order(x, y, n);
 	case 4:
 		return nsl_diff_first_deriv_fourth_order(x, y, n);
 	/*TODO: higher order */
 	default:
 		printf("nsl_diff_first_deriv() unsupported order %d\n", order);
 		return -1;
 	}
 }
 
 int nsl_diff_first_deriv_second_order(const double *x, double *y, const size_t n) {
 	if (n < 3)
 		return -1;
 
 	double dy=0, oldy=0, oldy2=0, xdata[3], ydata[3];
 	size_t i, j;
 	for (i=0; i < n; i++) {
 		if (i == 0) {
 			/* 3-point forward */
 			for (j=0; j < 3; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 			dy = nsl_sf_poly_interp_lagrange_2_deriv(x[0], xdata, ydata);
 		} else if (i == n-1) {
 			/* 3-point backward */
 			y[i] = nsl_sf_poly_interp_lagrange_2_deriv(x[i], xdata, ydata);
 			y[i-1] = oldy;
 		} else {
 			/* 3-point center */
 			for (j=0; j < 3; j++)
 				xdata[j]=x[i-1+j], ydata[j]=y[i-1+j];
 			dy = nsl_sf_poly_interp_lagrange_2_deriv(x[i], xdata, ydata);
 		}
 
 		if (i > 1)
 			y[i-2] = oldy2;
 		if (i > 0 && i < n-1)
 			oldy2 = oldy;
 
 		oldy = dy;
 	}
 
 	return 0;
 }
 
 int nsl_diff_first_deriv_fourth_order(const double *x, double *y, const size_t n) {
 	if (n < 5)
 		return -1;
 
 	double dy[5]={0}, xdata[5], ydata[5];
 	size_t i, j;
 	for (i=0; i < n; i++) {
 		if (i == 0)
 			for (j=0; j < 5; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 		else if (i > 1 && i < n-2)
 			for (j=0; j < 5; j++)
 				xdata[j]=x[i-2+j], ydata[j]=y[i-2+j];
 
 		/* 5-point rule */
 		dy[0] = nsl_sf_poly_interp_lagrange_4_deriv(x[i], xdata, ydata);
 
 		if (i == n-1)
 			for (j=0; j < 4; j++)
 				y[i-j] = dy[j];
 
 		if (i > 3)
 			y[i-4] = dy[4];
 		for (j=4; j > 0; j--)
 			if (i >= j-1)
 				dy[j] = dy[j-1];
 	}
 
 	return 0;
 }
 
 int nsl_diff_first_deriv_avg(const double *x, double *y, const size_t n) {
 	if (n < 1)
 		return -1;
 
 	size_t i;
 	double dy=0, oldy=0;
 	for (i=0; i<n; i++) {
 		if(i == 0)
 			dy = (y[1]-y[0])/(x[1]-x[0]);
 		else if (i == n-1)
 			y[i] = (y[i]-y[i-1])/(x[i]-x[i-1]);
 		else
 			dy = ( (y[i+1]-y[i])/(x[i+1]-x[i]) + (y[i]-y[i-1])/(x[i]-x[i-1]) )/2.;
 
 		if (i > 0)
 			y[i-1] = oldy;
 
 		oldy = dy;
 	}
 
 	return 0;
 }
 
 int nsl_diff_second_deriv(const double *x, double *y, const size_t n, int order) {
 	switch (order) {
 	case 1:
 		return nsl_diff_second_deriv_first_order(x, y, n);
 	case 2:
 		return nsl_diff_second_deriv_second_order(x, y, n);
 	case 3:
 		return nsl_diff_second_deriv_third_order(x, y, n);
 	/*TODO: higher order */
 	default:
 		printf("nsl_diff_second_deriv() unsupported order %d\n", order);
 		return -1;
 	}
 }
 
 int nsl_diff_second_deriv_first_order(const double *x, double *y, const size_t n) {
 	if (n < 3)
 		return -1;
 
 	double dy[3]={0}, xdata[3], ydata[3];
 	size_t i, j;
 	for (i=0; i<n; i++) {
 		if (i == 0)
 			for (j=0; j < 3; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 		else if (i > 1 && i < n-1)
 			for (j=0; j < 3; j++)
 				xdata[j]=x[i-1+j], ydata[j]=y[i-1+j];
 
 		/* 3-point rule */
 		dy[0] = nsl_sf_poly_interp_lagrange_2_deriv2(xdata, ydata);
 
 		if (i == n-1) {
 			y[i] = dy[0];
 			y[i-1] = dy[1];
 		}
 
 		if (i > 1)
 			y[i-2] = dy[2];
 		if (i > 0)
 			dy[2] = dy[1];
 
 		dy[1] = dy[0];
 	}
 
 	return 0;
 }
 
 int nsl_diff_second_deriv_second_order(const double *x, double *y, const size_t n) {
 	if (n < 4)
 		return -1;
 
 	double dy[4]={0}, xdata[4], ydata[4];
 	size_t i, j;
 	for (i=0; i<n; i++) {
 		if (i == 0) {
 			/* 4-point forward */
 			for (j=0; j < 4; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 			dy[0] = nsl_sf_poly_interp_lagrange_3_deriv2(x[i], xdata, ydata);
 		}
 		else if (i == n-1) {
 			/* 4-point backward */
 			for (j=0; j < 4; j++)
 				xdata[j]=x[i-3+j], ydata[j]=y[i-3+j];
 			y[i] = nsl_sf_poly_interp_lagrange_3_deriv2(x[i], xdata, ydata);
 			y[i-1] = dy[1];
 			y[i-2] = dy[2];
 		}
 		else {
 			/* 3-point center */
 			for (j=0; j < 3; j++)
 				xdata[j]=x[i-1+j], ydata[j]=y[i-1+j];
 			dy[0] = nsl_sf_poly_interp_lagrange_2_deriv2(xdata, ydata);
 		}
 
 		if (i > 2)
 			y[i-3] = dy[3];
 		for (j=3; j > 0; j--)
 			if (i >= j-1)
 				dy[j] = dy[j-1];
 	}
 
 	return 0;
 }
 
 int nsl_diff_second_deriv_third_order(const double *x, double *y, const size_t n) {
 	if (n < 5)
 		return -1;
 
 	double dy[5]={0}, xdata[5], ydata[5];
 	size_t i, j;
 	for (i=0; i < n; i++) {
 		if (i == 0)
 			for (j=0; j < 5; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 		else if (i > 2 && i < n-3)
 			for (j=0; j < 5; j++)
 				xdata[j]=x[i-2+j], ydata[j]=y[i-2+j];
 
 		/* 5-point rule */
 		dy[0] = nsl_sf_poly_interp_lagrange_4_deriv2(x[i], xdata, ydata);
 
 		if (i == n-1)
 			for (j=0; j < 4; j++)
 				y[i-j] = dy[j];
 
 		if (i > 3)
 			y[i-4] = dy[4];
 		for (j=4; j > 0; j--)
 			if (i >= j-1)
 				dy[j] = dy[j-1];
 	}
 
 	return 0;
 }
 
 int nsl_diff_third_deriv(const double *x, double *y, const size_t n, int order) {
 	switch (order) {
 	case 2:
 		return nsl_diff_third_deriv_second_order(x, y, n);
 	/*TODO: higher order */
 	default:
 		printf("nsl_diff_third_deriv() unsupported order %d\n", order);
 		return -1;
 	}
 }
 
 int nsl_diff_third_deriv_second_order(const double *x, double *y, const size_t n) {
 	if (n < 5)
 		return -1;
 
 	double dy[5]={0}, xdata[5], ydata[5];
 	size_t i, j;
 	for (i=0; i < n; i++) {
 		if (i == 0)
 			for (j=0; j < 5; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 		else if (i > 2 && i < n-3)
 			for (j=0; j < 5; j++)
 				xdata[j]=x[i-2+j], ydata[j]=y[i-2+j];
 
 		/* 5-point rule */
 		dy[0] = nsl_sf_poly_interp_lagrange_4_deriv3(x[i], xdata, ydata);
 
 		if (i == n-1)
 			for (j=0; j < 4; j++)
 				y[i-j] = dy[j];
 
 		if (i > 3)
 			y[i-4] = dy[4];
 		for (j=4; j > 0; j--)
 			if (i >= j-1)
 				dy[j] = dy[j-1];
 	}
 
 	return 0;
 }
 
 int nsl_diff_fourth_deriv(const double *x, double *y, const size_t n, int order) {
 	switch (order) {
 	case 1:
 		return nsl_diff_fourth_deriv_first_order(x, y, n);
 	case 3:
 		return nsl_diff_fourth_deriv_third_order(x, y, n);
 	/*TODO: higher order */
 	default:
 		printf("nsl_diff_fourth_deriv() unsupported order %d\n", order);
 		return -1;
 	}
 }
 
 int nsl_diff_fourth_deriv_first_order(const double *x, double *y, const size_t n) {
 	if (n < 5)
 		return -1;
 
 	double dy[5]={0}, xdata[5], ydata[5];
 	size_t i, j;
 	for (i=0; i < n; i++) {
 		if (i == 0)
 			for (j=0; j < 5; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 		else if (i > 2 && i < n-3)
 			for (j=0; j < 5; j++)
 				xdata[j]=x[i-2+j], ydata[j]=y[i-2+j];
 
 		/* 5-point rule */
 		dy[0] = nsl_sf_poly_interp_lagrange_4_deriv4(xdata, ydata);
 
 		if (i == n-1)
 			for (j=0; j < 4; j++)
 				y[i-j] = dy[j];
 
 		if (i > 3)
 			y[i-4] = dy[4];
 		for (j=4; j > 0; j--)
 			if (i >= j-1)
 				dy[j] = dy[j-1];
 	}
 
 	return 0;
 }
 
 int nsl_diff_fourth_deriv_third_order(const double *x, double *y, const size_t n) {
 	if (n < 7)
 		return -1;
 
-	double dy[7], xdata[7], ydata[7];
+	double dy[7]={0}, xdata[7], ydata[7];
 	size_t i, j;
 	for (i=0; i < n; i++) {
 		if (i == 0)
 			for (j=0; j < 7; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 		else if (i > 3 && i < n-4)
 			for (j=0; j < 7; j++)
 				xdata[j]=x[i-3+j], ydata[j]=y[i-3+j];
 
 		/* 7-point rule */
 		dy[0] = nsl_sf_poly_interp_lagrange_6_deriv4(x[i], xdata, ydata);
 
 		if (i == n-1)
 			for (j=0; j < 6; j++)
 				y[i-j] = dy[j];
 
 		if (i > 5)
 			y[i-6] = dy[6];
 		for (j=6; j > 0; j--)
 			if (i >= j-1)
 				dy[j] = dy[j-1];
 	}
 
 	return 0;
 }
 
 int nsl_diff_fifth_deriv(const double *x, double *y, const size_t n, int order) {
 	switch (order) {
 	case 2:
 		return nsl_diff_fifth_deriv_second_order(x, y, n);
 	/*TODO: higher order */
 	default:
 		printf("nsl_diff_fifth_deriv() unsupported order %d\n", order);
 		return -1;
 	}
 }
 
 int nsl_diff_fifth_deriv_second_order(const double *x, double *y, const size_t n) {
 	if (n < 7)
 		return -1;
 
 	double dy[7]={0}, xdata[7], ydata[7];
 	size_t i, j;
 	for (i=0; i < n; i++) {
 		if (i == 0)
 			for (j=0; j < 7; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 		else if (i > 3 && i < n-4)
 			for (j=0; j < 7; j++)
 				xdata[j]=x[i-3+j], ydata[j]=y[i-3+j];
 
 		/* 7-point rule */
 		dy[0] = nsl_sf_poly_interp_lagrange_6_deriv5(x[i], xdata, ydata);
 
 		if (i == n-1)
 			for (j=0; j < 6; j++)
 				y[i-j] = dy[j];
 
 		if (i > 5)
 			y[i-6] = dy[6];
 		for (j=6; j > 0; j--)
 			if (i >= j-1)
 				dy[j] = dy[j-1];
 	}
 
 	return 0;
 }
 
 int nsl_diff_sixth_deriv(const double *x, double *y, const size_t n, int order) {
 	switch (order) {
 	case 1:
 		return nsl_diff_sixth_deriv_first_order(x, y, n);
 	/*TODO: higher order */
 	default:
 		printf("nsl_diff_sixth_deriv() unsupported order %d\n", order);
 		return -1;
 	}
 }
 
 int nsl_diff_sixth_deriv_first_order(const double *x, double *y, const size_t n) {
 	if (n < 7)
 		return -1;
 
 	double dy[7]={0}, xdata[7], ydata[7];
 	size_t i, j;
 	for (i=0; i < n; i++) {
 		if (i == 0)
 			for (j=0; j < 7; j++)
 				xdata[j]=x[j], ydata[j]=y[j];
 		else if (i > 3 && i < n-4)
 			for (j=0; j < 7; j++)
 				xdata[j]=x[i-3+j], ydata[j]=y[i-3+j];
 
 		/* 7-point rule */
 		dy[0] = nsl_sf_poly_interp_lagrange_6_deriv6(xdata, ydata);
 
 		if (i == n-1)
 			for (j=0; j < 6; j++)
 				y[i-j] = dy[j];
 
 		if (i > 5)
 			y[i-6] = dy[6];
 		for (j=6; j > 0; j--)
 			if (i >= j-1)
 				dy[j] = dy[j-1];
 	}
 
 	return 0;
 }
 
diff --git a/src/backend/nsl/nsl_int.c b/src/backend/nsl/nsl_int.c
index 84d131944..c2d4b17a1 100644
--- a/src/backend/nsl/nsl_int.c
+++ b/src/backend/nsl/nsl_int.c
@@ -1,155 +1,155 @@
 /***************************************************************************
     File                 : nsl_int.c
     Project              : LabPlot
     Description          : NSL numerical integration functions
     --------------------------------------------------------------------
     Copyright            : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn)
  ***************************************************************************/
 
 /***************************************************************************
  *                                                                         *
  *  This program is free software; you can redistribute it and/or modify   *
  *  it under the terms of the GNU General Public License as published by   *
  *  the Free Software Foundation; either version 2 of the License, or      *
  *  (at your option) any later version.                                    *
  *                                                                         *
  *  This program is distributed in the hope that it will be useful,        *
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
  *  GNU General Public License for more details.                           *
  *                                                                         *
  *   You should have received a copy of the GNU General Public License     *
  *   along with this program; if not, write to the Free Software           *
  *   Foundation, Inc., 51 Franklin Street, Fifth Floor,                    *
  *   Boston, MA  02110-1301  USA                                           *
  *                                                                         *
  ***************************************************************************/
 
 /* TODO:
 	* absolute area for Simpson/Simpson-3/8 rules (needs more numerics)
 */
 
 #include "nsl_int.h"
 #include "nsl_common.h"
 #include "nsl_sf_poly.h"
 
 
 const char* nsl_int_method_name[] = {i18n("rectangle (1-point)"), i18n("trapezoid (2-point)"), i18n("Simpson's (3-point)"), i18n("Simpson's 3/8 (4-point)")};
 
 int nsl_int_rectangle(const double *x, double *y, const size_t n, int abs) {
 	if (n == 0)
 		return -1;
 
 	size_t i, j;
 	double s, sum=0, xdata[2];
 	for (i = 0; i < n-1; i++) {
 		for (j=0; j < 2; j++)
 			xdata[j] = x[i+j];
 		s = nsl_sf_poly_interp_lagrange_0_int(xdata, y[i]);
 		if (abs)
 			s = fabs(s);
 		y[i] = sum;
 		sum += s;
 	}
 	y[n-1] = sum;
 
 	return 0;
 }
 
 int nsl_int_trapezoid(const double *x, double *y, const size_t n, int abs) {
 	if (n < 2)
 		return -1;
 
 	size_t i, j;
 	double sum=0, xdata[2], ydata[2];
 	for (i = 0; i < n-1; i++) {
 		for (j=0; j < 2; j++)
 			xdata[j] = x[i+j], ydata[j] = y[i+j];
 		y[i] = sum;
 		if (abs)
 			sum += nsl_sf_poly_interp_lagrange_1_absint(xdata, ydata);
 		else
 			sum += nsl_sf_poly_interp_lagrange_1_int(xdata, ydata);
 	}
 	y[n-1] = sum;
 
 	return 0;
 }
 
 size_t nsl_int_simpson(double *x, double *y, const size_t n, int abs) {
 	if (n < 3)
 		return 0;
 	if (abs != 0) {
 		printf("absolute area Simpson rule not implemented yet.\n");
 		return 0;
 	}
 
 	size_t i, j, np=1;
 	double sum=0, xdata[3], ydata[3];
 	for (i = 0; i < n-2; i+=2) {
 		for (j=0; j < 3; j++)
 			xdata[j] = x[i+j], ydata[j] = y[i+j];
 
 		sum += nsl_sf_poly_interp_lagrange_2_int(xdata, ydata);
 		y[np] = sum;
 		x[np++] = (x[i]+x[i+1]+x[i+2])/3.;
 		/*printf("i/sum: %zu-%zu %g\n", i, i+2, sum);*/
 	}
 
 	/* handle possible last point: use trapezoid rule */
 	if (i == n-2) {
 		for (j=0; j < 2; j++)
 			xdata[j] = x[i+j], ydata[j] = y[i+j];
 		sum += nsl_sf_poly_interp_lagrange_1_int(xdata, ydata);
 		y[np] = sum;
 		x[np++] = x[i];
 	}
 
 	/* first point */
 	y[0]=0;
 
 	return np;
 }
 
 size_t nsl_int_simpson_3_8(double *x, double *y, const size_t n, int abs) {
 	if (n < 4) {
-		printf("minimum number of points is 4 (given %zu).\n", n);
+		printf("minimum number of points is 4 (given %d).\n", (int)n);
 		return 0;
 	}
 	if (abs != 0) {
 		printf("absolute area Simpson 3/8 rule not implemented yet.\n");
 		return 0;
 	}
 
 	size_t i, j, np=1;
 	double sum=0, xdata[4], ydata[4];
 	for (i = 0; i < n-3; i+=3) {
 		for (j=0; j < 4; j++)
 			xdata[j] = x[i+j], ydata[j] = y[i+j];
 
 		sum += nsl_sf_poly_interp_lagrange_3_int(xdata, ydata);
 		y[np] = sum;
 		x[np++] = (x[i]+x[i+1]+x[i+2]+x[i+3])/4.;
 		/*printf("i/sum: %zu-%zu %g\n", i, i+3, sum);*/
 	}
 
 	/* handle possible last point(s): use trapezoid (one point) or simpson rule (two points) */
 	if (i == n-2) {
 		for (j=0; j < 2; j++)
 			xdata[j] = x[i+j], ydata[j] = y[i+j];
 		sum += nsl_sf_poly_interp_lagrange_1_int(xdata, ydata);
 		y[np] = sum;
 		x[np++] = x[i];
 	} else if ( i == n-3) {
 		for (j=0; j < 3; j++)
 			xdata[j] = x[i+j], ydata[j] = y[i+j];
 		sum += nsl_sf_poly_interp_lagrange_2_int(xdata, ydata);
 		y[np] = sum;
 		x[np++] = (x[i]+x[i+1]+x[i+2])/3.;
 	}
 
 	/* first point */
 	y[0]=0;
 
 	return np;
 }
diff --git a/src/backend/nsl/nsl_smooth.c b/src/backend/nsl/nsl_smooth.c
index 64c053345..ac2bae9ed 100644
--- a/src/backend/nsl/nsl_smooth.c
+++ b/src/backend/nsl/nsl_smooth.c
@@ -1,547 +1,547 @@
 /***************************************************************************
     File                 : nsl_smooth.c
     Project              : LabPlot
     Description          : NSL smooth functions
     --------------------------------------------------------------------
     Copyright            : (C) 2010 by Knut Franke (knut.franke@gmx.de)
     Copyright            : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn)
  ***************************************************************************/
 
 /***************************************************************************
  *                                                                         *
  *  This program is free software; you can redistribute it and/or modify   *
  *  it under the terms of the GNU General Public License as published by   *
  *  the Free Software Foundation; either version 2 of the License, or      *
  *  (at your option) any later version.                                    *
  *                                                                         *
  *  This program is distributed in the hope that it will be useful,        *
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
  *  GNU General Public License for more details.                           *
  *                                                                         *
  *   You should have received a copy of the GNU General Public License     *
  *   along with this program; if not, write to the Free Software           *
  *   Foundation, Inc., 51 Franklin Street, Fifth Floor,                    *
  *   Boston, MA  02110-1301  USA                                           *
  *                                                                         *
  ***************************************************************************/
 
 #include "nsl_smooth.h"
 #include "nsl_common.h"
 #include "nsl_sf_kernel.h"
 #include "nsl_stats.h"
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_sf_gamma.h>   /* gsl_sf_choose */
 
 const char* nsl_smooth_type_name[] = { i18n("moving average (central)"), i18n("moving average (lagged)"), i18n("percentile"), i18n("Savitzky-Golay") };
 const char* nsl_smooth_pad_mode_name[] = { i18n("none"), i18n("interpolating"), i18n("mirror"), i18n("nearest"), i18n("constant"), i18n("periodic") };
 const char* nsl_smooth_weight_type_name[] = { i18n("uniform (rectangular)"), i18n("triangular"), i18n("binomial"), i18n("parabolic (Epanechnikov)"),
 		i18n("quartic (biweight)"), i18n("triweight"), i18n("tricube"), i18n("cosine")  };
 double nsl_smooth_pad_constant_lvalue = 0.0, nsl_smooth_pad_constant_rvalue = 0.0;
 
 int nsl_smooth_moving_average(double *data, size_t n, size_t points, nsl_smooth_weight_type weight, nsl_smooth_pad_mode mode) {
 	size_t i, j;
 	double *result = (double *)malloc(n*sizeof(double));
 	for (i = 0; i < n; i++)
 		result[i] = 0;
 
 	for (i = 0; i < n; i++) {
 		size_t np = points;
 		size_t half = (points-1)/2;
 		if (mode == nsl_smooth_pad_none) { /* reduce points */
 			half = GSL_MIN(GSL_MIN((points-1)/2, i), n-i-1);
 			np = 2 * half + 1;
 		}
 
 		/* weight */
 		double sum = 0.0;
 		double* w = (double*)malloc(np * sizeof(double));
 		switch (weight) {
 		case nsl_smooth_weight_uniform:
 			for (j = 0; j < np; j++)
 				w[j] = 1./np;
 			break;
 		case nsl_smooth_weight_triangular:
 			sum = gsl_pow_2((np + 1)/2);
 			for (j = 0; j < np; j++)
 				w[j] = GSL_MIN(j + 1, np - j)/sum;
 			break;
 		case nsl_smooth_weight_binomial:
 			sum = (np - 1)/2.;
 			for (j = 0; j < np; j++)
 				w[j] = gsl_sf_choose((unsigned int)(2 * sum), (unsigned int)((sum + fabs(j - sum))/pow(4., sum)));
 			break;
 		case nsl_smooth_weight_parabolic:
 			for (j = 0; j < np; j++) {
 				w[j] = nsl_sf_kernel_parabolic(2.*(j-(np-1)/2.)/(np+1));
 				sum += w[j];
 			}
 			for (j = 0; j < np; j++)
 				w[j] /= sum;
 			break;
 		case nsl_smooth_weight_quartic:
 			for (j = 0; j < np; j++) {
 				w[j]=nsl_sf_kernel_quartic(2.*(j-(np-1)/2.)/(np+1));
 				sum += w[j];
 			}
 			for (j = 0; j < np; j++)
 				w[j] /= sum;
 			break;
 		case nsl_smooth_weight_triweight:
 			for (j = 0; j < np; j++) {
 				w[j] = nsl_sf_kernel_triweight(2.*(j-(np-1)/2.)/(np+1));
 				sum += w[j];
 			}
 			for (j = 0 ; j < np; j++)
 				w[j] /= sum;
 			break;
 		case nsl_smooth_weight_tricube:
 			for (j = 0; j < np; j++) {
 				w[j] = nsl_sf_kernel_tricube(2.*(j-(np-1)/2.)/(np+1));
 				sum += w[j];
 			}
 			for (j = 0 ; j < np; j++)
 				w[j] /= sum;
 			break;
 		case nsl_smooth_weight_cosine:
 			for (j = 0; j < np; j++) {
 				w[j] = nsl_sf_kernel_cosine((j-(np-1)/2.)/((np+1)/2.));
 				sum += w[j];
 			}
 			for (j = 0 ; j < np; j++)
 				w[j] /= sum;
 			break;
 		}
 
 		/*printf("(%d) w:",i);
 		for(j=0;j<np;j++)
 			printf(" %g",w[j]);
 		printf(" (half=%d) index = ",half);*/
 
 		/* calculate weighted average */
 		for (j = 0; j < np; j++) {
 			int index = (int)(i-half+j);
 			switch(mode) {
 			case nsl_smooth_pad_none:
 				result[i] += w[j]*data[index];
 				break;
 			case nsl_smooth_pad_interp:
 				printf("not implemented yet\n");
 				break;
 			case nsl_smooth_pad_mirror:
 				index = abs((int)(i-half+j));
 				/*printf(" %d",GSL_MIN(index,2*(n-1)-index));*/
 				result[i] += w[j] * data[GSL_MIN(index, 2*((int)n-1)-index)];
 				break;
 			case nsl_smooth_pad_nearest:
 				/*printf(" %d",GSL_MIN(n-1,GSL_MAX(0,index)));*/
 				result[i] += w[j] * data[GSL_MIN((int)n-1, GSL_MAX(0, index))];
 				break;
 			case nsl_smooth_pad_constant:
 				if (index < 0)
 					result[i] += w[j]*nsl_smooth_pad_constant_lvalue;
 				else if (index > (int)n - 1)
 					result[i] += w[j]*nsl_smooth_pad_constant_rvalue;
 				else
 					result[i] += w[j]*data[index];
 				break;
 			case nsl_smooth_pad_periodic:
 				if (index < 0)
 					index = n+index;
 				else if (index > (int)n - 1)
 					index = index-n;
 				result[i] += w[j] * data[index];
 				break;
 			}
 		}
 		/*puts("");*/
 		free(w);
 	}
 	
 	for (i = 0; i < n; i++)
 		data[i] = result[i];
 	free(result);
 
 	return 0;
 }
 
 int nsl_smooth_moving_average_lagged(double *data, size_t n, size_t points, nsl_smooth_weight_type weight, nsl_smooth_pad_mode mode) {
 	size_t i, j;
 	double* result = (double *)malloc(n*sizeof(double));
 	for (i = 0; i < n; i++)
 		result[i] = 0;
 
 	for (i = 0; i < n; i++) {
 		size_t np = points;
 		size_t half = (points-1)/2;
 		if (mode == nsl_smooth_pad_none) { /* reduce points */
 			np = GSL_MIN(points, i+1);
 			half = np-1;
 		}
 
 		/* weight */
 		double sum = 0.0, *w = (double *)malloc(np*sizeof(double));
 		switch (weight) {
 		case nsl_smooth_weight_uniform:
 			for (j = 0; j < np; j++)
 				w[j] = 1./np;
 			break;
 		case nsl_smooth_weight_triangular:
 			sum = np*(np+1)/2;
 			for (j = 0; j < np; j++)
 				w[j] = (j+1)/sum;
 			break;
 		case nsl_smooth_weight_binomial:
 			for (j = 0; j < np; j++) {
 				w[j] = gsl_sf_choose(2*(np-1), j);
 				sum += w[j];
 			}
 			for (j = 0 ; j < np; j++)
 				w[j] /= sum;
 			break;
 		case nsl_smooth_weight_parabolic:
 			for (j = 0; j < np; j++) {
 				w[j] = nsl_sf_kernel_parabolic(1.-(1+j)/(double)np);
 				sum += w[j];
 			}
 			for (j = 0; j < np; j++)
 				w[j] /= sum;
 			break;
 		case nsl_smooth_weight_quartic:
 			for( j = 0; j < np; j++) {
 				w[j] = nsl_sf_kernel_quartic(1.-(1+j)/(double)np);
 				sum += w[j];
 			}
 			for (j = 0; j < np; j++)
 				w[j] /= sum;
 			break;
 		case nsl_smooth_weight_triweight:
 			for (j = 0; j < np; j++) {
 				w[j] = nsl_sf_kernel_triweight(1.-(1+j)/(double)np);
 				sum += w[j];
 			}
 			for (j = 0; j < np; j++)
 				w[j] /= sum;
 			break;
 		case nsl_smooth_weight_tricube:
 			for (j = 0; j < np; j++) {
 				w[j] = nsl_sf_kernel_tricube(1.-(1+j)/(double)np);
 				sum += w[j];
 			}
 			for (j = 0; j < np; j++)
 				w[j] /= sum;
 			break;
 		case nsl_smooth_weight_cosine:
 			for (j = 0; j < np; j++) {
 				w[j] = nsl_sf_kernel_cosine((np-1-j)/(double)np);
 				sum += w[j];
 			}
 			for (j = 0; j < np; j++)
 				w[j] /= sum;
 			break;
 		}
 
 		/*printf("(%d) w:",i);
 		for(j=0;j<np;j++)
 			printf(" %g",w[j]);
 		printf(" (half=%d) index = ",half);*/
 
 		/* calculate weighted average */
 		for (j = 0; j < np; j++) {
 			int index = (int)(i-np+1+j);
 			switch(mode) {
 			case nsl_smooth_pad_none:
 				result[i] += w[j]*data[i-half+j];
 				/*printf(" %d",index);*/
 				break;
 			case nsl_smooth_pad_interp:
 				printf("not implemented yet\n");
 				break;
 			case nsl_smooth_pad_mirror:
 				index=abs(index);
 				/*printf(" %d",index);*/
 				result[i] += w[j]*data[index];
 				break;
 			case nsl_smooth_pad_nearest:
 				/*printf(" %d",GSL_MAX(0,index));*/
 				result[i] += w[j]*data[index];
 				break;
 			case nsl_smooth_pad_constant:
 				if (index < 0)
 					result[i] += w[j]*nsl_smooth_pad_constant_lvalue;
 				else
 					result[i] += w[j]*data[index];
 
 				break;
 			case nsl_smooth_pad_periodic:
 				if (index < 0)
 					index += n;
 				/*printf(" %d",index);*/
 				result[i] += w[j]*data[index];
 				break;
 			}
 		}
 		/*puts("");*/
 		free(w);
 	}
 
 	for (i = 0; i < n; i++)
 		data[i] = result[i];
 	free(result);
 
 	return 0;
 }
 
 int nsl_smooth_percentile(double *data, size_t n, size_t points, double percentile, nsl_smooth_pad_mode mode) {
 	size_t i, j;
 	double *result = (double *)malloc(n * sizeof(double));
 
 	for (i = 0; i < n; i++) {
 		size_t np = points;
 		size_t half = (points-1)/2;
 		if (mode == nsl_smooth_pad_none) { /* reduce points */
 			half = GSL_MIN(GSL_MIN((points-1)/2, i), n-i-1);
 			np = 2*half+1;
 		}
 
 		double *values = (double *)malloc(np*sizeof(double));
 		for (j = 0; j < np; j++) {
 			int index = (int)(i-half+j);
 			switch(mode) {
 			case nsl_smooth_pad_none:
 				/*printf(" %d",index);*/
 				values[j] = data[index];
 				break;
 			case nsl_smooth_pad_interp:
 				printf("not implemented yet\n");
 				break;
 			case nsl_smooth_pad_mirror:
 				index=abs(index);
 				/*printf(" %d",GSL_MIN(index,2*(n-1)-index));*/
 				values[j] = data[GSL_MIN(index, 2*((int)n-1)-index)];
 				break;
 			case nsl_smooth_pad_nearest:
 				/*printf(" %d",GSL_MIN(n-1,GSL_MAX(0,index)));*/
 				values[j] = data[GSL_MIN((int)n-1, GSL_MAX(0, index))];
 				break;
 			case nsl_smooth_pad_constant:
 				if (index < 0)
 					values[j] = nsl_smooth_pad_constant_lvalue;
 				else if (index > (int)n-1)
 					values[j] = nsl_smooth_pad_constant_rvalue;
 				else
 					values[j] = data[index];
 				break;
 			case nsl_smooth_pad_periodic:
 				if (index < 0)
 					index = n+index;
 				else if (index > (int)n-1)
 					index = index-n;
 				/*printf(" %d",index);*/
 				values[j] = data[index];
 				break;
 			}
 		}
 		/*puts("");*/
 
 		/*using type 4 as default */
 		result[i] = nsl_stats_quantile(values, 1, np, percentile, nsl_stats_quantile_type4);
 		free(values);
 	}
 
 	for (i=0; i<n; i++)
 		data[i]=result[i];
 	free(result);
 
 	return 0;
 }
 
 /* taken from SciDAVis */
 int nsl_smooth_savgol_coeff(size_t points, int order, gsl_matrix *h) {
 	size_t i;
 	int j, error = 0;
 
 	/* compute Vandermonde matrix */
 	gsl_matrix *vandermonde = gsl_matrix_alloc(points, order+1);
 	for (i = 0; i < points; ++i) {
 		gsl_matrix_set(vandermonde, i, 0, 1.0);
 		for (j = 1; j <= order; ++j)
 			gsl_matrix_set(vandermonde, i, j, gsl_matrix_get(vandermonde, i, j-1) * i);
 	}
 
 	/* compute V^TV */
 	gsl_matrix *vtv = gsl_matrix_alloc(order+1, order+1);
 	error = gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, vandermonde, vandermonde, 0.0, vtv);
 
 	if (!error) {
 		/* compute (V^TV)^(-1) using LU decomposition */
 		gsl_permutation *p = gsl_permutation_alloc(order+1);
 		int signum;
 		error = gsl_linalg_LU_decomp(vtv, p, &signum);
 
 		if (!error) {
 			gsl_matrix *vtv_inv = gsl_matrix_alloc(order+1, order+1);
 			error = gsl_linalg_LU_invert(vtv, p, vtv_inv);
 			if (!error) {
 				/* compute (V^TV)^(-1)V^T */
 				gsl_matrix *vtv_inv_vt = gsl_matrix_alloc(order+1, points);
 				error = gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, vtv_inv, vandermonde, 0.0, vtv_inv_vt);
 
 				if (!error) {
 					/* finally, compute H = V(V^TV)^(-1)V^T */
 					error = gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, vandermonde, vtv_inv_vt, 0.0, h);
 				}
 				gsl_matrix_free(vtv_inv_vt);
 			}
 			gsl_matrix_free(vtv_inv);
 		}
 		gsl_permutation_free(p);
 	}
 	gsl_matrix_free(vtv);
 	gsl_matrix_free(vandermonde);
 
 	return error;
 }
 
 void nsl_smooth_pad_constant_set(double lvalue, double rvalue) {
 	nsl_smooth_pad_constant_lvalue = lvalue;
 	nsl_smooth_pad_constant_rvalue = rvalue;
 }
 
 int nsl_smooth_savgol(double *data, size_t n, size_t points, int order, nsl_smooth_pad_mode mode) {
 	size_t i, k;
 	int error = 0;
 	size_t half = (points-1)/2;	/* n//2 */
 
 	if (points > n) {
-		printf("Tried to smooth over more points (points=%zu) than given as input (%zu).", points, n);
+		printf("Tried to smooth over more points (points=%d) than given as input (%d).", (int)points, (int)n);
 		return -1;
 	}
 	if (order < 1 || (size_t)order > points-1) {
-		printf("The polynomial order must be between 1 and %zu (%d given).", points-1, order);
+		printf("The polynomial order must be between 1 and %d (%d given).", (int)(points-1), order);
 		return -2;
 	}
 
 	/* Savitzky-Golay coefficient matrix, y' = H y */
 	gsl_matrix *h = gsl_matrix_alloc(points, points);
 	error = nsl_smooth_savgol_coeff(points, order, h);
 	if (error) {
 		printf("Internal error in Savitzky-Golay algorithm:\n%s", gsl_strerror(error));
 		gsl_matrix_free(h);
 		return error;
 	}
 
 	double *result = (double *)malloc(n*sizeof(double));
 	for (i = 0; i < n; i++)
 		result[i] = 0;
 
 	/* left edge */
 	if(mode == nsl_smooth_pad_none) {
 		for (i = 0; i < half; i++) {
 			/*reduce points and order*/
 			size_t rpoints = 2*i+1, rorder = GSL_MIN((size_t)order, rpoints-GSL_MIN(rpoints, 2));
 
 			gsl_matrix *rh = gsl_matrix_alloc(rpoints, rpoints);
 			error = nsl_smooth_savgol_coeff(rpoints, rorder, rh);
 			if (error) {
 				printf("Internal error in Savitzky-Golay algorithm:\n%s",gsl_strerror(error));
 				gsl_matrix_free(rh);
 				free(result);
 				return error;
 			}
 			
 			for (k = 0; k < rpoints; k++)
 				result[i] += gsl_matrix_get(rh, i, k) * data[k];
 		}
 	} else {
 		for (i = 0; i < half; i++) {
 			for (k = 0; k < points; k++)
 				switch(mode) {
 				case nsl_smooth_pad_interp:
 					result[i] += gsl_matrix_get(h, i, k) * data[k];
 					break;
 				case nsl_smooth_pad_mirror:
 					result[i] += gsl_matrix_get(h, half, k) * data[abs((int)(k+i-half))];
 					break;
 				case nsl_smooth_pad_nearest:
 					result[i] += gsl_matrix_get(h, half, k) * data[i+k-GSL_MIN(half,i+k)];
 					break;
 				case nsl_smooth_pad_constant:
 					if (k<half-i)
 						result[i] += gsl_matrix_get(h, half, k) * nsl_smooth_pad_constant_lvalue;
 					else
 						result[i] += gsl_matrix_get(h, half, k) * data[i-half+k];
 					break;
 				case nsl_smooth_pad_periodic:
 					result[i] += gsl_matrix_get(h, half, k) * data[k<half-i?n+i+k-half:i-half+k];
 				case nsl_smooth_pad_none:
 					break;
 				}
 		}
 	}
 
 	/* central part: convolve with fixed row of h */
 	for (i = half; i < n-half; i++)
 		for (k = 0; k < points; k++)
 			result[i] += gsl_matrix_get(h, half, k) * data[i-half+k];
 
 	/* right edge */
 	if(mode == nsl_smooth_pad_none) {
 		for (i = n-half; i < n; i++) {
 			/*reduce points and order*/
 			size_t rpoints = 2*(n-1-i) + 1;
 			int rorder = (int)GSL_MIN((size_t)order, rpoints - GSL_MIN(2, rpoints));
 
 			gsl_matrix *rh = gsl_matrix_alloc(rpoints, rpoints);
 			error = nsl_smooth_savgol_coeff(rpoints, rorder, rh);
 			if (error) {
 				printf("Internal error in Savitzky-Golay algorithm:\n%s",gsl_strerror(error));
 				gsl_matrix_free(rh);
 				free(result);
 				return error;
 			}
 			
 			for (k = 0; k < rpoints; k++)
 				result[i] += gsl_matrix_get(rh, n-1-i, k) * data[n-rpoints+k];
 		}
 	} else {
 		for (i = n-half; i < n; i++) {
 			for (k = 0; k < points; k++)
 				switch(mode) {
 				case nsl_smooth_pad_interp:
 					result[i] += gsl_matrix_get(h, points-n+i, k) * data[n-points+k];
 					break;
 				case nsl_smooth_pad_mirror:
 					result[i] += gsl_matrix_get(h, half, k) * data[n-1-abs((int)(k+1+i-n-half))];
 					break;
 				case nsl_smooth_pad_nearest:
 					result[i] += gsl_matrix_get(h, half, k) * data[GSL_MIN(i-half+k,n-1)];
 					break;
 				case nsl_smooth_pad_constant:
 					if (k < n-i+half)
 						result[i] += gsl_matrix_get(h, half, k) * data[i-half+k];
 					else
 						result[i] += gsl_matrix_get(h, half, k) * nsl_smooth_pad_constant_rvalue;
 					break;
 				case nsl_smooth_pad_periodic:
 					result[i] += gsl_matrix_get(h, half, k) * data[(i-half+k) % n];
 				case nsl_smooth_pad_none:
 					break;
 				}
 		}
 	}
 
 	gsl_matrix_free(h);
 
 	for (i = 0; i < n; i++)
 		data[i] = result[i];
 	free(result);
 
 	return 0;
 }
 
 int nsl_smooth_savgol_default( double *data, size_t n, size_t points, int order) {
 	return nsl_smooth_savgol(data, n, points, order, nsl_smooth_pad_constant);
 }