[goffice] regression: improve accuracy a bit more.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [goffice] regression: improve accuracy a bit more.
- Date: Sun, 16 May 2010 13:06:34 +0000 (UTC)
commit c5bd02345bfbc79f8279a685831d89ad3150a1fe
Author: Morten Welinder <terra gnome org>
Date: Sun May 16 09:06:07 2010 -0400
regression: improve accuracy a bit more.
ChangeLog | 5 +
configure.in | 2 +
goffice/math/Makefile.am | 22 +++--
goffice/math/go-quad.c | 203 ++++++++++++++++++++++++++++++++++++++++++
goffice/math/go-quad.h | 52 +++++++++++
goffice/math/go-regression.c | 127 ++++++++++++++++++--------
goffice/math/goffice-math.h | 1 +
7 files changed, 364 insertions(+), 48 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 8687b68..64c0dda 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2010-05-16 Morten Welinder <terra gnome org>
+
+ * goffice/math/go-quad.c: quad precision routines.
+ * goffice/math/go-regression.c (refine): Use quad precision.
+
2010-05-13 Morten Welinder <terra gnome org>
* goffice/component/go-component-factory.c (go_component_type_service_read_xml):
diff --git a/configure.in b/configure.in
index 5e68813..412e93c 100644
--- a/configure.in
+++ b/configure.in
@@ -394,6 +394,8 @@ AC_CHECK_HEADERS(ieeefp.h ieee754.h)
dnl Check for implementations of passwd/group file entry functions
AC_CHECK_HEADERS(pwd.h grp.h)
+AC_CHECK_HEADERS(fpu_control.h)
+
dnl Check for some functions
AC_CHECK_FUNCS(random drand48 finite memmove mkdtemp uname times sysconf)
diff --git a/goffice/math/Makefile.am b/goffice/math/Makefile.am
index 70a8bd1..4ca683e 100644
--- a/goffice/math/Makefile.am
+++ b/goffice/math/Makefile.am
@@ -2,12 +2,13 @@ noinst_LTLIBRARIES = libgoffice-math.la
libgoffice_math_la_SOURCES = \
go-math.c \
- go-rangefunc.c \
- go-regression.c \
- go-cspline.c \
- go-complex.c \
+ go-rangefunc.c \
+ go-regression.c \
+ go-cspline.c \
+ go-complex.c \
go-fft.c \
- go-matrix3x3.c \
+ go-matrix3x3.c \
+ go-quad.c \
go-R.c \
go-distribution.c
@@ -15,12 +16,13 @@ libgoffice_math_ladir = $(goffice_include_dir)/math
libgoffice_math_la_HEADERS = \
goffice-math.h \
go-math.h \
- go-rangefunc.h \
- go-regression.h \
- go-cspline.h \
- go-complex.h \
+ go-rangefunc.h \
+ go-regression.h \
+ go-cspline.h \
+ go-complex.h \
go-fft.h \
- go-matrix3x3.h \
+ go-matrix3x3.h \
+ go-quad.h \
go-R.h \
go-distribution.h
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
new file mode 100644
index 0000000..eeb0aef
--- /dev/null
+++ b/goffice/math/go-quad.c
@@ -0,0 +1,203 @@
+/*
+ * go-quad.c: Extended precision routines.
+ *
+ * Authors:
+ * Morten Welinder <terra gnome org>
+ *
+ * This follows "A Floating-Point Technoque for Extending the Available
+ * Precision" by T. J. Dekker in _Numerische Mathematik_ 18.
+ * Springer Verlag 1971.
+ *
+ * Note: for this to work, the processor must evaluate with the right
+ * precision. For ix86 that means trouble as the default is to evaluate
+ * with long-double precision internally. We solve this by setting the
+ * relevant x87 control flag.
+ *
+ * Note: the compiler should not rearrange expressions.
+ */
+
+#include <goffice/goffice-config.h>
+#include "go-quad.h"
+#include <math.h>
+
+/* Normalize cpu id. */
+#if !defined(i386) && (defined(__i386__) || defined(__i386))
+#define i386 1
+#endif
+
+#if defined(i386) && defined(HAVE_FPU_CONTROL_H)
+#include <fpu_control.h>
+#endif
+
+#ifndef DOUBLE
+
+#define QUAD SUFFIX(GOQuad)
+
+#define DOUBLE double
+#define SUFFIX(_n) _n
+#define DOUBLE_MANT_DIG DBL_MANT_DIG
+
+#ifdef GOFFICE_WITH_LONG_DOUBLE
+#include "go-quad.c"
+#undef DOUBLE
+#undef SUFFIX
+#undef DOUBLE_MANT_DIG
+#define DOUBLE long double
+#define SUFFIX(_n) _n ## l
+#define DOUBLE_MANT_DIG LDBL_MANT_DIG
+#endif
+
+#endif
+
+gboolean
+SUFFIX(go_quad_functional) (void)
+{
+ if (FLT_RADIX != 2)
+ return FALSE;
+
+#ifdef i386
+ if (sizeof (DOUBLE) != sizeof (double))
+ return TRUE;
+
+#ifdef HAVE_FPU_CONTROL_H
+ return TRUE;
+#else
+ return FALSE;
+#endif
+#else
+ return TRUE;
+#endif
+}
+
+static DOUBLE SUFFIX(CST);
+
+void *
+SUFFIX(go_quad_start) (void)
+{
+ void *res = NULL;
+
+ if (!go_quad_functional ())
+ g_warning ("quad precision math may not be completely accurate.");
+
+#ifdef i386
+ if (sizeof (DOUBLE) == sizeof (double)) {
+#ifdef HAVE_FPU_CONTROL_H
+ fpu_control_t state, newstate;
+ fpu_control_t mask =
+ _FPU_EXTENDED | _FPU_DOUBLE | _FPU_SINGLE;
+ _FPU_GETCW (state);
+ newstate = (state & ~mask) | _FPU_DOUBLE;
+ _FPU_SETCW (newstate);
+
+ res = g_memdup (&state, sizeof (state));
+#else
+ /* Hope for the best. */
+#endif
+ }
+#endif
+
+ SUFFIX(CST) = 1 + SUFFIX(ldexp) (1.0, (DOUBLE_MANT_DIG + 1) / 2);
+
+ return res;
+}
+
+void
+SUFFIX(go_quad_end) (void *state)
+{
+ if (!state)
+ return;
+
+#ifdef i386
+#ifdef HAVE_FPU_CONTROL_H
+ _FPU_SETCW (*(fpu_control_t*)state);
+#endif
+#endif
+
+ g_free (state);
+}
+
+void
+SUFFIX(go_quad_init) (QUAD *res, DOUBLE h, DOUBLE l)
+{
+ res->h = h;
+ res->l = l;
+}
+
+DOUBLE
+SUFFIX(go_quad_value) (const QUAD *a)
+{
+ return a->h + a->l;
+}
+
+void
+SUFFIX(go_quad_add) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+ DOUBLE r = a->h + b->h;
+ DOUBLE s = SUFFIX(fabs) (a->h) > SUFFIX(fabs) (b->h)
+ ? a->h - r + b->h + b->l + a->l
+ : b->h - r + a->h + a->l + b->l;
+ res->h = r + s;
+ res->l = r - res->h + s;
+}
+
+void
+SUFFIX(go_quad_sub) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+ DOUBLE r = a->h - b->h;
+ DOUBLE s = SUFFIX(fabs) (a->h) > SUFFIX(fabs) (b->h)
+ ? +a->h - r - b->h - b->l + a->l
+ : -b->h - r + a->h + a->l - b->l;
+ res->h = r + s;
+ res->l = r - res->h + s;
+}
+
+void
+SUFFIX(go_quad_mul12) (QUAD *res, DOUBLE x, DOUBLE y)
+{
+ DOUBLE p1 = x * SUFFIX(CST);
+ DOUBLE hx = x - p1 + p1;
+ DOUBLE tx = x - hx;
+
+ DOUBLE p2 = y * SUFFIX(CST);
+ DOUBLE hy = y - p2 + p2;
+ DOUBLE ty = y - hy;
+
+ DOUBLE p = hx * hy;
+ DOUBLE q = hx * ty + tx * hy;
+ res->h = p + q;
+ res->l = p - res->h + q + tx * ty;
+}
+
+void
+SUFFIX(go_quad_mul) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+ QUAD c;
+ SUFFIX(go_quad_mul12) (&c, a->h, b->h);
+ c.l = a->h * b->l + a->l * b->h + c.l;
+ res->h = c.h + c.l;
+ res->l = c.h - res->h + c.l;
+}
+
+void
+SUFFIX(go_quad_div) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+ QUAD c, u;
+ c.h = a->h / b->h;
+ SUFFIX(go_quad_mul12) (&u, c.h, b->h);
+ c.l = (a->h - u.h - u.l + a->l - c.h * b->l) / b->h;
+ res->h = c.h + c.l;
+ res->l = c.h - res->h + c.l;
+}
+
+int
+SUFFIX(go_quad_sgn) (const QUAD *a)
+{
+ DOUBLE d = (a->h == 0 ? a->l : a->h);
+
+ if (d == 0)
+ return 0;
+ else if (d > 0)
+ return +1;
+ else
+ return -1;
+}
diff --git a/goffice/math/go-quad.h b/goffice/math/go-quad.h
new file mode 100644
index 0000000..53acf31
--- /dev/null
+++ b/goffice/math/go-quad.h
@@ -0,0 +1,52 @@
+#ifndef __GO_QUAD_H
+#define __GO_QUAD_H
+
+#include <glib.h>
+
+G_BEGIN_DECLS
+
+typedef struct {
+ double h;
+ double l;
+} GOQuad;
+
+gboolean go_quad_functional (void);
+void *go_quad_start (void);
+void go_quad_end (void *state);
+
+void go_quad_init (GOQuad *res, double h, double l);
+
+double go_quad_value (const GOQuad *a);
+void go_quad_add (GOQuad *res, const GOQuad *a, const GOQuad *b);
+void go_quad_sub (GOQuad *res, const GOQuad *a, const GOQuad *b);
+void go_quad_mul (GOQuad *res, const GOQuad *a, const GOQuad *b);
+void go_quad_div (GOQuad *res, const GOQuad *a, const GOQuad *b);
+int go_quad_sgn (const GOQuad *a);
+
+void go_quad_mul12 (GOQuad *res, double x, double y);
+
+#ifdef GOFFICE_WITH_LONG_DOUBLE
+typedef struct {
+ long double h;
+ long double l;
+} GOQuadl;
+
+gboolean go_quad_functionall (void);
+void *go_quad_startl (void);
+void go_quad_endl (void *state);
+
+void go_quad_initl (GOQuadl *res, long double h, long double l);
+
+long double go_quad_valuel (const GOQuadl *a);
+void go_quad_addl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
+void go_quad_subl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
+void go_quad_mull (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
+void go_quad_divl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
+int go_quad_sgnl (const GOQuadl *a);
+
+void go_quad_mul12l (GOQuadl *res, long double x, long double y);
+#endif
+
+G_END_DECLS
+
+#endif
diff --git a/goffice/math/go-regression.c b/goffice/math/go-regression.c
index 006dd9f..45464f8 100644
--- a/goffice/math/go-regression.c
+++ b/goffice/math/go-regression.c
@@ -17,20 +17,21 @@
#include <stdlib.h>
#include <stdio.h>
+/* ------------------------------------------------------------------------- */
+/* Self-inclusion magic. */
+
#ifndef DOUBLE
#define DEFINE_COMMON
#define DOUBLE double
#define DOUBLE_MANT_DIG DBL_MANT_DIG
#define SUFFIX(_n) _n
+#define FORMAT_f "f"
#define FORMAT_g "g"
-#define MATRIX DOUBLE **
-#define CONSTMATRIX DOUBLE *const *const
-
-#if 0
#ifdef GOFFICE_WITH_LONG_DOUBLE
+#if 0
/* We define this to cause certain "double" functions to be implemented in
terms of their "long double" counterparts. */
#define BOUNCE
@@ -49,6 +50,7 @@ general_linear_regressionl (long double *const *const xss, int xdim,
#undef DOUBLE
#undef DOUBLE_MANT_DIG
#undef SUFFIX
+#undef FORMAT_f
#undef FORMAT_g
#ifdef HAVE_SUNMATH_H
#include <sunmath.h>
@@ -56,13 +58,20 @@ general_linear_regressionl (long double *const *const xss, int xdim,
#define DOUBLE long double
#define DOUBLE_MANT_DIG LDBL_MANT_DIG
#define SUFFIX(_n) _n ## l
+#define FORMAT_f "Lf"
#define FORMAT_g "Lg"
#endif
#endif
+/* ------------------------------------------------------------------------- */
+
#ifdef DEFINE_COMMON
+#define MATRIX DOUBLE **
+#define CONSTMATRIX DOUBLE *const *const
+#define QUAD SUFFIX(GOQuad)
+
#undef DEBUG_NEAR_SINGULAR
#undef DEBUG_REFINEMENT
@@ -112,25 +121,25 @@ general_linear_regressionl (long double *const *const xss, int xdim,
int _i, _j, _d1, _d2; \
_d1 = (dim1); \
_d2 = (dim2); \
- for (_i = 0; _i < _d1; _i++) \
- { \
- for (_j = 0; _j < _d2; _j++) \
- g_printerr (" %19.10" FORMAT_g, (var)[_i][_j]); \
- g_printerr ("\n"); \
+ for (_j = 0; _j < _d2; _j++) { \
+ for (_i = 0; _i < _d1; _i++) { \
+ g_printerr (" %12.8" FORMAT_g, (var)[_i][_j]); \
} \
+ g_printerr ("\n"); \
+ } \
} while (0)
#endif
/*
- * ---> j
+ * ---> i
*
* | ********
* | ********
* | ******** A[i][j]
* v ********
* ********
- * i ********
+ * j ********
* ********
* ********
*
@@ -207,6 +216,10 @@ SUFFIX(QR) (CONSTMATRIX A, MATRIX Q, MATRIX R, int m, int n)
(void)SUFFIX(go_range_sumsq) (Q[k], n, &L);
L = SUFFIX(sqrt) (L);
+#if 0
+ PRINT_MATRIX (Q, m, n);
+ g_printerr ("L[%d] = %20.15" FORMAT_g "\n", k, L);
+#endif
if (L == 0)
return GO_REG_singular;
@@ -233,20 +246,28 @@ SUFFIX(QR) (CONSTMATRIX A, MATRIX Q, MATRIX R, int m, int n)
static void
SUFFIX(calc_residual) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
- const DOUBLE *res, DOUBLE *residual, DOUBLE *N2)
+ const DOUBLE *res, QUAD *residual, QUAD *N2)
{
int i, j;
- *N2 = 0;
+ SUFFIX(go_quad_init) (N2, 0, 0);
for (i = 0; i < n; i++) {
- DOUBLE d = 0;
- for (j = 0; j < dim; j++)
- d += A[j][i] * res[j];
- d = b[i] - d;
- (*N2) += d * d;
+ QUAD d;
+
+ SUFFIX(go_quad_init) (&d, b[i], 0);
+
+ for (j = 0; j < dim; j++) {
+ QUAD e;
+ SUFFIX(go_quad_mul12) (&e, A[j][i], res[j]);
+ SUFFIX(go_quad_sub) (&d, &d, &e);
+ }
+
if (residual)
residual[i] = d;
+
+ SUFFIX(go_quad_mul) (&d, &d, &d);
+ SUFFIX(go_quad_add) (N2, N2, &d);
}
}
@@ -254,50 +275,79 @@ static void
SUFFIX(refine) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
CONSTMATRIX Q, CONSTMATRIX R, DOUBLE *result)
{
- DOUBLE *residual = g_new (DOUBLE, n);
+ QUAD *residual = g_new (QUAD, n);
+ QUAD *newresidual = g_new (QUAD, n);
DOUBLE *delta = g_new (DOUBLE, dim);
DOUBLE *newresult = g_new (DOUBLE, dim);
int pass;
- DOUBLE best_N2;
+ QUAD best_N2;
+ void *state;
+
+ state = SUFFIX(go_quad_start) ();
SUFFIX(calc_residual) (A, b, dim, n, result, residual, &best_N2);
#ifdef DEBUG_REFINEMENT
- g_printerr ("-: Residual norm = %20.15" FORMAT_g "\n", best_N2);
+ g_printerr ("-: Residual norm = %20.15" FORMAT_g "\n",
+ SUFFIX(go_quad_value) (&best_N2));
#endif
- for (pass = 1; pass < 5; pass++) {
+ for (pass = 1; pass < 10; pass++) {
int i, j;
- DOUBLE N2;
+ QUAD dres, N2;
/* newresult = R^-1 Q^T residual */
for (i = dim - 1; i >= 0; i--) {
- DOUBLE acc = SUFFIX(dot_product) (Q[i], residual, n);
- for (j = i + 1; j < dim; j++)
- acc -= R[j][i] * delta[j];
- delta[i] = acc / R[i][i];
+ QUAD acc, d, Rii;
+
+ SUFFIX(go_quad_init) (&acc, 0, 0);
+
+ for (j = 0; j < n; j++) {
+ QUAD q, qr;
+ SUFFIX(go_quad_init) (&q, Q[i][j], 0);
+ SUFFIX(go_quad_mul) (&qr, &q, residual + j);
+ SUFFIX(go_quad_add) (&acc, &acc, &qr);
+ }
+ for (j = i + 1; j < dim; j++) {
+ QUAD Rd;
+ SUFFIX(go_quad_mul12) (&Rd, R[j][i], delta[j]);
+ SUFFIX(go_quad_sub) (&acc, &acc, &Rd);
+ }
+
+ SUFFIX(go_quad_init) (&Rii, R[i][i], 0);
+ SUFFIX(go_quad_div) (&d, &acc, &Rii);
+
+ delta[i] = SUFFIX(go_quad_value) (&d);
+ newresult[i] = delta[i] + result[i];
#ifdef DEBUG_REFINEMENT
- g_printerr ("d[%2d] = %20" FORMAT_g
- " %20" FORMAT_g "\n",
- i, delta[i], result[i]);
+ g_printerr ("d[%2d] = %20.15" FORMAT_f
+ " %20.15" FORMAT_f
+ " %20.15" FORMAT_f "\n",
+ i, delta[i], result[i], newresult[i]);
#endif
}
- for (i = 0; i < dim; i++)
- newresult[i] = delta[i] + result[i];
-
- SUFFIX(calc_residual) (A, b, dim, n, newresult, residual, &N2);
+ SUFFIX(calc_residual) (A, b, dim, n, newresult, newresidual, &N2);
+ SUFFIX (go_quad_sub) (&dres, &N2, &best_N2);
#ifdef DEBUG_REFINEMENT
- g_printerr ("%d: Residual norm = %20.15" FORMAT_g "\n", pass, N2);
+ g_printerr ("%d: Residual norm = %20.15" FORMAT_g
+ " (%.5" FORMAT_g ")\n",
+ pass,
+ SUFFIX(go_quad_value) (&N2),
+ SUFFIX(go_quad_value) (&dres));
#endif
- if (N2 >= best_N2)
+ if (SUFFIX(go_quad_sgn) (&dres) >= 0)
break;
best_N2 = N2;
COPY_VECTOR (result, newresult, dim);
+ COPY_VECTOR (residual, newresidual, n);
}
+ SUFFIX(go_quad_end) (state);
+
g_free (residual);
+ g_free (newresidual);
g_free (newresult);
g_free (delta);
}
@@ -661,6 +711,7 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
DOUBLE *inv = g_new (DOUBLE, xdim);
int err;
int k;
+ QUAD N2;
/* This should not fail since n >= 1. */
err = SUFFIX(go_range_average) (ys, n, &stat_->ybar);
@@ -680,8 +731,8 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
g_assert (err == 0);
}
- SUFFIX(calc_residual) (xss, ys, xdim, n, result, NULL,
- &stat_->ss_resid);
+ SUFFIX(calc_residual) (xss, ys, xdim, n, result, NULL, &N2);
+ stat_->ss_resid = SUFFIX(go_quad_value) (&N2);
stat_->sqr_r = (stat_->ss_total == 0)
? 1
diff --git a/goffice/math/goffice-math.h b/goffice/math/goffice-math.h
index 7762d3c..c307224 100644
--- a/goffice/math/goffice-math.h
+++ b/goffice/math/goffice-math.h
@@ -9,6 +9,7 @@
#include <goffice/math/go-matrix3x3.h>
#include <goffice/math/go-rangefunc.h>
#include <goffice/math/go-regression.h>
+#include <goffice/math/go-quad.h>
#include <goffice/math/go-R.h>
#endif
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]