New module to perform decimal floating point arithmetic for charts.

[pspp] / src / math / chart-geometry.c

diff --git a/src/math/chart-geometry.c b/src/math/chart-geometry.c
index eb0ad9282a854885768472319a899dbfba6c14ea..35380a6929b6c4529e0fe2d5579e03cad9b25f23 100644 (file)
--- 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;
+       }
+    }
+}