lib/arraystats: re-enable "Discont" algorithm (#5529)

Addresses heap buffer overflow in previously disabled 'Discont'
algorithm.

The immediate cause was too small allocation of variable `num` in
`AS_class_discont()`, which was fixed by an increase of its size
(broadly mirroring the original Fortran code).

This fix was complemented by some code related improvements:

- Code clean up and restructure of `AS_class_discont()`
- Re-enable "dis" algorithm in v.class and d.vect.thematic
- Remove realloc of classbreak, which wasn’t necessary and
  would need complicated solution to work

See: https://lists.osgeo.org/pipermail/grass-dev/2008-July/038951.html

Fixes: #5486

Author: nilason

display/d.vect.thematic/main.c

@@ -119,11 +119,10 @@ int main(int argc, char **argv)
     algo_opt->options = "int,std,qua,equ,dis";
     algo_opt->description = _("Algorithm to use for classification");
     desc = NULL;
-    G_asprintf(&desc, "int;%s;std;%s;qua;%s;equ;%s", _("simple intervals"),
-               _("standard deviations"), _("quantiles"),
-               _("equiprobable (normal distribution)"));
+    G_asprintf(&desc, "int;%s;std;%s;qua;%s;equ;%s;dis;%s",
+               _("simple intervals"), _("standard deviations"), _("quantiles"),
+               _("equiprobable (normal distribution)"), _("discontinuities"));
     algo_opt->descriptions = desc;
-    /*currently disabled because of bugs       "dis;discontinuities"); */
     algo_opt->guisection = _("Classes");
 
     nbclass_opt = G_define_option();

include/grass/defs/arraystats.h

@@ -2,17 +2,17 @@
 #define GRASS_ARRAYSTATSDEFS_H
 
 /* basic.c */
-void AS_eqdrt(double[], double[], int, int, double *);
-void AS_basic_stats(double *, int, struct GASTATS *);
+void AS_eqdrt(double[], double[], int, int, double *, double *, double *);
+void AS_basic_stats(const double[], int, struct GASTATS *);
 
 /* class.c */
 int AS_option_to_algorithm(const struct Option *);
-double AS_class_apply_algorithm(int, double *, int, int *, double *);
-int AS_class_interval(double *, int, int, double *);
-int AS_class_quant(double *, int, int, double *);
-double AS_class_discont(double *, int, int, double *);
-double AS_class_stdev(double *, int, int, double *);
-int AS_class_equiprob(double *, int, int *, double *);
-int AS_class_frequencies(double *, int, int, double *, int *);
+double AS_class_apply_algorithm(int, const double[], int, int *, double[]);
+int AS_class_interval(const double[], int, int, double[]);
+int AS_class_quant(const double[], int, int, double[]);
+double AS_class_discont(const double[], int, int, double[]);
+double AS_class_stdev(const double[], int, int, double[]);
+int AS_class_equiprob(const double[], int, int *, double[]);
+int AS_class_frequencies(const double[], int, int, double[], int[]);
 
 #endif

lib/arraystats/basic.c

@@ -1,8 +1,10 @@
 #include <math.h>
+
 #include <grass/arraystats.h>
+#include <grass/gis.h>
 
 /*provides basic univar stats */
-void AS_basic_stats(double *data, int count, struct GASTATS *stats)
+void AS_basic_stats(const double data[], int count, struct GASTATS *stats)
 {
     int i = 1;
     double sum = 0, sumsq = 0, sumabs = 0;
@@ -34,31 +36,40 @@ void AS_basic_stats(double *data, int count, struct GASTATS *stats)
     return;
 }
 
-void AS_eqdrt(double vectx[], double vecty[], int i1, int i2, double *vabc)
+void AS_eqdrt(double vectx[], double vecty[], int i1, int i2, double *a,
+              double *b, double *c)
 {
-    double bn = 0, bd = 0, x1 = 0, y1 = 0;
+    double x1 = vectx[i1];
+    double y1 = vecty[i1];
+    double x2 = vectx[i2];
+    double y2 = vecty[i2];
 
-    vabc[0] = 0;
-    vabc[1] = 0;
-    vabc[2] = 0;
     if (i1 == 0) {
-        x1 = 0;
-        y1 = 0;
+        x1 = 0.0;
+        y1 = 0.0;
     }
     else {
         x1 = vectx[i1];
         y1 = vecty[i1];
     }
-    bn = y1 - vecty[i2];
-    bd = x1 - vectx[i2];
-    if (bd != 0) {
-        vabc[1] = bn / bd;
-        vabc[0] = y1 - vabc[1] * x1;
-        return;
+
+    *a = 0.0;
+    *b = 0.0;
+    *c = 0.0;
+
+    double bn = y1 - y2;
+    double bd = x1 - x2;
+
+    if (fabs(bd) < GRASS_EPSILON) {
+        if (fabs(bn) < GRASS_EPSILON) {
+            G_debug(3, "Points are equal\n");
+        }
+        else {
+            *c = x1;
+        }
+    }
+    else {
+        *b = bn / bd;
+        *a = y1 - *b * x1;
     }
-    if (bn != 0)
-        vabc[2] = x1;
-    else
-        G_debug(3, "Points are equal\n");
-    return;
 }