[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [PATCH] New module to perform decimal floating point arithmetic for
|
From: |
John Darrington |
|
Subject: |
Re: [PATCH] New module to perform decimal floating point arithmetic for charts. |
|
Date: |
Tue, 13 Jan 2015 07:04:25 +0100 |
|
User-agent: |
Mutt/1.5.21 (2010-09-15) |
I'd appreciate a review of this patch before I push it.
On Tue, Jan 13, 2015 at 07:03:17AM +0100, John Darrington wrote:
This change adds a small module to perform floating point arithmetic
with a decimal base, and uses it to calculate the graticule marks for
charts.
Rationale: Graticule marks want to be displayed in decimal. However
not all decimal values can be represented in a double precision
floating point binary. This approach pushes the precision loss from
the value of the mark, to the position on the chart - we don't
particularly care if a tick mark is a fraction of a pixel out of place,
but we do care if 4.0 is displayed as 3.9999999999
---
src/math/automake.mk | 1 +
src/math/chart-geometry.c | 130 ++++-
src/math/chart-geometry.h | 7 +-
src/math/decimal.c | 572
++++++++++++++++++++
src/math/decimal.h | 115 ++++
src/math/histogram.c | 15 +-
src/output/cairo-chart.c | 48 +-
src/output/cairo-chart.h | 4 +-
src/output/charts/boxplot-cairo.c | 2 +-
src/output/charts/np-plot-cairo.c | 8 +-
src/output/charts/plot-hist-cairo.c | 21 +-
src/output/charts/roc-chart-cairo.c | 4 +-
src/output/charts/scatterplot-cairo.c | 4 +-
src/output/charts/scree-cairo.c | 4 +-
src/output/charts/spreadlevel-cairo.c | 4 +-
tests/automake.mk | 31 ++
.../math/chart-geometry-test.c | 72 +--
tests/math/chart-geometry.at | 35 ++
tests/math/chart-get-scale-test.c | 96 ++++
tests/math/decimal-test.c | 310 +++++++++++
tests/math/decimal.at | 7 +
21 files changed, 1394 insertions(+), 96 deletions(-)
create mode 100644 src/math/decimal.c
create mode 100644 src/math/decimal.h
copy src/math/chart-geometry.c => tests/math/chart-geometry-test.c (50%)
create mode 100644 tests/math/chart-geometry.at
create mode 100644 tests/math/chart-get-scale-test.c
create mode 100644 tests/math/decimal-test.c
create mode 100644 tests/math/decimal.at
diff --git a/src/math/automake.mk b/src/math/automake.mk
index cdcc2c8..f40aac5 100644
--- a/src/math/automake.mk
+++ b/src/math/automake.mk
@@ -17,6 +17,7 @@ src_math_libpspp_math_la_SOURCES = \
src/math/covariance.h \
src/math/correlation.c \
src/math/correlation.h \
+ src/math/decimal.c src/math/decimal.h \
src/math/extrema.c src/math/extrema.h \
src/math/histogram.c src/math/histogram.h \
src/math/interaction.c src/math/interaction.h \
diff --git a/src/math/chart-geometry.c b/src/math/chart-geometry.c
index eb0ad92..35380a6 100644
--- a/src/math/chart-geometry.c
+++ b/src/math/chart-geometry.c
@@ -1,5 +1,5 @@
/* PSPP - a program for statistical analysis.
- Copyright (C) 2004 Free Software Foundation, Inc.
+ Copyright (C) 2004, 2015 Free Software Foundation, Inc.
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
@@ -17,41 +17,135 @@
#include <config.h>
#include <math.h>
#include <float.h>
+#include <assert.h>
#include "chart-geometry.h"
+#include "decimal.h"
+#include <stdlib.h>
-static const double standard_ticks[] = {1, 2, 5, 10};
-
+static const double standard_tick[] = {1, 2, 5, 10};
/* Adjust tick to be a sensible value
ie: ... 0.1,0.2,0.5, 1,2,5, 10,20,50 ... */
-double
-chart_rounded_tick (double tick)
+void
+chart_rounded_tick (double tick, struct decimal *result)
{
int i;
- double diff = DBL_MAX;
- double t = tick;
-
- double factor;
+ struct decimal ddif = {1, 1000};
/* Avoid arithmetic problems with very small values */
if (fabs (tick) < DBL_EPSILON)
- return 0;
+ {
+ result->ordinate = 0;
+ result->mantissa = 0;
+ return;
+ }
- factor = pow (10,ceil (log10 (standard_ticks[0] / tick)));
+ struct decimal dt;
+ decimal_from_double (&dt, tick);
+
+ double expd = dec_log10 (&dt) - 1;
- for (i = 3 ; i >= 0 ; --i)
+ for (i = 0 ; i < 4 ; ++i)
{
- const double d = fabs (tick - standard_ticks[i] / factor);
+ struct decimal candidate;
+ struct decimal delta;
- if ( d < diff )
+ decimal_init (&candidate, standard_tick[i], expd);
+
+ delta = dt;
+ decimal_subtract (&delta, &candidate);
+ delta.ordinate = llabs (delta.ordinate);
+
+ if (decimal_cmp (&delta, &ddif) < 0)
{
- diff = d;
- t = standard_ticks[i] / factor ;
+ ddif = delta;
+ *result = candidate;
}
}
-
- return t;
}
+/*
+ Find a set {LOWER, INTERVAL, N_TICKS} such that:
+
+ LOWER <= LOWDBL,
+ LOWER + INTERVAL > LOWDBL,
+
+ LOWER + N_TICKS * INTERVAL < HIGHDBL
+ LOWER + (N_TICKS + 1) * INTERVAL >= HIGHDBL
+
+ INTERVAL = X * 10^N
+ where: N is integer
+ and X is an element of {1, 2, 5}
+
+ In other words:
+
+ INTERVAL
+ > <
+ |....+....+....+. .+....|
+ LOWER 1 2 3 N_TICKS
+ ^LOWDBL ^HIGHDBL
+*/
+void
+chart_get_scale (double highdbl, double lowdbl,
+ struct decimal *lower,
+ struct decimal *interval,
+ int *n_ticks)
+{
+ int i;
+ double fitness = DBL_MAX;
+ *n_ticks = 0;
+ struct decimal high;
+ struct decimal low;
+
+ assert (highdbl >= lowdbl);
+
+ decimal_from_double (&high, highdbl);
+ decimal_from_double (&low, lowdbl);
+
+ struct decimal diff = high;
+ decimal_subtract (&diff, &low);
+
+ double expd = dec_log10 (&diff) - 2;
+
+ /* Find the most pleasing interval */
+ for (i = 1; i < 4; ++i)
+ {
+ struct decimal clbound = low;
+ struct decimal cubound = high;
+ struct decimal candidate;
+ decimal_init (&candidate, standard_tick[i], expd);
+
+ decimal_divide (&clbound, &candidate);
+ int fl = decimal_floor (&clbound);
+ decimal_int_multiply (&candidate, fl);
+ clbound = candidate;
+
+
+ decimal_init (&candidate, standard_tick[i], expd);
+ decimal_divide (&cubound, &candidate);
+ int fu = decimal_ceil (&cubound);
+ decimal_int_multiply (&candidate, fu);
+
+ cubound = candidate;
+
+ decimal_init (&candidate, standard_tick[i], expd);
+ decimal_subtract (&cubound, &clbound);
+ decimal_divide (&cubound, &candidate);
+
+
+ ord_t n_int = decimal_floor (&cubound);
+
+ /* We prefer to have between 5 and 10 tick marks on a scale */
+ double f = 1 - exp (-0.2 * fabs (n_int - 7.5) / 7.5);
+
+ if (f < fitness)
+ {
+ fitness = f;
+ *lower = clbound;
+ *interval = candidate;
+ *n_ticks = n_int;
+ }
+ }
+}
diff --git a/src/math/chart-geometry.h b/src/math/chart-geometry.h
index ef481fe..c543086 100644
--- a/src/math/chart-geometry.h
+++ b/src/math/chart-geometry.h
@@ -18,6 +18,11 @@
#ifndef CHART_GEOMETRY_H
#define CHART_GEOMETRY_H
-double chart_rounded_tick(double tick);
+struct decimal;
+void chart_rounded_tick(double tick, struct decimal *);
+
+void chart_get_scale (double high, double low,
+ struct decimal *lower, struct decimal *interval, int
*n_ticks);
+
#endif
diff --git a/src/math/decimal.c b/src/math/decimal.c
new file mode 100644
index 0000000..ba33e76
--- /dev/null
+++ b/src/math/decimal.c
@@ -0,0 +1,572 @@
+/* PSPP - a program for statistical analysis.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+
+ 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 3 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, see <http://www.gnu.org/licenses/>.
*/
+
+
+#include <config.h>
+#include <stdio.h>
+#include <stdint.h>
+#include <string.h>
+#include <limits.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <stdbool.h>
+#include <math.h>
+#include "libpspp/i18n.h"
+
+#include "decimal.h"
+
+int dec_warning;
+
+static bool
+down (struct decimal *dec)
+{
+ if (dec->ordinate % 10 == 0 && dec->mantissa < MANT_MAX - 1)
+ {
+ dec->ordinate /= 10;
+ dec->mantissa++;
+ return true;
+ }
+
+ return false;
+}
+
+static bool
+up (struct decimal *dec)
+{
+ if (llabs (dec->ordinate) < ORD_MAX / 10 && dec->mantissa >
MANT_MIN)
+ {
+ dec->ordinate *= 10;
+ dec->mantissa--;
+ return true;
+ }
+ return false;
+}
+
+
+/* Reduce the absolute value of the ordinate to the smallest possible,
+ without loosing precision */
+static void
+reduce (struct decimal *dec)
+{
+ if (dec->ordinate == 0)
+ {
+ dec->mantissa = 0;
+ return;
+ }
+
+ while (dec->ordinate % 10 == 0)
+ {
+ if (! down (dec))
+ break;
+ }
+}
+
+/* Attempt to make the mantissas of BM and SM equal.
+ Prerequisite: the mantissa SM must be no greater than that of BM.
+ */
+static void
+normalisebs (struct decimal *sm, struct decimal *bm)
+{
+ while (sm->mantissa < bm->mantissa)
+ {
+ if (down (sm))
+ ;
+ else if (up (bm))
+ ;
+ else
+ {
+ dec_warning = DEC_PREC;
+ break;
+ }
+ }
+
+ while (sm->mantissa < bm->mantissa)
+ {
+ sm->ordinate /= 10;
+ sm->mantissa++;
+ }
+}
+
+
+/* arrange d1 and d2 such that thier mantissas are equal */
+void
+normalise (struct decimal *d1, struct decimal *d2)
+{
+ normalisebs (d1, d2);
+ normalisebs (d2, d1);
+}
+
+
+
+/* Return log base 10 of D */
+mant_t
+dec_log10 (const struct decimal *d_)
+{
+ struct decimal d = *d_;
+
+ while (llabs (d.ordinate) > 0)
+ {
+ d.ordinate /= 10;
+ d.mantissa++;
+ }
+
+ return d.mantissa;
+}
+
+
+
+/* Return the smallest integer >= d */
+static ord_t
+decimal_ceil_pos (const struct decimal *d)
+{
+ mant_t m = d->mantissa;
+ ord_t o = d->ordinate;
+
+ assert (d->ordinate >= 0);
+
+ while (m > 0)
+ {
+ o *= 10;
+ m--;
+ }
+
+ while (m < 0)
+ {
+ bool flag = o % 10;
+ o /= 10;
+ if (flag)
+ o++;
+ m++;
+ }
+
+ return o;
+}
+
+/* Return the largest integer <= d */
+static ord_t
+decimal_floor_pos (const struct decimal *d)
+{
+ mant_t m = d->mantissa;
+ ord_t o = d->ordinate;
+
+ assert (d->ordinate >= 0);
+
+ while (m > 0)
+ {
+ m--;
+ o *= 10;
+ }
+
+ while (m < 0)
+ {
+ m++;
+ o /= 10;
+ }
+
+
+ return o;
+}
+
+/* Return the smallest integer which is no less than D.
+ (round towards minus infinity) */
+ord_t
+decimal_floor (const struct decimal *d)
+{
+ if (d->ordinate >= 0)
+ return decimal_floor_pos (d);
+ else
+ {
+ struct decimal dd = *d;
+ dd.ordinate = llabs (dd.ordinate);
+ return -decimal_ceil_pos (&dd);
+ }
+}
+
+/* Return the largest integer which is no greater than D.
+ (round towards plus infinity) */
+ord_t
+decimal_ceil (const struct decimal *d)
+{
+ if (d->ordinate >= 0)
+ return decimal_ceil_pos (d);
+ else
+ {
+ struct decimal dd = *d;
+ dd.ordinate = llabs (dd.ordinate);
+ return -decimal_floor_pos (&dd);
+ }
+}
+
+/* Add SRC onto DEST */
+void
+decimal_add (struct decimal *dest, const struct decimal *src_)
+{
+ struct decimal src = *src_;
+
+ src.ordinate = -src.ordinate;
+
+ decimal_subtract (dest, &src);
+}
+
+/* Subtract SRC from DEST */
+void
+decimal_subtract (struct decimal *dest, const struct decimal *src_)
+{
+ struct decimal src = *src_;
+
+ normalise (dest, &src);
+
+ bool dest_neg = dest->ordinate < 0;
+ bool src_neg = src.ordinate < 0;
+
+ bool expected_neg = dest_neg * src_neg;
+
+ if (dest->ordinate == src.ordinate)
+ {
+ expected_neg = 0;
+ }
+ else if (llabs (src.ordinate) > llabs (dest->ordinate))
+ {
+ if (dest_neg == src_neg)
+ expected_neg = !expected_neg;
+ }
+
+ dest->ordinate -= src.ordinate;
+
+ bool result_neg = dest->ordinate < 0;
+
+ if (expected_neg != result_neg)
+ {
+ /* The operation has resulted in an overflow.
+ To resolve this, undo the operation,
+ reduce the precision and try again */
+
+ dest->ordinate += src.ordinate;
+
+ dest->ordinate /= 10;
+ src.ordinate /= 10;
+
+ dest->mantissa ++;
+ src.mantissa ++;
+
+ dest->ordinate -= src.ordinate;
+ }
+
+ reduce (dest);
+
+}
+
+/* Initialise DEC with ordinate ORD and mantissa MANT */
+void
+decimal_init (struct decimal *dec, ord_t ord, mant_t mant)
+{
+ dec->ordinate = ord;
+ dec->mantissa = mant;
+ reduce (dec);
+}
+
+
+/*
+ Compare D1 and D2.
+
+ Returns zero if equal, +1 if D1 > D2 and -1 if D1 < D2
+*/
+int
+decimal_cmp (const struct decimal *d1, const struct decimal *d2)
+{
+ struct decimal td1 = *d1;
+ struct decimal td2 = *d2;
+
+ normalise (&td1, &td2);
+
+ if (td1.ordinate < td2.ordinate)
+ return -1;
+
+ return (td1.ordinate > td2.ordinate);
+}
+
+
+/* Multiply DEST by SRC */
+void
+decimal_multiply (struct decimal *dest, const struct decimal *src)
+{
+ dest->ordinate *= src->ordinate;
+ dest->mantissa += src->mantissa;
+}
+
+
+/* Multiply DEST by M */
+void
+decimal_int_multiply (struct decimal *dest, ord_t m)
+{
+ dest->ordinate *= m;
+ reduce (dest);
+}
+
+/* Divide DEST by M */
+void
+decimal_int_divide (struct decimal *dest, ord_t m)
+{
+ while (dest->ordinate % m)
+ {
+ if (labs (dest->ordinate) > ORD_MAX / 10)
+ {
+ dec_warning = DEC_PREC;
+ break;
+ }
+ up (dest);
+ }
+ dest->ordinate /= m;
+}
+
+/* Divide DEST by SRC */
+void
+decimal_divide (struct decimal *dest, const struct decimal *src)
+{
+ while (dest->ordinate % src->ordinate)
+ {
+ if (labs (dest->ordinate) > ORD_MAX / 10)
+ {
+ dec_warning = DEC_PREC;
+ break;
+ }
+ up (dest);
+ }
+
+ dest->ordinate /= src->ordinate;
+ dest->mantissa -= src->mantissa;
+}
+
+/* Print the value of DEC to F. Probably useful only for debugging */
+void
+decimal_show (const struct decimal *dec, FILE *f)
+{
+ fprintf (f, PR_ORD " x 10^" PR_MANT "\n", dec->ordinate, dec->mantissa);
+}
+
+
+/* Reverse the characters in string S which has length LEN */
+static void
+reverse (char *s, int len)
+{
+ int i;
+ for (i = 0; i < len / 2; ++i)
+ {
+ char temp = s[len - i - 1];
+ s[len - i - 1] = s[i];
+ s[i] = temp;
+ }
+}
+
+/* Return a string representation of DEC on the heap.
+ The caller is responsible for freeing the string */
+char *
+decimal_to_string (const struct decimal *dec)
+{
+ int cap = 16;
+ int len = 0;
+ char *s = calloc (cap, 1);
+ ord_t ordinate = dec->ordinate;
+
+ while (len < dec->mantissa)
+ {
+ s[len++] = '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+
+ while (ordinate)
+ {
+ s[len++] = labs (ordinate % 10) + '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ ordinate /= 10;
+ }
+
+ if (ordinate < 0)
+ ordinate = -ordinate;
+
+ while (len < -dec->mantissa)
+ {
+ s[len++] = '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+
+ if (dec->mantissa < 0 )
+ {
+ if (len <= -dec->mantissa)
+ {
+ s[len++] = get_system_decimal ();
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ s[len++] = '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+ else
+ {
+ int i;
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ for (i = len - 1 ; i >= -dec->mantissa ; --i)
+ s[i + 1] = s[i];
+ s[i + 1] = get_system_decimal ();
+ len++;
+ }
+ }
+
+ if (dec->ordinate < 0)
+ {
+ s[len++] = '-';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+
+
+ reverse (s, len);
+
+ {
+ int abs_len = len;
+ if (dec->ordinate < 0)
+ abs_len--;
+
+ while (abs_len++ <= dec->mantissa)
+ {
+ s[len++] = '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+ }
+
+ return s;
+}
+
+
+/* Initialise DECIMAL from INPUT.
+ INPUT should be a convential decimal representation.
+ */
+void
+decimal_init_from_string (struct decimal *decimal, const char *input)
+{
+ ord_t ordinate = 0;
+
+ int point = -1;
+ int lsd = -1;
+ int fsd = -1;
+ int i = 0;
+ int len = 0;
+ int sign = 1;
+
+ const char *p;
+
+ for (p = input; *p ; p++)
+ {
+ if (*p == '-')
+ {
+ sign = -1;
+ }
+ else if (*p == get_system_decimal ())
+ {
+ assert (point == -1);
+ point = i;
+ }
+ else if (*p > '0' && *p <= '9')
+ {
+ lsd = i;
+ if (fsd == -1)
+ fsd = i;
+ }
+ else if (*p == '0')
+ /* ignore */
+ ;
+ else
+ {
+ fprintf (stderr, "Error: invalid character %c\n", *p);
+ return;
+ }
+
+ i++;
+ }
+ len = i;
+
+ if (point == -1)
+ point = len;
+
+ mant_t mantissa = 0;
+ if (fsd >= 0)
+ {
+ mant_t m = 1;
+ for (i = lsd ; i >= fsd ; --i)
+ {
+ if (input[i] != get_system_decimal ())
+ {
+ if (ordinate > ORD_MAX - m * (input[i] - '0'))
+ {
+ fprintf (stderr, "Overflow reading value %s\n", input);
+ break;
+ }
+ ordinate += m * (input[i] - '0');
+ m *= 10;
+ }
+ }
+
+ if (lsd > point)
+ mantissa = point - lsd;
+ else
+ mantissa = point - lsd - 1;
+ }
+
+ decimal->ordinate = ordinate * sign;
+ decimal->mantissa = mantissa;
+}
+
+
+
+/* Initialise DEC from the binary fp value X */
+void
+decimal_from_double (struct decimal *dec, double x)
+{
+ dec->mantissa = 0;
+
+ while (trunc (x) != x)
+ {
+ if (fabs (x) > ORD_MAX / 10.0)
+ {
+ dec_warning = DEC_PREC;
+ break;
+ }
+ x *= 10.0;
+ dec->mantissa--;
+ }
+
+ dec->ordinate = x;
+}
+
+/* Return a binary floating point value
+ approximating DEC */
+double
+decimal_to_double (const struct decimal *dec)
+{
+ double x = dec->ordinate;
+ int mult = dec->mantissa;
+
+ while (mult < 0)
+ {
+ x /= 10.0;
+ mult++;
+ }
+
+ while (mult > 0)
+ {
+ x *= 10.0;
+ mult--;
+ }
+
+ return x;
+}
diff --git a/src/math/decimal.h b/src/math/decimal.h
new file mode 100644
index 0000000..2a3fd60
--- /dev/null
+++ b/src/math/decimal.h
@@ -0,0 +1,115 @@
+/* PSPP - a program for statistical analysis.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+
+ 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 3 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, see <http://www.gnu.org/licenses/>.
*/
+
+
+#ifndef DECIMAL_H
+#define DECIMAL_H
+
+/* This module provides a rudimentary floating point implementation using
a decimal
+ base. It can be used for floating point calculations where it is
desirable that
+ the result is representable in decimal base.
+
+ Any of the functions may set the static variable dec_warning to
non-zero if a
+ loss of precision or other issue occurs.
+
+ This does not purport to be efficient, either in time or space.
+ */
+
+#include <stdio.h>
+#include <stdbool.h>
+
+#include <limits.h>
+
+#define DEC_PREC 1 /* operation resulted in a loss of precision */
+extern int dec_warning;
+
+
+#define ORDINATE_LONG
+#define MANTISSA_LONG
+
+#ifdef ORDINATE_SHORT
+typedef short ord_t;
+static const short ORD_MAX = SHRT_MAX;
+#define PR_ORD "%d"
+#endif
+
+#ifdef ORDINATE_INT
+typedef int ord_t;
+static const int ORD_MAX = INT_MAX;
+#define PR_ORD "%d"
+#endif
+
+#ifdef ORDINATE_LONG
+typedef long ord_t;
+static const long ORD_MAX = LONG_MAX;
+#define PR_ORD "%ld"
+#endif
+
+
+
+#ifdef MANTISSA_SHORT
+typedef short mant_t;
+static const short MANT_MAX = SHRT_MAX;
+#define PR_MANT "%d"
+#endif
+
+#ifdef MANTISSA_INT
+typedef int mant_t;
+static const int MANT_MAX = INT_MAX;
+#define PR_MANT "%d"
+#endif
+
+#ifdef MANTISSA_LONG
+typedef long mant_t;
+static const long MANT_MAX = LONG_MAX;
+#define PR_MANT "%ld"
+#endif
+
+
+
+#define MANT_MIN (-MANT_MAX - 1)
+#define ORD_MIN (-ORD_MAX - 1)
+
+struct decimal
+{
+ ord_t ordinate;
+ mant_t mantissa;
+};
+
+void normalise (struct decimal *d1, struct decimal *d2);
+void decimal_init (struct decimal *dec, ord_t ord, mant_t mant);
+void decimal_init_from_string (struct decimal *dec, const char *s);
+int decimal_cmp (const struct decimal *d1, const struct decimal *d2);
+void decimal_multiply (struct decimal *dest, const struct decimal *src);
+void decimal_int_multiply (struct decimal *dest, ord_t m);
+void decimal_int_divide (struct decimal *dest, ord_t m);
+void decimal_divide (struct decimal *dest, const struct decimal *src);
+void decimal_show (const struct decimal *dec, FILE *f);
+char *decimal_to_string (const struct decimal *dec);
+
+void decimal_add (struct decimal *dest, const struct decimal *);
+void decimal_subtract (struct decimal *dest, const struct decimal *);
+ord_t decimal_ceil (const struct decimal *d);
+ord_t decimal_floor (const struct decimal *d);
+mant_t dec_log10 (const struct decimal *d);
+
+
+void decimal_from_double (struct decimal *dec, double x);
+double decimal_to_double (const struct decimal *dec);
+
+
+
+#endif
diff --git a/src/math/histogram.c b/src/math/histogram.c
index b000672..3c324ef 100644
--- a/src/math/histogram.c
+++ b/src/math/histogram.c
@@ -17,6 +17,7 @@
#include <config.h>
#include "math/histogram.h"
+#include "math/decimal.h"
#include <gsl/gsl_histogram.h>
#include <math.h>
@@ -211,8 +212,9 @@ adjust_bin_ranges (double bin_width, double min,
double max, double *adj_min, do
struct histogram *
-histogram_create (double bin_width, double min, double max)
+histogram_create (double bin_width_in, double min, double max)
{
+ struct decimal bin_width;
const int MAX_BINS = 25;
struct histogram *h;
struct statistic *stat;
@@ -225,16 +227,17 @@ histogram_create (double bin_width, double min,
double max)
return NULL;
}
- assert (bin_width > 0);
+ assert (bin_width_in > 0);
- bin_width = chart_rounded_tick (bin_width);
- bins = adjust_bin_ranges (bin_width, min, max, &adjusted_min,
&adjusted_max);
+ chart_rounded_tick (bin_width_in, &bin_width);
+ bins = adjust_bin_ranges (decimal_to_double (&bin_width),
+ min, max, &adjusted_min, &adjusted_max);
/* Force the number of bins to lie in a sensible range. */
if (bins > MAX_BINS)
{
- bin_width = chart_rounded_tick ((max - min) / (double) (MAX_BINS -
1));
- bins = adjust_bin_ranges (bin_width,
+ chart_rounded_tick ((max - min) / (double) (MAX_BINS - 1),
&bin_width);
+ bins = adjust_bin_ranges (decimal_to_double (&bin_width),
min, max, &adjusted_min, &adjusted_max);
}
diff --git a/src/output/cairo-chart.c b/src/output/cairo-chart.c
index 8886b26..01c5d5c 100644
--- a/src/output/cairo-chart.c
+++ b/src/output/cairo-chart.c
@@ -1,5 +1,5 @@
/* PSPP - a program for statistical analysis.
- Copyright (C) 2004, 2009, 2010, 2011, 2014 Free Software Foundation,
Inc.
+ Copyright (C) 2004, 2009, 2010, 2011, 2014, 2015 Free Software
Foundation, Inc.
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
@@ -17,6 +17,8 @@
#include <config.h>
#include "output/cairo-chart.h"
+#include "math/decimal.h"
+#include "math/chart-geometry.h"
#include <assert.h>
#include <cairo/cairo.h>
@@ -347,45 +349,57 @@ xrchart_write_title (cairo_t *cr, const struct
xrchart_geometry *geom,
static void
xrchart_write_scale (cairo_t *cr, struct xrchart_geometry *geom,
- double smin, double smax, int ticks, enum tick_orientation
orient)
+ double smin, double smax, enum tick_orientation orient)
{
int s;
+ int ticks;
- const double tick_interval =
- chart_rounded_tick ((smax - smin) / (double) ticks);
+ struct decimal dinterval;
+ struct decimal dlower;
+ struct decimal dupper;
- int upper = ceil (smax / tick_interval);
- int lower = floor (smin / tick_interval);
+ chart_get_scale (smax, smin, &dlower, &dinterval, &ticks);
- geom->axis[orient].max = tick_interval * upper;
- geom->axis[orient].min = tick_interval * lower;
+ dupper = dinterval;
+ decimal_int_multiply (&dupper, ticks);
+ decimal_add (&dupper, &dlower);
+ double tick_interval = decimal_to_double (&dinterval);
+
+ geom->axis[orient].max = decimal_to_double (&dupper);
+ geom->axis[orient].min = decimal_to_double (&dlower);
+
geom->axis[orient].scale = (fabs (geom->axis[orient].data_max -
geom->axis[orient].data_min)
- / fabs (geom->axis[orient].max - geom->axis[orient].min));
+ / fabs (geom->axis[orient].max -
geom->axis[orient].min));
+
+ struct decimal pos = dlower;
- for (s = 0 ; s < upper - lower; ++s)
+ for (s = 0 ; s < ticks; ++s)
{
- double pos = (s + lower) * tick_interval;
+ char *str = decimal_to_string (&pos);
draw_tick (cr, geom, orient, false,
- s * tick_interval * geom->axis[orient].scale, "%.*g",
- DBL_DIG + 1, pos);
+ s * tick_interval * geom->axis[orient].scale,
+ "%s", str);
+ free (str);
+
+ decimal_add (&pos, &dinterval);
}
}
/* Set the scale for the ordinate */
void
xrchart_write_yscale (cairo_t *cr, struct xrchart_geometry *geom,
- double smin, double smax, int ticks)
+ double smin, double smax)
{
- xrchart_write_scale (cr, geom, smin, smax, ticks, SCALE_ORDINATE);
+ xrchart_write_scale (cr, geom, smin, smax, SCALE_ORDINATE);
}
/* Set the scale for the abscissa */
void
xrchart_write_xscale (cairo_t *cr, struct xrchart_geometry *geom,
- double smin, double smax, int ticks)
+ double smin, double smax)
{
- xrchart_write_scale (cr, geom, smin, smax, ticks, SCALE_ABSCISSA);
+ xrchart_write_scale (cr, geom, smin, smax, SCALE_ABSCISSA);
}
diff --git a/src/output/cairo-chart.h b/src/output/cairo-chart.h
index edf4ed1..27000cb 100644
--- a/src/output/cairo-chart.h
+++ b/src/output/cairo-chart.h
@@ -118,12 +118,12 @@ void xrchart_write_title (cairo_t *, const struct
xrchart_geometry *,
/* Set the scale for the abscissa */
void xrchart_write_xscale (cairo_t *, struct xrchart_geometry *,
- double min, double max, int ticks);
+ double min, double max);
/* Set the scale for the ordinate */
void xrchart_write_yscale (cairo_t *, struct xrchart_geometry *,
- double smin, double smax, int ticks);
+ double smin, double smax);
void xrchart_write_xlabel (cairo_t *, const struct xrchart_geometry *,
const char *label) ;
diff --git a/src/output/charts/boxplot-cairo.c
b/src/output/charts/boxplot-cairo.c
index 946b7d4..a35b277 100644
--- a/src/output/charts/boxplot-cairo.c
+++ b/src/output/charts/boxplot-cairo.c
@@ -143,7 +143,7 @@ xrchart_draw_boxplot (const struct chart_item
*chart_item, cairo_t *cr,
double box_width;
size_t i;
- xrchart_write_yscale (cr, geom, boxplot->y_min, boxplot->y_max, 5);
+ xrchart_write_yscale (cr, geom, boxplot->y_min, boxplot->y_max);
xrchart_write_title (cr, geom, "%s", chart_item->title);
box_width = (geom->axis[SCALE_ABSCISSA].data_max -
geom->axis[SCALE_ABSCISSA].data_min) / boxplot->n_boxes / 2.0;
diff --git a/src/output/charts/np-plot-cairo.c
b/src/output/charts/np-plot-cairo.c
index 334284b..e9a1e07 100644
--- a/src/output/charts/np-plot-cairo.c
+++ b/src/output/charts/np-plot-cairo.c
@@ -39,8 +39,8 @@ np_plot_chart_draw (const struct chart_item *chart_item,
cairo_t *cr,
xrchart_write_ylabel (cr, geom, _("Expected Normal"));
xrchart_write_xscale (cr, geom,
npp->x_lower - npp->slack,
- npp->x_upper + npp->slack, 5);
- xrchart_write_yscale (cr, geom, npp->y_first, npp->y_last, 5);
+ npp->x_upper + npp->slack);
+ xrchart_write_yscale (cr, geom, npp->y_first, npp->y_last);
data = casereader_clone (npp->data);
for (; (c = casereader_read (data)) != NULL; case_unref (c))
@@ -64,8 +64,8 @@ dnp_plot_chart_draw (const struct chart_item
*chart_item, cairo_t *cr,
xrchart_write_title (cr, geom, _("Detrended Normal Q-Q Plot of %s"),
chart_item->title);
xrchart_write_xlabel (cr, geom, _("Observed Value"));
xrchart_write_ylabel (cr, geom, _("Dev from Normal"));
- xrchart_write_xscale (cr, geom, dnpp->y_min, dnpp->y_max, 5);
- xrchart_write_yscale (cr, geom, dnpp->dns_min, dnpp->dns_max, 5);
+ xrchart_write_xscale (cr, geom, dnpp->y_min, dnpp->y_max);
+ xrchart_write_yscale (cr, geom, dnpp->dns_min, dnpp->dns_max);
data = casereader_clone (dnpp->data);
for (; (c = casereader_read (data)) != NULL; case_unref (c))
diff --git a/src/output/charts/plot-hist-cairo.c
b/src/output/charts/plot-hist-cairo.c
index 93133c2..006b392 100644
--- a/src/output/charts/plot-hist-cairo.c
+++ b/src/output/charts/plot-hist-cairo.c
@@ -15,7 +15,7 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
-
+#include "math/decimal.h"
#include "output/charts/plot-hist.h"
#include <float.h>
@@ -103,9 +103,20 @@ hist_draw_bar (cairo_t *cr, const struct
xrchart_geometry *geom,
cairo_stroke (cr);
if (label)
- draw_tick (cr, geom, SCALE_ABSCISSA, bins > 10,
- x_pos + width / 2.0, "%.*g",
- DBL_DIG + 1, (upper + lower) / 2.0);
+ {
+ struct decimal decupper;
+ struct decimal declower;
+ struct decimal middle;
+ decimal_from_double (&declower, lower);
+ decimal_from_double (&decupper, upper);
+ middle = declower;
+ decimal_add (&middle, &decupper);
+ decimal_int_divide (&middle, 2);
+ char *str = decimal_to_string (&middle);
+ draw_tick (cr, geom, SCALE_ABSCISSA, bins > 10,
+ x_pos + width / 2.0, "%s", str);
+ free (str);
+ }
}
void
@@ -129,7 +140,7 @@ xrchart_draw_histogram (const struct chart_item
*chart_item, cairo_t *cr,
bins = gsl_histogram_bins (h->gsl_hist);
- xrchart_write_yscale (cr, geom, 0, gsl_histogram_max_val (h->gsl_hist),
5);
+ xrchart_write_yscale (cr, geom, 0, gsl_histogram_max_val (h->gsl_hist));
for (i = 0; i < bins; i++)
{
diff --git a/src/output/charts/roc-chart-cairo.c
b/src/output/charts/roc-chart-cairo.c
index f6da350..2ac494e 100644
--- a/src/output/charts/roc-chart-cairo.c
+++ b/src/output/charts/roc-chart-cairo.c
@@ -37,8 +37,8 @@ xrchart_draw_roc (const struct chart_item *chart_item,
cairo_t *cr,
xrchart_write_xlabel (cr, geom, _("1 - Specificity"));
xrchart_write_ylabel (cr, geom, _("Sensitivity"));
- xrchart_write_xscale (cr, geom, 0, 1, 5);
- xrchart_write_yscale (cr, geom, 0, 1, 5);
+ xrchart_write_xscale (cr, geom, 0, 1);
+ xrchart_write_yscale (cr, geom, 0, 1);
if ( rc->reference )
{
diff --git a/src/output/charts/scatterplot-cairo.c
b/src/output/charts/scatterplot-cairo.c
index b555a12..f25b2cd 100644
--- a/src/output/charts/scatterplot-cairo.c
+++ b/src/output/charts/scatterplot-cairo.c
@@ -51,8 +51,8 @@ xrchart_draw_scatterplot (const struct chart_item
*chart_item, cairo_t *cr,
xrchart_write_xscale (cr, geom,
spc->x_min,
- spc->x_max, 5);
- xrchart_write_yscale (cr, geom, spc->y_min, spc->y_max, 5);
+ spc->x_max);
+ xrchart_write_yscale (cr, geom, spc->y_min, spc->y_max);
xrchart_write_title (cr, geom, _("Scatterplot %s"), chart_item->title);
xrchart_write_xlabel (cr, geom, var_to_string(spc->xvar));
xrchart_write_ylabel (cr, geom, var_to_string(spc->yvar));
diff --git a/src/output/charts/scree-cairo.c
b/src/output/charts/scree-cairo.c
index f9765c6..db0ce13 100644
--- a/src/output/charts/scree-cairo.c
+++ b/src/output/charts/scree-cairo.c
@@ -44,8 +44,8 @@ xrchart_draw_scree (const struct chart_item *chart_item,
cairo_t *cr,
else
max = fabs (min);
- xrchart_write_yscale (cr, geom, 0, max, max);
- xrchart_write_xscale (cr, geom, 0, rc->eval->size + 1, rc->eval->size +
1);
+ xrchart_write_yscale (cr, geom, 0, max);
+ xrchart_write_xscale (cr, geom, 0, rc->eval->size + 1);
xrchart_vector_start (cr, geom, "");
for (i = 0 ; i < rc->eval->size; ++i)
diff --git a/src/output/charts/spreadlevel-cairo.c
b/src/output/charts/spreadlevel-cairo.c
index eeb3c81..13b7a7a 100644
--- a/src/output/charts/spreadlevel-cairo.c
+++ b/src/output/charts/spreadlevel-cairo.c
@@ -39,8 +39,8 @@ xrchart_draw_spreadlevel (const struct chart_item
*chart_item, cairo_t *cr,
xrchart_write_ylabel (cr, geom, _("Spread"));
- xrchart_write_xscale (cr, geom, sl->x_lower, sl->x_upper, 5);
- xrchart_write_yscale (cr, geom, sl->y_lower, sl->y_upper, 5);
+ xrchart_write_xscale (cr, geom, sl->x_lower, sl->x_upper);
+ xrchart_write_yscale (cr, geom, sl->y_lower, sl->y_upper);
for (i = 0 ; i < sl->n_data; ++i)
{
diff --git a/tests/automake.mk b/tests/automake.mk
index 31a7877..4cb4283 100644
--- a/tests/automake.mk
+++ b/tests/automake.mk
@@ -31,6 +31,9 @@ check_PROGRAMS += \
tests/libpspp/tower-test \
tests/libpspp/u8-istream-test \
tests/libpspp/zip-test \
+ tests/math/chart-geometry-test \
+ tests/math/chart-get-scale-test \
+ tests/math/decimal-test \
tests/output/render-test \
tests/ui/syntax-gen-test
@@ -209,6 +212,32 @@ tests_libpspp_zip_test_LDADD = \
src/libpspp-core.la \
gl/libgl.la
+check_PROGRAMS += tests/math/chart-geometry-test
+tests_math_chart_geometry_test_SOURCES = tests/math/chart-geometry-test.c
+tests_math_chart_geometry_test_LDADD = \
+ src/math/libpspp-math.la \
+ src/libpspp/liblibpspp.la \
+ src/libpspp-core.la \
+ gl/libgl.la
+
+check_PROGRAMS += tests/math/chart-get-scale-test
+tests_math_chart_get_scale_test_SOURCES =
tests/math/chart-get-scale-test.c
+tests_math_chart_get_scale_test_LDADD = \
+ src/math/libpspp-math.la \
+ src/libpspp/liblibpspp.la \
+ src/libpspp-core.la \
+ gl/libgl.la
+
+
+check_PROGRAMS += tests/math/decimal-test
+tests_math_decimal_test_SOURCES = tests/math/decimal-test.c
+tests_math_decimal_test_LDADD = \
+ src/math/libpspp-math.la \
+ src/libpspp/liblibpspp.la \
+ src/libpspp-core.la \
+ gl/libgl.la
+
+
check_PROGRAMS += tests/output/render-test
tests_output_render_test_SOURCES = tests/output/render-test.c
tests_output_render_test_LDADD = \
@@ -367,6 +396,8 @@ TESTSUITE_AT = \
tests/libpspp/tower.at \
tests/libpspp/u8-istream.at \
tests/libpspp/zip.at \
+ tests/math/chart-geometry.at \
+ tests/math/decimal.at \
tests/math/moments.at \
tests/math/randist.at \
tests/output/ascii.at \
diff --git a/src/math/chart-geometry.c b/tests/math/chart-geometry-test.c
similarity index 50%
copy from src/math/chart-geometry.c
copy to tests/math/chart-geometry-test.c
index eb0ad92..a713cdb 100644
--- a/src/math/chart-geometry.c
+++ b/tests/math/chart-geometry-test.c
@@ -1,5 +1,5 @@
/* PSPP - a program for statistical analysis.
- Copyright (C) 2004 Free Software Foundation, Inc.
+ Copyright (C) 2015 Free Software Foundation, Inc.
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
@@ -15,43 +15,47 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
-#include <math.h>
-#include <float.h>
-
-#include "chart-geometry.h"
-
-static const double standard_ticks[] = {1, 2, 5, 10};
-
-
-/* Adjust tick to be a sensible value
- ie: ... 0.1,0.2,0.5, 1,2,5, 10,20,50 ... */
-double
-chart_rounded_tick (double tick)
+#include <stdlib.h>
+#include "math/chart-geometry.h"
+#include "math/decimal.h"
+
+
+const double in[20] =
+ {
+ 0.00648687,
+ 728815,
+ 8.14431e-07,
+ 77611.4,
+ 3.33497,
+ 180.426,
+ 0.676168,
+ 2.00744e+08,
+ 14099.3,
+ 19.5186,
+ 1.17473e-07,
+ 166337,
+ 0.00163644,
+ 1.94724e-09,
+ 2.31564e-06,
+ 3.10674e+06,
+ 5.10314e-05,
+ 1.95101,
+ 1.40884e+09,
+ 78217.6
+ };
+
+int
+main ()
{
int i;
-
- double diff = DBL_MAX;
- double t = tick;
-
- double factor;
-
- /* Avoid arithmetic problems with very small values */
- if (fabs (tick) < DBL_EPSILON)
- return 0;
-
- factor = pow (10,ceil (log10 (standard_ticks[0] / tick)));
-
- for (i = 3 ; i >= 0 ; --i)
+ for (i = 0; i < 20; ++i)
{
- const double d = fabs (tick - standard_ticks[i] / factor);
-
- if ( d < diff )
- {
- diff = d;
- t = standard_ticks[i] / factor ;
- }
+ struct decimal dout;
+ chart_rounded_tick (in[i], &dout);
+
+ printf ("%g %s\n", in[i], decimal_to_string (&dout));
}
- return t;
+ return 0;
}
diff --git a/tests/math/chart-geometry.at b/tests/math/chart-geometry.at
new file mode 100644
index 0000000..14e1584
--- /dev/null
+++ b/tests/math/chart-geometry.at
@@ -0,0 +1,35 @@
+AT_BANNER([Chart Geometry])
+
+AT_SETUP([Chart Rounding])
+
+AT_CHECK([../../math/chart-geometry-test], [0], [dnl
+0.00648687 0.005
+728815 500000
+8.14431e-07 0.000001
+77611.4 100000
+3.33497 2
+180.426 200
+0.676168 0.5
+2.00744e+08 200000000
+14099.3 10000
+19.5186 20
+1.17473e-07 0.0000001
+166337 200000
+0.00163644 0.002
+1.94724e-09 0.000000002
+2.31564e-06 0.000002
+3.10674e+06 2000000
+5.10314e-05 0.00005
+1.95101 2
+1.40884e+09 1000000000
+78217.6 100000
+])
+
+AT_CLEANUP
+
+
+AT_SETUP([Chart Scale])
+
+AT_CHECK([../../math/chart-get-scale-test], [0], [ignore])
+
+AT_CLEANUP
diff --git a/tests/math/chart-get-scale-test.c
b/tests/math/chart-get-scale-test.c
new file mode 100644
index 0000000..11a572d
--- /dev/null
+++ b/tests/math/chart-get-scale-test.c
@@ -0,0 +1,96 @@
+/* PSPP - a program for statistical analysis.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+
+ 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 3 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, see <http://www.gnu.org/licenses/>.
*/
+
+#include <config.h>
+
+#include <stdlib.h>
+#include <stdbool.h>
+#include <assert.h>
+#include <string.h>
+
+#include "math/decimal.h"
+#include <limits.h>
+#include <float.h>
+#include <math.h>
+
+void
+dump_scale (const struct decimal *low, const struct decimal *interval,
int n_ticks)
+{
+ int i;
+ struct decimal tick = *low;
+ for (i = 0; i <= n_ticks; ++i)
+ {
+ printf ("Tick %d: %s (%g)\n", i, decimal_to_string (&tick),
decimal_to_double (&tick));
+ decimal_add (&tick, interval);
+ }
+}
+
+
+void
+test_range (double low, double high)
+{
+ int n_ticks = 0;
+ struct decimal interval;
+ struct decimal lower;
+
+
+ chart_get_scale (high, low,
+ &lower, &interval, &n_ticks);
+
+ assert (n_ticks > 0);
+ assert (n_ticks < 12);
+
+ // dump_scale (&lower, &interval, n_ticks);
+
+ assert (decimal_to_double (&lower) <= low);
+
+ {
+ struct decimal l = lower;
+ decimal_add (&l, &interval);
+ assert (decimal_to_double (&l) > low);
+ }
+
+ {
+ struct decimal i = interval;
+
+ decimal_int_multiply (&i, n_ticks - 1);
+
+ decimal_add (&i, &lower);
+
+ /* i now contains the upper bound minus one tick */
+ assert (decimal_to_double (&i) < high);
+
+ decimal_add (&i, &interval);
+
+ assert (decimal_to_double (&i) >= high);
+ }
+
+}
+
+
+int
+main (int argc, char **argv)
+{
+ test_range (0.2, 11);
+ test_range (-0.2, 11);
+ test_range (-10, 0.2);
+ test_range (-10, -0.2);
+
+ test_range (102, 50030);
+ test_range (0.00102, 0.0050030);
+
+ return 0;
+}
diff --git a/tests/math/decimal-test.c b/tests/math/decimal-test.c
new file mode 100644
index 0000000..aa4d244
--- /dev/null
+++ b/tests/math/decimal-test.c
@@ -0,0 +1,310 @@
+/* PSPP - a program for statistical analysis.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+
+ 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 3 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, see <http://www.gnu.org/licenses/>.
*/
+
+#include <config.h>
+
+#include <stdlib.h>
+#include <stdbool.h>
+#include <assert.h>
+#include <string.h>
+
+#include "math/decimal.h"
+#include <limits.h>
+#include <float.h>
+#include <math.h>
+
+/* Canonicalise a string holding the decimal representation of a number.
+ For example, leading zeros left of the decimal point are removed, as
are
+ trailing zeros to the right.
+
+ This function is used purely for testing, and need not and is not
intended
+ to be efficient.
+ */
+char *
+canonicalise_string (const char *s)
+{
+ char *out;
+ char *dot = NULL;
+ bool negative = false;
+ char *p;
+ char *temp = malloc (strlen (s) + 3);
+ char *last_leading_zero = NULL;
+
+ /* Strip leading - if present */
+ if (*s == '-')
+ {
+ negative = true;
+ s++;
+ }
+
+ strcpy (temp, "00");
+ strcat (temp, s);
+
+ char *first_trailing_zero = NULL;
+ char *significant_digit = NULL;
+ for (p = temp; *p; p++)
+ {
+ if (*p == '0' && dot == NULL && significant_digit == NULL)
+ last_leading_zero = p;
+
+ if (*p == '0' && first_trailing_zero == NULL)
+ first_trailing_zero = p;
+
+ if (*p == '.')
+ {
+ dot = p;
+ first_trailing_zero = NULL;
+ }
+
+ if (*p >= '1' && *p <= '9')
+ {
+ significant_digit = p;
+ first_trailing_zero = NULL;
+ }
+ }
+
+ if (first_trailing_zero && dot)
+ *first_trailing_zero = '\0';
+
+ if (last_leading_zero)
+ {
+ /* Strip leading zeros */
+ out = last_leading_zero + 1;
+
+ /* But if we now start with . put a zero back again */
+ if (dot == last_leading_zero + 1)
+ out--;
+ }
+
+
+ if (negative)
+ {
+ out--;
+ out[0] = '-';
+ }
+
+ if (!significant_digit)
+ {
+ *out = '0';
+ *(out+1) = '\0';
+ }
+
+
+ return out;
+}
+
+
+/* Tests both the decimal_to_string function, and the
decimal_input_from_string
+ function */
+void
+test_run (const char *input)
+ {
+ struct decimal test;
+ struct decimal number;
+ decimal_init_from_string (&number, input);
+
+ char *s = decimal_to_string (&number);
+ char *canon = canonicalise_string (input);
+ if (0 != strcmp (canon, s))
+ {
+ fprintf (stdout, "\"%s\" does not match \n\"%s\"\n", canon, s);
+ exit (1);
+ }
+
+ decimal_init_from_string (&test, s);
+ assert (0 == decimal_cmp (&test, &number));
+
+ free (s);
+ }
+
+
+void
+test_can (const char *in, const char *soll)
+{
+ char *ist = canonicalise_string (in);
+ if (0 != strcmp (soll, ist))
+ {
+ printf ("\"%s\" canonicalises to \"%s\" (should be \"%s\")\n", in,
ist, soll);
+ }
+}
+
+
+void
+dump_scale (const struct decimal *low, const struct decimal *interval,
int n_ticks)
+{
+ int i;
+ struct decimal tick = *interval;
+ printf ("Lowest: %s\n", decimal_to_string (low));
+ for (i = 0; i <= n_ticks; ++i)
+ {
+ printf ("Tick %d: %s (%g)\n", i, decimal_to_string (&tick),
decimal_to_double (&tick));
+ decimal_add (&tick, interval);
+ }
+}
+
+
+
+void
+test_ceil (double x)
+{
+ struct decimal dx;
+ decimal_from_double (&dx, x);
+ int act = decimal_ceil (&dx);
+ int expected = ceil (x);
+
+ assert (act == expected);
+}
+
+void
+test_floor (double x)
+{
+ struct decimal dx;
+ decimal_from_double (&dx, x);
+ int act = decimal_floor (&dx);
+ int expected = floor (x);
+
+ assert (act == expected);
+}
+
+
+void
+test_addition (const struct decimal *one_, const struct decimal *two)
+{
+ struct decimal one = *one_;
+ double d1 = decimal_to_double (&one);
+ double d2 = decimal_to_double (two);
+
+ decimal_add (&one, two);
+
+ double dsum = decimal_to_double (&one);
+
+ char sdsum1[256];
+ char sdsum2[256];
+
+ snprintf (sdsum1, 256, "%s", decimal_to_string (&one));
+ snprintf (sdsum2, 256, "%g", dsum);
+
+ assert (strcmp (sdsum1, sdsum2) == 0);
+}
+
+void
+test_multiplication (const struct decimal *one_, const struct decimal
*two)
+{
+ struct decimal one = *one_;
+ double d1 = decimal_to_double (&one);
+ double d2 = decimal_to_double (two);
+
+ decimal_multiply (&one, two);
+
+ double dprod = decimal_to_double (&one);
+
+ assert (dprod == d1 * d2);
+}
+
+
+int
+main (int argc, char **argv)
+{
+ /* Test that our canonicalise function works for all corner cases we
+ can think of. */
+
+ test_can ("500", "500");
+ test_can ("5", "5");
+ test_can ("-3", "-3");
+ test_can ("-3.001", "-3.001");
+ test_can ("-03.001", "-3.001");
+ test_can ("-.0301", "-0.0301");
+ test_can ("0314.09", "314.09");
+ test_can ("0314.090", "314.09");
+ test_can ("0314.0900340", "314.090034");
+ test_can ("0.0", "0");
+ test_can ("0.", "0");
+ test_can (".0", "0");
+ test_can ("-.1", "-0.1");
+ test_can (".090", "0.09");
+ test_can ("03410.098700", "3410.0987");
+ test_can ("-03410.098700", "-3410.0987");
+
+ /* Test the conversion functions */
+ test_run ("-90000");
+ test_run ("-3");
+ test_run ("50001");
+ test_run ("500");
+ test_run ("350");
+ test_run ("050");
+ test_run ("4");
+ test_run ("0");
+ test_run (".45");
+ test_run ("-.45");
+ test_run ("666666666");
+ test_run ("6000000000");
+ test_run ("0.000000005");
+ test_run ("0.00000000000000000000000000000000000000005");
+ test_run ("0.0234");
+ test_run ("0.234");
+ test_run ("-0123.45600");
+
+ test_ceil (5.21);
+ test_ceil (-4.32);
+ test_ceil (0);
+ test_ceil (0.0009);
+
+ test_floor (4.09);
+ test_floor (-4.09);
+ test_floor (0);
+ test_floor (0.004);
+
+
+ {
+ struct decimal high = {2, 0};
+ struct decimal low = {2, -1};
+
+ test_addition (&high, &low);
+ }
+
+
+ {
+ struct decimal high = {10, 0};
+ struct decimal low = {2, -1};
+
+ test_addition (&high, &low);
+ }
+
+
+ {
+ struct decimal high = {10, 0};
+ struct decimal low = {-2, -1};
+
+ test_addition (&high, &low);
+ }
+
+ {
+ struct decimal high = {12, -5};
+ struct decimal low = {-2, -1};
+
+ test_addition (&high, &low);
+ }
+
+ {
+ struct decimal high = {-112, -1};
+ struct decimal low = {2, -1};
+
+ test_addition (&high, &low);
+ }
+
+
+ return 0;
+}
diff --git a/tests/math/decimal.at b/tests/math/decimal.at
new file mode 100644
index 0000000..b3e71d1
--- /dev/null
+++ b/tests/math/decimal.at
@@ -0,0 +1,7 @@
+AT_BANNER([Decimal floating point])
+
+AT_SETUP([Decimal float])
+
+AT_CHECK([../../math/decimal-test], [0], [ignore])
+
+AT_CLEANUP
--
1.7.10.4
--
PGP Public key ID: 1024D/2DE827B3
fingerprint = 8797 A26D 0854 2EAB 0285 A290 8A67 719C 2DE8 27B3
See http://sks-keyservers.net or any PGP keyserver for public key.
signature.asc
Description: Digital signature