From: John Darrington Date: Tue, 6 Jan 2015 06:25:21 +0000 (+0100) Subject: New module to perform decimal floating point arithmetic for charts. X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?p=pspp;a=commitdiff_plain;h=aad0ae4913ecd01ccc954f8828623ad5da35cd1b New module to perform decimal floating point arithmetic for charts. 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 Reviewed-by: Ben Pfaff --- diff --git a/src/math/automake.mk b/src/math/automake.mk index cdcc2c8cb7..f40aac57e1 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 eb0ad9282a..35380a6929 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 #include #include +#include #include "chart-geometry.h" +#include "decimal.h" +#include -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 ef481feb40..c543086a39 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 0000000000..0655834225 --- /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 . */ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#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 M */ +void +decimal_int_multiply (struct decimal *dest, ord_t m) +{ + if (m != 0) + while (llabs (dest->ordinate) > llabs (ORD_MAX / m)) + { + dest->ordinate /= 10; + dest->mantissa++; + } + + 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 0000000000..b5740f60ac --- /dev/null +++ b/src/math/decimal.h @@ -0,0 +1,114 @@ +/* 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 . */ + + +#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 +#include + +#include + +#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_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 b00067210f..3c324eff24 100644 --- a/src/math/histogram.c +++ b/src/math/histogram.c @@ -17,6 +17,7 @@ #include #include "math/histogram.h" +#include "math/decimal.h" #include #include @@ -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 8886b2642d..01c5d5c044 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 #include "output/cairo-chart.h" +#include "math/decimal.h" +#include "math/chart-geometry.h" #include #include @@ -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 edf4ed1148..27000cb56f 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 946b7d443c..a35b277b8c 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 334284ba8f..e9a1e07033 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 93133c2e9e..006b39225e 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 . */ #include - +#include "math/decimal.h" #include "output/charts/plot-hist.h" #include @@ -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 f6da350667..2ac494e3cd 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 b555a12088..f25b2cde4e 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 f9765c69f2..db0ce1370f 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 eeb3c813be..13b7a7a3b5 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 31a7877d07..4cb4283e9d 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/tests/math/chart-geometry-test.c b/tests/math/chart-geometry-test.c new file mode 100644 index 0000000000..a713cdb0d5 --- /dev/null +++ b/tests/math/chart-geometry-test.c @@ -0,0 +1,61 @@ +/* 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 . */ + +#include +#include +#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; + for (i = 0; i < 20; ++i) + { + struct decimal dout; + chart_rounded_tick (in[i], &dout); + + printf ("%g %s\n", in[i], decimal_to_string (&dout)); + } + + return 0; +} + diff --git a/tests/math/chart-geometry.at b/tests/math/chart-geometry.at new file mode 100644 index 0000000000..14e15844ed --- /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 0000000000..11a572deeb --- /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 . */ + +#include + +#include +#include +#include +#include + +#include "math/decimal.h" +#include +#include +#include + +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 0000000000..3e1e64eaaf --- /dev/null +++ b/tests/math/decimal-test.c @@ -0,0 +1,347 @@ +/* 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 . */ + +#include + +#include +#include +#include +#include + +#include "math/decimal.h" +#include +#include +#include + +/* 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 *d, int m) +{ + char b1[256]; + char b2[256]; + struct decimal dest = *d; + double x = decimal_to_double (&dest); + + decimal_int_multiply (&dest, m); + + double y = decimal_to_double (&dest); + + snprintf (b1, 256, "%g", m * x); + snprintf (b2, 256, "%g", y); + assert (0 == strcmp (b1, b2)); +} + + + +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); + } + + + { + struct decimal m = {10, 0}; + + test_multiplication (&m, 11); + } + + + { + struct decimal m = {ORD_MAX - 2, 0}; + + test_multiplication (&m, 11); + } + + + { + struct decimal m = {34, 0}; + + test_multiplication (&m, 0); + } + + { + struct decimal m = {34, -20}; + + test_multiplication (&m, 33); + } + + { + struct decimal m = {304, 2}; + + test_multiplication (&m, -33); + } + + return 0; +} diff --git a/tests/math/decimal.at b/tests/math/decimal.at new file mode 100644 index 0000000000..b3e71d1fab --- /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