LCOV - code coverage report
Current view: top level - gdk - gdk_aggr.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 1246 1835 67.9 %
Date: 2020-06-29 20:00:14 Functions: 44 45 97.8 %

          Line data    Source code
       1             : /*
       2             :  * This Source Code Form is subject to the terms of the Mozilla Public
       3             :  * License, v. 2.0.  If a copy of the MPL was not distributed with this
       4             :  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
       5             :  *
       6             :  * Copyright 1997 - July 2008 CWI, August 2008 - 2020 MonetDB B.V.
       7             :  */
       8             : 
       9             : #include "monetdb_config.h"
      10             : #include "gdk.h"
      11             : #include "gdk_private.h"
      12             : #include "gdk_calc_private.h"
      13             : #include <math.h>
      14             : 
      15             : /* grouped aggregates
      16             :  *
      17             :  * The following functions take two to four input BATs and produce a
      18             :  * single output BAT.
      19             :  *
      20             :  * The input BATs are
      21             :  * - b, a dense-headed BAT with the values to work on in the tail;
      22             :  * - g, a dense-headed BAT, aligned with b, with group ids (OID) in
      23             :  *   the tail;
      24             :  * - e, optional but recommended, a dense-headed BAT with the list of
      25             :  *   group ids in the head(!) (the tail is completely ignored);
      26             :  * - s, optional, a dense-headed bat with a list of candidate ids in
      27             :  *   the tail.
      28             :  *
      29             :  * The tail values of s refer to the head of b and g.  Only entries at
      30             :  * the specified ids are taken into account for the grouped
      31             :  * aggregates.  All other values are ignored.  s is compatible with
      32             :  * the result of BATselect().
      33             :  *
      34             :  * If e is not specified, we need to do an extra scan over g to find
      35             :  * out the range of the group ids that are used.  e is defined in such
      36             :  * a way that it can be either the extents or the histo result from
      37             :  * BATgroups().
      38             :  *
      39             :  * All functions calculate grouped aggregates.  There are as many
      40             :  * groups as there are entries in e.  If e is not specified, the
      41             :  * number of groups is equal to the difference between the maximum and
      42             :  * minimum values in g.
      43             :  *
      44             :  * If a group is empty, the result for that group is nil.
      45             :  *
      46             :  * If there is overflow during the calculation of an aggregate, the
      47             :  * whole operation fails if abort_on_error is set to non-zero,
      48             :  * otherwise the result of the group in which the overflow occurred is
      49             :  * nil.
      50             :  *
      51             :  * If skip_nils is non-zero, a nil value in b is ignored, otherwise a
      52             :  * nil in b results in a nil result for the group.
      53             :  */
      54             : 
      55             : /* helper function
      56             :  *
      57             :  * This function finds the minimum and maximum group id (and the
      58             :  * number of groups) and initializes the variables for candidates
      59             :  * selection.
      60             :  */
      61             : const char *
      62       22969 : BATgroupaggrinit(BAT *b, BAT *g, BAT *e, BAT *s,
      63             :                  /* outputs: */
      64             :                  oid *minp, oid *maxp, BUN *ngrpp,
      65             :                  struct canditer *ci, BUN *ncand)
      66             : {
      67       22969 :         oid min, max;
      68       22969 :         BUN i, ngrp;
      69       22969 :         const oid *restrict gids;
      70             : 
      71       22969 :         if (b == NULL)
      72             :                 return "b must exist";
      73       22969 :         if (g) {
      74       15865 :                 if (BATcount(b) != BATcount(g) ||
      75        8726 :                     (BATcount(b) != 0 && b->hseqbase != g->hseqbase))
      76             :                         return "b and g must be aligned";
      77       15863 :                 assert(BATttype(g) == TYPE_oid);
      78             :         }
      79       22967 :         if (g == NULL) {
      80             :                 min = 0;
      81             :                 max = 0;
      82             :                 ngrp = 1;
      83       15863 :         } else if (e == NULL) {
      84             :                 /* we need to find out the min and max of g */
      85           0 :                 PROPrec *prop;
      86             : 
      87           0 :                 prop = BATgetprop(g, GDK_MAX_VALUE);
      88           0 :                 if (prop) {
      89           0 :                         assert(prop->v.vtype == TYPE_oid);
      90           0 :                         min = 0; /* just assume it starts at 0 */
      91           0 :                         max = prop->v.val.oval;
      92             :                 } else {
      93           0 :                         min = oid_nil;  /* note that oid_nil > 0! (unsigned) */
      94           0 :                         max = 0;
      95           0 :                         if (BATtdense(g)) {
      96           0 :                                 min = g->tseqbase;
      97           0 :                                 max = g->tseqbase + BATcount(g) - 1;
      98           0 :                         } else if (g->tsorted) {
      99           0 :                                 gids = (const oid *) Tloc(g, 0);
     100             :                                 /* find first non-nil */
     101           0 :                                 for (i = 0, ngrp = BATcount(g); i < ngrp; i++) {
     102           0 :                                         if (!is_oid_nil(gids[i])) {
     103             :                                                 min = gids[i];
     104             :                                                 break;
     105             :                                         }
     106             :                                 }
     107           0 :                                 if (!is_oid_nil(min)) {
     108             :                                         /* found a non-nil, max must be last
     109             :                                          * value (and there is one!) */
     110           0 :                                         max = gids[BUNlast(g) - 1];
     111             :                                 }
     112             :                         } else {
     113             :                                 /* we'll do a complete scan */
     114           0 :                                 gids = (const oid *) Tloc(g, 0);
     115           0 :                                 for (i = 0, ngrp = BATcount(g); i < ngrp; i++) {
     116           0 :                                         if (!is_oid_nil(gids[i])) {
     117           0 :                                                 if (gids[i] < min)
     118             :                                                         min = gids[i];
     119           0 :                                                 if (gids[i] > max)
     120             :                                                         max = gids[i];
     121             :                                         }
     122             :                                 }
     123             :                                 /* note: max < min is possible if all groups
     124             :                                  * are nil (or BATcount(g)==0) */
     125             :                         }
     126             :                 }
     127           0 :                 ngrp = max < min ? 0 : max - min + 1;
     128             :         } else {
     129       15863 :                 ngrp = BATcount(e);
     130       15863 :                 min = e->hseqbase;
     131       15863 :                 max = e->hseqbase + ngrp - 1;
     132             :         }
     133       22967 :         *minp = min;
     134       22967 :         *maxp = max;
     135       22967 :         *ngrpp = ngrp;
     136             : 
     137       22967 :         *ncand = canditer_init(ci, b, s);
     138             : 
     139       22969 :         return NULL;
     140             : }
     141             : 
     142             : /* ---------------------------------------------------------------------- */
     143             : /* sum */
     144             : 
     145             : static inline bool
     146        5297 : samesign(double x, double y)
     147             : {
     148        5297 :         return (x >= 0) == (y >= 0);
     149             : }
     150             : 
     151             : /* Add two values, producing the sum and the remainder due to limited
     152             :  * range of floating point arithmetic.  This function depends on the
     153             :  * fact that the sum returns INFINITY in *hi of the correct sign
     154             :  * (i.e. isinf() returns TRUE) in case of overflow. */
     155             : static inline void
     156      273899 : twosum(volatile double *hi, volatile double *lo, double x, double y)
     157             : {
     158      273899 :         volatile double yr;
     159             : 
     160      273899 :         assert(fabs(x) >= fabs(y));
     161             : 
     162      273899 :         *hi = x + y;
     163      273899 :         yr = *hi - x;
     164      273899 :         *lo = y - yr;
     165      273899 : }
     166             : 
     167             : static inline void
     168       79122 : exchange(double *x, double *y)
     169             : {
     170       79122 :         double t = *x;
     171       79122 :         *x = *y;
     172       79122 :         *y = t;
     173       79122 : }
     174             : 
     175             : /* this function was adapted from https://bugs.python.org/file10357/msum4.py */
     176             : BUN
     177        1301 : dofsum(const void *restrict values, oid seqb,
     178             :        struct canditer *restrict ci, BUN ncand,
     179             :        void *restrict results, BUN ngrp, int tp1, int tp2,
     180             :        const oid *restrict gids,
     181             :        oid min, oid max, bool skip_nils, bool abort_on_error,
     182             :        bool nil_if_empty)
     183             : {
     184        1301 :         struct pergroup {
     185             :                 int npartials;
     186             :                 int maxpartials;
     187             :                 bool valseen;
     188             : #ifdef INFINITES_ALLOWED
     189             :                 float infs;
     190             : #else
     191             :                 int infs;
     192             : #endif
     193             :                 double *partials;
     194             :         } *pergroup;
     195        1301 :         BUN listi;
     196        1301 :         int parti;
     197        1301 :         int i;
     198        1301 :         BUN grp;
     199        1301 :         double x, y;
     200        1301 :         volatile double lo, hi;
     201        1301 :         double twopow = pow((double) FLT_RADIX, (double) (DBL_MAX_EXP - 1));
     202        1301 :         BUN nils = 0;
     203        1301 :         volatile flt f;
     204             : 
     205             :         /* we only deal with the two floating point types */
     206        1301 :         assert(tp1 == TYPE_flt || tp1 == TYPE_dbl);
     207        1301 :         assert(tp2 == TYPE_flt || tp2 == TYPE_dbl);
     208             :         /* if input is dbl, then so it output */
     209        1301 :         assert(tp2 == TYPE_flt || tp1 == TYPE_dbl);
     210             :         /* if no gids, then we have a single group */
     211        1301 :         assert(ngrp == 1 || gids != NULL);
     212        1301 :         if (gids == NULL || ngrp == 1) {
     213        1270 :                 min = max = 0;
     214        1270 :                 ngrp = 1;
     215        1270 :                 gids = NULL;
     216             :         }
     217        1301 :         pergroup = GDKmalloc(ngrp * sizeof(*pergroup));
     218        1301 :         if (pergroup == NULL)
     219             :                 return BUN_NONE;
     220       53635 :         for (grp = 0; grp < ngrp; grp++) {
     221       52334 :                 pergroup[grp].npartials = 0;
     222       52334 :                 pergroup[grp].valseen = false;
     223       52334 :                 pergroup[grp].maxpartials = 2;
     224       52334 :                 pergroup[grp].infs = 0;
     225       52334 :                 pergroup[grp].partials = GDKmalloc(pergroup[grp].maxpartials * sizeof(double));
     226       52334 :                 if (pergroup[grp].partials == NULL) {
     227           0 :                         while (grp > 0)
     228           0 :                                 GDKfree(pergroup[--grp].partials);
     229           0 :                         GDKfree(pergroup);
     230           0 :                         return BUN_NONE;
     231             :                 }
     232             :         }
     233      149115 :         while (ncand > 0) {
     234      147814 :                 ncand--;
     235      147814 :                 listi = canditer_next(ci) - seqb;
     236      147795 :                 grp = gids ? gids[listi] : 0;
     237      147795 :                 if (grp < min || grp > max)
     238           0 :                         continue;
     239      147795 :                 if (pergroup[grp].partials == NULL)
     240           0 :                         continue;
     241      147795 :                 if (tp1 == TYPE_flt && !is_flt_nil(((const flt *) values)[listi]))
     242         426 :                         x = ((const flt *) values)[listi];
     243      147369 :                 else if (tp1 == TYPE_dbl && !is_dbl_nil(((const dbl *) values)[listi]))
     244             :                         x = ((const dbl *) values)[listi];
     245             :                 else {
     246             :                         /* it's a nil */
     247         373 :                         if (!skip_nils) {
     248           0 :                                 if (tp2 == TYPE_flt)
     249           0 :                                         ((flt *) results)[grp] = flt_nil;
     250             :                                 else
     251           0 :                                         ((dbl *) results)[grp] = dbl_nil;
     252           0 :                                 GDKfree(pergroup[grp].partials);
     253           0 :                                 pergroup[grp].partials = NULL;
     254           0 :                                 if (++nils == ngrp)
     255             :                                         break;
     256             :                         }
     257         373 :                         continue;
     258             :                 }
     259      147422 :                 pergroup[grp].valseen = true;
     260             : #ifdef INFINITES_ALLOWED
     261             :                 if (isinf(x)) {
     262             :                         pergroup[grp].infs += x;
     263             :                         continue;
     264             :                 }
     265             : #endif
     266      147422 :                 i = 0;
     267      401775 :                 for (parti = 0; parti < pergroup[grp].npartials; parti++) {
     268      254334 :                         y = pergroup[grp].partials[parti];
     269      254334 :                         if (fabs(x) < fabs(y))
     270       79114 :                                 exchange(&x, &y);
     271      254334 :                         twosum(&hi, &lo, x, y);
     272      253920 :                         if (isinf(hi)) {
     273          58 :                                 int sign = hi > 0 ? 1 : -1;
     274          58 :                                 hi = x - twopow * sign;
     275          58 :                                 x = hi - twopow * sign;
     276          58 :                                 pergroup[grp].infs += sign;
     277          58 :                                 if (fabs(x) < fabs(y))
     278           8 :                                         exchange(&x, &y);
     279          58 :                                 twosum(&hi, &lo, x, y);
     280             :                         }
     281      254353 :                         if (lo != 0)
     282      179801 :                                 pergroup[grp].partials[i++] = lo;
     283      254353 :                         x = hi;
     284             :                 }
     285      147441 :                 if (x != 0) {
     286      140048 :                         if (i == pergroup[grp].maxpartials) {
     287        6179 :                                 double *temp;
     288        6179 :                                 pergroup[grp].maxpartials += pergroup[grp].maxpartials;
     289        6179 :                                 temp = GDKrealloc(pergroup[grp].partials, pergroup[grp].maxpartials * sizeof(double));
     290        6179 :                                 if (temp == NULL)
     291           0 :                                         goto bailout;
     292        6179 :                                 pergroup[grp].partials = temp;
     293             :                         }
     294      140048 :                         pergroup[grp].partials[i++] = x;
     295             :                 }
     296      147441 :                 pergroup[grp].npartials = i;
     297             :         }
     298             : 
     299       54127 :         for (grp = 0; grp < ngrp; grp++) {
     300       52854 :                 if (pergroup[grp].partials == NULL)
     301           0 :                         continue;
     302       52854 :                 if (!pergroup[grp].valseen) {
     303         104 :                         if (tp2 == TYPE_flt)
     304          11 :                                 ((flt *) results)[grp] = nil_if_empty ? flt_nil : 0;
     305             :                         else
     306          93 :                                 ((dbl *) results)[grp] = nil_if_empty ? dbl_nil : 0;
     307         104 :                         nils += nil_if_empty;
     308         104 :                         GDKfree(pergroup[grp].partials);
     309          74 :                         pergroup[grp].partials = NULL;
     310          74 :                         continue;
     311             :                 }
     312             : #ifdef INFINITES_ALLOWED
     313             :                 if (isinf(pergroup[grp].infs) || isnan(pergroup[grp].infs)) {
     314             :                         if (abort_on_error) {
     315             :                                 goto overflow;
     316             :                         }
     317             :                         if (tp2 == TYPE_flt)
     318             :                                 ((flt *) results)[grp] = flt_nil;
     319             :                         else
     320             :                                 ((dbl *) results)[grp] = dbl_nil;
     321             :                         nils++;
     322             :                         GDKfree(pergroup[grp].partials);
     323             :                         pergroup[grp].partials = NULL;
     324             :                         continue;
     325             :                 }
     326             : #endif
     327             : 
     328       52750 :                 if ((pergroup[grp].infs == 1 || pergroup[grp].infs == -1) &&
     329          42 :                     pergroup[grp].npartials > 0 &&
     330          38 :                     !samesign(pergroup[grp].infs, pergroup[grp].partials[pergroup[grp].npartials - 1])) {
     331          34 :                         twosum(&hi, &lo, pergroup[grp].infs * twopow, pergroup[grp].partials[pergroup[grp].npartials - 1] / 2);
     332          34 :                         if (isinf(2 * hi)) {
     333          22 :                                 y = 2 * lo;
     334          22 :                                 x = hi + y;
     335          22 :                                 x -= hi;
     336          22 :                                 if (x == y &&
     337          18 :                                     pergroup[grp].npartials > 1 &&
     338          12 :                                     samesign(lo, pergroup[grp].partials[pergroup[grp].npartials - 2])) {
     339           6 :                                         GDKfree(pergroup[grp].partials);
     340           6 :                                         pergroup[grp].partials = NULL;
     341           6 :                                         x = 2 * (hi + y);
     342           6 :                                         if (tp2 == TYPE_flt) {
     343           0 :                                                 f = (flt) x;
     344           0 :                                                 if (isinf(f) ||
     345           0 :                                                     isnan(f) ||
     346           0 :                                                     is_flt_nil(f)) {
     347           0 :                                                         if (abort_on_error)
     348           0 :                                                                 goto overflow;
     349           0 :                                                         ((flt *) results)[grp] = flt_nil;
     350           0 :                                                         nils++;
     351             :                                                 } else
     352           0 :                                                         ((flt *) results)[grp] = f;
     353           6 :                                         } else if (is_dbl_nil(x)) {
     354           0 :                                                 if (abort_on_error)
     355           0 :                                                         goto overflow;
     356           0 :                                                 ((dbl *) results)[grp] = dbl_nil;
     357           0 :                                                 nils++;
     358             :                                         } else
     359           6 :                                                 ((dbl *) results)[grp] = x;
     360           6 :                                         continue;
     361             :                                 }
     362             :                         } else {
     363          12 :                                 if (lo) {
     364           2 :                                         if (pergroup[grp].npartials == pergroup[grp].maxpartials) {
     365           0 :                                                 double *temp;
     366             :                                                 /* we need space for one more */
     367           0 :                                                 pergroup[grp].maxpartials++;
     368           0 :                                                 temp = GDKrealloc(pergroup[grp].partials, pergroup[grp].maxpartials * sizeof(double));
     369           0 :                                                 if (temp == NULL)
     370           0 :                                                         goto bailout;
     371           0 :                                                 pergroup[grp].partials = temp;
     372             :                                         }
     373           2 :                                         pergroup[grp].partials[pergroup[grp].npartials - 1] = 2 * lo;
     374           2 :                                         pergroup[grp].partials[pergroup[grp].npartials++] = 2 * hi;
     375             :                                 } else {
     376          10 :                                         pergroup[grp].partials[pergroup[grp].npartials - 1] = 2 * hi;
     377             :                                 }
     378          12 :                                 pergroup[grp].infs = 0;
     379             :                         }
     380             :                 }
     381             : 
     382       52744 :                 if (pergroup[grp].infs != 0)
     383          28 :                         goto overflow;
     384             : 
     385       52716 :                 if (pergroup[grp].npartials == 0) {
     386        6284 :                         GDKfree(pergroup[grp].partials);
     387        6284 :                         pergroup[grp].partials = NULL;
     388        6284 :                         if (tp2 == TYPE_flt)
     389           0 :                                 ((flt *) results)[grp] = 0;
     390             :                         else
     391        6284 :                                 ((dbl *) results)[grp] = 0;
     392        6284 :                         continue;
     393             :                 }
     394             : 
     395             :                 /* accumulate into hi */
     396       46432 :                 hi = pergroup[grp].partials[--pergroup[grp].npartials];
     397       46514 :                 while (pergroup[grp].npartials > 0) {
     398       19946 :                         twosum(&hi, &lo, hi, pergroup[grp].partials[--pergroup[grp].npartials]);
     399       19817 :                         if (lo) {
     400       19735 :                                 pergroup[grp].partials[pergroup[grp].npartials++] = lo;
     401       19735 :                                 break;
     402             :                         }
     403             :                 }
     404             : 
     405       46303 :                 if (pergroup[grp].npartials >= 2 &&
     406        5247 :                     samesign(pergroup[grp].partials[pergroup[grp].npartials - 1], pergroup[grp].partials[pergroup[grp].npartials - 2]) &&
     407        2495 :                     hi + 2 * pergroup[grp].partials[pergroup[grp].npartials - 1] - hi == 2 * pergroup[grp].partials[pergroup[grp].npartials - 1]) {
     408        1299 :                         hi += 2 * pergroup[grp].partials[pergroup[grp].npartials - 1];
     409        1299 :                         pergroup[grp].partials[pergroup[grp].npartials - 1] = -pergroup[grp].partials[pergroup[grp].npartials - 1];
     410             :                 }
     411             : 
     412       46303 :                 GDKfree(pergroup[grp].partials);
     413       46462 :                 pergroup[grp].partials = NULL;
     414       46462 :                 if (tp2 == TYPE_flt) {
     415         115 :                         f = (flt) hi;
     416         115 :                         if (isinf(f) || isnan(f) || is_flt_nil(f)) {
     417           0 :                                 if (abort_on_error)
     418           0 :                                         goto overflow;
     419           0 :                                 ((flt *) results)[grp] = flt_nil;
     420           0 :                                 nils++;
     421             :                         } else
     422         115 :                                 ((flt *) results)[grp] = f;
     423       46347 :                 } else if (is_dbl_nil(hi)) {
     424           0 :                         if (abort_on_error)
     425           0 :                                 goto overflow;
     426           0 :                         ((dbl *) results)[grp] = dbl_nil;
     427           0 :                         nils++;
     428             :                 } else
     429       46347 :                         ((dbl *) results)[grp] = hi;
     430             :         }
     431        1273 :         GDKfree(pergroup);
     432        1273 :         return nils;
     433             : 
     434          28 :   overflow:
     435          28 :         GDKerror("22003!overflow in calculation.\n");
     436             :   bailout:
     437          56 :         for (grp = 0; grp < ngrp; grp++)
     438          28 :                 GDKfree(pergroup[grp].partials);
     439          28 :         GDKfree(pergroup);
     440          28 :         return BUN_NONE;
     441             : }
     442             : 
     443             : #define AGGR_SUM(TYPE1, TYPE2)                                          \
     444             :         do {                                                            \
     445             :                 TYPE1 x;                                                \
     446             :                 const TYPE1 *restrict vals = (const TYPE1 *) values;    \
     447             :                 if (ngrp == 1 && ci->tpe == cand_dense) {            \
     448             :                         /* single group, no candidate list */           \
     449             :                         TYPE2 sum;                                      \
     450             :                         *algo = "no candidates, no groups";           \
     451             :                         sum = 0;                                        \
     452             :                         if (nonil) {                                    \
     453             :                                 *seen = ncand > 0;                   \
     454             :                                 for (i = 0; i < ncand && nils == 0; i++) { \
     455             :                                         x = vals[ci->seq + i - seqb];        \
     456             :                                         ADD_WITH_CHECK(x, sum,          \
     457             :                                                        TYPE2, sum,      \
     458             :                                                        GDK_##TYPE2##_max, \
     459             :                                                        goto overflow);  \
     460             :                                 }                                       \
     461             :                         } else {                                        \
     462             :                                 bool seenval = false;                   \
     463             :                                 for (i = 0; i < ncand && nils == 0; i++) { \
     464             :                                         x = vals[ci->seq + i - seqb];        \
     465             :                                         if (is_##TYPE1##_nil(x)) {      \
     466             :                                                 if (!skip_nils) {       \
     467             :                                                         sum = TYPE2##_nil; \
     468             :                                                         nils = 1;       \
     469             :                                                 }                       \
     470             :                                         } else {                        \
     471             :                                                 ADD_WITH_CHECK(x, sum,  \
     472             :                                                                TYPE2, sum, \
     473             :                                                                GDK_##TYPE2##_max, \
     474             :                                                                goto overflow); \
     475             :                                                 seenval = true;         \
     476             :                                         }                               \
     477             :                                 }                                       \
     478             :                                 *seen = seenval;                        \
     479             :                         }                                               \
     480             :                         if (*seen)                                      \
     481             :                                 *sums = sum;                            \
     482             :                 } else if (ngrp == 1) {                                 \
     483             :                         /* single group, with candidate list */         \
     484             :                         TYPE2 sum;                                      \
     485             :                         bool seenval = false;                           \
     486             :                         *algo = "with candidates, no groups";         \
     487             :                         sum = 0;                                        \
     488             :                         for (i = 0; i < ncand && nils == 0; i++) {   \
     489             :                                 x = vals[canditer_next(ci) - seqb];     \
     490             :                                 if (is_##TYPE1##_nil(x)) {              \
     491             :                                         if (!skip_nils) {               \
     492             :                                                 sum = TYPE2##_nil;      \
     493             :                                                 nils = 1;               \
     494             :                                         }                               \
     495             :                                 } else {                                \
     496             :                                         ADD_WITH_CHECK(x, sum,          \
     497             :                                                        TYPE2, sum,      \
     498             :                                                        GDK_##TYPE2##_max, \
     499             :                                                        goto overflow);  \
     500             :                                         seenval = true;                 \
     501             :                                 }                                       \
     502             :                         }                                               \
     503             :                         if (seenval)                                    \
     504             :                                 *sums = sum;                            \
     505             :                 } else if (ci->tpe == cand_dense) {                  \
     506             :                         /* multiple groups, no candidate list */        \
     507             :                         *algo = "no candidates, with groups";         \
     508             :                         for (i = 0; i < ncand; i++) {                        \
     509             :                                 if (gids == NULL ||                     \
     510             :                                     (gids[i] >= min && gids[i] <= max)) { \
     511             :                                         gid = gids ? gids[i] - min : (oid) i; \
     512             :                                         x = vals[ci->seq + i - seqb];        \
     513             :                                         if (is_##TYPE1##_nil(x)) {      \
     514             :                                                 if (!skip_nils) {       \
     515             :                                                         sums[gid] = TYPE2##_nil; \
     516             :                                                         nils++;         \
     517             :                                                 }                       \
     518             :                                         } else {                        \
     519             :                                                 if (nil_if_empty &&     \
     520             :                                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     521             :                                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
     522             :                                                         sums[gid] = 0;  \
     523             :                                                 }                       \
     524             :                                                 if (!is_##TYPE2##_nil(sums[gid])) { \
     525             :                                                         ADD_WITH_CHECK( \
     526             :                                                                 x,      \
     527             :                                                                 sums[gid], \
     528             :                                                                 TYPE2,  \
     529             :                                                                 sums[gid], \
     530             :                                                                 GDK_##TYPE2##_max, \
     531             :                                                                 goto overflow); \
     532             :                                                 }                       \
     533             :                                         }                               \
     534             :                                 }                                       \
     535             :                         }                                               \
     536             :                 } else {                                                \
     537             :                         /* multiple groups, with candidate list */      \
     538             :                         *algo = "with candidates, with groups";               \
     539             :                         while (ncand > 0) {                          \
     540             :                                 ncand--;                                \
     541             :                                 i = canditer_next(ci) - seqb;           \
     542             :                                 if (gids == NULL ||                     \
     543             :                                     (gids[i] >= min && gids[i] <= max)) { \
     544             :                                         gid = gids ? gids[i] - min : (oid) i; \
     545             :                                         x = vals[i];                    \
     546             :                                         if (is_##TYPE1##_nil(x)) {      \
     547             :                                                 if (!skip_nils) {       \
     548             :                                                         sums[gid] = TYPE2##_nil; \
     549             :                                                         nils++;         \
     550             :                                                 }                       \
     551             :                                         } else {                        \
     552             :                                                 if (nil_if_empty &&     \
     553             :                                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     554             :                                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
     555             :                                                         sums[gid] = 0;  \
     556             :                                                 }                       \
     557             :                                                 if (!is_##TYPE2##_nil(sums[gid])) { \
     558             :                                                         ADD_WITH_CHECK( \
     559             :                                                                 x,      \
     560             :                                                                 sums[gid], \
     561             :                                                                 TYPE2,  \
     562             :                                                                 sums[gid], \
     563             :                                                                 GDK_##TYPE2##_max, \
     564             :                                                                 goto overflow); \
     565             :                                                 }                       \
     566             :                                         }                               \
     567             :                                 }                                       \
     568             :                         }                                               \
     569             :                 }                                                       \
     570             :         } while (0)
     571             : 
     572             : #define AGGR_SUM_NOOVL(TYPE1, TYPE2)                                    \
     573             :         do {                                                            \
     574             :                 TYPE1 x;                                                \
     575             :                 const TYPE1 *restrict vals = (const TYPE1 *) values;    \
     576             :                 if (ngrp == 1 && ci->tpe == cand_dense) {            \
     577             :                         /* single group, no candidate list */           \
     578             :                         TYPE2 sum;                                      \
     579             :                         sum = 0;                                        \
     580             :                         if (nonil) {                                    \
     581             :                                 *algo = "no candidates, no groups, no nils, no overflow"; \
     582             :                                 *seen = ncand > 0;                   \
     583             :                                 for (i = 0; i < ncand && nils == 0; i++) { \
     584             :                                         sum += vals[ci->seq + i - seqb]; \
     585             :                                 }                                       \
     586             :                         } else {                                        \
     587             :                                 bool seenval = false;                   \
     588             :                                 *algo = "no candidates, no groups, no overflow"; \
     589             :                                 for (i = 0; i < ncand && nils == 0; i++) { \
     590             :                                         x = vals[ci->seq + i - seqb];        \
     591             :                                         if (is_##TYPE1##_nil(x)) {      \
     592             :                                                 if (!skip_nils) {       \
     593             :                                                         sum = TYPE2##_nil; \
     594             :                                                         nils = 1;       \
     595             :                                                 }                       \
     596             :                                         } else {                        \
     597             :                                                 sum += x;               \
     598             :                                                 seenval = true;         \
     599             :                                         }                               \
     600             :                                 }                                       \
     601             :                                 *seen = seenval;                        \
     602             :                         }                                               \
     603             :                         if (*seen)                                      \
     604             :                                 *sums = sum;                            \
     605             :                 } else if (ngrp == 1) {                                 \
     606             :                         /* single group, with candidate list */         \
     607             :                         TYPE2 sum;                                      \
     608             :                         bool seenval = false;                           \
     609             :                         *algo = "with candidates, no groups, no overflow"; \
     610             :                         sum = 0;                                        \
     611             :                         for (i = 0; i < ncand && nils == 0; i++) {   \
     612             :                                 x = vals[canditer_next(ci) - seqb];     \
     613             :                                 if (is_##TYPE1##_nil(x)) {              \
     614             :                                         if (!skip_nils) {               \
     615             :                                                 sum = TYPE2##_nil;      \
     616             :                                                 nils = 1;               \
     617             :                                         }                               \
     618             :                                 } else {                                \
     619             :                                         sum += x;                       \
     620             :                                         seenval = true;                 \
     621             :                                 }                                       \
     622             :                         }                                               \
     623             :                         if (seenval)                                    \
     624             :                                 *sums = sum;                            \
     625             :                 } else if (ci->tpe == cand_dense) {                  \
     626             :                         /* multiple groups, no candidate list */        \
     627             :                         if (nonil) {                                    \
     628             :                                 *algo = "no candidates, with groups, no nils, no overflow"; \
     629             :                                 for (i = 0; i < ncand; i++) {                \
     630             :                                         if (gids == NULL ||             \
     631             :                                             (gids[i] >= min && gids[i] <= max)) { \
     632             :                                                 gid = gids ? gids[i] - min : (oid) i; \
     633             :                                                 x = vals[ci->seq + i - seqb]; \
     634             :                                                 if (nil_if_empty &&     \
     635             :                                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     636             :                                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
     637             :                                                         sums[gid] = 0;  \
     638             :                                                 }                       \
     639             :                                                 sums[gid] += x;         \
     640             :                                         }                               \
     641             :                                 }                                       \
     642             :                         } else {                                        \
     643             :                                 *algo = "no candidates, with groups, no overflow"; \
     644             :                                 for (i = 0; i < ncand; i++) {                \
     645             :                                         if (gids == NULL ||             \
     646             :                                             (gids[i] >= min && gids[i] <= max)) { \
     647             :                                                 gid = gids ? gids[i] - min : (oid) i; \
     648             :                                                 x = vals[ci->seq + i - seqb]; \
     649             :                                                 if (is_##TYPE1##_nil(x)) { \
     650             :                                                         if (!skip_nils) { \
     651             :                                                                 sums[gid] = TYPE2##_nil; \
     652             :                                                                 nils++; \
     653             :                                                         }               \
     654             :                                                 } else {                \
     655             :                                                         if (nil_if_empty && \
     656             :                                                             !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     657             :                                                                 seen[gid >> 5] |= 1U << (gid & 0x1F); \
     658             :                                                                 sums[gid] = 0; \
     659             :                                                         }               \
     660             :                                                         if (!is_##TYPE2##_nil(sums[gid])) { \
     661             :                                                                 sums[gid] += x; \
     662             :                                                         }               \
     663             :                                                 }                       \
     664             :                                         }                               \
     665             :                                 }                                       \
     666             :                         }                                               \
     667             :                 } else {                                                \
     668             :                         /* multiple groups, with candidate list */      \
     669             :                         *algo = "with candidates, with groups, no overflow"; \
     670             :                         while (ncand > 0) {                          \
     671             :                                 ncand--;                                \
     672             :                                 i = canditer_next(ci) - seqb;           \
     673             :                                 if (gids == NULL ||                     \
     674             :                                     (gids[i] >= min && gids[i] <= max)) { \
     675             :                                         gid = gids ? gids[i] - min : (oid) i; \
     676             :                                         x = vals[i];                    \
     677             :                                         if (is_##TYPE1##_nil(x)) {      \
     678             :                                                 if (!skip_nils) {       \
     679             :                                                         sums[gid] = TYPE2##_nil; \
     680             :                                                         nils++;         \
     681             :                                                 }                       \
     682             :                                         } else {                        \
     683             :                                                 if (nil_if_empty &&     \
     684             :                                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     685             :                                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
     686             :                                                         sums[gid] = 0;  \
     687             :                                                 }                       \
     688             :                                                 if (!is_##TYPE2##_nil(sums[gid])) { \
     689             :                                                         sums[gid] += x; \
     690             :                                                 }                       \
     691             :                                         }                               \
     692             :                                 }                                       \
     693             :                         }                                               \
     694             :                 }                                                       \
     695             :         } while (0)
     696             : 
     697             : static BUN
     698        9128 : dosum(const void *restrict values, bool nonil, oid seqb,
     699             :       struct canditer *restrict ci, BUN ncand,
     700             :       void *restrict results, BUN ngrp, int tp1, int tp2,
     701             :       const oid *restrict gids,
     702             :       oid min, oid max, bool skip_nils, bool abort_on_error,
     703             :       bool nil_if_empty, const char *func, const char **algo)
     704             : {
     705        9128 :         BUN nils = 0;
     706        9128 :         BUN i;
     707        9128 :         oid gid;
     708        9128 :         unsigned int *restrict seen = NULL; /* bitmask for groups that we've seen */
     709             : 
     710        9128 :         switch (tp2) {
     711          12 :         case TYPE_flt:
     712          12 :                 if (tp1 != TYPE_flt)
     713           0 :                         goto unsupported;
     714             :                 /* fall through */
     715             :         case TYPE_dbl:
     716        1125 :                 if (tp1 != TYPE_flt && tp1 != TYPE_dbl)
     717           0 :                         goto unsupported;
     718        1125 :                 *algo = "floating sum";
     719        1125 :                 return dofsum(values, seqb, ci, ncand, results, ngrp, tp1, tp2,
     720             :                               gids, min, max, skip_nils, abort_on_error,
     721             :                               nil_if_empty);
     722             :         }
     723             : 
     724             :         /* allocate bitmap for seen group ids */
     725        8003 :         seen = GDKzalloc(((ngrp + 31) / 32) * sizeof(int));
     726        8003 :         if (seen == NULL) {
     727             :                 return BUN_NONE;
     728             :         }
     729             : 
     730        8003 :         switch (tp2) {
     731           0 :         case TYPE_bte: {
     732           0 :                 bte *restrict sums = (bte *) results;
     733           0 :                 switch (tp1) {
     734           0 :                 case TYPE_bte:
     735           0 :                         AGGR_SUM(bte, bte);
     736             :                         break;
     737           0 :                 default:
     738           0 :                         goto unsupported;
     739             :                 }
     740             :                 break;
     741             :         }
     742           0 :         case TYPE_sht: {
     743           0 :                 sht *restrict sums = (sht *) results;
     744           0 :                 switch (tp1) {
     745           0 :                 case TYPE_bte:
     746           0 :                         if (ncand < ((BUN) 1 << ((sizeof(sht) - sizeof(bte)) << 3)))
     747           0 :                                 AGGR_SUM_NOOVL(bte, sht);
     748             :                         else
     749           0 :                                 AGGR_SUM(bte, sht);
     750             :                         break;
     751           0 :                 case TYPE_sht:
     752           0 :                         AGGR_SUM(sht, sht);
     753             :                         break;
     754           0 :                 default:
     755           0 :                         goto unsupported;
     756             :                 }
     757             :                 break;
     758             :         }
     759          10 :         case TYPE_int: {
     760          10 :                 int *restrict sums = (int *) results;
     761          10 :                 switch (tp1) {
     762           0 :                 case TYPE_bte:
     763           0 :                         if (ncand < ((BUN) 1 << ((sizeof(int) - sizeof(bte)) << 3)))
     764           0 :                                 AGGR_SUM_NOOVL(bte, int);
     765             :                         else
     766           0 :                                 AGGR_SUM(bte, int);
     767             :                         break;
     768           0 :                 case TYPE_sht:
     769           0 :                         if (ncand < ((BUN) 1 << ((sizeof(int) - sizeof(sht)) << 3)))
     770           0 :                                 AGGR_SUM_NOOVL(sht, int);
     771             :                         else
     772           0 :                                 AGGR_SUM(sht, int);
     773             :                         break;
     774          10 :                 case TYPE_int:
     775          26 :                         AGGR_SUM(int, int);
     776             :                         break;
     777           0 :                 default:
     778           0 :                         goto unsupported;
     779             :                 }
     780             :                 break;
     781             :         }
     782        6237 :         case TYPE_lng: {
     783        6237 :                 lng *restrict sums = (lng *) results;
     784        6237 :                 switch (tp1) {
     785             : #if SIZEOF_BUN == 4
     786             :                 case TYPE_bte:
     787             :                         AGGR_SUM_NOOVL(bte, lng);
     788             :                         break;
     789             :                 case TYPE_sht:
     790             :                         AGGR_SUM_NOOVL(sht, lng);
     791             :                         break;
     792             :                 case TYPE_int:
     793             :                         AGGR_SUM_NOOVL(int, lng);
     794             :                         break;
     795             : #else
     796           0 :                 case TYPE_bte:
     797           0 :                         if (ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(bte)) << 3)))
     798           0 :                                 AGGR_SUM_NOOVL(bte, lng);
     799             :                         else
     800           0 :                                 AGGR_SUM(bte, lng);
     801             :                         break;
     802           0 :                 case TYPE_sht:
     803           0 :                         if (ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(sht)) << 3)))
     804           0 :                                 AGGR_SUM_NOOVL(sht, lng);
     805             :                         else
     806           0 :                                 AGGR_SUM(sht, lng);
     807             :                         break;
     808           0 :                 case TYPE_int:
     809           0 :                         if (ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(int)) << 3)))
     810           0 :                                 AGGR_SUM_NOOVL(int, lng);
     811             :                         else
     812           0 :                                 AGGR_SUM(int, lng);
     813             :                         break;
     814             : #endif
     815        6237 :                 case TYPE_lng:
     816     7259960 :                         AGGR_SUM(lng, lng);
     817             :                         break;
     818           0 :                 default:
     819           0 :                         goto unsupported;
     820             :                 }
     821             :                 break;
     822             :         }
     823             : #ifdef HAVE_HGE
     824        1756 :         case TYPE_hge: {
     825        1756 :                 hge *sums = (hge *) results;
     826        1756 :                 switch (tp1) {
     827          94 :                 case TYPE_bte:
     828     2972960 :                         AGGR_SUM_NOOVL(bte, hge);
     829             :                         break;
     830           4 :                 case TYPE_sht:
     831           9 :                         AGGR_SUM_NOOVL(sht, hge);
     832             :                         break;
     833         948 :                 case TYPE_int:
     834    53398700 :                         AGGR_SUM_NOOVL(int, hge);
     835             :                         break;
     836          46 :                 case TYPE_lng:
     837     5163900 :                         AGGR_SUM_NOOVL(lng, hge);
     838             :                         break;
     839         664 :                 case TYPE_hge:
     840    14781100 :                         AGGR_SUM(hge, hge);
     841             :                         break;
     842           0 :                 default:
     843           0 :                         goto unsupported;
     844             :                 }
     845             :                 break;
     846             :         }
     847             : #endif
     848           0 :         default:
     849           0 :                 goto unsupported;
     850             :         }
     851             : 
     852        8003 :         if (nils == 0 && nil_if_empty) {
     853             :                 /* figure out whether there were any empty groups
     854             :                  * (that result in a nil value) */
     855        8002 :                 if (ngrp & 0x1F) {
     856             :                         /* fill last slot */
     857        7983 :                         seen[ngrp >> 5] |= ~0U << (ngrp & 0x1F);
     858             :                 }
     859      196630 :                 for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
     860      188954 :                         if (seen[i] != ~0U) {
     861             :                                 nils = 1;
     862             :                                 break;
     863             :                         }
     864             :                 }
     865             :         }
     866        8003 :         GDKfree(seen);
     867             : 
     868        8003 :         return nils;
     869             : 
     870           0 :   unsupported:
     871           0 :         GDKfree(seen);
     872           0 :         GDKerror("%s: type combination (sum(%s)->%s) not supported.\n",
     873             :                  func, ATOMname(tp1), ATOMname(tp2));
     874           0 :         return BUN_NONE;
     875             : 
     876           0 :   overflow:
     877           0 :         GDKfree(seen);
     878           0 :         GDKerror("22003!overflow in calculation.\n");
     879           0 :         return BUN_NONE;
     880             : }
     881             : 
     882             : /* calculate group sums with optional candidates list */
     883             : BAT *
     884        4779 : BATgroupsum(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
     885             : {
     886        4779 :         const oid *restrict gids;
     887        4779 :         oid min, max;
     888        4779 :         BUN ngrp;
     889        4779 :         BUN nils;
     890        4779 :         BAT *bn;
     891        4779 :         struct canditer ci;
     892        4779 :         BUN ncand;
     893        4779 :         const char *err;
     894        4779 :         const char *algo = NULL;
     895        4779 :         lng t0 = 0;
     896             : 
     897        4779 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
     898             : 
     899        4779 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
     900           0 :                 GDKerror("%s\n", err);
     901           0 :                 return NULL;
     902             :         }
     903        4779 :         if (g == NULL) {
     904           0 :                 GDKerror("b and g must be aligned\n");
     905           0 :                 return NULL;
     906             :         }
     907             : 
     908        4779 :         if (BATcount(b) == 0 || ngrp == 0) {
     909             :                 /* trivial: no sums, so return bat aligned with g with
     910             :                  * nil in the tail */
     911        1691 :                 return BATconstant(ngrp == 0 ? 0 : min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
     912             :         }
     913             : 
     914        3088 :         if ((e == NULL ||
     915        3088 :              (BATcount(e) == BATcount(b) && e->hseqbase == b->hseqbase)) &&
     916         786 :             (BATtdense(g) || (g->tkey && g->tnonil))) {
     917             :                 /* trivial: singleton groups, so all results are equal
     918             :                  * to the inputs (but possibly a different type) */
     919         786 :                 return BATconvert(b, s, NULL, tp, abort_on_error);
     920             :         }
     921             : 
     922        2302 :         bn = BATconstant(min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
     923        2302 :         if (bn == NULL) {
     924             :                 return NULL;
     925             :         }
     926             : 
     927        2302 :         if (BATtdense(g))
     928             :                 gids = NULL;
     929             :         else
     930        2168 :                 gids = (const oid *) Tloc(g, 0);
     931             : 
     932        4604 :         nils = dosum(Tloc(b, 0), b->tnonil, b->hseqbase, &ci, ncand,
     933        2302 :                      Tloc(bn, 0), ngrp, b->ttype, tp, gids, min, max,
     934             :                      skip_nils, abort_on_error, true, __func__, &algo);
     935             : 
     936        2302 :         if (nils < BUN_NONE) {
     937        2302 :                 BATsetcount(bn, ngrp);
     938        2302 :                 bn->tkey = BATcount(bn) <= 1;
     939        2302 :                 bn->tsorted = BATcount(bn) <= 1;
     940        2302 :                 bn->trevsorted = BATcount(bn) <= 1;
     941        2302 :                 bn->tnil = nils != 0;
     942        2302 :                 bn->tnonil = nils == 0;
     943             :         } else {
     944           0 :                 BBPunfix(bn->batCacheid);
     945           0 :                 bn = NULL;
     946             :         }
     947             : 
     948        2302 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
     949             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
     950             :                   "; start " OIDFMT ", count " BUNFMT " (%s -- " LLFMT " usec)\n",
     951             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
     952             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
     953             :                   ci.seq, ncand, algo ? algo : "", GDKusec() - t0);
     954             :         return bn;
     955             : }
     956             : 
     957             : gdk_return
     958        7029 : BATsum(void *res, int tp, BAT *b, BAT *s, bool skip_nils, bool abort_on_error, bool nil_if_empty)
     959             : {
     960        7029 :         oid min, max;
     961        7029 :         BUN ngrp;
     962        7029 :         BUN nils;
     963        7029 :         struct canditer ci;
     964        7029 :         BUN ncand;
     965        7029 :         const char *err;
     966        7029 :         const char *algo = NULL;
     967        7029 :         lng t0 = 0;
     968             : 
     969        7029 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
     970             : 
     971        7029 :         if ((err = BATgroupaggrinit(b, NULL, NULL, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
     972           0 :                 GDKerror("%s\n", err);
     973           0 :                 return GDK_FAIL;
     974             :         }
     975        7029 :         switch (tp) {
     976           0 :         case TYPE_bte:
     977           0 :                 * (bte *) res = nil_if_empty ? bte_nil : 0;
     978           0 :                 break;
     979           0 :         case TYPE_sht:
     980           0 :                 * (sht *) res = nil_if_empty ? sht_nil : 0;
     981           0 :                 break;
     982          10 :         case TYPE_int:
     983          10 :                 * (int *) res = nil_if_empty ? int_nil : 0;
     984          10 :                 break;
     985        5236 :         case TYPE_lng:
     986        5236 :                 * (lng *) res = nil_if_empty ? lng_nil : 0;
     987        5236 :                 break;
     988             : #ifdef HAVE_HGE
     989         678 :         case TYPE_hge:
     990         678 :                 * (hge *) res = nil_if_empty ? hge_nil : 0;
     991         678 :                 break;
     992             : #endif
     993        1105 :         case TYPE_flt:
     994             :         case TYPE_dbl:
     995        1105 :                 switch (b->ttype) {
     996           1 :                 case TYPE_bte:
     997             :                 case TYPE_sht:
     998             :                 case TYPE_int:
     999             :                 case TYPE_lng:
    1000             : #ifdef HAVE_HGE
    1001             :                 case TYPE_hge:
    1002             : #endif
    1003             :                 {
    1004             :                         /* special case for summing integer types into
    1005             :                          * a floating point: We calculate the average
    1006             :                          * (which is done exactly), and multiply the
    1007             :                          * result by the count to get the sum.  Note
    1008             :                          * that if we just summed into a floating
    1009             :                          * point number, we could loose too much
    1010             :                          * accuracy, and if we summed into lng first,
    1011             :                          * we could get unnecessary overflow. */
    1012           1 :                         dbl avg;
    1013           1 :                         BUN cnt;
    1014             : 
    1015           1 :                         if (BATcalcavg(b, s, &avg, &cnt, 0) != GDK_SUCCEED)
    1016             :                                 return GDK_FAIL;
    1017           1 :                         if (cnt == 0) {
    1018           0 :                                 avg = nil_if_empty ? dbl_nil : 0;
    1019             :                         }
    1020           1 :                         if (cnt < BATcount(b) && !skip_nils) {
    1021             :                                 /* there were nils */
    1022           0 :                                 avg = dbl_nil;
    1023             :                         }
    1024           1 :                         if (tp == TYPE_flt) {
    1025           0 :                                 if (is_dbl_nil(avg))
    1026           0 :                                         *(flt *) res = flt_nil;
    1027           0 :                                 else if (cnt > 0 &&
    1028           0 :                                          GDK_flt_max / cnt < fabs(avg)) {
    1029           0 :                                         if (abort_on_error) {
    1030           0 :                                                 GDKerror("22003!overflow in calculation.\n");
    1031           0 :                                                 return GDK_FAIL;
    1032             :                                         }
    1033           0 :                                         *(flt *) res = flt_nil;
    1034             :                                 } else {
    1035           0 :                                         *(flt *) res = (flt) avg * cnt;
    1036             :                                 }
    1037             :                         } else {
    1038           1 :                                 if (is_dbl_nil(avg)) {
    1039           0 :                                         *(dbl *) res = dbl_nil;
    1040           1 :                                 } else if (cnt > 0 &&
    1041           1 :                                            GDK_dbl_max / cnt < fabs(avg)) {
    1042           0 :                                         if (abort_on_error) {
    1043           0 :                                                 GDKerror("22003!overflow in calculation.\n");
    1044           0 :                                                 return GDK_FAIL;
    1045             :                                         }
    1046           0 :                                         *(dbl *) res = dbl_nil;
    1047             :                                 } else {
    1048           1 :                                         *(dbl *) res = avg * cnt;
    1049             :                                 }
    1050             :                         }
    1051             :                         return GDK_SUCCEED;
    1052             :                 }
    1053             :                 default:
    1054        1104 :                         break;
    1055             :                 }
    1056        1104 :                 if (b->ttype == TYPE_dbl)
    1057        1092 :                         * (dbl *) res = nil_if_empty ? dbl_nil : 0;
    1058             :                 else
    1059          12 :                         * (flt *) res = nil_if_empty ? flt_nil : 0;
    1060             :                 break;
    1061           0 :         default:
    1062           0 :                 GDKerror("type combination (sum(%s)->%s) not supported.\n",
    1063             :                          ATOMname(b->ttype), ATOMname(tp));
    1064           0 :                 return GDK_FAIL;
    1065             :         }
    1066        7028 :         if (BATcount(b) == 0)
    1067             :                 return GDK_SUCCEED;
    1068       13651 :         nils = dosum(Tloc(b, 0), b->tnonil, b->hseqbase, &ci, ncand,
    1069        6825 :                      res, true, b->ttype, tp, &min, min, max,
    1070             :                      skip_nils, abort_on_error, nil_if_empty, __func__, &algo);
    1071        6826 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",s=" ALGOOPTBATFMT "; "
    1072             :                   "start " OIDFMT ", count " BUNFMT " (%s -- " LLFMT " usec)\n",
    1073             :                   ALGOBATPAR(b), ALGOOPTBATPAR(s),
    1074             :                   ci.seq, ncand, algo ? algo : "", GDKusec() - t0);
    1075        6826 :         return nils < BUN_NONE ? GDK_SUCCEED : GDK_FAIL;
    1076             : }
    1077             : 
    1078             : /* ---------------------------------------------------------------------- */
    1079             : /* product */
    1080             : 
    1081             : #define AGGR_PROD(TYPE1, TYPE2, TYPE3)                                  \
    1082             :         do {                                                            \
    1083             :                 const TYPE1 *restrict vals = (const TYPE1 *) values;    \
    1084             :                 gid = 0;        /* doesn't change if gidincr == false */ \
    1085             :                 while (ncand > 0) {                                  \
    1086             :                         ncand--;                                        \
    1087             :                         i = canditer_next(ci) - seqb;                   \
    1088             :                         if (gids == NULL || !gidincr ||                 \
    1089             :                             (gids[i] >= min && gids[i] <= max)) { \
    1090             :                                 if (gidincr) {                          \
    1091             :                                         if (gids)                       \
    1092             :                                                 gid = gids[i] - min;    \
    1093             :                                         else                            \
    1094             :                                                 gid = (oid) i;          \
    1095             :                                 }                                       \
    1096             :                                 if (is_##TYPE1##_nil(vals[i])) {        \
    1097             :                                         if (!skip_nils) {               \
    1098             :                                                 prods[gid] = TYPE2##_nil; \
    1099             :                                                 nils++;                 \
    1100             :                                         }                               \
    1101             :                                 } else {                                \
    1102             :                                         if (nil_if_empty &&             \
    1103             :                                             !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
    1104             :                                                 seen[gid >> 5] |= 1U << (gid & 0x1F); \
    1105             :                                                 prods[gid] = 1;         \
    1106             :                                         }                               \
    1107             :                                         if (!is_##TYPE2##_nil(prods[gid])) { \
    1108             :                                                 MUL4_WITH_CHECK(        \
    1109             :                                                         vals[i],        \
    1110             :                                                         prods[gid],     \
    1111             :                                                         TYPE2, prods[gid], \
    1112             :                                                         GDK_##TYPE2##_max, \
    1113             :                                                         TYPE3,          \
    1114             :                                                         goto overflow); \
    1115             :                                         }                               \
    1116             :                                 }                                       \
    1117             :                         }                                               \
    1118             :                 }                                                       \
    1119             :         } while (0)
    1120             : 
    1121             : #ifdef HAVE_HGE
    1122             : #define AGGR_PROD_HGE(TYPE)                                             \
    1123             :         do {                                                            \
    1124             :                 const TYPE *vals = (const TYPE *) values;               \
    1125             :                 gid = 0;        /* doesn't change if gidincr == false */ \
    1126             :                 while (ncand > 0) {                                  \
    1127             :                         ncand--;                                        \
    1128             :                         i = canditer_next(ci) - seqb;                   \
    1129             :                         if (gids == NULL || !gidincr ||                 \
    1130             :                             (gids[i] >= min && gids[i] <= max)) { \
    1131             :                                 if (gidincr) {                          \
    1132             :                                         if (gids)                       \
    1133             :                                                 gid = gids[i] - min;    \
    1134             :                                         else                            \
    1135             :                                                 gid = (oid) i;          \
    1136             :                                 }                                       \
    1137             :                                 if (nil_if_empty &&                     \
    1138             :                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
    1139             :                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
    1140             :                                         prods[gid] = 1;                 \
    1141             :                                 }                                       \
    1142             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    1143             :                                         if (!skip_nils) {               \
    1144             :                                                 prods[gid] = hge_nil;   \
    1145             :                                                 nils++;                 \
    1146             :                                         }                               \
    1147             :                                 } else if (!is_hge_nil(prods[gid])) {   \
    1148             :                                         HGEMUL_CHECK(vals[i],           \
    1149             :                                                      prods[gid],        \
    1150             :                                                      prods[gid],        \
    1151             :                                                      GDK_hge_max,       \
    1152             :                                                      goto overflow);    \
    1153             :                                 }                                       \
    1154             :                         }                                               \
    1155             :                 }                                                       \
    1156             :         } while (0)
    1157             : #else
    1158             : #define AGGR_PROD_LNG(TYPE)                                             \
    1159             :         do {                                                            \
    1160             :                 const TYPE *restrict vals = (const TYPE *) values;      \
    1161             :                 gid = 0;        /* doesn't change if gidincr == false */ \
    1162             :                 while (ncand > 0) {                                  \
    1163             :                         ncand--;                                        \
    1164             :                         i = canditer_next(ci) - seqb;                   \
    1165             :                         if (gids == NULL || !gidincr ||                 \
    1166             :                             (gids[i] >= min && gids[i] <= max)) { \
    1167             :                                 if (gidincr) {                          \
    1168             :                                         if (gids)                       \
    1169             :                                                 gid = gids[i] - min;    \
    1170             :                                         else                            \
    1171             :                                                 gid = (oid) i;          \
    1172             :                                 }                                       \
    1173             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    1174             :                                         if (!skip_nils) {               \
    1175             :                                                 prods[gid] = lng_nil;   \
    1176             :                                                 nils++;                 \
    1177             :                                         }                               \
    1178             :                                 } else {                                \
    1179             :                                         if (nil_if_empty &&             \
    1180             :                                             !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
    1181             :                                                 seen[gid >> 5] |= 1U << (gid & 0x1F); \
    1182             :                                                 prods[gid] = 1;         \
    1183             :                                         }                               \
    1184             :                                         if (!is_lng_nil(prods[gid])) {  \
    1185             :                                                 LNGMUL_CHECK(           \
    1186             :                                                         vals[i],        \
    1187             :                                                         prods[gid],     \
    1188             :                                                         prods[gid],     \
    1189             :                                                         GDK_lng_max,    \
    1190             :                                                         goto overflow); \
    1191             :                                         }                               \
    1192             :                                 }                                       \
    1193             :                         }                                               \
    1194             :                 }                                                       \
    1195             :         } while (0)
    1196             : #endif
    1197             : 
    1198             : #define AGGR_PROD_FLOAT(TYPE1, TYPE2)                                   \
    1199             :         do {                                                            \
    1200             :                 const TYPE1 *restrict vals = (const TYPE1 *) values;    \
    1201             :                 gid = 0;        /* doesn't change if gidincr == false */ \
    1202             :                 while (ncand > 0) {                                  \
    1203             :                         ncand--;                                        \
    1204             :                         i = canditer_next(ci) - seqb;                   \
    1205             :                         if (gids == NULL || !gidincr ||                 \
    1206             :                             (gids[i] >= min && gids[i] <= max)) { \
    1207             :                                 if (gidincr) {                          \
    1208             :                                         if (gids)                       \
    1209             :                                                 gid = gids[i] - min;    \
    1210             :                                         else                            \
    1211             :                                                 gid = (oid) i;          \
    1212             :                                 }                                       \
    1213             :                                 if (is_##TYPE1##_nil(vals[i])) {        \
    1214             :                                         if (!skip_nils) {               \
    1215             :                                                 prods[gid] = TYPE2##_nil; \
    1216             :                                                 nils++;                 \
    1217             :                                         }                               \
    1218             :                                 } else {                                \
    1219             :                                         if (nil_if_empty &&             \
    1220             :                                             !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
    1221             :                                                 seen[gid >> 5] |= 1U << (gid & 0x1F); \
    1222             :                                                 prods[gid] = 1;         \
    1223             :                                         }                               \
    1224             :                                         if (!is_##TYPE2##_nil(prods[gid])) { \
    1225             :                                                 if (ABSOLUTE(vals[i]) > 1 && \
    1226             :                                                     GDK_##TYPE2##_max / ABSOLUTE(vals[i]) < ABSOLUTE(prods[gid])) { \
    1227             :                                                         if (abort_on_error) \
    1228             :                                                                 goto overflow; \
    1229             :                                                         prods[gid] = TYPE2##_nil; \
    1230             :                                                         nils++;         \
    1231             :                                                 } else {                \
    1232             :                                                         prods[gid] *= vals[i]; \
    1233             :                                                 }                       \
    1234             :                                         }                               \
    1235             :                                 }                                       \
    1236             :                         }                                               \
    1237             :                 }                                                       \
    1238             :         } while (0)
    1239             : 
    1240             : static BUN
    1241          31 : doprod(const void *restrict values, oid seqb, struct canditer *restrict ci, BUN ncand,
    1242             :        void *restrict results, BUN ngrp, int tp1, int tp2,
    1243             :        const oid *restrict gids, bool gidincr, oid min, oid max,
    1244             :        bool skip_nils, bool abort_on_error, bool nil_if_empty, const char *func)
    1245             : {
    1246          31 :         BUN nils = 0;
    1247          31 :         BUN i;
    1248          31 :         oid gid;
    1249          31 :         unsigned int *restrict seen; /* bitmask for groups that we've seen */
    1250             : 
    1251             :         /* allocate bitmap for seen group ids */
    1252          31 :         seen = GDKzalloc(((ngrp + 31) / 32) * sizeof(int));
    1253          31 :         if (seen == NULL) {
    1254             :                 return BUN_NONE;
    1255             :         }
    1256             : 
    1257          31 :         switch (tp2) {
    1258           0 :         case TYPE_bte: {
    1259           0 :                 bte *restrict prods = (bte *) results;
    1260           0 :                 switch (tp1) {
    1261             :                 case TYPE_bte:
    1262           0 :                         AGGR_PROD(bte, bte, sht);
    1263             :                         break;
    1264           0 :                 default:
    1265           0 :                         goto unsupported;
    1266             :                 }
    1267             :                 break;
    1268             :         }
    1269           0 :         case TYPE_sht: {
    1270           0 :                 sht *restrict prods = (sht *) results;
    1271           0 :                 switch (tp1) {
    1272             :                 case TYPE_bte:
    1273           0 :                         AGGR_PROD(bte, sht, int);
    1274             :                         break;
    1275             :                 case TYPE_sht:
    1276           0 :                         AGGR_PROD(sht, sht, int);
    1277             :                         break;
    1278           0 :                 default:
    1279           0 :                         goto unsupported;
    1280             :                 }
    1281             :                 break;
    1282             :         }
    1283           0 :         case TYPE_int: {
    1284           0 :                 int *restrict prods = (int *) results;
    1285           0 :                 switch (tp1) {
    1286             :                 case TYPE_bte:
    1287           0 :                         AGGR_PROD(bte, int, lng);
    1288             :                         break;
    1289             :                 case TYPE_sht:
    1290           0 :                         AGGR_PROD(sht, int, lng);
    1291             :                         break;
    1292             :                 case TYPE_int:
    1293           0 :                         AGGR_PROD(int, int, lng);
    1294             :                         break;
    1295           0 :                 default:
    1296           0 :                         goto unsupported;
    1297             :                 }
    1298             :                 break;
    1299             :         }
    1300             : #ifdef HAVE_HGE
    1301           0 :         case TYPE_lng: {
    1302           0 :                 lng *prods = (lng *) results;
    1303           0 :                 switch (tp1) {
    1304             :                 case TYPE_bte:
    1305           0 :                         AGGR_PROD(bte, lng, hge);
    1306             :                         break;
    1307             :                 case TYPE_sht:
    1308           0 :                         AGGR_PROD(sht, lng, hge);
    1309             :                         break;
    1310             :                 case TYPE_int:
    1311           0 :                         AGGR_PROD(int, lng, hge);
    1312             :                         break;
    1313             :                 case TYPE_lng:
    1314           0 :                         AGGR_PROD(lng, lng, hge);
    1315             :                         break;
    1316           0 :                 default:
    1317           0 :                         goto unsupported;
    1318             :                 }
    1319             :                 break;
    1320             :         }
    1321          25 :         case TYPE_hge: {
    1322          25 :                 hge *prods = (hge *) results;
    1323          25 :                 switch (tp1) {
    1324             :                 case TYPE_bte:
    1325          10 :                         AGGR_PROD_HGE(bte);
    1326             :                         break;
    1327             :                 case TYPE_sht:
    1328          10 :                         AGGR_PROD_HGE(sht);
    1329             :                         break;
    1330             :                 case TYPE_int:
    1331          10 :                         AGGR_PROD_HGE(int);
    1332             :                         break;
    1333             :                 case TYPE_lng:
    1334          20 :                         AGGR_PROD_HGE(lng);
    1335             :                         break;
    1336             :                 case TYPE_hge:
    1337          25 :                         AGGR_PROD_HGE(hge);
    1338             :                         break;
    1339           0 :                 default:
    1340           0 :                         goto unsupported;
    1341             :                 }
    1342             :                 break;
    1343             :         }
    1344             : #else
    1345             :         case TYPE_lng: {
    1346             :                 lng *restrict prods = (lng *) results;
    1347             :                 switch (tp1) {
    1348             :                 case TYPE_bte:
    1349             :                         AGGR_PROD_LNG(bte);
    1350             :                         break;
    1351             :                 case TYPE_sht:
    1352             :                         AGGR_PROD_LNG(sht);
    1353             :                         break;
    1354             :                 case TYPE_int:
    1355             :                         AGGR_PROD_LNG(int);
    1356             :                         break;
    1357             :                 case TYPE_lng:
    1358             :                         AGGR_PROD_LNG(lng);
    1359             :                         break;
    1360             :                 default:
    1361             :                         goto unsupported;
    1362             :                 }
    1363             :                 break;
    1364             :         }
    1365             : #endif
    1366           0 :         case TYPE_flt: {
    1367           0 :                 flt *restrict prods = (flt *) results;
    1368           0 :                 switch (tp1) {
    1369             :                 case TYPE_bte:
    1370           0 :                         AGGR_PROD_FLOAT(bte, flt);
    1371             :                         break;
    1372             :                 case TYPE_sht:
    1373           0 :                         AGGR_PROD_FLOAT(sht, flt);
    1374             :                         break;
    1375             :                 case TYPE_int:
    1376           0 :                         AGGR_PROD_FLOAT(int, flt);
    1377             :                         break;
    1378             :                 case TYPE_lng:
    1379           0 :                         AGGR_PROD_FLOAT(lng, flt);
    1380             :                         break;
    1381             : #ifdef HAVE_HGE
    1382             :                 case TYPE_hge:
    1383           0 :                         AGGR_PROD_FLOAT(hge, flt);
    1384             :                         break;
    1385             : #endif
    1386             :                 case TYPE_flt:
    1387           0 :                         AGGR_PROD_FLOAT(flt, flt);
    1388             :                         break;
    1389           0 :                 default:
    1390           0 :                         goto unsupported;
    1391             :                 }
    1392             :                 break;
    1393             :         }
    1394           6 :         case TYPE_dbl: {
    1395           6 :                 dbl *restrict prods = (dbl *) results;
    1396           6 :                 switch (tp1) {
    1397             :                 case TYPE_bte:
    1398           0 :                         AGGR_PROD_FLOAT(bte, dbl);
    1399             :                         break;
    1400             :                 case TYPE_sht:
    1401           0 :                         AGGR_PROD_FLOAT(sht, dbl);
    1402             :                         break;
    1403             :                 case TYPE_int:
    1404           0 :                         AGGR_PROD_FLOAT(int, dbl);
    1405             :                         break;
    1406             :                 case TYPE_lng:
    1407           0 :                         AGGR_PROD_FLOAT(lng, dbl);
    1408             :                         break;
    1409             : #ifdef HAVE_HGE
    1410             :                 case TYPE_hge:
    1411           0 :                         AGGR_PROD_FLOAT(hge, dbl);
    1412             :                         break;
    1413             : #endif
    1414             :                 case TYPE_flt:
    1415           0 :                         AGGR_PROD_FLOAT(flt, dbl);
    1416             :                         break;
    1417             :                 case TYPE_dbl:
    1418          24 :                         AGGR_PROD_FLOAT(dbl, dbl);
    1419             :                         break;
    1420           0 :                 default:
    1421           0 :                         goto unsupported;
    1422             :                 }
    1423             :                 break;
    1424             :         }
    1425           0 :         default:
    1426           0 :                 goto unsupported;
    1427             :         }
    1428             : 
    1429          31 :         if (nils == 0 && nil_if_empty) {
    1430             :                 /* figure out whether there were any empty groups
    1431             :                  * (that result in a nil value) */
    1432          31 :                 if (ngrp & 0x1F) {
    1433             :                         /* fill last slot */
    1434          31 :                         seen[ngrp >> 5] |= ~0U << (ngrp & 0x1F);
    1435             :                 }
    1436          62 :                 for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
    1437          31 :                         if (seen[i] != ~0U) {
    1438             :                                 nils = 1;
    1439             :                                 break;
    1440             :                         }
    1441             :                 }
    1442             :         }
    1443          31 :         GDKfree(seen);
    1444             : 
    1445          31 :         return nils;
    1446             : 
    1447           0 :   unsupported:
    1448           0 :         GDKfree(seen);
    1449           0 :         GDKerror("%s: type combination (mul(%s)->%s) not supported.\n",
    1450             :                  func, ATOMname(tp1), ATOMname(tp2));
    1451           0 :         return BUN_NONE;
    1452             : 
    1453           0 :   overflow:
    1454           0 :         GDKfree(seen);
    1455           0 :         GDKerror("22003!overflow in calculation.\n");
    1456           0 :         return BUN_NONE;
    1457             : }
    1458             : 
    1459             : /* calculate group products with optional candidates list */
    1460             : BAT *
    1461          16 : BATgroupprod(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    1462             : {
    1463          16 :         const oid *restrict gids;
    1464          16 :         oid min, max;
    1465          16 :         BUN ngrp;
    1466          16 :         BUN nils;
    1467          16 :         BAT *bn;
    1468          16 :         struct canditer ci;
    1469          16 :         BUN ncand;
    1470          16 :         const char *err;
    1471          16 :         lng t0 = 0;
    1472             : 
    1473          16 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    1474             : 
    1475          16 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    1476           0 :                 GDKerror("%s\n", err);
    1477           0 :                 return NULL;
    1478             :         }
    1479          16 :         if (g == NULL) {
    1480           0 :                 GDKerror("b and g must be aligned\n");
    1481           0 :                 return NULL;
    1482             :         }
    1483             : 
    1484          16 :         if (BATcount(b) == 0 || ngrp == 0) {
    1485             :                 /* trivial: no products, so return bat aligned with g
    1486             :                  * with nil in the tail */
    1487           9 :                 return BATconstant(ngrp == 0 ? 0 : min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
    1488             :         }
    1489             : 
    1490           7 :         if ((e == NULL ||
    1491           7 :              (BATcount(e) == BATcount(b) && e->hseqbase == b->hseqbase)) &&
    1492           6 :             (BATtdense(g) || (g->tkey && g->tnonil))) {
    1493             :                 /* trivial: singleton groups, so all results are equal
    1494             :                  * to the inputs (but possibly a different type) */
    1495           6 :                 return BATconvert(b, s, NULL, tp, abort_on_error);
    1496             :         }
    1497             : 
    1498           1 :         bn = BATconstant(min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
    1499           1 :         if (bn == NULL) {
    1500             :                 return NULL;
    1501             :         }
    1502             : 
    1503           1 :         if (BATtdense(g))
    1504             :                 gids = NULL;
    1505             :         else
    1506           1 :                 gids = (const oid *) Tloc(g, 0);
    1507             : 
    1508           2 :         nils = doprod(Tloc(b, 0), b->hseqbase, &ci, ncand, Tloc(bn, 0), ngrp,
    1509           1 :                       b->ttype, tp, gids, true, min, max, skip_nils,
    1510             :                       abort_on_error, true, __func__);
    1511             : 
    1512           1 :         if (nils < BUN_NONE) {
    1513           1 :                 BATsetcount(bn, ngrp);
    1514           1 :                 bn->tkey = BATcount(bn) <= 1;
    1515           1 :                 bn->tsorted = BATcount(bn) <= 1;
    1516           1 :                 bn->trevsorted = BATcount(bn) <= 1;
    1517           1 :                 bn->tnil = nils != 0;
    1518           1 :                 bn->tnonil = nils == 0;
    1519             :         } else {
    1520           0 :                 BBPunfix(bn->batCacheid);
    1521           0 :                 bn = NULL;
    1522             :         }
    1523             : 
    1524           1 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    1525             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    1526             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    1527             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    1528             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    1529             :                   ci.seq, ncand, GDKusec() - t0);
    1530             : 
    1531             :         return bn;
    1532             : }
    1533             : 
    1534             : gdk_return
    1535          30 : BATprod(void *res, int tp, BAT *b, BAT *s, bool skip_nils, bool abort_on_error, bool nil_if_empty)
    1536             : {
    1537          30 :         oid min, max;
    1538          30 :         BUN ngrp;
    1539          30 :         BUN nils;
    1540          30 :         struct canditer ci;
    1541          30 :         BUN ncand;
    1542          30 :         const char *err;
    1543          30 :         lng t0 = 0;
    1544             : 
    1545          30 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    1546             : 
    1547          30 :         if ((err = BATgroupaggrinit(b, NULL, NULL, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    1548           0 :                 GDKerror("%s\n", err);
    1549           0 :                 return GDK_FAIL;
    1550             :         }
    1551          30 :         switch (tp) {
    1552           0 :         case TYPE_bte:
    1553           0 :                 * (bte *) res = nil_if_empty ? bte_nil : (bte) 1;
    1554           0 :                 break;
    1555           0 :         case TYPE_sht:
    1556           0 :                 * (sht *) res = nil_if_empty ? sht_nil : (sht) 1;
    1557           0 :                 break;
    1558           0 :         case TYPE_int:
    1559           0 :                 * (int *) res = nil_if_empty ? int_nil : (int) 1;
    1560           0 :                 break;
    1561           0 :         case TYPE_lng:
    1562           0 :                 * (lng *) res = nil_if_empty ? lng_nil : (lng) 1;
    1563           0 :                 break;
    1564             : #ifdef HAVE_HGE
    1565          25 :         case TYPE_hge:
    1566          25 :                 * (hge *) res = nil_if_empty ? hge_nil : (hge) 1;
    1567          25 :                 break;
    1568             : #endif
    1569           0 :         case TYPE_flt:
    1570           0 :                 * (flt *) res = nil_if_empty ? flt_nil : (flt) 1;
    1571           0 :                 break;
    1572           5 :         case TYPE_dbl:
    1573           5 :                 * (dbl *) res = nil_if_empty ? dbl_nil : (dbl) 1;
    1574           5 :                 break;
    1575           0 :         default:
    1576           0 :                 GDKerror("type combination (prod(%s)->%s) not supported.\n",
    1577             :                          ATOMname(b->ttype), ATOMname(tp));
    1578           0 :                 return GDK_FAIL;
    1579             :         }
    1580          30 :         if (BATcount(b) == 0)
    1581             :                 return GDK_SUCCEED;
    1582          60 :         nils = doprod(Tloc(b, 0), b->hseqbase, &ci, ncand, res, true,
    1583          30 :                       b->ttype, tp, &min, false, min, max,
    1584             :                       skip_nils, abort_on_error, nil_if_empty, __func__);
    1585          30 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",s=" ALGOOPTBATFMT "; "
    1586             :                   "start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    1587             :                   ALGOBATPAR(b), ALGOOPTBATPAR(s),
    1588             :                   ci.seq, ncand, GDKusec() - t0);
    1589          30 :         return nils < BUN_NONE ? GDK_SUCCEED : GDK_FAIL;
    1590             : }
    1591             : 
    1592             : /* ---------------------------------------------------------------------- */
    1593             : /* average */
    1594             : 
    1595             : #define AGGR_AVG(TYPE)                                                  \
    1596             :         do {                                                            \
    1597             :                 const TYPE *restrict vals = (const TYPE *) Tloc(b, 0);  \
    1598             :                 TYPE *restrict avgs = GDKzalloc(ngrp * sizeof(TYPE));   \
    1599             :                 if (avgs == NULL)                                       \
    1600             :                         goto alloc_fail;                                \
    1601             :                 while (ncand > 0) {                                  \
    1602             :                         ncand--;                                        \
    1603             :                         i = canditer_next(&ci) - b->hseqbase;            \
    1604             :                         if (gids == NULL ||                             \
    1605             :                             (gids[i] >= min && gids[i] <= max)) { \
    1606             :                                 if (gids)                               \
    1607             :                                         gid = gids[i] - min;            \
    1608             :                                 else                                    \
    1609             :                                         gid = (oid) i;                  \
    1610             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    1611             :                                         if (!skip_nils)                 \
    1612             :                                                 cnts[gid] = lng_nil;    \
    1613             :                                 } else if (!is_lng_nil(cnts[gid])) {    \
    1614             :                                         AVERAGE_ITER(TYPE, vals[i],     \
    1615             :                                                      avgs[gid],         \
    1616             :                                                      rems[gid],         \
    1617             :                                                      cnts[gid]);        \
    1618             :                                 }                                       \
    1619             :                         }                                               \
    1620             :                 }                                                       \
    1621             :                 for (i = 0; i < ngrp; i++) {                         \
    1622             :                         if (cnts[i] == 0 || is_lng_nil(cnts[i])) {      \
    1623             :                                 dbls[i] = dbl_nil;                      \
    1624             :                                 cnts[i] = 0;                            \
    1625             :                                 nils++;                                 \
    1626             :                         } else {                                        \
    1627             :                                 dbls[i] = avgs[i] + (dbl) rems[i] / cnts[i]; \
    1628             :                         }                                               \
    1629             :                 }                                                       \
    1630             :                 GDKfree(avgs);                                          \
    1631             :         } while (0)
    1632             : 
    1633             : #define AGGR_AVG_FLOAT(TYPE)                                            \
    1634             :         do {                                                            \
    1635             :                 const TYPE *restrict vals = (const TYPE *) Tloc(b, 0);  \
    1636             :                 for (i = 0; i < ngrp; i++)                           \
    1637             :                         dbls[i] = 0;                                    \
    1638             :                 while (ncand > 0) {                                  \
    1639             :                         ncand--;                                        \
    1640             :                         i = canditer_next(&ci) - b->hseqbase;            \
    1641             :                         if (gids == NULL ||                             \
    1642             :                             (gids[i] >= min && gids[i] <= max)) { \
    1643             :                                 if (gids)                               \
    1644             :                                         gid = gids[i] - min;            \
    1645             :                                 else                                    \
    1646             :                                         gid = (oid) i;                  \
    1647             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    1648             :                                         if (!skip_nils)                 \
    1649             :                                                 cnts[gid] = lng_nil;    \
    1650             :                                 } else if (!is_lng_nil(cnts[gid])) {    \
    1651             :                                         AVERAGE_ITER_FLOAT(TYPE, vals[i], \
    1652             :                                                            dbls[gid],   \
    1653             :                                                            cnts[gid]);  \
    1654             :                                 }                                       \
    1655             :                         }                                               \
    1656             :                 }                                                       \
    1657             :                 for (i = 0; i < ngrp; i++) {                         \
    1658             :                         if (cnts[i] == 0 || is_lng_nil(cnts[i])) {      \
    1659             :                                 dbls[i] = dbl_nil;                      \
    1660             :                                 cnts[i] = 0;                            \
    1661             :                                 nils++;                                 \
    1662             :                         }                                               \
    1663             :                 }                                                       \
    1664             :         } while (0)
    1665             : 
    1666             : /* calculate group averages with optional candidates list */
    1667             : gdk_return
    1668         337 : BATgroupavg(BAT **bnp, BAT **cntsp, BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error, int scale)
    1669             : {
    1670         337 :         const oid *restrict gids;
    1671         337 :         oid gid;
    1672         337 :         oid min, max;
    1673         337 :         BUN i, ngrp;
    1674         337 :         BUN nils = 0;
    1675         337 :         BUN *restrict rems = NULL;
    1676         337 :         lng *restrict cnts = NULL;
    1677         337 :         dbl *restrict dbls;
    1678         337 :         BAT *bn = NULL, *cn = NULL;
    1679         337 :         struct canditer ci;
    1680         337 :         BUN ncand;
    1681         337 :         const char *err;
    1682         337 :         lng t0 = 0;
    1683             : 
    1684         337 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    1685             : 
    1686         337 :         assert(tp == TYPE_dbl);
    1687         337 :         (void) tp;              /* compatibility (with other BATgroup*
    1688             :                                  * functions) argument */
    1689             : 
    1690         337 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    1691           0 :                 GDKerror("%s\n", err);
    1692           0 :                 return GDK_FAIL;
    1693             :         }
    1694         337 :         if (g == NULL) {
    1695           0 :                 GDKerror("b and g must be aligned\n");
    1696           0 :                 return GDK_FAIL;
    1697             :         }
    1698             : 
    1699         337 :         if (BATcount(b) == 0 || ngrp == 0) {
    1700             :                 /* trivial: no averages, so return bat aligned with g
    1701             :                  * with nil in the tail */
    1702          22 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
    1703          22 :                 if (bn == NULL) {
    1704             :                         return GDK_FAIL;
    1705             :                 }
    1706          22 :                 if (cntsp) {
    1707          14 :                         lng zero = 0;
    1708          14 :                         if ((cn = BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT)) == NULL) {
    1709           0 :                                 BBPreclaim(bn);
    1710           0 :                                 return GDK_FAIL;
    1711             :                         }
    1712          14 :                         *cntsp = cn;
    1713             :                 }
    1714          22 :                 *bnp = bn;
    1715          22 :                 return GDK_SUCCEED;
    1716             :         }
    1717             : 
    1718         315 :         if ((!skip_nils || cntsp == NULL || b->tnonil) &&
    1719         257 :             (e == NULL ||
    1720         257 :              (BATcount(e) == BATcount(b) && e->hseqbase == b->hseqbase)) &&
    1721         125 :             (BATtdense(g) || (g->tkey && g->tnonil))) {
    1722             :                 /* trivial: singleton groups, so all results are equal
    1723             :                  * to the inputs (but possibly a different type) */
    1724         125 :                 if ((bn = BATconvert(b, s, NULL, TYPE_dbl, abort_on_error)) == NULL)
    1725             :                         return GDK_FAIL;
    1726         125 :                 if (cntsp) {
    1727          22 :                         lng one = 1;
    1728          22 :                         if ((cn = BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &one, ngrp, TRANSIENT)) == NULL) {
    1729           0 :                                 BBPreclaim(bn);
    1730           0 :                                 return GDK_FAIL;
    1731             :                         }
    1732          22 :                         *cntsp = cn;
    1733             :                 }
    1734         125 :                 *bnp = bn;
    1735         125 :                 return GDK_SUCCEED;
    1736             :         }
    1737             : 
    1738             :         /* allocate temporary space to do per group calculations */
    1739         190 :         switch (b->ttype) {
    1740          74 :         case TYPE_bte:
    1741             :         case TYPE_sht:
    1742             :         case TYPE_int:
    1743             :         case TYPE_lng:
    1744             : #ifdef HAVE_HGE
    1745             :         case TYPE_hge:
    1746             : #endif
    1747          74 :                 rems = GDKzalloc(ngrp * sizeof(BUN));
    1748          74 :                 if (rems == NULL)
    1749           0 :                         goto alloc_fail;
    1750             :                 break;
    1751             :         default:
    1752             :                 break;
    1753             :         }
    1754         190 :         if (cntsp) {
    1755         128 :                 if ((cn = COLnew(min, TYPE_lng, ngrp, TRANSIENT)) == NULL)
    1756           0 :                         goto alloc_fail;
    1757         128 :                 cnts = (lng *) Tloc(cn, 0);
    1758         128 :                 memset(cnts, 0, ngrp * sizeof(lng));
    1759             :         } else {
    1760          62 :                 cnts = GDKzalloc(ngrp * sizeof(lng));
    1761          62 :                 if (cnts == NULL)
    1762           0 :                         goto alloc_fail;
    1763             :         }
    1764             : 
    1765         190 :         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    1766         190 :         if (bn == NULL)
    1767           0 :                 goto alloc_fail;
    1768         190 :         dbls = (dbl *) Tloc(bn, 0);
    1769             : 
    1770         190 :         if (BATtdense(g))
    1771             :                 gids = NULL;
    1772             :         else
    1773         161 :                 gids = (const oid *) Tloc(g, 0);
    1774             : 
    1775         190 :         switch (b->ttype) {
    1776           0 :         case TYPE_bte:
    1777           0 :                 AGGR_AVG(bte);
    1778           0 :                 break;
    1779           4 :         case TYPE_sht:
    1780          31 :                 AGGR_AVG(sht);
    1781           4 :                 break;
    1782          57 :         case TYPE_int:
    1783    15228700 :                 AGGR_AVG(int);
    1784          57 :                 break;
    1785          13 :         case TYPE_lng:
    1786          90 :                 AGGR_AVG(lng);
    1787          13 :                 break;
    1788             : #ifdef HAVE_HGE
    1789           0 :         case TYPE_hge:
    1790           0 :                 AGGR_AVG(hge);
    1791           0 :                 break;
    1792             : #endif
    1793           6 :         case TYPE_flt:
    1794        2959 :                 AGGR_AVG_FLOAT(flt);
    1795             :                 break;
    1796         110 :         case TYPE_dbl:
    1797     1419470 :                 AGGR_AVG_FLOAT(dbl);
    1798             :                 break;
    1799           0 :         default:
    1800           0 :                 GDKfree(rems);
    1801           0 :                 if (cn)
    1802           0 :                         BBPreclaim(cn);
    1803             :                 else
    1804           0 :                         GDKfree(cnts);
    1805           0 :                 BBPunfix(bn->batCacheid);
    1806           0 :                 GDKerror("type (%s) not supported.\n", ATOMname(b->ttype));
    1807           0 :                 return GDK_FAIL;
    1808             :         }
    1809         190 :         GDKfree(rems);
    1810         190 :         if (cn == NULL)
    1811          62 :                 GDKfree(cnts);
    1812             :         else {
    1813         128 :                 BATsetcount(cn, ngrp);
    1814         128 :                 cn->tkey = BATcount(cn) <= 1;
    1815         128 :                 cn->tsorted = BATcount(cn) <= 1;
    1816         128 :                 cn->trevsorted = BATcount(cn) <= 1;
    1817         128 :                 cn->tnil = false;
    1818         128 :                 cn->tnonil = true;
    1819         128 :                 *cntsp = cn;
    1820             :         }
    1821         190 :         if (scale != 0) {
    1822           0 :                 dbl fac = pow(10.0, (double) scale);
    1823           0 :                 for (i = 0; i < ngrp; i++) {
    1824           0 :                         if (!is_dbl_nil(dbls[i]))
    1825           0 :                                 dbls[i] *= fac;
    1826             :                 }
    1827             :         }
    1828         190 :         BATsetcount(bn, ngrp);
    1829         190 :         bn->tkey = BATcount(bn) <= 1;
    1830         190 :         bn->tsorted = BATcount(bn) <= 1;
    1831         190 :         bn->trevsorted = BATcount(bn) <= 1;
    1832         190 :         bn->tnil = nils != 0;
    1833         190 :         bn->tnonil = nils == 0;
    1834         190 :         *bnp = bn;
    1835         190 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    1836             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    1837             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    1838             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    1839             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    1840             :                   ci.seq, ncand, GDKusec() - t0);
    1841             :         return GDK_SUCCEED;
    1842             : 
    1843           0 :   alloc_fail:
    1844           0 :         if (bn)
    1845           0 :                 BBPunfix(bn->batCacheid);
    1846           0 :         GDKfree(rems);
    1847           0 :         if (cntsp) {
    1848           0 :                 BBPreclaim(*cntsp);
    1849           0 :                 *cntsp = NULL;
    1850           0 :         } else if (cnts) {
    1851           0 :                 GDKfree(cnts);
    1852             :         }
    1853             :         return GDK_FAIL;
    1854             : }
    1855             : 
    1856             : #define AVERAGE_TYPE_LNG_HGE(TYPE,lng_hge)                              \
    1857             :         do {                                                            \
    1858             :                 TYPE x, a;                                              \
    1859             :                                                                         \
    1860             :                 /* first try to calculate the sum of all values into a */ \
    1861             :                 /* lng/hge */                                           \
    1862             :                 while (ncand > 0) {                                  \
    1863             :                         ncand--;                                        \
    1864             :                         i = canditer_next(&ci) - b->hseqbase;            \
    1865             :                         x = ((const TYPE *) src)[i];                    \
    1866             :                         if (is_##TYPE##_nil(x))                         \
    1867             :                                 continue;                               \
    1868             :                         ADD_WITH_CHECK(x, sum,                          \
    1869             :                                        lng_hge, sum,                    \
    1870             :                                        GDK_##lng_hge##_max,             \
    1871             :                                        goto overflow##TYPE);            \
    1872             :                         /* don't count value until after overflow check */ \
    1873             :                         n++;                                            \
    1874             :                 }                                                       \
    1875             :                 /* the sum fit, so now we can calculate the average */  \
    1876             :                 *avg = n > 0 ? (dbl) sum / n : dbl_nil;                      \
    1877             :                 if (0) {                                                \
    1878             :                   overflow##TYPE:                                       \
    1879             :                         /* we get here if sum(x[0],...,x[i]) doesn't */ \
    1880             :                         /* fit in a lng/hge but sum(x[0],...,x[i-1]) did */ \
    1881             :                         /* the variable sum contains that sum */        \
    1882             :                         /* the rest of the calculation is done */       \
    1883             :                         /* according to the loop invariant described */ \
    1884             :                         /* in the below loop */                         \
    1885             :                         /* note that n necessarily is > 0 (else no */        \
    1886             :                         /* overflow possible) */                        \
    1887             :                         assert(n > 0);                                       \
    1888             :                         if (sum >= 0) {                                      \
    1889             :                                 a = (TYPE) (sum / (lng_hge) n); /* this fits */ \
    1890             :                                 r = (BUN) (sum % (SBUN) n);             \
    1891             :                         } else {                                        \
    1892             :                                 sum = -sum;                             \
    1893             :                                 a = - (TYPE) (sum / (lng_hge) n); /* this fits */ \
    1894             :                                 r = (BUN) (sum % (SBUN) n);             \
    1895             :                                 if (r) {                                \
    1896             :                                         a--;                            \
    1897             :                                         r = n - r;                      \
    1898             :                                 }                                       \
    1899             :                         }                                               \
    1900             :                         while (ncand > 0) {                          \
    1901             :                                 /* loop invariant: */                   \
    1902             :                                 /* a + r/n == average(x[0],...,x[n]); */ \
    1903             :                                 /* 0 <= r < n */                  \
    1904             :                                 ncand--;                                \
    1905             :                                 i = canditer_next(&ci) - b->hseqbase;    \
    1906             :                                 x = ((const TYPE *) src)[i];            \
    1907             :                                 if (is_##TYPE##_nil(x))                 \
    1908             :                                         continue;                       \
    1909             :                                 AVERAGE_ITER(TYPE, x, a, r, n);         \
    1910             :                         }                                               \
    1911             :                         *avg = a + (dbl) r / n;                         \
    1912             :                 }                                                       \
    1913             :         } while (0)
    1914             : 
    1915             : #ifdef HAVE_HGE
    1916             : #define AVERAGE_TYPE(TYPE) AVERAGE_TYPE_LNG_HGE(TYPE,hge)
    1917             : #else
    1918             : #define AVERAGE_TYPE(TYPE) AVERAGE_TYPE_LNG_HGE(TYPE,lng)
    1919             : #endif
    1920             : 
    1921             : #define AVERAGE_FLOATTYPE(TYPE)                                 \
    1922             :         do {                                                    \
    1923             :                 double a = 0;                                   \
    1924             :                 TYPE x;                                         \
    1925             :                 while (ncand > 0) {                          \
    1926             :                         ncand--;                                \
    1927             :                         i = canditer_next(&ci) - b->hseqbase;    \
    1928             :                         x = ((const TYPE *) src)[i];            \
    1929             :                         if (is_##TYPE##_nil(x))                 \
    1930             :                                 continue;                       \
    1931             :                         AVERAGE_ITER_FLOAT(TYPE, x, a, n);      \
    1932             :                 }                                               \
    1933             :                 *avg = n > 0 ? a : dbl_nil;                  \
    1934             :         } while (0)
    1935             : 
    1936             : gdk_return
    1937        4104 : BATcalcavg(BAT *b, BAT *s, dbl *avg, BUN *vals, int scale)
    1938             : {
    1939        4104 :         BUN n = 0, r = 0, i = 0;
    1940             : #ifdef HAVE_HGE
    1941        4104 :         hge sum = 0;
    1942             : #else
    1943             :         lng sum = 0;
    1944             : #endif
    1945        4104 :         struct canditer ci;
    1946        4104 :         BUN ncand;
    1947        4104 :         const void *restrict src;
    1948             :         /* these two needed for ADD_WITH_CHECK macro */
    1949        4104 :         bool abort_on_error = true;
    1950        4104 :         BUN nils = 0;
    1951             : 
    1952        4104 :         ncand = canditer_init(&ci, b, s);
    1953             : 
    1954        4102 :         src = Tloc(b, 0);
    1955             : 
    1956        4102 :         switch (b->ttype) {
    1957             :         case TYPE_bte:
    1958          10 :                 AVERAGE_TYPE(bte);
    1959             :                 break;
    1960             :         case TYPE_sht:
    1961           0 :                 AVERAGE_TYPE(sht);
    1962             :                 break;
    1963             :         case TYPE_int:
    1964     2390480 :                 AVERAGE_TYPE(int);
    1965             :                 break;
    1966             :         case TYPE_lng:
    1967    10000000 :                 AVERAGE_TYPE(lng);
    1968             :                 break;
    1969             : #ifdef HAVE_HGE
    1970             :         case TYPE_hge:
    1971           0 :                 AVERAGE_TYPE(hge);
    1972             :                 break;
    1973             : #endif
    1974             :         case TYPE_flt:
    1975   100000000 :                 AVERAGE_FLOATTYPE(flt);
    1976           8 :                 break;
    1977             :         case TYPE_dbl:
    1978    10366900 :                 AVERAGE_FLOATTYPE(dbl);
    1979         150 :                 break;
    1980           0 :         default:
    1981           0 :                 GDKerror("average of type %s unsupported.\n",
    1982             :                          ATOMname(b->ttype));
    1983           0 :                 return GDK_FAIL;
    1984             :         }
    1985        4103 :         if (scale != 0 && !is_dbl_nil(*avg))
    1986           0 :                 *avg *= pow(10.0, (double) scale);
    1987        4103 :         if (vals)
    1988        4103 :                 *vals = n;
    1989             :         return GDK_SUCCEED;
    1990             : }
    1991             : 
    1992             : /* ---------------------------------------------------------------------- */
    1993             : /* count */
    1994             : 
    1995             : #define AGGR_COUNT(TYPE)                                                \
    1996             :         do {                                                            \
    1997             :                 const TYPE *restrict vals = (const TYPE *) Tloc(b, 0);  \
    1998             :                 while (ncand > 0) {                                  \
    1999             :                         ncand--;                                        \
    2000             :                         i = canditer_next(&ci) - b->hseqbase;            \
    2001             :                         if (gids == NULL ||                             \
    2002             :                             (gids[i] >= min && gids[i] <= max)) { \
    2003             :                                 if (gids)                               \
    2004             :                                         gid = gids[i] - min;            \
    2005             :                                 else                                    \
    2006             :                                         gid = (oid) i;                  \
    2007             :                                 if (!is_##TYPE##_nil(vals[i])) {        \
    2008             :                                         cnts[gid]++;                    \
    2009             :                                 }                                       \
    2010             :                         }                                               \
    2011             :                 }                                                       \
    2012             :         } while (0)
    2013             : 
    2014             : /* calculate group counts with optional candidates list */
    2015             : BAT *
    2016        9861 : BATgroupcount(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    2017             : {
    2018        9861 :         const oid *restrict gids;
    2019        9861 :         oid gid;
    2020        9861 :         oid min, max;
    2021        9861 :         BUN i, ngrp;
    2022        9861 :         lng *restrict cnts;
    2023        9861 :         BAT *bn = NULL;
    2024        9861 :         int t;
    2025        9861 :         const void *nil;
    2026        9861 :         int (*atomcmp)(const void *, const void *);
    2027        9861 :         BATiter bi;
    2028        9861 :         struct canditer ci;
    2029        9861 :         BUN ncand;
    2030        9861 :         const char *err;
    2031        9861 :         lng t0 = 0;
    2032             : 
    2033        9861 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    2034             : 
    2035        9861 :         assert(tp == TYPE_lng);
    2036        9861 :         (void) tp;              /* compatibility (with other BATgroup* */
    2037        9861 :         (void) abort_on_error;  /* functions) argument */
    2038             : 
    2039        9861 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    2040           0 :                 GDKerror("%s\n", err);
    2041           0 :                 return NULL;
    2042             :         }
    2043        9858 :         if (g == NULL) {
    2044           0 :                 GDKerror("b and g must be aligned\n");
    2045           0 :                 return NULL;
    2046             :         }
    2047             : 
    2048        9858 :         if (BATcount(b) == 0 || ngrp == 0) {
    2049             :                 /* trivial: no counts, so return bat aligned with g
    2050             :                  * with zero in the tail */
    2051        5344 :                 lng zero = 0;
    2052        5344 :                 return BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT);
    2053             :         }
    2054             : 
    2055        4514 :         bn = COLnew(min, TYPE_lng, ngrp, TRANSIENT);
    2056        4517 :         if (bn == NULL)
    2057             :                 return NULL;
    2058        4517 :         cnts = (lng *) Tloc(bn, 0);
    2059        4517 :         memset(cnts, 0, ngrp * sizeof(lng));
    2060             : 
    2061        4517 :         if (BATtdense(g))
    2062             :                 gids = NULL;
    2063             :         else
    2064        4095 :                 gids = (const oid *) Tloc(g, 0);
    2065             : 
    2066        4517 :         if (!skip_nils || b->tnonil) {
    2067             :                 /* if nils are nothing special, or if there are no
    2068             :                  * nils, we don't need to look at the values at all */
    2069        4462 :                 if (gids) {
    2070     4685140 :                         while (ncand > 0) {
    2071     4681070 :                                 ncand--;
    2072     4681070 :                                 i = canditer_next(&ci) - b->hseqbase;
    2073     4681070 :                                 if (gids[i] >= min && gids[i] <= max)
    2074     4681070 :                                         cnts[gids[i] - min]++;
    2075             :                         }
    2076             :                 } else {
    2077       35174 :                         while (ncand > 0) {
    2078       34782 :                                 ncand--;
    2079       34782 :                                 i = canditer_next(&ci) - b->hseqbase;
    2080       34782 :                                 cnts[i] = 1;
    2081             :                         }
    2082             :                 }
    2083             :         } else {
    2084          55 :                 t = b->ttype;
    2085          55 :                 nil = ATOMnilptr(t);
    2086          55 :                 atomcmp = ATOMcompare(t);
    2087          55 :                 t = ATOMbasetype(t);
    2088             : 
    2089          55 :                 switch (t) {
    2090           2 :                 case TYPE_bte:
    2091           7 :                         AGGR_COUNT(bte);
    2092             :                         break;
    2093           0 :                 case TYPE_sht:
    2094           0 :                         AGGR_COUNT(sht);
    2095             :                         break;
    2096          43 :                 case TYPE_int:
    2097       30828 :                         AGGR_COUNT(int);
    2098             :                         break;
    2099           0 :                 case TYPE_lng:
    2100           0 :                         AGGR_COUNT(lng);
    2101             :                         break;
    2102             : #ifdef HAVE_HGE
    2103           0 :                 case TYPE_hge:
    2104           0 :                         AGGR_COUNT(hge);
    2105             :                         break;
    2106             : #endif
    2107           0 :                 case TYPE_flt:
    2108           0 :                         AGGR_COUNT(flt);
    2109             :                         break;
    2110           5 :                 case TYPE_dbl:
    2111         150 :                         AGGR_COUNT(dbl);
    2112             :                         break;
    2113           5 :                 default:
    2114           5 :                         bi = bat_iterator(b);
    2115             : 
    2116         117 :                         while (ncand > 0) {
    2117         113 :                                 ncand--;
    2118         113 :                                 i = canditer_next(&ci) - b->hseqbase;
    2119         112 :                                 if (gids == NULL ||
    2120         110 :                                     (gids[i] >= min && gids[i] <= max)) {
    2121         111 :                                         if (gids)
    2122         109 :                                                 gid = gids[i] - min;
    2123             :                                         else
    2124             :                                                 gid = (oid) i;
    2125         111 :                                         if ((*atomcmp)(BUNtail(bi, i), nil) != 0) {
    2126          53 :                                                 cnts[gid]++;
    2127             :                                         }
    2128             :                                 }
    2129             :                         }
    2130             :                         break;
    2131             :                 }
    2132             :         }
    2133        4516 :         BATsetcount(bn, ngrp);
    2134        4516 :         bn->tkey = BATcount(bn) <= 1;
    2135        4516 :         bn->tsorted = BATcount(bn) <= 1;
    2136        4516 :         bn->trevsorted = BATcount(bn) <= 1;
    2137        4516 :         bn->tnil = false;
    2138        4516 :         bn->tnonil = true;
    2139        4516 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    2140             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    2141             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    2142             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    2143             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    2144             :                   ci.seq, ncand, GDKusec() - t0);
    2145             :         return bn;
    2146             : }
    2147             : 
    2148             : /* calculate group sizes (number of TRUE values) with optional
    2149             :  * candidates list */
    2150             : BAT *
    2151           0 : BATgroupsize(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    2152             : {
    2153           0 :         const oid *restrict gids;
    2154           0 :         oid min, max;
    2155           0 :         BUN i, ngrp;
    2156           0 :         const bit *restrict bits;
    2157           0 :         lng *restrict cnts;
    2158           0 :         BAT *bn = NULL;
    2159           0 :         struct canditer ci;
    2160           0 :         BUN ncand;
    2161           0 :         const char *err;
    2162           0 :         lng t0 = 0;
    2163             : 
    2164           0 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    2165             : 
    2166           0 :         assert(tp == TYPE_lng);
    2167           0 :         assert(b->ttype == TYPE_bit);
    2168             :         /* compatibility arguments */
    2169           0 :         (void) tp;
    2170           0 :         (void) abort_on_error;
    2171           0 :         (void) skip_nils;
    2172             : 
    2173           0 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    2174           0 :                 GDKerror("%s\n", err);
    2175           0 :                 return NULL;
    2176             :         }
    2177           0 :         if (g == NULL) {
    2178           0 :                 GDKerror("b and g must be aligned\n");
    2179           0 :                 return NULL;
    2180             :         }
    2181             : 
    2182           0 :         if (BATcount(b) == 0 || ngrp == 0) {
    2183             :                 /* trivial: no products, so return bat aligned with g
    2184             :                  * with zero in the tail */
    2185           0 :                 lng zero = 0;
    2186           0 :                 return BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT);
    2187             :         }
    2188             : 
    2189           0 :         bn = COLnew(min, TYPE_lng, ngrp, TRANSIENT);
    2190           0 :         if (bn == NULL)
    2191             :                 return NULL;
    2192           0 :         cnts = (lng *) Tloc(bn, 0);
    2193           0 :         memset(cnts, 0, ngrp * sizeof(lng));
    2194             : 
    2195           0 :         if (BATtdense(g))
    2196             :                 gids = NULL;
    2197             :         else
    2198           0 :                 gids = (const oid *) Tloc(g, 0);
    2199             : 
    2200           0 :         bits = (const bit *) Tloc(b, 0);
    2201             : 
    2202           0 :         while (ncand > 0) {
    2203           0 :                 ncand--;
    2204           0 :                 i = canditer_next(&ci) - b->hseqbase;
    2205           0 :                 if (bits[i] == 1 &&
    2206           0 :                     (gids == NULL || (gids[i] >= min && gids[i] <= max))) {
    2207           0 :                         cnts[gids ? gids[i] - min : (oid) i]++;
    2208             :                 }
    2209             :         }
    2210           0 :         BATsetcount(bn, ngrp);
    2211           0 :         bn->tkey = BATcount(bn) <= 1;
    2212           0 :         bn->tsorted = BATcount(bn) <= 1;
    2213           0 :         bn->trevsorted = BATcount(bn) <= 1;
    2214           0 :         bn->tnil = false;
    2215           0 :         bn->tnonil = true;
    2216           0 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    2217             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    2218             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    2219             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    2220             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    2221             :                   ci.seq, ncand, GDKusec() - t0);
    2222             :         return bn;
    2223             : }
    2224             : 
    2225             : /* ---------------------------------------------------------------------- */
    2226             : /* min and max */
    2227             : 
    2228             : #define AGGR_CMP(TYPE, OP)                                              \
    2229             :         do {                                                            \
    2230             :                 const TYPE *restrict vals = (const TYPE *) Tloc(b, 0);  \
    2231             :                 if (ngrp == cnt) {                                      \
    2232             :                         /* single element groups */                     \
    2233             :                         while (ncand > 0) {                          \
    2234             :                                 ncand--;                                \
    2235             :                                 i = canditer_next(ci) - b->hseqbase; \
    2236             :                                 if (!skip_nils ||                       \
    2237             :                                     !is_##TYPE##_nil(vals[i])) {        \
    2238             :                                         oids[i] = i + b->hseqbase;   \
    2239             :                                         nils--;                         \
    2240             :                                 }                                       \
    2241             :                         }                                               \
    2242             :                 } else {                                                \
    2243             :                         gid = 0; /* in case gids == NULL */             \
    2244             :                         while (ncand > 0) {                          \
    2245             :                                 ncand--;                                \
    2246             :                                 i = canditer_next(ci) - b->hseqbase; \
    2247             :                                 if (gids == NULL ||                     \
    2248             :                                     (gids[i] >= min && gids[i] <= max)) { \
    2249             :                                         if (gids)                       \
    2250             :                                                 gid = gids[i] - min;    \
    2251             :                                         if (!skip_nils || !is_##TYPE##_nil(vals[i])) { \
    2252             :                                                 if (is_oid_nil(oids[gid])) { \
    2253             :                                                         oids[gid] = i + b->hseqbase; \
    2254             :                                                         nils--;         \
    2255             :                                                 } else if (!is_##TYPE##_nil(vals[oids[gid] - b->hseqbase]) && \
    2256             :                                                            (is_##TYPE##_nil(vals[i]) || \
    2257             :                                                             OP(vals[i], vals[oids[gid] - b->hseqbase]))) \
    2258             :                                                         oids[gid] = i + b->hseqbase; \
    2259             :                                         }                               \
    2260             :                                 }                                       \
    2261             :                         }                                               \
    2262             :                 }                                                       \
    2263             :         } while (0)
    2264             : 
    2265             : /* calculate group minimums with optional candidates list
    2266             :  *
    2267             :  * note that this functions returns *positions* of where the minimum
    2268             :  * values occur */
    2269             : static BUN
    2270         773 : do_groupmin(oid *restrict oids, BAT *b, const oid *restrict gids, BUN ngrp,
    2271             :             oid min, oid max, struct canditer *restrict ci, BUN ncand,
    2272             :             BUN cnt, bool skip_nils, bool gdense)
    2273             : {
    2274         773 :         oid gid;
    2275         773 :         BUN i, nils;
    2276         773 :         int t;
    2277         773 :         const void *nil;
    2278         773 :         int (*atomcmp)(const void *, const void *);
    2279         773 :         BATiter bi;
    2280             : 
    2281         773 :         nils = ngrp;
    2282       16828 :         for (i = 0; i < ngrp; i++)
    2283       16055 :                 oids[i] = oid_nil;
    2284         773 :         if (cnt == 0)
    2285             :                 return nils;
    2286             : 
    2287         773 :         t = b->ttype;
    2288         773 :         nil = ATOMnilptr(t);
    2289         773 :         atomcmp = ATOMcompare(t);
    2290         773 :         t = ATOMbasetype(t);
    2291         773 :         switch (t) {
    2292          28 :         case TYPE_bte:
    2293       28723 :                 AGGR_CMP(bte, LT);
    2294             :                 break;
    2295          26 :         case TYPE_sht:
    2296        2046 :                 AGGR_CMP(sht, LT);
    2297             :                 break;
    2298         506 :         case TYPE_int:
    2299      388480 :                 AGGR_CMP(int, LT);
    2300             :                 break;
    2301          57 :         case TYPE_lng:
    2302       53892 :                 AGGR_CMP(lng, LT);
    2303             :                 break;
    2304             : #ifdef HAVE_HGE
    2305           2 :         case TYPE_hge:
    2306           8 :                 AGGR_CMP(hge, LT);
    2307             :                 break;
    2308             : #endif
    2309           2 :         case TYPE_flt:
    2310          10 :                 AGGR_CMP(flt, LT);
    2311             :                 break;
    2312          43 :         case TYPE_dbl:
    2313         530 :                 AGGR_CMP(dbl, LT);
    2314             :                 break;
    2315           0 :         case TYPE_void:
    2316           0 :                 if (!gdense && gids == NULL) {
    2317           0 :                         oids[0] = canditer_next(ci);
    2318           0 :                         nils--;
    2319           0 :                         break;
    2320             :                 }
    2321             :                 /* fall through */
    2322             :         default:
    2323         109 :                 assert(b->ttype != TYPE_oid);
    2324         109 :                 bi = bat_iterator(b);
    2325             : 
    2326         109 :                 if (gdense) {
    2327             :                         /* single element groups */
    2328           0 :                         while (ncand > 0) {
    2329           0 :                                 ncand--;
    2330           0 :                                 i = canditer_next(ci) - b->hseqbase;
    2331           0 :                                 if (!skip_nils ||
    2332           0 :                                     (*atomcmp)(BUNtail(bi, i), nil) != 0) {
    2333           0 :                                         oids[i] = i + b->hseqbase;
    2334           0 :                                         nils--;
    2335             :                                 }
    2336             :                         }
    2337             :                 } else {
    2338             :                         gid = 0; /* in case gids == NULL */
    2339      647466 :                         while (ncand > 0) {
    2340      647357 :                                 ncand--;
    2341      647357 :                                 i = canditer_next(ci) - b->hseqbase;
    2342      647214 :                                 if (gids == NULL ||
    2343          21 :                                     (gids[i] >= min && gids[i] <= max)) {
    2344      647214 :                                         const void *v = BUNtail(bi, i);
    2345      647214 :                                         if (gids)
    2346         532 :                                                 gid = gids[i] - min;
    2347     1300300 :                                         if (!skip_nils ||
    2348      646715 :                                             (*atomcmp)(v, nil) != 0) {
    2349      650791 :                                                 if (is_oid_nil(oids[gid])) {
    2350         105 :                                                         oids[gid] = i + b->hseqbase;
    2351         105 :                                                         nils--;
    2352      650686 :                                                 } else if (t != TYPE_void) {
    2353      650748 :                                                         const void *g = BUNtail(bi, (BUN) (oids[gid] - b->hseqbase));
    2354     1290500 :                                                         if ((*atomcmp)(g, nil) != 0 &&
    2355     1290060 :                                                             ((*atomcmp)(v, nil) == 0 ||
    2356      639756 :                                                              LT((*atomcmp)(v, g), 0)))
    2357         130 :                                                                 oids[gid] = i + b->hseqbase;
    2358             :                                                 }
    2359             :                                         }
    2360             :                                 }
    2361             :                         }
    2362             :                 }
    2363             :                 break;
    2364             :         }
    2365             : 
    2366             :         return nils;
    2367             : }
    2368             : 
    2369             : /* calculate group maximums with optional candidates list
    2370             :  *
    2371             :  * note that this functions returns *positions* of where the maximum
    2372             :  * values occur */
    2373             : static BUN
    2374         780 : do_groupmax(oid *restrict oids, BAT *b, const oid *restrict gids, BUN ngrp,
    2375             :             oid min, oid max, struct canditer *restrict ci, BUN ncand,
    2376             :             BUN cnt, bool skip_nils, bool gdense)
    2377             : {
    2378         780 :         oid gid;
    2379         780 :         BUN i, nils;
    2380         780 :         int t;
    2381         780 :         const void *nil;
    2382         780 :         int (*atomcmp)(const void *, const void *);
    2383         780 :         BATiter bi;
    2384             : 
    2385         780 :         nils = ngrp;
    2386      198613 :         for (i = 0; i < ngrp; i++)
    2387      197833 :                 oids[i] = oid_nil;
    2388         780 :         if (cnt == 0)
    2389             :                 return nils;
    2390             : 
    2391         780 :         t = b->ttype;
    2392         780 :         nil = ATOMnilptr(t);
    2393         780 :         atomcmp = ATOMcompare(t);
    2394         780 :         t = ATOMbasetype(t);
    2395         780 :         switch (t) {
    2396          32 :         case TYPE_bte:
    2397       28736 :                 AGGR_CMP(bte, GT);
    2398             :                 break;
    2399          26 :         case TYPE_sht:
    2400        2046 :                 AGGR_CMP(sht, GT);
    2401             :                 break;
    2402         501 :         case TYPE_int:
    2403      284488 :                 AGGR_CMP(int, GT);
    2404             :                 break;
    2405          65 :         case TYPE_lng:
    2406       59721 :                 AGGR_CMP(lng, GT);
    2407             :                 break;
    2408             : #ifdef HAVE_HGE
    2409          10 :         case TYPE_hge:
    2410     5654390 :                 AGGR_CMP(hge, GT);
    2411             :                 break;
    2412             : #endif
    2413           2 :         case TYPE_flt:
    2414          10 :                 AGGR_CMP(flt, GT);
    2415             :                 break;
    2416          34 :         case TYPE_dbl:
    2417         407 :                 AGGR_CMP(dbl, GT);
    2418             :                 break;
    2419           0 :         case TYPE_void:
    2420           0 :                 if (!gdense && gids == NULL) {
    2421           0 :                         oids[0] = canditer_last(ci);
    2422           0 :                         nils--;
    2423           0 :                         break;
    2424             :                 }
    2425             :                 /* fall through */
    2426             :         default:
    2427         110 :                 assert(b->ttype != TYPE_oid);
    2428         110 :                 bi = bat_iterator(b);
    2429             : 
    2430         110 :                 if (gdense) {
    2431             :                         /* single element groups */
    2432           6 :                         while (ncand > 0) {
    2433           3 :                                 ncand--;
    2434           3 :                                 i = canditer_next(ci) - b->hseqbase;
    2435           6 :                                 if (!skip_nils ||
    2436           3 :                                     (*atomcmp)(BUNtail(bi, i), nil) != 0) {
    2437           3 :                                         oids[i] = i + b->hseqbase;
    2438           3 :                                         nils--;
    2439             :                                 }
    2440             :                         }
    2441             :                 } else {
    2442             :                         gid = 0; /* in case gids == NULL */
    2443      661966 :                         while (ncand > 0) {
    2444      661859 :                                 ncand--;
    2445      661859 :                                 i = canditer_next(ci) - b->hseqbase;
    2446      661919 :                                 if (gids == NULL ||
    2447          26 :                                     (gids[i] >= min && gids[i] <= max)) {
    2448      661919 :                                         const void *v = BUNtail(bi, i);
    2449      661919 :                                         if (gids)
    2450          28 :                                                 gid = gids[i] - min;
    2451     1325180 :                                         if (!skip_nils ||
    2452      661916 :                                             (*atomcmp)(v, nil) != 0) {
    2453      660838 :                                                 if (is_oid_nil(oids[gid])) {
    2454         104 :                                                         oids[gid] = i + b->hseqbase;
    2455         104 :                                                         nils--;
    2456             :                                                 } else {
    2457      660734 :                                                         const void *g = BUNtail(bi, (BUN) (oids[gid] - b->hseqbase));
    2458     1320230 :                                                         if (t == TYPE_void ||
    2459     1319820 :                                                             ((*atomcmp)(g, nil) != 0 &&
    2460     1318850 :                                                              ((*atomcmp)(v, nil) == 0 ||
    2461      659082 :                                                               GT((*atomcmp)(v, g), 0))))
    2462         895 :                                                                 oids[gid] = i + b->hseqbase;
    2463             :                                                 }
    2464             :                                         }
    2465             :                                 }
    2466             :                         }
    2467             :                 }
    2468             :                 break;
    2469             :         }
    2470             : 
    2471             :         return nils;
    2472             : }
    2473             : 
    2474             : static BAT *
    2475         689 : BATgroupminmax(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils,
    2476             :                bool abort_on_error,
    2477             :                BUN (*minmax)(oid *restrict, BAT *, const oid *restrict, BUN,
    2478             :                              oid, oid, struct canditer *restrict, BUN,
    2479             :                              BUN, bool, bool),
    2480             :                const char *name)
    2481             : {
    2482         689 :         const oid *restrict gids;
    2483         689 :         oid min, max;
    2484         689 :         BUN ngrp;
    2485         689 :         oid *restrict oids;
    2486         689 :         BAT *bn = NULL;
    2487         689 :         BUN nils;
    2488         689 :         struct canditer ci;
    2489         689 :         BUN ncand;
    2490         689 :         const char *err;
    2491         689 :         lng t0 = 0;
    2492             : 
    2493         689 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    2494             : 
    2495         689 :         assert(tp == TYPE_oid);
    2496         689 :         (void) tp;              /* compatibility (with other BATgroup* */
    2497         689 :         (void) abort_on_error;  /* functions) argument */
    2498             : 
    2499         689 :         if (!ATOMlinear(b->ttype)) {
    2500           0 :                 GDKerror("%s: cannot determine minimum on "
    2501             :                          "non-linear type %s\n", name, ATOMname(b->ttype));
    2502           0 :                 return NULL;
    2503             :         }
    2504             : 
    2505         689 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    2506           0 :                 GDKerror("%s: %s\n", name, err);
    2507           0 :                 return NULL;
    2508             :         }
    2509             : 
    2510         689 :         if (BATcount(b) == 0 || ngrp == 0) {
    2511             :                 /* trivial: no minimums, so return bat aligned with g
    2512             :                  * with nil in the tail */
    2513          68 :                 return BATconstant(ngrp == 0 ? 0 : min, TYPE_oid, &oid_nil, ngrp, TRANSIENT);
    2514             :         }
    2515             : 
    2516         621 :         bn = COLnew(min, TYPE_oid, ngrp, TRANSIENT);
    2517         621 :         if (bn == NULL)
    2518             :                 return NULL;
    2519         621 :         oids = (oid *) Tloc(bn, 0);
    2520             : 
    2521         621 :         if (g == NULL || BATtdense(g))
    2522             :                 gids = NULL;
    2523             :         else
    2524         283 :                 gids = (const oid *) Tloc(g, 0);
    2525             : 
    2526         621 :         nils = (*minmax)(oids, b, gids, ngrp, min, max, &ci, ncand,
    2527         621 :                          BATcount(b), skip_nils, g && BATtdense(g));
    2528             : 
    2529         621 :         BATsetcount(bn, ngrp);
    2530             : 
    2531         621 :         bn->tkey = BATcount(bn) <= 1;
    2532         621 :         bn->tsorted = BATcount(bn) <= 1;
    2533         621 :         bn->trevsorted = BATcount(bn) <= 1;
    2534         621 :         bn->tnil = nils != 0;
    2535         621 :         bn->tnonil = nils == 0;
    2536         621 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    2537             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    2538             :                   "; start " OIDFMT ", count " BUNFMT " (%s -- " LLFMT " usec)\n",
    2539             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    2540             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    2541             :                   ci.seq, ncand, name, GDKusec() - t0);
    2542             :         return bn;
    2543             : }
    2544             : 
    2545             : BAT *
    2546         345 : BATgroupmin(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    2547             :             bool skip_nils, bool abort_on_error)
    2548             : {
    2549         345 :         return BATgroupminmax(b, g, e, s, tp, skip_nils, abort_on_error,
    2550             :                               do_groupmin, __func__);
    2551             : }
    2552             : 
    2553             : /* return pointer to smallest non-nil value in b, or pointer to nil if
    2554             :  * there is no such value (no values at all, or only nil) */
    2555             : void *
    2556         706 : BATmin_skipnil(BAT *b, void *aggr, bit skipnil)
    2557             : {
    2558         706 :         PROPrec *prop;
    2559         706 :         const void *res;
    2560         706 :         size_t s;
    2561         706 :         BATiter bi;
    2562         706 :         lng t0 = 0;
    2563             : 
    2564         706 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    2565             : 
    2566         706 :         if (!ATOMlinear(b->ttype)) {
    2567             :                 /* there is no such thing as a smallest value if you
    2568             :                  * can't compare values */
    2569           0 :                 GDKerror("non-linear type");
    2570           0 :                 return NULL;
    2571             :         }
    2572         706 :         if (BATcount(b) == 0) {
    2573         218 :                 res = ATOMnilptr(b->ttype);
    2574         488 :         } else if ((prop = BATgetprop(b, GDK_MIN_VALUE)) != NULL) {
    2575          27 :                 res = VALptr(&prop->v);
    2576             :         } else {
    2577         461 :                 oid pos;
    2578         461 :                 BAT *pb = NULL;
    2579             : 
    2580         461 :                 if (BATcheckorderidx(b) ||
    2581         461 :                     (VIEWtparent(b) &&
    2582         128 :                      (pb = BBPdescriptor(VIEWtparent(b))) != NULL &&
    2583         128 :                      pb->theap.base == b->theap.base &&
    2584          57 :                      BATcount(pb) == BATcount(b) &&
    2585          58 :                      pb->hseqbase == b->hseqbase &&
    2586          29 :                      BATcheckorderidx(pb))) {
    2587           0 :                         const oid *ords = (const oid *) (pb ? pb->torderidx->base : b->torderidx->base) + ORDERIDXOFF;
    2588           0 :                         BUN r;
    2589           0 :                         if (!b->tnonil) {
    2590           0 :                                 r = binsearch(ords, 0, b->ttype, Tloc(b, 0),
    2591           0 :                                               b->tvheap ? b->tvheap->base : NULL,
    2592           0 :                                               b->twidth, 0, BATcount(b),
    2593           0 :                                               ATOMnilptr(b->ttype), 1, 1);
    2594           0 :                                 if (r == 0) {
    2595           0 :                                         b->tnonil = true;
    2596           0 :                                         b->batDirtydesc = true;
    2597             :                                 }
    2598             :                         } else {
    2599             :                                 r = 0;
    2600             :                         }
    2601           0 :                         if (r == BATcount(b)) {
    2602             :                                 /* no non-nil values */
    2603           0 :                                 pos = oid_nil;
    2604             :                         } else {
    2605           0 :                                 pos = ords[r];
    2606             :                         }
    2607         461 :                 } else if ((VIEWtparent(b) == 0 ||
    2608         490 :                             BATcount(b) == BATcount(BBPdescriptor(VIEWtparent(b)))) &&
    2609         362 :                            BATcheckimprints(b)) {
    2610           0 :                         Imprints *imprints = VIEWtparent(b) ? BBPdescriptor(VIEWtparent(b))->timprints : b->timprints;
    2611           0 :                         int i;
    2612             : 
    2613           0 :                         pos = oid_nil;
    2614             :                         /* find first non-empty bin */
    2615           0 :                         for (i = 0; i < imprints->bits; i++) {
    2616           0 :                                 if (imprints->stats[i + 128]) {
    2617           0 :                                         pos = imprints->stats[i] + b->hseqbase;
    2618           0 :                                         break;
    2619             :                                 }
    2620             :                         }
    2621             :                 } else {
    2622         461 :                         struct canditer ci;
    2623         461 :                         BUN ncand = canditer_init(&ci, b, NULL);
    2624         461 :                         (void) do_groupmin(&pos, b, NULL, 1, 0, 0, &ci, ncand,
    2625             :                                            BATcount(b), skipnil, false);
    2626             :                 }
    2627         461 :                 if (is_oid_nil(pos)) {
    2628          37 :                         res = ATOMnilptr(b->ttype);
    2629             :                 } else {
    2630         424 :                         bi = bat_iterator(b);
    2631         424 :                         res = BUNtail(bi, pos - b->hseqbase);
    2632         424 :                         BATsetprop(b, GDK_MIN_VALUE, b->ttype, res);
    2633             :                 }
    2634             :         }
    2635         706 :         if (aggr == NULL) {
    2636         488 :                 s = ATOMlen(b->ttype, res);
    2637         488 :                 aggr = GDKmalloc(s);
    2638             :         } else {
    2639         436 :                 s = ATOMsize(ATOMtype(b->ttype));
    2640             :         }
    2641         706 :         if (aggr != NULL)       /* else: malloc error */
    2642         706 :                 memcpy(aggr, res, s);
    2643         706 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
    2644             :                   ALGOBATPAR(b), skipnil, GDKusec() - t0);
    2645             :         return aggr;
    2646             : }
    2647             : 
    2648             : void *
    2649         477 : BATmin(BAT *b, void *aggr)
    2650             : {
    2651         477 :         return BATmin_skipnil(b, aggr, 1);
    2652             : }
    2653             : 
    2654             : BAT *
    2655         344 : BATgroupmax(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    2656             :             bool skip_nils, bool abort_on_error)
    2657             : {
    2658         344 :         return BATgroupminmax(b, g, e, s, tp, skip_nils, abort_on_error,
    2659             :                               do_groupmax, __func__);
    2660             : }
    2661             : 
    2662             : void *
    2663        1037 : BATmax_skipnil(BAT *b, void *aggr, bit skipnil)
    2664             : {
    2665        1037 :         PROPrec *prop;
    2666        1037 :         const void *res;
    2667        1037 :         size_t s;
    2668        1037 :         BATiter bi;
    2669        1037 :         lng t0 = 0;
    2670             : 
    2671        1037 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    2672             : 
    2673        1037 :         if (!ATOMlinear(b->ttype)) {
    2674           0 :                 GDKerror("non-linear type");
    2675           0 :                 return NULL;
    2676             :         }
    2677        1037 :         if (BATcount(b) == 0) {
    2678         216 :                 res = ATOMnilptr(b->ttype);
    2679         821 :         } else if ((prop = BATgetprop(b, GDK_MAX_VALUE)) != NULL) {
    2680         349 :                 res = VALptr(&prop->v);
    2681             :         } else {
    2682         472 :                 oid pos;
    2683         472 :                 BAT *pb = NULL;
    2684             : 
    2685         472 :                 if (BATcheckorderidx(b) ||
    2686         472 :                     (VIEWtparent(b) &&
    2687         134 :                      (pb = BBPdescriptor(VIEWtparent(b))) != NULL &&
    2688         134 :                      pb->theap.base == b->theap.base &&
    2689          68 :                      BATcount(pb) == BATcount(b) &&
    2690          82 :                      pb->hseqbase == b->hseqbase &&
    2691          41 :                      BATcheckorderidx(pb))) {
    2692           1 :                         const oid *ords = (const oid *) (pb ? pb->torderidx->base : b->torderidx->base) + ORDERIDXOFF;
    2693             : 
    2694           1 :                         pos = ords[BATcount(b) - 1];
    2695             :                         /* nils are first, ie !skipnil, check for nils */
    2696           1 :                         if (!skipnil) {
    2697           0 :                                 BUN z = ords[0];
    2698             : 
    2699           0 :                                 bi = bat_iterator(b);
    2700           0 :                                 res = BUNtail(bi, z - b->hseqbase);
    2701             : 
    2702           0 :                                 if (ATOMcmp(b->ttype, res, ATOMnilptr(b->ttype)) == 0)
    2703           0 :                                         pos = z;
    2704             :                         }
    2705         472 :                 } else if ((VIEWtparent(b) == 0 ||
    2706         513 :                             BATcount(b) == BATcount(BBPdescriptor(VIEWtparent(b)))) &&
    2707         379 :                            BATcheckimprints(b)) {
    2708           0 :                         Imprints *imprints = VIEWtparent(b) ? BBPdescriptor(VIEWtparent(b))->timprints : b->timprints;
    2709           0 :                         int i;
    2710             : 
    2711           0 :                         pos = oid_nil;
    2712             :                         /* find last non-empty bin */
    2713           0 :                         for (i = imprints->bits - 1; i >= 0; i--) {
    2714           0 :                                 if (imprints->stats[i + 128]) {
    2715           0 :                                         pos = imprints->stats[i + 64] + b->hseqbase;
    2716           0 :                                         break;
    2717             :                                 }
    2718             :                         }
    2719             :                 } else {
    2720         472 :                         struct canditer ci;
    2721         472 :                         BUN ncand = canditer_init(&ci, b, NULL);
    2722         472 :                         (void) do_groupmax(&pos, b, NULL, 1, 0, 0, &ci, ncand,
    2723             :                                            BATcount(b), skipnil, false);
    2724             :                 }
    2725         472 :                 if (is_oid_nil(pos)) {
    2726          34 :                         res = ATOMnilptr(b->ttype);
    2727             :                 } else {
    2728         438 :                         bi = bat_iterator(b);
    2729         438 :                         res = BUNtail(bi, pos - b->hseqbase);
    2730         438 :                         if (b->tnonil)
    2731         384 :                                 BATsetprop(b, GDK_MAX_VALUE, b->ttype, res);
    2732             :                 }
    2733             :         }
    2734        1037 :         if (aggr == NULL) {
    2735         483 :                 s = ATOMlen(b->ttype, res);
    2736         483 :                 aggr = GDKmalloc(s);
    2737             :         } else {
    2738        1108 :                 s = ATOMsize(ATOMtype(b->ttype));
    2739             :         }
    2740        1037 :         if (aggr != NULL)       /* else: malloc error */
    2741        1037 :                 memcpy(aggr, res, s);
    2742        1037 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
    2743             :                   ALGOBATPAR(b), skipnil, GDKusec() - t0);
    2744             :         return aggr;
    2745             : }
    2746             : 
    2747             : void *
    2748         789 : BATmax(BAT *b, void *aggr)
    2749             : {
    2750         789 :         return BATmax_skipnil(b, aggr, 1);
    2751             : }
    2752             : 
    2753             : 
    2754             : /* ---------------------------------------------------------------------- */
    2755             : /* quantiles/median */
    2756             : 
    2757             : #if SIZEOF_OID == SIZEOF_INT
    2758             : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_int(indir, offset, (const int *) vals, lo, hi, (int) (v), ordering, last)
    2759             : #endif
    2760             : #if SIZEOF_OID == SIZEOF_LNG
    2761             : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_lng(indir, offset, (const lng *) vals, lo, hi, (lng) (v), ordering, last)
    2762             : #endif
    2763             : 
    2764             : #define DO_QUANTILE_AVG(TPE)                                            \
    2765             :         do {                                                            \
    2766             :                 TPE low = *(TPE*) BUNtail(bi, r + (BUN) hi);            \
    2767             :                 TPE high = *(TPE*) BUNtail(bi, r + (BUN) lo);           \
    2768             :                 if (is_##TPE##_nil(low) || is_##TPE##_nil(high)) {      \
    2769             :                         val = dbl_nil;                                  \
    2770             :                         nils++;                                         \
    2771             :                 } else {                                                \
    2772             :                         val = (f - lo) * low + (lo + 1 - f) * high;     \
    2773             :                 }                                                       \
    2774             :         } while (0)
    2775             : 
    2776             : static BAT *
    2777          61 : doBATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
    2778             :                    bool skip_nils, bool abort_on_error, bool average)
    2779             : {
    2780          61 :         bool freeb = false, freeg = false;
    2781          61 :         oid min, max;
    2782          61 :         BUN ngrp;
    2783          61 :         BUN nils = 0;
    2784          61 :         BAT *bn = NULL;
    2785          61 :         struct canditer ci;
    2786          61 :         BUN ncand;
    2787          61 :         BAT *t1, *t2;
    2788          61 :         BATiter bi;
    2789          61 :         const void *v;
    2790          61 :         const void *nil = ATOMnilptr(tp);
    2791          61 :         const void *dnil = nil;
    2792          61 :         dbl val;                /* only used for average */
    2793          61 :         int (*atomcmp)(const void *, const void *) = ATOMcompare(tp);
    2794          61 :         const char *err;
    2795          61 :         (void) abort_on_error;
    2796          61 :         lng t0 = 0;
    2797             : 
    2798          61 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    2799             : 
    2800          61 :         if (average) {
    2801          16 :                 switch (ATOMbasetype(b->ttype)) {
    2802             :                 case TYPE_bte:
    2803             :                 case TYPE_sht:
    2804             :                 case TYPE_int:
    2805             :                 case TYPE_lng:
    2806             : #ifdef HAVE_HGE
    2807             :                 case TYPE_hge:
    2808             : #endif
    2809             :                 case TYPE_flt:
    2810             :                 case TYPE_dbl:
    2811             :                         break;
    2812           0 :                 default:
    2813           0 :                         GDKerror("incompatible type\n");
    2814           0 :                         return NULL;
    2815             :                 }
    2816             :                 dnil = &dbl_nil;
    2817             :         }
    2818          61 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    2819           0 :                 GDKerror("%s\n", err);
    2820           0 :                 return NULL;
    2821             :         }
    2822          61 :         assert(tp == b->ttype);
    2823          61 :         if (!ATOMlinear(tp)) {
    2824           0 :                 GDKerror("cannot determine quantile on "
    2825             :                          "non-linear type %s\n", ATOMname(tp));
    2826           0 :                 return NULL;
    2827             :         }
    2828          61 :         if (quantile < 0 || quantile > 1) {
    2829           0 :                 GDKerror("cannot determine quantile for "
    2830             :                          "p=%f (p has to be in [0,1])\n", quantile);
    2831           0 :                 return NULL;
    2832             :         }
    2833             : 
    2834          61 :         if (BATcount(b) == 0 || ngrp == 0 || is_dbl_nil(quantile)) {
    2835             :                 /* trivial: no values, thus also no quantiles,
    2836             :                  * so return bat aligned with e with nil in the tail
    2837             :                  * The same happens for a NULL quantile */
    2838          22 :                 return BATconstant(ngrp == 0 ? 0 : min, average ? TYPE_dbl : tp, dnil, ngrp, TRANSIENT);
    2839             :         }
    2840             : 
    2841          50 :         if (s) {
    2842             :                 /* there is a candidate list, replace b (and g, if
    2843             :                  * given) with just the values we're interested in */
    2844           0 :                 b = BATproject(s, b);
    2845           0 :                 if (b == NULL)
    2846             :                         return NULL;
    2847           0 :                 freeb = true;
    2848           0 :                 if (g) {
    2849           0 :                         g = BATproject(s, g);
    2850           0 :                         if (g == NULL)
    2851           0 :                                 goto bunins_failed;
    2852             :                         freeg = true;
    2853             :                 }
    2854             :         }
    2855             : 
    2856             :         /* we want to sort b so that we can figure out the quantile, but
    2857             :          * if g is given, sort g and subsort b so that we can get the
    2858             :          * quantile for each group */
    2859          50 :         if (g) {
    2860          16 :                 const oid *restrict grps;
    2861          16 :                 oid prev;
    2862          16 :                 BUN p, q, r;
    2863             : 
    2864          16 :                 if (BATtdense(g)) {
    2865             :                         /* singleton groups, so calculating quantile is
    2866             :                          * easy */
    2867           1 :                         if (average)
    2868           0 :                                 bn = BATconvert(b, NULL, NULL, TYPE_dbl, abort_on_error);
    2869             :                         else
    2870           1 :                                 bn = COLcopy(b, tp, false, TRANSIENT);
    2871           1 :                         BAThseqbase(bn, g->tseqbase); /* deals with NULL */
    2872           1 :                         if (freeb)
    2873           0 :                                 BBPunfix(b->batCacheid);
    2874           1 :                         if (freeg)
    2875           0 :                                 BBPunfix(g->batCacheid);
    2876           1 :                         return bn;
    2877             :                 }
    2878          15 :                 if (BATsort(&t1, &t2, NULL, g, NULL, NULL, false, false, false) != GDK_SUCCEED)
    2879           0 :                         goto bunins_failed;
    2880          15 :                 if (freeg)
    2881           0 :                         BBPunfix(g->batCacheid);
    2882          15 :                 g = t1;
    2883          15 :                 freeg = true;
    2884             : 
    2885          15 :                 if (BATsort(&t1, NULL, NULL, b, t2, g, false, false, false) != GDK_SUCCEED) {
    2886           0 :                         BBPunfix(t2->batCacheid);
    2887           0 :                         goto bunins_failed;
    2888             :                 }
    2889          15 :                 if (freeb)
    2890           0 :                         BBPunfix(b->batCacheid);
    2891          15 :                 b = t1;
    2892          15 :                 freeb = true;
    2893          15 :                 BBPunfix(t2->batCacheid);
    2894             : 
    2895          15 :                 if (average)
    2896           2 :                         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    2897             :                 else
    2898          13 :                         bn = COLnew(min, tp, ngrp, TRANSIENT);
    2899          15 :                 if (bn == NULL)
    2900           0 :                         goto bunins_failed;
    2901             : 
    2902          15 :                 bi = bat_iterator(b);
    2903             : 
    2904          15 :                 grps = (const oid *) Tloc(g, 0);
    2905             :                 /* for each group (r and p are the beginning and end
    2906             :                  * of the current group, respectively) */
    2907          62 :                 for (r = 0, q = BATcount(g); r < q; r = p) {
    2908          47 :                         BUN qindex;
    2909          47 :                         prev = grps[r];
    2910             :                         /* search for end of current group (grps is
    2911             :                          * sorted so we can use binary search) */
    2912          47 :                         p = binsearch_oid(NULL, 0, grps, r, q - 1, prev, 1, 1);
    2913          47 :                         if (skip_nils && !b->tnonil) {
    2914             :                                 /* within group, locate start of non-nils */
    2915          15 :                                 r = binsearch(NULL, 0, tp, Tloc(b, 0),
    2916          15 :                                               b->tvheap ? b->tvheap->base : NULL,
    2917          15 :                                               b->twidth, r, p, nil,
    2918             :                                               1, 1);
    2919             :                         }
    2920          47 :                         if (r == p) {
    2921             :                                 /* no non-nils found */
    2922           5 :                                 v = dnil;
    2923           5 :                                 nils++;
    2924          42 :                         } else if (average) {
    2925           4 :                                 double f = (p - r - 1) * quantile;
    2926           4 :                                 double lo = floor(f);
    2927           4 :                                 double hi = ceil(f);
    2928           4 :                                 switch (ATOMbasetype(tp)) {
    2929           4 :                                 case TYPE_bte:
    2930           4 :                                         DO_QUANTILE_AVG(bte);
    2931             :                                         break;
    2932           0 :                                 case TYPE_sht:
    2933           0 :                                         DO_QUANTILE_AVG(sht);
    2934             :                                         break;
    2935           0 :                                 case TYPE_int:
    2936           0 :                                         DO_QUANTILE_AVG(int);
    2937             :                                         break;
    2938           0 :                                 case TYPE_lng:
    2939           0 :                                         DO_QUANTILE_AVG(lng);
    2940             :                                         break;
    2941             : #ifdef HAVE_HGE
    2942           0 :                                 case TYPE_hge:
    2943           0 :                                         DO_QUANTILE_AVG(hge);
    2944             :                                         break;
    2945             : #endif
    2946           0 :                                 case TYPE_flt:
    2947           0 :                                         DO_QUANTILE_AVG(flt);
    2948             :                                         break;
    2949           0 :                                 case TYPE_dbl:
    2950           0 :                                         DO_QUANTILE_AVG(dbl);
    2951             :                                         break;
    2952             :                                 }
    2953             :                                 v = &val;
    2954             :                         } else {
    2955             :                                 /* round *down* to nearest integer */
    2956          38 :                                 double f = (p - r - 1) * quantile;
    2957          38 :                                 qindex = r + p - (BUN) (p + 0.5 - f);
    2958             :                                 /* be a little paranoid about the index */
    2959          38 :                                 assert(qindex >= r && qindex <  p);
    2960          38 :                                 v = BUNtail(bi, qindex);
    2961          38 :                                 if (!skip_nils && !b->tnonil)
    2962           0 :                                         nils += (*atomcmp)(v, dnil) == 0;
    2963             :                         }
    2964          94 :                         if (bunfastapp_nocheck(bn, BUNlast(bn), v, Tsize(bn)) != GDK_SUCCEED)
    2965           0 :                                 goto bunins_failed;
    2966             :                 }
    2967          15 :                 nils += ngrp - BATcount(bn);
    2968          15 :                 while (BATcount(bn) < ngrp) {
    2969           0 :                         if (bunfastapp_nocheck(bn, BUNlast(bn), dnil, Tsize(bn)) != GDK_SUCCEED)
    2970           0 :                                 goto bunins_failed;
    2971             :                 }
    2972          15 :                 bn->theap.dirty = true;
    2973          15 :                 BBPunfix(g->batCacheid);
    2974             :         } else {
    2975          34 :                 BUN index, r, p = BATcount(b);
    2976          34 :                 BAT *pb = NULL;
    2977          34 :                 const oid *ords;
    2978             : 
    2979          62 :                 bn = COLnew(0, average ? TYPE_dbl : tp, 1, TRANSIENT);
    2980          34 :                 if (bn == NULL)
    2981           0 :                         goto bunins_failed;
    2982             : 
    2983          34 :                 t1 = NULL;
    2984             : 
    2985          34 :                 if (BATcheckorderidx(b) ||
    2986          34 :                     (VIEWtparent(b) &&
    2987           4 :                      (pb = BBPdescriptor(VIEWtparent(b))) != NULL &&
    2988           4 :                      pb->theap.base == b->theap.base &&
    2989           4 :                      BATcount(pb) == BATcount(b) &&
    2990           8 :                      pb->hseqbase == b->hseqbase &&
    2991           4 :                      BATcheckorderidx(pb))) {
    2992           0 :                         ords = (const oid *) (pb ? pb->torderidx->base : b->torderidx->base) + ORDERIDXOFF;
    2993             :                 } else {
    2994          34 :                         if (BATsort(NULL, &t1, NULL, b, NULL, g, false, false, false) != GDK_SUCCEED)
    2995           0 :                                 goto bunins_failed;
    2996          34 :                         if (BATtdense(t1))
    2997             :                                 ords = NULL;
    2998             :                         else
    2999          21 :                                 ords = (const oid *) Tloc(t1, 0);
    3000             :                 }
    3001             : 
    3002          34 :                 if (skip_nils && !b->tnonil)
    3003           6 :                         r = binsearch(ords, 0, tp, Tloc(b, 0),
    3004           6 :                                       b->tvheap ? b->tvheap->base : NULL,
    3005           6 :                                       b->twidth, 0, p,
    3006             :                                       nil, 1, 1);
    3007             :                 else
    3008             :                         r = 0;
    3009             : 
    3010          34 :                 bi = bat_iterator(b);
    3011          34 :                 if (r == p) {
    3012             :                         /* no non-nil values, so quantile is nil */
    3013             :                         v = dnil;
    3014             :                         nils++;
    3015          31 :                 } else if (average) {
    3016           6 :                         double f = (p - r - 1) * quantile;
    3017           6 :                         double lo = floor(f);
    3018           6 :                         double hi = ceil(f);
    3019           6 :                         switch (ATOMbasetype(tp)) {
    3020           5 :                         case TYPE_bte:
    3021           5 :                                 DO_QUANTILE_AVG(bte);
    3022             :                                 break;
    3023           1 :                         case TYPE_sht:
    3024           1 :                                 DO_QUANTILE_AVG(sht);
    3025             :                                 break;
    3026           0 :                         case TYPE_int:
    3027           0 :                                 DO_QUANTILE_AVG(int);
    3028             :                                 break;
    3029           0 :                         case TYPE_lng:
    3030           0 :                                 DO_QUANTILE_AVG(lng);
    3031             :                                 break;
    3032             : #ifdef HAVE_HGE
    3033           0 :                         case TYPE_hge:
    3034           0 :                                 DO_QUANTILE_AVG(hge);
    3035             :                                 break;
    3036             : #endif
    3037           0 :                         case TYPE_flt:
    3038           0 :                                 DO_QUANTILE_AVG(flt);
    3039             :                                 break;
    3040           0 :                         case TYPE_dbl:
    3041           0 :                                 DO_QUANTILE_AVG(dbl);
    3042             :                                 break;
    3043             :                         }
    3044             :                         v = &val;
    3045             :                 } else {
    3046          25 :                         double f;
    3047             :                         /* round (p-r-1)*quantile *down* to nearest
    3048             :                          * integer (i.e., 1.49 and 1.5 are rounded to
    3049             :                          * 1, 1.51 is rounded to 2) */
    3050          25 :                         f = (p - r - 1) * quantile;
    3051          25 :                         index = r + p - (BUN) (p + 0.5 - f);
    3052          25 :                         if (ords)
    3053          19 :                                 index = ords[index] - b->hseqbase;
    3054             :                         else
    3055           6 :                                 index = index + t1->tseqbase;
    3056          25 :                         v = BUNtail(bi, index);
    3057          25 :                         nils += (*atomcmp)(v, dnil) == 0;
    3058             :                 }
    3059          34 :                 if (t1)
    3060          34 :                         BBPunfix(t1->batCacheid);
    3061          34 :                 if (BUNappend(bn, v, false) != GDK_SUCCEED)
    3062           0 :                         goto bunins_failed;
    3063             :         }
    3064             : 
    3065          49 :         if (freeb)
    3066          15 :                 BBPunfix(b->batCacheid);
    3067             : 
    3068          49 :         bn->tkey = BATcount(bn) <= 1;
    3069          49 :         bn->tsorted = BATcount(bn) <= 1;
    3070          49 :         bn->trevsorted = BATcount(bn) <= 1;
    3071          49 :         bn->tnil = nils != 0;
    3072          49 :         bn->tnonil = nils == 0;
    3073          49 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    3074             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
    3075             :                   ",quantile=%g,average=%s -> " ALGOOPTBATFMT
    3076             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    3077             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    3078             :                   ALGOOPTBATPAR(s), quantile, average ? "true" : "false",
    3079             :                   ALGOOPTBATPAR(bn), ci.seq, ncand, GDKusec() - t0);
    3080             :         return bn;
    3081             : 
    3082           0 :   bunins_failed:
    3083           0 :         if (freeb)
    3084           0 :                 BBPunfix(b->batCacheid);
    3085           0 :         if (freeg)
    3086           0 :                 BBPunfix(g->batCacheid);
    3087           0 :         if (bn)
    3088           0 :                 BBPunfix(bn->batCacheid);
    3089             :         return NULL;
    3090             : }
    3091             : 
    3092             : BAT *
    3093          24 : BATgroupmedian(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    3094             :                bool skip_nils, bool abort_on_error)
    3095             : {
    3096          24 :         return doBATgroupquantile(b, g, e, s, tp, 0.5,
    3097             :                                   skip_nils, abort_on_error, false);
    3098             : }
    3099             : 
    3100             : BAT *
    3101          29 : BATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
    3102             :                  bool skip_nils, bool abort_on_error)
    3103             : {
    3104          29 :         return doBATgroupquantile(b, g, e, s, tp, quantile,
    3105             :                                   skip_nils, abort_on_error, false);
    3106             : }
    3107             : 
    3108             : BAT *
    3109           5 : BATgroupmedian_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    3110             :                    bool skip_nils, bool abort_on_error)
    3111             : {
    3112           5 :         return doBATgroupquantile(b, g, e, s, tp, 0.5,
    3113             :                                   skip_nils, abort_on_error, true);
    3114             : }
    3115             : 
    3116             : BAT *
    3117           3 : BATgroupquantile_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
    3118             :                      bool skip_nils, bool abort_on_error)
    3119             : {
    3120           3 :         return doBATgroupquantile(b, g, e, s, tp, quantile,
    3121             :                                   skip_nils, abort_on_error, true);
    3122             : }
    3123             : 
    3124             : /* ---------------------------------------------------------------------- */
    3125             : /* standard deviation (both biased and non-biased) */
    3126             : 
    3127             : #define AGGR_STDEV_SINGLE(TYPE)                         \
    3128             :         do {                                            \
    3129             :                 TYPE x;                                 \
    3130             :                 for (i = 0; i < cnt; i++) {          \
    3131             :                         x = ((const TYPE *) values)[i]; \
    3132             :                         if (is_##TYPE##_nil(x))         \
    3133             :                                 continue;               \
    3134             :                         n++;                            \
    3135             :                         delta = (dbl) x - mean;         \
    3136             :                         mean += delta / n;              \
    3137             :                         m2 += delta * ((dbl) x - mean); \
    3138             :                         if (isinf(m2))                  \
    3139             :                                 goto overflow;          \
    3140             :                 }                                       \
    3141             :         } while (0)
    3142             : 
    3143             : static dbl
    3144          24 : calcvariance(dbl *restrict avgp, const void *restrict values, BUN cnt, int tp, bool issample)
    3145             : {
    3146          24 :         BUN n = 0, i;
    3147          24 :         dbl mean = 0;
    3148          24 :         dbl m2 = 0;
    3149          24 :         dbl delta;
    3150             : 
    3151          24 :         switch (tp) {
    3152             :         case TYPE_bte:
    3153          14 :                 AGGR_STDEV_SINGLE(bte);
    3154             :                 break;
    3155             :         case TYPE_sht:
    3156          14 :                 AGGR_STDEV_SINGLE(sht);
    3157             :                 break;
    3158             :         case TYPE_int:
    3159          33 :                 AGGR_STDEV_SINGLE(int);
    3160             :                 break;
    3161             :         case TYPE_lng:
    3162          28 :                 AGGR_STDEV_SINGLE(lng);
    3163             :                 break;
    3164             : #ifdef HAVE_HGE
    3165             :         case TYPE_hge:
    3166           0 :                 AGGR_STDEV_SINGLE(hge);
    3167             :                 break;
    3168             : #endif
    3169             :         case TYPE_flt:
    3170           0 :                 AGGR_STDEV_SINGLE(flt);
    3171             :                 break;
    3172             :         case TYPE_dbl:
    3173          13 :                 AGGR_STDEV_SINGLE(dbl);
    3174             :                 break;
    3175           0 :         default:
    3176           0 :                 GDKerror("type (%s) not supported.\n", ATOMname(tp));
    3177           0 :                 return dbl_nil;
    3178             :         }
    3179          22 :         if (n <= (BUN) issample) {
    3180           7 :                 if (avgp)
    3181           0 :                         *avgp = dbl_nil;
    3182           7 :                 return dbl_nil;
    3183             :         }
    3184          15 :         if (avgp)
    3185           0 :                 *avgp = mean;
    3186          15 :         return m2 / (n - issample);
    3187           2 :   overflow:
    3188           2 :         GDKerror("22003!overflow in calculation.\n");
    3189           2 :         return dbl_nil;
    3190             : }
    3191             : 
    3192             : dbl
    3193          13 : BATcalcstdev_population(dbl *avgp, BAT *b)
    3194             : {
    3195          13 :         lng t0 = 0;
    3196             : 
    3197          13 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3198          26 :         dbl v = calcvariance(avgp, (const void *) Tloc(b, 0),
    3199          13 :                              BATcount(b), b->ttype, false);
    3200          13 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
    3201             :                   ALGOBATPAR(b), GDKusec() - t0);
    3202          13 :         return is_dbl_nil(v) ? dbl_nil : sqrt(v);
    3203             : }
    3204             : 
    3205             : dbl
    3206           7 : BATcalcstdev_sample(dbl *avgp, BAT *b)
    3207             : {
    3208           7 :         lng t0 = 0;
    3209             : 
    3210           7 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3211          14 :         dbl v = calcvariance(avgp, (const void *) Tloc(b, 0),
    3212           7 :                              BATcount(b), b->ttype, true);
    3213           7 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
    3214             :                   ALGOBATPAR(b), GDKusec() - t0);
    3215           7 :         return is_dbl_nil(v) ? dbl_nil : sqrt(v);
    3216             : }
    3217             : 
    3218             : dbl
    3219           2 : BATcalcvariance_population(dbl *avgp, BAT *b)
    3220             : {
    3221           2 :         lng t0 = 0;
    3222             : 
    3223           2 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3224           4 :         dbl v = calcvariance(avgp, (const void *) Tloc(b, 0),
    3225           2 :                              BATcount(b), b->ttype, false);
    3226           2 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
    3227             :                   ALGOBATPAR(b), GDKusec() - t0);
    3228           2 :         return v;
    3229             : }
    3230             : 
    3231             : dbl
    3232           2 : BATcalcvariance_sample(dbl *avgp, BAT *b)
    3233             : {
    3234           2 :         lng t0 = 0;
    3235             : 
    3236           2 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3237           4 :         dbl v = calcvariance(avgp, (const void *) Tloc(b, 0),
    3238           2 :                              BATcount(b), b->ttype, true);
    3239           2 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
    3240             :                   ALGOBATPAR(b), GDKusec() - t0);
    3241           2 :         return v;
    3242             : }
    3243             : 
    3244             : #define AGGR_COVARIANCE_SINGLE(TYPE)                                    \
    3245             :         do {                                                            \
    3246             :                 TYPE x, y;                                              \
    3247             :                 for (i = 0; i < cnt; i++) {                          \
    3248             :                         x = ((const TYPE *) v1)[i];                     \
    3249             :                         y = ((const TYPE *) v2)[i];                     \
    3250             :                         if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y))   \
    3251             :                                 continue;                               \
    3252             :                         n++;                                            \
    3253             :                         delta1 = (dbl) x - mean1;                       \
    3254             :                         mean1 += delta1 / n;                            \
    3255             :                         delta2 = (dbl) y - mean2;                       \
    3256             :                         mean2 += delta2 / n;                            \
    3257             :                         m2 += delta1 * ((dbl) y - mean2);               \
    3258             :                         if (isinf(m2))                  \
    3259             :                                 goto overflow;          \
    3260             :                 }                                                       \
    3261             :         } while (0)
    3262             : 
    3263             : static dbl
    3264          10 : calccovariance(const void *v1, const void *v2, BUN cnt, int tp, bool issample)
    3265             : {
    3266          10 :         BUN n = 0, i;
    3267          10 :         dbl mean1 = 0, mean2 = 0, m2 = 0, delta1, delta2;
    3268             : 
    3269          10 :         switch (tp) {
    3270             :         case TYPE_bte:
    3271           0 :                 AGGR_COVARIANCE_SINGLE(bte);
    3272             :                 break;
    3273             :         case TYPE_sht:
    3274           0 :                 AGGR_COVARIANCE_SINGLE(sht);
    3275             :                 break;
    3276             :         case TYPE_int:
    3277          83 :                 AGGR_COVARIANCE_SINGLE(int);
    3278             :                 break;
    3279             :         case TYPE_lng:
    3280           0 :                 AGGR_COVARIANCE_SINGLE(lng);
    3281             :                 break;
    3282             : #ifdef HAVE_HGE
    3283             :         case TYPE_hge:
    3284           0 :                 AGGR_COVARIANCE_SINGLE(hge);
    3285             :                 break;
    3286             : #endif
    3287             :         case TYPE_flt:
    3288           0 :                 AGGR_COVARIANCE_SINGLE(flt);
    3289             :                 break;
    3290             :         case TYPE_dbl:
    3291           2 :                 AGGR_COVARIANCE_SINGLE(dbl);
    3292             :                 break;
    3293           0 :         default:
    3294           0 :                 GDKerror("type (%s) not supported.\n", ATOMname(tp));
    3295           0 :                 return dbl_nil;
    3296             :         }
    3297           9 :         if (n <= (BUN) issample)
    3298           0 :                 return dbl_nil;
    3299           9 :         return m2 / (n - issample);
    3300           1 :   overflow:
    3301           1 :         GDKerror("22003!overflow in calculation.\n");
    3302           1 :         return dbl_nil;
    3303             : }
    3304             : 
    3305             : dbl
    3306           6 : BATcalccovariance_population(BAT *b1, BAT *b2)
    3307             : {
    3308           6 :         lng t0 = 0;
    3309             : 
    3310           6 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3311          12 :         dbl v = calccovariance(Tloc(b1, 0), Tloc(b2, 0),
    3312           6 :                                BATcount(b1), b1->ttype, false);
    3313           6 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
    3314             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
    3315           6 :         return v;
    3316             : }
    3317             : 
    3318             : dbl
    3319           4 : BATcalccovariance_sample(BAT *b1, BAT *b2)
    3320             : {
    3321           4 :         lng t0 = 0;
    3322             : 
    3323           4 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3324           8 :         dbl v = calccovariance(Tloc(b1, 0), Tloc(b2, 0),
    3325           4 :                               BATcount(b1), b1->ttype, true);
    3326           4 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
    3327             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
    3328           4 :         return v;
    3329             : }
    3330             : 
    3331             : #define AGGR_CORRELATION_SINGLE(TYPE)                                   \
    3332             :         do {                                                            \
    3333             :                 TYPE x, y;                                              \
    3334             :                 for (i = 0; i < cnt; i++) {                          \
    3335             :                         x = ((const TYPE *) v1)[i];                     \
    3336             :                         y = ((const TYPE *) v2)[i];                     \
    3337             :                         if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y))   \
    3338             :                                 continue;                               \
    3339             :                         n++;                                            \
    3340             :                         delta1 = (dbl) x - mean1;                       \
    3341             :                         mean1 += delta1 / n;                            \
    3342             :                         delta2 = (dbl) y - mean2;                       \
    3343             :                         mean2 += delta2 / n;                            \
    3344             :                         aux = (dbl) y - mean2;                          \
    3345             :                         up += delta1 * aux;                             \
    3346             :                         down1 += delta1 * ((dbl) x - mean1);            \
    3347             :                         down2 += delta2 * aux;                          \
    3348             :                         if (isinf(up) || isinf(down1) || isinf(down2))          \
    3349             :                                 goto overflow;          \
    3350             :                 }                                                       \
    3351             :         } while (0)
    3352             : 
    3353             : dbl
    3354          17 : BATcalccorrelation(BAT *b1, BAT *b2)
    3355             : {
    3356          17 :         BUN n = 0, i, cnt = BATcount(b1);
    3357          17 :         dbl mean1 = 0, mean2 = 0, up = 0, down1 = 0, down2 = 0, delta1, delta2, aux;
    3358          17 :         const void *v1 = (const void *) Tloc(b1, 0), *v2 = (const void *) Tloc(b2, 0);
    3359          17 :         int tp = b1->ttype;
    3360          17 :         lng t0 = 0;
    3361             : 
    3362          17 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3363             : 
    3364          17 :         switch (tp) {
    3365             :         case TYPE_bte:
    3366           0 :                 AGGR_CORRELATION_SINGLE(bte);
    3367             :                 break;
    3368             :         case TYPE_sht:
    3369           0 :                 AGGR_CORRELATION_SINGLE(sht);
    3370             :                 break;
    3371             :         case TYPE_int:
    3372          77 :                 AGGR_CORRELATION_SINGLE(int);
    3373             :                 break;
    3374             :         case TYPE_lng:
    3375           0 :                 AGGR_CORRELATION_SINGLE(lng);
    3376             :                 break;
    3377             : #ifdef HAVE_HGE
    3378             :         case TYPE_hge:
    3379           0 :                 AGGR_CORRELATION_SINGLE(hge);
    3380             :                 break;
    3381             : #endif
    3382             :         case TYPE_flt:
    3383           0 :                 AGGR_CORRELATION_SINGLE(flt);
    3384             :                 break;
    3385             :         case TYPE_dbl:
    3386          11 :                 AGGR_CORRELATION_SINGLE(dbl);
    3387             :                 break;
    3388           0 :         default:
    3389           0 :                 GDKerror("type (%s) not supported.\n",
    3390             :                          ATOMname(tp));
    3391           0 :                 return dbl_nil;
    3392             :         }
    3393          16 :         if (n != 0 && down1 != 0 && down2 != 0)
    3394           9 :                 aux = (up / n) / (sqrt(down1 / n) * sqrt(down2 / n));
    3395             :         else
    3396           7 :                 aux = dbl_nil;
    3397          16 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
    3398             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
    3399             :         return aux;
    3400           1 :   overflow:
    3401           1 :         GDKerror("22003!overflow in calculation.\n");
    3402           1 :         return dbl_nil;
    3403             : }
    3404             : 
    3405             : #define AGGR_STDEV(TYPE)                                                \
    3406             :         do {                                                            \
    3407             :                 const TYPE *restrict vals = (const TYPE *) Tloc(b, 0);  \
    3408             :                 while (ncand > 0) {                                  \
    3409             :                         ncand--;                                        \
    3410             :                         i = canditer_next(&ci) - b->hseqbase;            \
    3411             :                         if (gids == NULL ||                             \
    3412             :                             (gids[i] >= min && gids[i] <= max)) { \
    3413             :                                 if (gids)                               \
    3414             :                                         gid = gids[i] - min;            \
    3415             :                                 else                                    \
    3416             :                                         gid = (oid) i;                  \
    3417             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    3418             :                                         if (!skip_nils)                 \
    3419             :                                                 cnts[gid] = BUN_NONE;   \
    3420             :                                 } else if (cnts[gid] != BUN_NONE) {     \
    3421             :                                         cnts[gid]++;                    \
    3422             :                                         delta[gid] = (dbl) vals[i] - mean[gid]; \
    3423             :                                         mean[gid] += delta[gid] / cnts[gid]; \
    3424             :                                         m2[gid] += delta[gid] * ((dbl) vals[i] - mean[gid]); \
    3425             :                                 }                                       \
    3426             :                         }                                               \
    3427             :                 }                                                       \
    3428             :                 for (i = 0; i < ngrp; i++) {                         \
    3429             :                         if (cnts[i] == 0 || cnts[i] == BUN_NONE) {      \
    3430             :                                 dbls[i] = dbl_nil;                      \
    3431             :                                 mean[i] = dbl_nil;                      \
    3432             :                                 nils++;                                 \
    3433             :                         } else if (cnts[i] == 1) {                      \
    3434             :                                 dbls[i] = issample ? dbl_nil : 0;       \
    3435             :                                 nils2++;                                \
    3436             :                         } else if (isinf(m2[i])) {                      \
    3437             :                                 goto overflow;          \
    3438             :                         } else if (variance) {                          \
    3439             :                                 dbls[i] = m2[i] / (cnts[i] - issample); \
    3440             :                         } else {                                        \
    3441             :                                 dbls[i] = sqrt(m2[i] / (cnts[i] - issample)); \
    3442             :                         }                                               \
    3443             :                 }                                                       \
    3444             :         } while (0)
    3445             : 
    3446             : /* Calculate group standard deviation (population (i.e. biased) or
    3447             :  * sample (i.e. non-biased)) with optional candidates list.
    3448             :  *
    3449             :  * Note that this helper function is prepared to return two BATs: one
    3450             :  * (as return value) with the standard deviation per group, and one
    3451             :  * (as return argument) with the average per group.  This isn't
    3452             :  * currently used since it doesn't fit into the mold of grouped
    3453             :  * aggregates. */
    3454             : static BAT *
    3455          16 : dogroupstdev(BAT **avgb, BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    3456             :              bool skip_nils, bool issample, bool variance, const char *func)
    3457             : {
    3458          16 :         const oid *restrict gids;
    3459          16 :         oid gid;
    3460          16 :         oid min, max;
    3461          16 :         BUN i, ngrp;
    3462          16 :         BUN nils = 0, nils2 = 0;
    3463          16 :         BUN *restrict cnts = NULL;
    3464          16 :         dbl *restrict dbls, *restrict mean, *restrict delta, *restrict m2;
    3465          16 :         BAT *bn = NULL, *an = NULL;
    3466          16 :         struct canditer ci;
    3467          16 :         BUN ncand;
    3468          16 :         const char *err;
    3469          16 :         lng t0 = 0;
    3470             : 
    3471          16 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3472             : 
    3473          16 :         assert(tp == TYPE_dbl);
    3474          16 :         (void) tp;              /* compatibility (with other BATgroup*
    3475             :                                  * functions) argument */
    3476             : 
    3477          16 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    3478           0 :                 GDKerror("%s: %s\n", func, err);
    3479           0 :                 return NULL;
    3480             :         }
    3481          16 :         if (g == NULL) {
    3482           0 :                 GDKerror("%s: b and g must be aligned\n", func);
    3483           0 :                 return NULL;
    3484             :         }
    3485             : 
    3486          16 :         if (BATcount(b) == 0 || ngrp == 0) {
    3487             :                 /* trivial: no products, so return bat aligned with g
    3488             :                  * with nil in the tail */
    3489           2 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
    3490           2 :                 goto doreturn;
    3491             :         }
    3492             : 
    3493          14 :         if ((e == NULL ||
    3494          14 :              (BATcount(e) == BATcount(b) && e->hseqbase == b->hseqbase)) &&
    3495           6 :             (BATtdense(g) || (g->tkey && g->tnonil)) &&
    3496           3 :             (issample || b->tnonil)) {
    3497             :                 /* trivial: singleton groups, so all results are equal
    3498             :                  * to zero (population) or nil (sample) */
    3499           5 :                 dbl v = issample ? dbl_nil : 0;
    3500           5 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
    3501           5 :                 goto doreturn;
    3502             :         }
    3503             : 
    3504           9 :         delta = GDKmalloc(ngrp * sizeof(dbl));
    3505           9 :         m2 = GDKmalloc(ngrp * sizeof(dbl));
    3506           9 :         cnts = GDKzalloc(ngrp * sizeof(BUN));
    3507           9 :         if (avgb) {
    3508           0 :                 an = COLnew(0, TYPE_dbl, ngrp, TRANSIENT);
    3509           0 :                 *avgb = an;
    3510           0 :                 if (an == NULL) {
    3511           0 :                         mean = NULL;
    3512           0 :                         goto alloc_fail;
    3513             :                 }
    3514           0 :                 mean = (dbl *) Tloc(an, 0);
    3515             :         } else {
    3516           9 :                 mean = GDKmalloc(ngrp * sizeof(dbl));
    3517             :         }
    3518           9 :         if (mean == NULL || delta == NULL || m2 == NULL || cnts == NULL)
    3519           0 :                 goto alloc_fail;
    3520             : 
    3521           9 :         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    3522           9 :         if (bn == NULL)
    3523           0 :                 goto alloc_fail;
    3524           9 :         dbls = (dbl *) Tloc(bn, 0);
    3525             : 
    3526     1080040 :         for (i = 0; i < ngrp; i++) {
    3527     1080030 :                 mean[i] = 0;
    3528     1080030 :                 delta[i] = 0;
    3529     1080030 :                 m2[i] = 0;
    3530             :         }
    3531             : 
    3532           9 :         if (BATtdense(g))
    3533             :                 gids = NULL;
    3534             :         else
    3535           8 :                 gids = (const oid *) Tloc(g, 0);
    3536             : 
    3537           9 :         switch (b->ttype) {
    3538           0 :         case TYPE_bte:
    3539           0 :                 AGGR_STDEV(bte);
    3540             :                 break;
    3541           0 :         case TYPE_sht:
    3542           0 :                 AGGR_STDEV(sht);
    3543             :                 break;
    3544           4 :         case TYPE_int:
    3545     5760020 :                 AGGR_STDEV(int);
    3546             :                 break;
    3547           0 :         case TYPE_lng:
    3548           0 :                 AGGR_STDEV(lng);
    3549             :                 break;
    3550             : #ifdef HAVE_HGE
    3551           0 :         case TYPE_hge:
    3552           0 :                 AGGR_STDEV(hge);
    3553             :                 break;
    3554             : #endif
    3555           0 :         case TYPE_flt:
    3556           0 :                 AGGR_STDEV(flt);
    3557             :                 break;
    3558           5 :         case TYPE_dbl:
    3559          77 :                 AGGR_STDEV(dbl);
    3560             :                 break;
    3561           0 :         default:
    3562           0 :                 if (an)
    3563           0 :                         BBPreclaim(an);
    3564             :                 else
    3565           0 :                         GDKfree(mean);
    3566           0 :                 GDKfree(delta);
    3567           0 :                 GDKfree(m2);
    3568           0 :                 GDKfree(cnts);
    3569           0 :                 BBPunfix(bn->batCacheid);
    3570           0 :                 GDKerror("%s: type (%s) not supported.\n",
    3571             :                          func, ATOMname(b->ttype));
    3572           0 :                 return NULL;
    3573             :         }
    3574           7 :         if (an) {
    3575           0 :                 BATsetcount(an, ngrp);
    3576           0 :                 an->tkey = ngrp <= 1;
    3577           0 :                 an->tsorted = ngrp <= 1;
    3578           0 :                 an->trevsorted = ngrp <= 1;
    3579           0 :                 an->tnil = nils != 0;
    3580           0 :                 an->tnonil = nils == 0;
    3581             :         } else {
    3582           7 :                 GDKfree(mean);
    3583             :         }
    3584           7 :         if (issample)
    3585           3 :                 nils += nils2;
    3586           7 :         GDKfree(delta);
    3587           7 :         GDKfree(m2);
    3588           7 :         GDKfree(cnts);
    3589           7 :         BATsetcount(bn, ngrp);
    3590           7 :         bn->tkey = ngrp <= 1;
    3591           7 :         bn->tsorted = ngrp <= 1;
    3592           7 :         bn->trevsorted = ngrp <= 1;
    3593           7 :         bn->tnil = nils != 0;
    3594           7 :         bn->tnonil = nils == 0;
    3595          14 :   doreturn:
    3596          14 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOBATFMT ",e=" ALGOOPTBATFMT
    3597             :                   ",s=" ALGOOPTBATFMT
    3598             :                   ",skip_nils=%s,issample=%s,variance=%s -> " ALGOOPTBATFMT
    3599             :                   ",avgb=" ALGOOPTBATFMT " (%s -- " LLFMT " usec)\n",
    3600             :                   ALGOBATPAR(b), ALGOBATPAR(g), ALGOOPTBATPAR(e),
    3601             :                   ALGOOPTBATPAR(s),
    3602             :                   skip_nils ? "true" : "false",
    3603             :                   issample ? "true" : "false",
    3604             :                   variance ? "true" : "false",
    3605             :                   ALGOOPTBATPAR(bn), ALGOOPTBATPAR(an),
    3606             :                   func, GDKusec() - t0);
    3607             :         return bn;
    3608           2 :   overflow:
    3609           2 :         GDKerror("22003!overflow in calculation.\n");
    3610           2 :   alloc_fail:
    3611           2 :         if (an)
    3612           0 :                 BBPreclaim(an);
    3613             :         else
    3614           2 :                 GDKfree(mean);
    3615           2 :         BBPreclaim(bn);
    3616           2 :         GDKfree(delta);
    3617           2 :         GDKfree(m2);
    3618           2 :         GDKfree(cnts);
    3619           2 :         return NULL;
    3620             : }
    3621             : 
    3622             : BAT *
    3623           7 : BATgroupstdev_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    3624             :                      bool skip_nils, bool abort_on_error)
    3625             : {
    3626           7 :         (void) abort_on_error;
    3627           7 :         return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, false,
    3628             :                             __func__);
    3629             : }
    3630             : 
    3631             : BAT *
    3632           5 : BATgroupstdev_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    3633             :                          bool skip_nils, bool abort_on_error)
    3634             : {
    3635           5 :         (void) abort_on_error;
    3636           5 :         return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, false,
    3637             :                             __func__);
    3638             : }
    3639             : 
    3640             : BAT *
    3641           1 : BATgroupvariance_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    3642             :                         bool skip_nils, bool abort_on_error)
    3643             : {
    3644           1 :         (void) abort_on_error;
    3645           1 :         return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, true,
    3646             :                             __func__);
    3647             : }
    3648             : 
    3649             : BAT *
    3650           3 : BATgroupvariance_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    3651             :                             bool skip_nils, bool abort_on_error)
    3652             : {
    3653           3 :         (void) abort_on_error;
    3654           3 :         return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, true,
    3655             :                             __func__);
    3656             : }
    3657             : 
    3658             : #define AGGR_COVARIANCE(TYPE)                                           \
    3659             :         do {                                                            \
    3660             :                 const TYPE *vals1 = (const TYPE *) Tloc(b1, 0);         \
    3661             :                 const TYPE *vals2 = (const TYPE *) Tloc(b2, 0);         \
    3662             :                 while (ncand > 0) {                                  \
    3663             :                         ncand--;                                        \
    3664             :                         i = canditer_next(&ci) - b1->hseqbase;           \
    3665             :                         if (gids == NULL ||                             \
    3666             :                             (gids[i] >= min && gids[i] <= max)) { \
    3667             :                                 if (gids)                               \
    3668             :                                         gid = gids[i] - min;            \
    3669             :                                 else                                    \
    3670             :                                         gid = (oid) i;                  \
    3671             :                                 if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
    3672             :                                         if (!skip_nils)                 \
    3673             :                                                 cnts[gid] = BUN_NONE;   \
    3674             :                                 } else if (cnts[gid] != BUN_NONE) {     \
    3675             :                                         cnts[gid]++;                    \
    3676             :                                         delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
    3677             :                                         mean1[gid] += delta1[gid] / cnts[gid]; \
    3678             :                                         delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
    3679             :                                         mean2[gid] += delta2[gid] / cnts[gid]; \
    3680             :                                         m2[gid] += delta1[gid] * ((dbl) vals2[i] - mean2[gid]); \
    3681             :                                 }                                       \
    3682             :                         }                                               \
    3683             :                 }                                                       \
    3684             :                 for (i = 0; i < ngrp; i++) {                         \
    3685             :                         if (cnts[i] == 0 || cnts[i] == BUN_NONE) {      \
    3686             :                                 dbls[i] = dbl_nil;                      \
    3687             :                                 nils++;                                 \
    3688             :                         } else if (cnts[i] == 1) {                      \
    3689             :                                 dbls[i] = issample ? dbl_nil : 0;       \
    3690             :                                 nils2++;                                \
    3691             :                         } else if (isinf(m2[i])) {              \
    3692             :                                 goto overflow;          \
    3693             :                         } else {                                        \
    3694             :                                 dbls[i] = m2[i] / (cnts[i] - issample); \
    3695             :                         }                                               \
    3696             :                 }                                                       \
    3697             :         } while (0)
    3698             : 
    3699             : static BAT *
    3700          60 : dogroupcovariance(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp,
    3701             :                   bool skip_nils, bool issample, const char *func)
    3702             : {
    3703          60 :         const oid *restrict gids;
    3704          60 :         oid gid, min, max;
    3705          60 :         BUN i, ngrp, nils = 0, nils2 = 0, ncand;
    3706          60 :         BUN *restrict cnts = NULL;
    3707          60 :         dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict m2;
    3708          60 :         BAT *bn = NULL;
    3709          60 :         struct canditer ci;
    3710          60 :         const char *err;
    3711          60 :         lng t0 = 0;
    3712             : 
    3713          60 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3714             : 
    3715         180 :         assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
    3716          60 :         (void) tp;
    3717             : 
    3718          60 :         if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    3719           0 :                 GDKerror("%s: %s\n", func, err);
    3720           0 :                 return NULL;
    3721             :         }
    3722          60 :         if (g == NULL) {
    3723           0 :                 GDKerror("%s: b1, b2 and g must be aligned\n", func);
    3724           0 :                 return NULL;
    3725             :         }
    3726             : 
    3727          60 :         if (BATcount(b1) == 0 || ngrp == 0) {
    3728           0 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
    3729           0 :                 goto doreturn;
    3730             :         }
    3731             : 
    3732          60 :         if ((e == NULL ||
    3733          60 :              (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
    3734          20 :             (BATtdense(g) || (g->tkey && g->tnonil)) &&
    3735          10 :             (issample || (b1->tnonil && b2->tnonil))) {
    3736             :                 /* trivial: singleton groups, so all results are equal
    3737             :                  * to zero (population) or nil (sample) */
    3738          13 :                 dbl v = issample ? dbl_nil : 0;
    3739          13 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
    3740          13 :                 goto doreturn;
    3741             :         }
    3742             : 
    3743          47 :         delta1 = GDKmalloc(ngrp * sizeof(dbl));
    3744          47 :         delta2 = GDKmalloc(ngrp * sizeof(dbl));
    3745          47 :         m2 = GDKmalloc(ngrp * sizeof(dbl));
    3746          47 :         cnts = GDKzalloc(ngrp * sizeof(BUN));
    3747          47 :         mean1 = GDKmalloc(ngrp * sizeof(dbl));
    3748          47 :         mean2 = GDKmalloc(ngrp * sizeof(dbl));
    3749             : 
    3750          47 :         if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || m2 == NULL || cnts == NULL)
    3751           0 :                 goto alloc_fail;
    3752             : 
    3753         357 :         for (i = 0; i < ngrp; i++) {
    3754         310 :                 m2[i] = 0;
    3755         310 :                 mean1[i] = 0;
    3756         310 :                 mean2[i] = 0;
    3757             :         }
    3758             : 
    3759          47 :         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    3760          47 :         if (bn == NULL)
    3761           0 :                 goto alloc_fail;
    3762          47 :         dbls = (dbl *) Tloc(bn, 0);
    3763             : 
    3764          47 :         if (!g || BATtdense(g))
    3765             :                 gids = NULL;
    3766             :         else
    3767          47 :                 gids = (const oid *) Tloc(g, 0);
    3768             : 
    3769          47 :         switch (b1->ttype) {
    3770          14 :         case TYPE_bte:
    3771         246 :                 AGGR_COVARIANCE(bte);
    3772             :                 break;
    3773           0 :         case TYPE_sht:
    3774           0 :                 AGGR_COVARIANCE(sht);
    3775             :                 break;
    3776          33 :         case TYPE_int:
    3777         563 :                 AGGR_COVARIANCE(int);
    3778             :                 break;
    3779           0 :         case TYPE_lng:
    3780           0 :                 AGGR_COVARIANCE(lng);
    3781             :                 break;
    3782             : #ifdef HAVE_HGE
    3783           0 :         case TYPE_hge:
    3784           0 :                 AGGR_COVARIANCE(hge);
    3785             :                 break;
    3786             : #endif
    3787           0 :         case TYPE_flt:
    3788           0 :                 AGGR_COVARIANCE(flt);
    3789             :                 break;
    3790           0 :         case TYPE_dbl:
    3791           0 :                 AGGR_COVARIANCE(dbl);
    3792             :                 break;
    3793           0 :         default:
    3794           0 :                 BBPreclaim(bn);
    3795           0 :                 GDKfree(mean1);
    3796           0 :                 GDKfree(mean2);
    3797           0 :                 GDKfree(delta1);
    3798           0 :                 GDKfree(delta2);
    3799           0 :                 GDKfree(m2);
    3800           0 :                 GDKfree(cnts);
    3801           0 :                 GDKerror("%s: type (%s) not supported.\n", func, ATOMname(b1->ttype));
    3802           0 :                 return NULL;
    3803             :         }
    3804          47 :         GDKfree(mean1);
    3805          47 :         GDKfree(mean2);
    3806             : 
    3807          47 :         if (issample)
    3808          20 :                 nils += nils2;
    3809          47 :         GDKfree(delta1);
    3810          47 :         GDKfree(delta2);
    3811          47 :         GDKfree(m2);
    3812          47 :         GDKfree(cnts);
    3813          47 :         BATsetcount(bn, ngrp);
    3814          46 :         bn->tkey = ngrp <= 1;
    3815          46 :         bn->tsorted = ngrp <= 1;
    3816          46 :         bn->trevsorted = ngrp <= 1;
    3817          46 :         bn->tnil = nils != 0;
    3818          46 :         bn->tnonil = nils == 0;
    3819          59 :   doreturn:
    3820          59 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
    3821             :                   ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
    3822             :                   ",skip_nils=%s,issample=%s -> " ALGOOPTBATFMT
    3823             :                   " (%s -- " LLFMT " usec)\n",
    3824             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
    3825             :                   ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
    3826             :                   skip_nils ? "true" : "false",
    3827             :                   issample ? "true" : "false",
    3828             :                   ALGOOPTBATPAR(bn),
    3829             :                   func, GDKusec() - t0);
    3830             :         return bn;
    3831           0 :   overflow:
    3832           0 :         GDKerror("22003!overflow in calculation.\n");
    3833           0 :   alloc_fail:
    3834           0 :         BBPreclaim(bn);
    3835           0 :         GDKfree(mean1);
    3836           0 :         GDKfree(mean2);
    3837           0 :         GDKfree(delta1);
    3838           0 :         GDKfree(delta2);
    3839           0 :         GDKfree(m2);
    3840           0 :         GDKfree(cnts);
    3841           0 :         return NULL;
    3842             : }
    3843             : 
    3844             : BAT *
    3845          30 : BATgroupcovariance_sample(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    3846             : {
    3847          30 :         (void) abort_on_error;
    3848          30 :         return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, true,
    3849             :                                  __func__);
    3850             : }
    3851             : 
    3852             : BAT *
    3853          30 : BATgroupcovariance_population(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    3854             : {
    3855          30 :         (void) abort_on_error;
    3856          30 :         return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, false,
    3857             :                                  __func__);
    3858             : }
    3859             : 
    3860             : #define AGGR_CORRELATION(TYPE)                                          \
    3861             :         do {                                                            \
    3862             :                 const TYPE *vals1 = (const TYPE *) Tloc(b1, 0);         \
    3863             :                 const TYPE *vals2 = (const TYPE *) Tloc(b2, 0);         \
    3864             :                 while (ncand > 0) {                                  \
    3865             :                         ncand--;                                        \
    3866             :                         i = canditer_next(&ci) - b1->hseqbase;           \
    3867             :                         if (gids == NULL ||                             \
    3868             :                             (gids[i] >= min && gids[i] <= max)) { \
    3869             :                                 if (gids)                               \
    3870             :                                         gid = gids[i] - min;            \
    3871             :                                 else                                    \
    3872             :                                         gid = (oid) i;                  \
    3873             :                                 if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
    3874             :                                         if (!skip_nils)                 \
    3875             :                                                 cnts[gid] = BUN_NONE;   \
    3876             :                                 } else if (cnts[gid] != BUN_NONE) {     \
    3877             :                                         cnts[gid]++;                    \
    3878             :                                         delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
    3879             :                                         mean1[gid] += delta1[gid] / cnts[gid]; \
    3880             :                                         delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
    3881             :                                         mean2[gid] += delta2[gid] / cnts[gid]; \
    3882             :                                         aux = (dbl) vals2[i] - mean2[gid]; \
    3883             :                                         up[gid] += delta1[gid] * aux;   \
    3884             :                                         down1[gid] += delta1[gid] * ((dbl) vals1[i] - mean1[gid]); \
    3885             :                                         down2[gid] += delta2[gid] * aux; \
    3886             :                                 }                                       \
    3887             :                         }                                               \
    3888             :                 }                                                       \
    3889             :                 for (i = 0; i < ngrp; i++) {                         \
    3890             :                         if (cnts[i] <= 1 || cnts[i] == BUN_NONE || down1[i] == 0 || down2[i] == 0) { \
    3891             :                                 dbls[i] = dbl_nil;                      \
    3892             :                                 nils++;                                 \
    3893             :                         } else if (isinf(up[i]) || isinf(down1[i]) || isinf(down2[i])) {        \
    3894             :                                 goto overflow;          \
    3895             :                         } else {                                        \
    3896             :                                 dbls[i] = (up[i] / cnts[i]) / (sqrt(down1[i] / cnts[i]) * sqrt(down2[i] / cnts[i])); \
    3897             :                                 assert(!is_dbl_nil(dbls[i]));           \
    3898             :                         }                                               \
    3899             :                 }                                                       \
    3900             :         } while (0)
    3901             : 
    3902             : BAT *
    3903          46 : BATgroupcorrelation(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    3904             : {
    3905          46 :         const oid *restrict gids;
    3906          46 :         oid gid, min, max;
    3907          46 :         BUN i, ngrp, nils = 0, ncand;
    3908          46 :         BUN *restrict cnts = NULL;
    3909          46 :         dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict up, *restrict down1, *restrict down2, aux;
    3910          46 :         BAT *bn = NULL;
    3911          46 :         struct canditer ci;
    3912          46 :         const char *err;
    3913          46 :         lng t0 = 0;
    3914             : 
    3915          46 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3916             : 
    3917         138 :         assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
    3918          46 :         (void) tp;
    3919          46 :         (void) abort_on_error;
    3920             : 
    3921          46 :         if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    3922           0 :                 GDKerror("%s\n", err);
    3923           0 :                 return NULL;
    3924             :         }
    3925          46 :         if (g == NULL) {
    3926           0 :                 GDKerror("b1, b2 and g must be aligned\n");
    3927           0 :                 return NULL;
    3928             :         }
    3929             : 
    3930          46 :         if (BATcount(b1) == 0 || ngrp == 0) {
    3931           1 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
    3932           1 :                 goto doreturn;
    3933             :         }
    3934             : 
    3935          45 :         if ((e == NULL ||
    3936          45 :              (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
    3937          11 :             (BATtdense(g) || (g->tkey && g->tnonil))) {
    3938          11 :                 dbl v = dbl_nil;
    3939          11 :                 bn = BATconstant(min, TYPE_dbl, &v, ngrp, TRANSIENT);
    3940          11 :                 goto doreturn;
    3941             :         }
    3942             : 
    3943          34 :         delta1 = GDKmalloc(ngrp * sizeof(dbl));
    3944          34 :         delta2 = GDKmalloc(ngrp * sizeof(dbl));
    3945          34 :         up = GDKmalloc(ngrp * sizeof(dbl));
    3946          34 :         down1 = GDKmalloc(ngrp * sizeof(dbl));
    3947          34 :         down2 = GDKmalloc(ngrp * sizeof(dbl));
    3948          33 :         cnts = GDKzalloc(ngrp * sizeof(BUN));
    3949          34 :         mean1 = GDKmalloc(ngrp * sizeof(dbl));
    3950          34 :         mean2 = GDKmalloc(ngrp * sizeof(dbl));
    3951             : 
    3952          34 :         if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || up == NULL || down1 == NULL || down2 == NULL || cnts == NULL)
    3953           0 :                 goto alloc_fail;
    3954             : 
    3955         239 :         for (i = 0; i < ngrp; i++) {
    3956         205 :                 up[i] = 0;
    3957         205 :                 down1[i] = 0;
    3958         205 :                 down2[i] = 0;
    3959         205 :                 mean1[i] = 0;
    3960         205 :                 mean2[i] = 0;
    3961             :         }
    3962             : 
    3963          34 :         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    3964          34 :         if (bn == NULL)
    3965           0 :                 goto alloc_fail;
    3966          34 :         dbls = (dbl *) Tloc(bn, 0);
    3967             : 
    3968          34 :         if (!g || BATtdense(g))
    3969             :                 gids = NULL;
    3970             :         else
    3971          34 :                 gids = (const oid *) Tloc(g, 0);
    3972             : 
    3973          34 :         switch (b1->ttype) {
    3974           6 :         case TYPE_bte:
    3975         102 :                 AGGR_CORRELATION(bte);
    3976             :                 break;
    3977           0 :         case TYPE_sht:
    3978           0 :                 AGGR_CORRELATION(sht);
    3979             :                 break;
    3980          23 :         case TYPE_int:
    3981         562 :                 AGGR_CORRELATION(int);
    3982             :                 break;
    3983           4 :         case TYPE_lng:
    3984         312 :                 AGGR_CORRELATION(lng);
    3985             :                 break;
    3986             : #ifdef HAVE_HGE
    3987           0 :         case TYPE_hge:
    3988           0 :                 AGGR_CORRELATION(hge);
    3989             :                 break;
    3990             : #endif
    3991           0 :         case TYPE_flt:
    3992           0 :                 AGGR_CORRELATION(flt);
    3993             :                 break;
    3994           1 :         case TYPE_dbl:
    3995           3 :                 AGGR_CORRELATION(dbl);
    3996             :                 break;
    3997           0 :         default:
    3998           0 :                 BBPreclaim(bn);
    3999           0 :                 GDKfree(mean1);
    4000           0 :                 GDKfree(mean2);
    4001           0 :                 GDKfree(delta1);
    4002           0 :                 GDKfree(delta2);
    4003           0 :                 GDKfree(up);
    4004           0 :                 GDKfree(down1);
    4005           0 :                 GDKfree(down2);
    4006           0 :                 GDKfree(cnts);
    4007           0 :                 GDKerror("type (%s) not supported.\n", ATOMname(b1->ttype));
    4008           0 :                 return NULL;
    4009             :         }
    4010          33 :         GDKfree(mean1);
    4011          33 :         GDKfree(mean2);
    4012          33 :         GDKfree(delta1);
    4013          33 :         GDKfree(delta2);
    4014          33 :         GDKfree(up);
    4015          33 :         GDKfree(down1);
    4016          33 :         GDKfree(down2);
    4017          33 :         GDKfree(cnts);
    4018          33 :         BATsetcount(bn, ngrp);
    4019          33 :         bn->tkey = ngrp <= 1;
    4020          33 :         bn->tsorted = ngrp <= 1;
    4021          33 :         bn->trevsorted = ngrp <= 1;
    4022          33 :         bn->tnil = nils != 0;
    4023          33 :         bn->tnonil = nils == 0;
    4024          45 :   doreturn:
    4025          45 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
    4026             :                   ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
    4027             :                   ",skip_nils=%s -> " ALGOOPTBATFMT
    4028             :                   " (" LLFMT " usec)\n",
    4029             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
    4030             :                   ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
    4031             :                   skip_nils ? "true" : "false",
    4032             :                   ALGOOPTBATPAR(bn),
    4033             :                   GDKusec() - t0);
    4034             :         return bn;
    4035           1 :   overflow:
    4036           1 :         GDKerror("22003!overflow in calculation.\n");
    4037           1 :   alloc_fail:
    4038           1 :         BBPreclaim(bn);
    4039           1 :         GDKfree(mean1);
    4040           1 :         GDKfree(mean2);
    4041           1 :         GDKfree(delta1);
    4042           1 :         GDKfree(delta2);
    4043           1 :         GDKfree(up);
    4044           1 :         GDKfree(down1);
    4045           1 :         GDKfree(down2);
    4046           1 :         GDKfree(cnts);
    4047           1 :         return NULL;
    4048             : }

Generated by: LCOV version 1.14