LCOV - code coverage report
Current view: top level - gdk - gdk_aggr.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 1406 2312 60.8 %
Date: 2021-10-13 02:24:04 Functions: 51 54 94.4 %

          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 - 2021 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             : 
      14             : /* grouped aggregates
      15             :  *
      16             :  * The following functions take two to four input BATs and produce a
      17             :  * single output BAT.
      18             :  *
      19             :  * The input BATs are
      20             :  * - b, a dense-headed BAT with the values to work on in the tail;
      21             :  * - g, a dense-headed BAT, aligned with b, with group ids (OID) in
      22             :  *   the tail;
      23             :  * - e, optional but recommended, a dense-headed BAT with the list of
      24             :  *   group ids in the head(!) (the tail is completely ignored);
      25             :  * - s, optional, a dense-headed bat with a list of candidate ids in
      26             :  *   the tail.
      27             :  *
      28             :  * The tail values of s refer to the head of b and g.  Only entries at
      29             :  * the specified ids are taken into account for the grouped
      30             :  * aggregates.  All other values are ignored.  s is compatible with
      31             :  * the result of BATselect().
      32             :  *
      33             :  * If e is not specified, we need to do an extra scan over g to find
      34             :  * out the range of the group ids that are used.  e is defined in such
      35             :  * a way that it can be either the extents or the histo result from
      36             :  * BATgroups().
      37             :  *
      38             :  * All functions calculate grouped aggregates.  There are as many
      39             :  * groups as there are entries in e.  If e is not specified, the
      40             :  * number of groups is equal to the difference between the maximum and
      41             :  * minimum values in g.
      42             :  *
      43             :  * If a group is empty, the result for that group is nil.
      44             :  *
      45             :  * If there is overflow during the calculation of an aggregate, the
      46             :  * whole operation fails if abort_on_error is set to non-zero,
      47             :  * otherwise the result of the group in which the overflow occurred is
      48             :  * nil.
      49             :  *
      50             :  * If skip_nils is non-zero, a nil value in b is ignored, otherwise a
      51             :  * nil in b results in a nil result for the group.
      52             :  */
      53             : 
      54             : /* helper function
      55             :  *
      56             :  * This function finds the minimum and maximum group id (and the
      57             :  * number of groups) and initializes the variables for candidates
      58             :  * selection.
      59             :  */
      60             : const char *
      61       31696 : BATgroupaggrinit(BAT *b, BAT *g, BAT *e, BAT *s,
      62             :                  /* outputs: */
      63             :                  oid *minp, oid *maxp, BUN *ngrpp,
      64             :                  struct canditer *ci, BUN *ncand)
      65             : {
      66             :         oid min, max;
      67             :         BUN i, ngrp;
      68             :         const oid *restrict gids;
      69             : 
      70       31696 :         if (b == NULL)
      71             :                 return "b must exist";
      72       31696 :         *ncand = canditer_init(ci, b, s);
      73       31696 :         if (g) {
      74       23216 :                 if (ci->ncand != BATcount(g) ||
      75       14369 :                     (ci->ncand != 0 && ci->seq != g->hseqbase))
      76             :                         return "b with s and g must be aligned";
      77       23216 :                 assert(BATttype(g) == TYPE_oid);
      78             :         }
      79       31696 :         if (g == NULL) {
      80             :                 min = 0;
      81             :                 max = 0;
      82             :                 ngrp = 1;
      83       23216 :         } else if (e == NULL) {
      84           0 :                 if (g->tmaxpos != BUN_NONE) {
      85             :                         min = 0;
      86           0 :                         max = BUNtoid(g, g->tmaxpos);
      87             :                 } else {
      88           0 :                         min = oid_nil;  /* note that oid_nil > 0! (unsigned) */
      89             :                         max = 0;
      90           0 :                         if (BATtdense(g)) {
      91             :                                 min = g->tseqbase;
      92           0 :                                 max = g->tseqbase + BATcount(g) - 1;
      93           0 :                         } else if (g->tsorted) {
      94           0 :                                 gids = (const oid *) Tloc(g, 0);
      95             :                                 /* find first non-nil */
      96           0 :                                 for (i = 0, ngrp = BATcount(g); i < ngrp; i++) {
      97           0 :                                         if (!is_oid_nil(gids[i])) {
      98             :                                                 min = gids[i];
      99             :                                                 break;
     100             :                                         }
     101             :                                 }
     102           0 :                                 if (!is_oid_nil(min)) {
     103             :                                         /* found a non-nil, max must be last
     104             :                                          * value (and there is one!) */
     105           0 :                                         max = gids[BUNlast(g) - 1];
     106             :                                 }
     107             :                         } else {
     108             :                                 /* we'll do a complete scan */
     109           0 :                                 gids = (const oid *) Tloc(g, 0);
     110           0 :                                 for (i = 0, ngrp = BATcount(g); i < ngrp; i++) {
     111           0 :                                         if (!is_oid_nil(gids[i])) {
     112             :                                                 if (gids[i] < min)
     113             :                                                         min = gids[i];
     114             :                                                 if (gids[i] > max)
     115             :                                                         max = gids[i];
     116             :                                         }
     117             :                                 }
     118             :                                 /* note: max < min is possible if all groups
     119             :                                  * are nil (or BATcount(g)==0) */
     120             :                         }
     121             :                 }
     122           0 :                 ngrp = max < min ? 0 : max - min + 1;
     123             :         } else {
     124       23216 :                 ngrp = BATcount(e);
     125       23216 :                 min = e->hseqbase;
     126       23216 :                 max = e->hseqbase + ngrp - 1;
     127             :         }
     128       31696 :         *minp = min;
     129       31696 :         *maxp = max;
     130       31696 :         *ngrpp = ngrp;
     131             : 
     132       31696 :         return NULL;
     133             : }
     134             : 
     135             : /* ---------------------------------------------------------------------- */
     136             : /* sum */
     137             : 
     138             : static inline bool
     139             : samesign(double x, double y)
     140             : {
     141         326 :         return (x >= 0) == (y >= 0);
     142             : }
     143             : 
     144             : /* Add two values, producing the sum and the remainder due to limited
     145             :  * range of floating point arithmetic.  This function depends on the
     146             :  * fact that the sum returns INFINITY in *hi of the correct sign
     147             :  * (i.e. isinf() returns TRUE) in case of overflow. */
     148             : static inline void
     149       21427 : twosum(volatile double *hi, volatile double *lo, double x, double y)
     150             : {
     151             :         volatile double yr;
     152             : 
     153       21427 :         assert(fabs(x) >= fabs(y));
     154             : 
     155       21427 :         *hi = x + y;
     156       21427 :         yr = *hi - x;
     157       21427 :         *lo = y - yr;
     158       21427 : }
     159             : 
     160             : static inline void
     161             : exchange(double *x, double *y)
     162             : {
     163             :         double t = *x;
     164             :         *x = *y;
     165             :         *y = t;
     166        9943 : }
     167             : 
     168             : /* this function was adapted from https://bugs.python.org/file10357/msum4.py */
     169             : BUN
     170        1100 : dofsum(const void *restrict values, oid seqb,
     171             :        struct canditer *restrict ci, BUN ncand,
     172             :        void *restrict results, BUN ngrp, int tp1, int tp2,
     173             :        const oid *restrict gids,
     174             :        oid min, oid max, bool skip_nils, bool abort_on_error,
     175             :        bool nil_if_empty)
     176             : {
     177             :         struct pergroup {
     178             :                 int npartials;
     179             :                 int maxpartials;
     180             :                 bool valseen;
     181             : #ifdef INFINITES_ALLOWED
     182             :                 float infs;
     183             : #else
     184             :                 int infs;
     185             : #endif
     186             :                 double *partials;
     187             :         } *pergroup;
     188             :         BUN listi;
     189             :         int parti;
     190             :         int i;
     191             :         BUN grp;
     192             :         double x, y;
     193             :         volatile double lo, hi;
     194             :         double twopow = pow((double) FLT_RADIX, (double) (DBL_MAX_EXP - 1));
     195             :         BUN nils = 0;
     196             :         volatile flt f;
     197             : 
     198             :         lng timeoffset = 0;
     199        1100 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
     200        1100 :         if (qry_ctx != NULL) {
     201        1100 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
     202             :         }
     203             : 
     204             :         /* we only deal with the two floating point types */
     205        1100 :         assert(tp1 == TYPE_flt || tp1 == TYPE_dbl);
     206        1100 :         assert(tp2 == TYPE_flt || tp2 == TYPE_dbl);
     207             :         /* if input is dbl, then so it output */
     208        1100 :         assert(tp2 == TYPE_flt || tp1 == TYPE_dbl);
     209             :         /* if no gids, then we have a single group */
     210        1100 :         assert(ngrp == 1 || gids != NULL);
     211        1100 :         if (gids == NULL || ngrp == 1) {
     212             :                 min = max = 0;
     213             :                 ngrp = 1;
     214             :                 gids = NULL;
     215             :         }
     216        1100 :         pergroup = GDKmalloc(ngrp * sizeof(*pergroup));
     217        1100 :         if (pergroup == NULL)
     218             :                 return BUN_NONE;
     219       10720 :         for (grp = 0; grp < ngrp; grp++) {
     220        9620 :                 pergroup[grp].npartials = 0;
     221        9620 :                 pergroup[grp].valseen = false;
     222        9620 :                 pergroup[grp].maxpartials = 2;
     223        9620 :                 pergroup[grp].infs = 0;
     224        9620 :                 pergroup[grp].partials = GDKmalloc(pergroup[grp].maxpartials * sizeof(double));
     225        9620 :                 if (pergroup[grp].partials == NULL) {
     226           0 :                         while (grp > 0)
     227           0 :                                 GDKfree(pergroup[--grp].partials);
     228           0 :                         GDKfree(pergroup);
     229           0 :                         return BUN_NONE;
     230             :                 }
     231             :         }
     232       25498 :         TIMEOUT_LOOP(ncand, timeoffset) {
     233       22198 :                 listi = canditer_next(ci) - seqb;
     234       22198 :                 grp = gids ? gids[listi] : 0;
     235       22198 :                 if (grp < min || grp > max)
     236           0 :                         continue;
     237       22198 :                 if (pergroup[grp].partials == NULL)
     238           0 :                         continue;
     239       22198 :                 if (tp1 == TYPE_flt && !is_flt_nil(((const flt *) values)[listi]))
     240          16 :                         x = ((const flt *) values)[listi];
     241       22182 :                 else if (tp1 == TYPE_dbl && !is_dbl_nil(((const dbl *) values)[listi]))
     242             :                         x = ((const dbl *) values)[listi];
     243             :                 else {
     244             :                         /* it's a nil */
     245          90 :                         if (!skip_nils) {
     246           0 :                                 if (tp2 == TYPE_flt)
     247           0 :                                         ((flt *) results)[grp] = flt_nil;
     248             :                                 else
     249           0 :                                         ((dbl *) results)[grp] = dbl_nil;
     250           0 :                                 GDKfree(pergroup[grp].partials);
     251           0 :                                 pergroup[grp].partials = NULL;
     252           0 :                                 if (++nils == ngrp)
     253           0 :                                         TIMEOUT_LOOP_BREAK;
     254             :                         }
     255          90 :                         continue;
     256             :                 }
     257       22108 :                 pergroup[grp].valseen = true;
     258             : #ifdef INFINITES_ALLOWED
     259             :                 if (isinf(x)) {
     260             :                         pergroup[grp].infs += x;
     261             :                         continue;
     262             :                 }
     263             : #endif
     264             :                 i = 0;
     265       42757 :                 for (parti = 0; parti < pergroup[grp].npartials; parti++) {
     266       20649 :                         y = pergroup[grp].partials[parti];
     267       20649 :                         if (fabs(x) < fabs(y))
     268             :                                 exchange(&x, &y);
     269       20649 :                         twosum(&hi, &lo, x, y);
     270       20649 :                         if (isinf(hi)) {
     271          58 :                                 int sign = hi > 0 ? 1 : -1;
     272          58 :                                 hi = x - twopow * sign;
     273          58 :                                 x = hi - twopow * sign;
     274          58 :                                 pergroup[grp].infs += sign;
     275          58 :                                 if (fabs(x) < fabs(y))
     276             :                                         exchange(&x, &y);
     277          58 :                                 twosum(&hi, &lo, x, y);
     278             :                         }
     279       20649 :                         if (lo != 0)
     280        9152 :                                 pergroup[grp].partials[i++] = lo;
     281       20649 :                         x = hi;
     282             :                 }
     283       22108 :                 if (x != 0) {
     284       22039 :                         if (i == pergroup[grp].maxpartials) {
     285             :                                 double *temp;
     286         281 :                                 pergroup[grp].maxpartials += pergroup[grp].maxpartials;
     287         281 :                                 temp = GDKrealloc(pergroup[grp].partials, pergroup[grp].maxpartials * sizeof(double));
     288         281 :                                 if (temp == NULL)
     289           0 :                                         goto bailout;
     290         281 :                                 pergroup[grp].partials = temp;
     291             :                         }
     292       22039 :                         pergroup[grp].partials[i++] = x;
     293             :                 }
     294       22108 :                 pergroup[grp].npartials = i;
     295             :         }
     296        1100 :         TIMEOUT_CHECK(timeoffset, GOTO_LABEL_TIMEOUT_HANDLER(bailout));
     297       10692 :         for (grp = 0; grp < ngrp; grp++) {
     298        9620 :                 if (pergroup[grp].partials == NULL)
     299           0 :                         continue;
     300        9620 :                 if (!pergroup[grp].valseen) {
     301          31 :                         if (tp2 == TYPE_flt)
     302           5 :                                 ((flt *) results)[grp] = nil_if_empty ? flt_nil : 0;
     303             :                         else
     304          26 :                                 ((dbl *) results)[grp] = nil_if_empty ? dbl_nil : 0;
     305          31 :                         nils += nil_if_empty;
     306          31 :                         GDKfree(pergroup[grp].partials);
     307          31 :                         pergroup[grp].partials = NULL;
     308          31 :                         continue;
     309             :                 }
     310             : #ifdef INFINITES_ALLOWED
     311             :                 if (isinf(pergroup[grp].infs) || isnan(pergroup[grp].infs)) {
     312             :                         if (abort_on_error) {
     313             :                                 goto overflow;
     314             :                         }
     315             :                         if (tp2 == TYPE_flt)
     316             :                                 ((flt *) results)[grp] = flt_nil;
     317             :                         else
     318             :                                 ((dbl *) results)[grp] = dbl_nil;
     319             :                         nils++;
     320             :                         GDKfree(pergroup[grp].partials);
     321             :                         pergroup[grp].partials = NULL;
     322             :                         continue;
     323             :                 }
     324             : #endif
     325             : 
     326        9589 :                 if ((pergroup[grp].infs == 1 || pergroup[grp].infs == -1) &&
     327          42 :                     pergroup[grp].npartials > 0 &&
     328          38 :                     !samesign(pergroup[grp].infs, pergroup[grp].partials[pergroup[grp].npartials - 1])) {
     329          34 :                         twosum(&hi, &lo, pergroup[grp].infs * twopow, pergroup[grp].partials[pergroup[grp].npartials - 1] / 2);
     330          34 :                         if (isinf(2 * hi)) {
     331          22 :                                 y = 2 * lo;
     332          22 :                                 x = hi + y;
     333          22 :                                 x -= hi;
     334          22 :                                 if (x == y &&
     335          18 :                                     pergroup[grp].npartials > 1 &&
     336          12 :                                     samesign(lo, pergroup[grp].partials[pergroup[grp].npartials - 2])) {
     337           6 :                                         GDKfree(pergroup[grp].partials);
     338           6 :                                         pergroup[grp].partials = NULL;
     339           6 :                                         x = 2 * (hi + y);
     340           6 :                                         if (tp2 == TYPE_flt) {
     341           0 :                                                 f = (flt) x;
     342           0 :                                                 if (isinf(f) ||
     343           0 :                                                     isnan(f) ||
     344           0 :                                                     is_flt_nil(f)) {
     345           0 :                                                         if (abort_on_error)
     346           0 :                                                                 goto overflow;
     347           0 :                                                         ((flt *) results)[grp] = flt_nil;
     348           0 :                                                         nils++;
     349             :                                                 } else
     350           0 :                                                         ((flt *) results)[grp] = f;
     351           6 :                                         } else if (is_dbl_nil(x)) {
     352           0 :                                                 if (abort_on_error)
     353           0 :                                                         goto overflow;
     354           0 :                                                 ((dbl *) results)[grp] = dbl_nil;
     355           0 :                                                 nils++;
     356             :                                         } else
     357           6 :                                                 ((dbl *) results)[grp] = x;
     358           6 :                                         continue;
     359             :                                 }
     360             :                         } else {
     361          12 :                                 if (lo) {
     362           2 :                                         if (pergroup[grp].npartials == pergroup[grp].maxpartials) {
     363             :                                                 double *temp;
     364             :                                                 /* we need space for one more */
     365           0 :                                                 pergroup[grp].maxpartials++;
     366           0 :                                                 temp = GDKrealloc(pergroup[grp].partials, pergroup[grp].maxpartials * sizeof(double));
     367           0 :                                                 if (temp == NULL)
     368           0 :                                                         goto bailout;
     369           0 :                                                 pergroup[grp].partials = temp;
     370             :                                         }
     371           2 :                                         pergroup[grp].partials[pergroup[grp].npartials - 1] = 2 * lo;
     372           2 :                                         pergroup[grp].partials[pergroup[grp].npartials++] = 2 * hi;
     373             :                                 } else {
     374          10 :                                         pergroup[grp].partials[pergroup[grp].npartials - 1] = 2 * hi;
     375             :                                 }
     376          12 :                                 pergroup[grp].infs = 0;
     377             :                         }
     378             :                 }
     379             : 
     380        9583 :                 if (pergroup[grp].infs != 0)
     381          28 :                         goto overflow;
     382             : 
     383        9555 :                 if (pergroup[grp].npartials == 0) {
     384          12 :                         GDKfree(pergroup[grp].partials);
     385          12 :                         pergroup[grp].partials = NULL;
     386          12 :                         if (tp2 == TYPE_flt)
     387           0 :                                 ((flt *) results)[grp] = 0;
     388             :                         else
     389          12 :                                 ((dbl *) results)[grp] = 0;
     390          12 :                         continue;
     391             :                 }
     392             : 
     393             :                 /* accumulate into hi */
     394        9543 :                 hi = pergroup[grp].partials[--pergroup[grp].npartials];
     395        9543 :                 while (pergroup[grp].npartials > 0) {
     396         686 :                         twosum(&hi, &lo, hi, pergroup[grp].partials[--pergroup[grp].npartials]);
     397         686 :                         if (lo) {
     398         686 :                                 pergroup[grp].partials[pergroup[grp].npartials++] = lo;
     399         686 :                                 break;
     400             :                         }
     401             :                 }
     402             : 
     403        9543 :                 if (pergroup[grp].npartials >= 2 &&
     404         276 :                     samesign(pergroup[grp].partials[pergroup[grp].npartials - 1], pergroup[grp].partials[pergroup[grp].npartials - 2]) &&
     405          41 :                     hi + 2 * pergroup[grp].partials[pergroup[grp].npartials - 1] - hi == 2 * pergroup[grp].partials[pergroup[grp].npartials - 1]) {
     406          29 :                         hi += 2 * pergroup[grp].partials[pergroup[grp].npartials - 1];
     407          29 :                         pergroup[grp].partials[pergroup[grp].npartials - 1] = -pergroup[grp].partials[pergroup[grp].npartials - 1];
     408             :                 }
     409             : 
     410        9543 :                 GDKfree(pergroup[grp].partials);
     411        9543 :                 pergroup[grp].partials = NULL;
     412        9543 :                 if (tp2 == TYPE_flt) {
     413          11 :                         f = (flt) hi;
     414          11 :                         if (isinf(f) || isnan(f) || is_flt_nil(f)) {
     415           0 :                                 if (abort_on_error)
     416           0 :                                         goto overflow;
     417           0 :                                 ((flt *) results)[grp] = flt_nil;
     418           0 :                                 nils++;
     419             :                         } else
     420          11 :                                 ((flt *) results)[grp] = f;
     421        9532 :                 } else if (is_dbl_nil(hi)) {
     422           0 :                         if (abort_on_error)
     423           0 :                                 goto overflow;
     424           0 :                         ((dbl *) results)[grp] = dbl_nil;
     425           0 :                         nils++;
     426             :                 } else
     427        9532 :                         ((dbl *) results)[grp] = hi;
     428             :         }
     429        1072 :         GDKfree(pergroup);
     430        1072 :         return nils;
     431             : 
     432          28 :   overflow:
     433          28 :         GDKerror("22003!overflow in sum aggregate.\n");
     434             :   bailout:
     435          56 :         for (grp = 0; grp < ngrp; grp++)
     436          28 :                 GDKfree(pergroup[grp].partials);
     437          28 :         GDKfree(pergroup);
     438          28 :         return BUN_NONE;
     439             : }
     440             : 
     441             : #define AGGR_SUM(TYPE1, TYPE2)                                          \
     442             :         do {                                                            \
     443             :                 TYPE1 x;                                                \
     444             :                 const TYPE1 *restrict vals = (const TYPE1 *) values;    \
     445             :                 if (ngrp == 1 && ci->tpe == cand_dense) {            \
     446             :                         /* single group, no candidate list */           \
     447             :                         TYPE2 sum;                                      \
     448             :                         *algo = "sum: no candidates, no groups";      \
     449             :                         sum = 0;                                        \
     450             :                         if (nonil) {                                    \
     451             :                                 *seen = ncand > 0;                   \
     452             :                                 TIMEOUT_LOOP_IDX(i, ncand, timeoffset) { \
     453             :                                         x = vals[ci->seq + i - seqb];        \
     454             :                                         ADD_WITH_CHECK(x, sum,          \
     455             :                                                        TYPE2, sum,      \
     456             :                                                        GDK_##TYPE2##_max, \
     457             :                                                        goto overflow);  \
     458             :                                 }                                       \
     459             :                                 TIMEOUT_CHECK(timeoffset,               \
     460             :                                               GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     461             :                         } else {                                        \
     462             :                                 bool seenval = false;                   \
     463             :                                 TIMEOUT_LOOP_IDX(i, ncand, timeoffset) { \
     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             :                                                         TIMEOUT_LOOP_BREAK; \
     470             :                                                 }                       \
     471             :                                         } else {                        \
     472             :                                                 ADD_WITH_CHECK(x, sum,  \
     473             :                                                                TYPE2, sum, \
     474             :                                                                GDK_##TYPE2##_max, \
     475             :                                                                goto overflow); \
     476             :                                                 seenval = true;         \
     477             :                                         }                               \
     478             :                                 }                                       \
     479             :                                 TIMEOUT_CHECK(timeoffset,               \
     480             :                                               GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     481             :                                 *seen = seenval;                        \
     482             :                         }                                               \
     483             :                         if (*seen)                                      \
     484             :                                 *sums = sum;                            \
     485             :                 } else if (ngrp == 1) {                                 \
     486             :                         /* single group, with candidate list */         \
     487             :                         TYPE2 sum;                                      \
     488             :                         bool seenval = false;                           \
     489             :                         *algo = "sum: with candidates, no groups";    \
     490             :                         sum = 0;                                        \
     491             :                         TIMEOUT_LOOP_IDX(i, ncand, timeoffset) {        \
     492             :                                 x = vals[canditer_next(ci) - seqb];     \
     493             :                                 if (is_##TYPE1##_nil(x)) {              \
     494             :                                         if (!skip_nils) {               \
     495             :                                                 sum = TYPE2##_nil;      \
     496             :                                                 nils = 1;               \
     497             :                                                 TIMEOUT_LOOP_BREAK;     \
     498             :                                         }                               \
     499             :                                 } else {                                \
     500             :                                         ADD_WITH_CHECK(x, sum,          \
     501             :                                                        TYPE2, sum,      \
     502             :                                                        GDK_##TYPE2##_max, \
     503             :                                                        goto overflow);  \
     504             :                                         seenval = true;                 \
     505             :                                 }                                       \
     506             :                         }                                               \
     507             :                         TIMEOUT_CHECK(timeoffset,                       \
     508             :                                       GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     509             :                         if (seenval)                                    \
     510             :                                 *sums = sum;                            \
     511             :                 } else if (ci->tpe == cand_dense) {                  \
     512             :                         /* multiple groups, no candidate list */        \
     513             :                         *algo = "sum: no candidates, with groups";    \
     514             :                         TIMEOUT_LOOP_IDX(i, ncand, timeoffset) {        \
     515             :                                 if (gids == NULL ||                     \
     516             :                                     (gids[i] >= min && gids[i] <= max)) { \
     517             :                                         gid = gids ? gids[i] - min : (oid) i; \
     518             :                                         x = vals[ci->seq + i - seqb];        \
     519             :                                         if (is_##TYPE1##_nil(x)) {      \
     520             :                                                 if (!skip_nils) {       \
     521             :                                                         sums[gid] = TYPE2##_nil; \
     522             :                                                         nils++;         \
     523             :                                                 }                       \
     524             :                                         } else {                        \
     525             :                                                 if (nil_if_empty &&     \
     526             :                                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     527             :                                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
     528             :                                                         sums[gid] = 0;  \
     529             :                                                 }                       \
     530             :                                                 if (!is_##TYPE2##_nil(sums[gid])) { \
     531             :                                                         ADD_WITH_CHECK( \
     532             :                                                                 x,      \
     533             :                                                                 sums[gid], \
     534             :                                                                 TYPE2,  \
     535             :                                                                 sums[gid], \
     536             :                                                                 GDK_##TYPE2##_max, \
     537             :                                                                 goto overflow); \
     538             :                                                 }                       \
     539             :                                         }                               \
     540             :                                 }                                       \
     541             :                         }                                               \
     542             :                         TIMEOUT_CHECK(timeoffset,                       \
     543             :                                       GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     544             :                 } else {                                                \
     545             :                         /* multiple groups, with candidate list */      \
     546             :                         *algo = "sum: with candidates, with groups";  \
     547             :                         TIMEOUT_LOOP(ncand, timeoffset) {               \
     548             :                                 i = canditer_next(ci) - seqb;           \
     549             :                                 if (gids == NULL ||                     \
     550             :                                     (gids[i] >= min && gids[i] <= max)) { \
     551             :                                         gid = gids ? gids[i] - min : (oid) i; \
     552             :                                         x = vals[i];                    \
     553             :                                         if (is_##TYPE1##_nil(x)) {      \
     554             :                                                 if (!skip_nils) {       \
     555             :                                                         sums[gid] = TYPE2##_nil; \
     556             :                                                         nils++;         \
     557             :                                                 }                       \
     558             :                                         } else {                        \
     559             :                                                 if (nil_if_empty &&     \
     560             :                                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     561             :                                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
     562             :                                                         sums[gid] = 0;  \
     563             :                                                 }                       \
     564             :                                                 if (!is_##TYPE2##_nil(sums[gid])) { \
     565             :                                                         ADD_WITH_CHECK( \
     566             :                                                                 x,      \
     567             :                                                                 sums[gid], \
     568             :                                                                 TYPE2,  \
     569             :                                                                 sums[gid], \
     570             :                                                                 GDK_##TYPE2##_max, \
     571             :                                                                 goto overflow); \
     572             :                                                 }                       \
     573             :                                         }                               \
     574             :                                 }                                       \
     575             :                         }                                               \
     576             :                         TIMEOUT_CHECK(timeoffset,                       \
     577             :                                       GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     578             :                 }                                                       \
     579             :         } while (0)
     580             : 
     581             : #define AGGR_SUM_NOOVL(TYPE1, TYPE2)                                    \
     582             :         do {                                                            \
     583             :                 TYPE1 x;                                                \
     584             :                 const TYPE1 *restrict vals = (const TYPE1 *) values;    \
     585             :                 if (ngrp == 1 && ci->tpe == cand_dense) {            \
     586             :                         /* single group, no candidate list */           \
     587             :                         TYPE2 sum;                                      \
     588             :                         sum = 0;                                        \
     589             :                         if (nonil) {                                    \
     590             :                                 *algo = "sum: no candidates, no groups, no nils, no overflow"; \
     591             :                                 *seen = ncand > 0;                   \
     592             :                                 TIMEOUT_LOOP_IDX(i, ncand, timeoffset) { \
     593             :                                         sum += vals[ci->seq + i - seqb]; \
     594             :                                 }                                       \
     595             :                                 TIMEOUT_CHECK(timeoffset,               \
     596             :                                               GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     597             :                         } else {                                        \
     598             :                                 bool seenval = false;                   \
     599             :                                 *algo = "sum: no candidates, no groups, no overflow"; \
     600             :                                 TIMEOUT_LOOP_IDX(i, ncand, timeoffset) { \
     601             :                                         x = vals[ci->seq + i - seqb];        \
     602             :                                         if (is_##TYPE1##_nil(x)) {      \
     603             :                                                 if (!skip_nils) {       \
     604             :                                                         sum = TYPE2##_nil; \
     605             :                                                         nils = 1;       \
     606             :                                                         TIMEOUT_LOOP_BREAK; \
     607             :                                                 }                       \
     608             :                                         } else {                        \
     609             :                                                 sum += x;               \
     610             :                                                 seenval = true;         \
     611             :                                         }                               \
     612             :                                 }                                       \
     613             :                                 TIMEOUT_CHECK(timeoffset,               \
     614             :                                               GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     615             :                                 *seen = seenval;                        \
     616             :                         }                                               \
     617             :                         if (*seen)                                      \
     618             :                                 *sums = sum;                            \
     619             :                 } else if (ngrp == 1) {                                 \
     620             :                         /* single group, with candidate list */         \
     621             :                         TYPE2 sum;                                      \
     622             :                         bool seenval = false;                           \
     623             :                         *algo = "sum: with candidates, no groups, no overflow"; \
     624             :                         sum = 0;                                        \
     625             :                         TIMEOUT_LOOP_IDX(i, ncand, timeoffset) {        \
     626             :                                 x = vals[canditer_next(ci) - seqb];     \
     627             :                                 if (is_##TYPE1##_nil(x)) {              \
     628             :                                         if (!skip_nils) {               \
     629             :                                                 sum = TYPE2##_nil;      \
     630             :                                                 nils = 1;               \
     631             :                                                 TIMEOUT_LOOP_BREAK;     \
     632             :                                         }                               \
     633             :                                 } else {                                \
     634             :                                         sum += x;                       \
     635             :                                         seenval = true;                 \
     636             :                                 }                                       \
     637             :                         }                                               \
     638             :                         TIMEOUT_CHECK(timeoffset,                       \
     639             :                                       GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     640             :                         if (seenval)                                    \
     641             :                                 *sums = sum;                            \
     642             :                 } else if (ci->tpe == cand_dense) {                  \
     643             :                         /* multiple groups, no candidate list */        \
     644             :                         if (nonil) {                                    \
     645             :                                 *algo = "sum: no candidates, with groups, no nils, no overflow"; \
     646             :                                 TIMEOUT_LOOP_IDX(i, ncand, timeoffset) { \
     647             :                                         if (gids == NULL ||             \
     648             :                                             (gids[i] >= min && gids[i] <= max)) { \
     649             :                                                 gid = gids ? gids[i] - min : (oid) i; \
     650             :                                                 x = vals[ci->seq + i - seqb]; \
     651             :                                                 if (nil_if_empty &&     \
     652             :                                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     653             :                                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
     654             :                                                         sums[gid] = 0;  \
     655             :                                                 }                       \
     656             :                                                 sums[gid] += x;         \
     657             :                                         }                               \
     658             :                                 }                                       \
     659             :                                 TIMEOUT_CHECK(timeoffset,               \
     660             :                                               GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     661             :                         } else {                                        \
     662             :                                 *algo = "sum: no candidates, with groups, no overflow"; \
     663             :                                 TIMEOUT_LOOP_IDX(i, ncand, timeoffset) { \
     664             :                                         if (gids == NULL ||             \
     665             :                                             (gids[i] >= min && gids[i] <= max)) { \
     666             :                                                 gid = gids ? gids[i] - min : (oid) i; \
     667             :                                                 x = vals[ci->seq + i - seqb]; \
     668             :                                                 if (is_##TYPE1##_nil(x)) { \
     669             :                                                         if (!skip_nils) { \
     670             :                                                                 sums[gid] = TYPE2##_nil; \
     671             :                                                                 nils++; \
     672             :                                                         }               \
     673             :                                                 } else {                \
     674             :                                                         if (nil_if_empty && \
     675             :                                                             !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     676             :                                                                 seen[gid >> 5] |= 1U << (gid & 0x1F); \
     677             :                                                                 sums[gid] = 0; \
     678             :                                                         }               \
     679             :                                                         if (!is_##TYPE2##_nil(sums[gid])) { \
     680             :                                                                 sums[gid] += x; \
     681             :                                                         }               \
     682             :                                                 }                       \
     683             :                                         }                               \
     684             :                                 }                                       \
     685             :                                 TIMEOUT_CHECK(timeoffset,               \
     686             :                                               GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     687             :                         }                                               \
     688             :                 } else {                                                \
     689             :                         /* multiple groups, with candidate list */      \
     690             :                         *algo = "sum: with candidates, with groups, no overflow"; \
     691             :                         TIMEOUT_LOOP(ncand, timeoffset) {               \
     692             :                                 i = canditer_next(ci) - seqb;           \
     693             :                                 if (gids == NULL ||                     \
     694             :                                     (gids[i] >= min && gids[i] <= max)) { \
     695             :                                         gid = gids ? gids[i] - min : (oid) i; \
     696             :                                         x = vals[i];                    \
     697             :                                         if (is_##TYPE1##_nil(x)) {      \
     698             :                                                 if (!skip_nils) {       \
     699             :                                                         sums[gid] = TYPE2##_nil; \
     700             :                                                         nils++;         \
     701             :                                                 }                       \
     702             :                                         } else {                        \
     703             :                                                 if (nil_if_empty &&     \
     704             :                                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
     705             :                                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
     706             :                                                         sums[gid] = 0;  \
     707             :                                                 }                       \
     708             :                                                 if (!is_##TYPE2##_nil(sums[gid])) { \
     709             :                                                         sums[gid] += x; \
     710             :                                                 }                       \
     711             :                                         }                               \
     712             :                                 }                                       \
     713             :                         }                                               \
     714             :                         TIMEOUT_CHECK(timeoffset,                       \
     715             :                                       GOTO_LABEL_TIMEOUT_HANDLER(bailout)); \
     716             :                 }                                                       \
     717             :         } while (0)
     718             : 
     719             : static BUN
     720       11042 : dosum(const void *restrict values, bool nonil, oid seqb,
     721             :       struct canditer *restrict ci, BUN ncand,
     722             :       void *restrict results, BUN ngrp, int tp1, int tp2,
     723             :       const oid *restrict gids,
     724             :       oid min, oid max, bool skip_nils, bool abort_on_error,
     725             :       bool nil_if_empty, const char *func, const char **algo)
     726             : {
     727             :         BUN nils = 0;
     728             :         BUN i;
     729             :         oid gid;
     730             :         unsigned int *restrict seen = NULL; /* bitmask for groups that we've seen */
     731             : 
     732             :         lng timeoffset = 0;
     733       11042 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
     734       11042 :         if (qry_ctx != NULL) {
     735       11042 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
     736             :         }
     737             : 
     738       11042 :         switch (tp2) {
     739          12 :         case TYPE_flt:
     740          12 :                 if (tp1 != TYPE_flt)
     741           0 :                         goto unsupported;
     742             :                 /* fall through */
     743             :         case TYPE_dbl:
     744        1095 :                 if (tp1 != TYPE_flt && tp1 != TYPE_dbl)
     745           0 :                         goto unsupported;
     746        1095 :                 *algo = "sum: floating point";
     747        1095 :                 return dofsum(values, seqb, ci, ncand, results, ngrp, tp1, tp2,
     748             :                               gids, min, max, skip_nils, abort_on_error,
     749             :                               nil_if_empty);
     750             :         }
     751             : 
     752             :         /* allocate bitmap for seen group ids */
     753        9947 :         seen = GDKzalloc(((ngrp + 31) / 32) * sizeof(int));
     754        9947 :         if (seen == NULL) {
     755             :                 return BUN_NONE;
     756             :         }
     757             : 
     758        9947 :         switch (tp2) {
     759           0 :         case TYPE_bte: {
     760             :                 bte *restrict sums = (bte *) results;
     761             :                 switch (tp1) {
     762           0 :                 case TYPE_bte:
     763           0 :                         AGGR_SUM(bte, bte);
     764             :                         break;
     765           0 :                 default:
     766           0 :                         goto unsupported;
     767             :                 }
     768             :                 break;
     769             :         }
     770           0 :         case TYPE_sht: {
     771             :                 sht *restrict sums = (sht *) results;
     772             :                 switch (tp1) {
     773           0 :                 case TYPE_bte:
     774           0 :                         if (ncand < ((BUN) 1 << ((sizeof(sht) - sizeof(bte)) << 3)))
     775           0 :                                 AGGR_SUM_NOOVL(bte, sht);
     776             :                         else
     777           0 :                                 AGGR_SUM(bte, sht);
     778             :                         break;
     779           0 :                 case TYPE_sht:
     780           0 :                         AGGR_SUM(sht, sht);
     781             :                         break;
     782           0 :                 default:
     783           0 :                         goto unsupported;
     784             :                 }
     785             :                 break;
     786             :         }
     787          12 :         case TYPE_int: {
     788             :                 int *restrict sums = (int *) results;
     789             :                 switch (tp1) {
     790           0 :                 case TYPE_bte:
     791           0 :                         if (ncand < ((BUN) 1 << ((sizeof(int) - sizeof(bte)) << 3)))
     792           0 :                                 AGGR_SUM_NOOVL(bte, int);
     793             :                         else
     794           0 :                                 AGGR_SUM(bte, int);
     795             :                         break;
     796           0 :                 case TYPE_sht:
     797           0 :                         if (ncand < ((BUN) 1 << ((sizeof(int) - sizeof(sht)) << 3)))
     798           0 :                                 AGGR_SUM_NOOVL(sht, int);
     799             :                         else
     800           0 :                                 AGGR_SUM(sht, int);
     801             :                         break;
     802          12 :                 case TYPE_int:
     803          58 :                         AGGR_SUM(int, int);
     804             :                         break;
     805           0 :                 default:
     806           0 :                         goto unsupported;
     807             :                 }
     808             :                 break;
     809             :         }
     810        7002 :         case TYPE_lng: {
     811             :                 lng *restrict sums = (lng *) results;
     812             :                 switch (tp1) {
     813             : #if SIZEOF_BUN == 4
     814             :                 case TYPE_bte:
     815             :                         AGGR_SUM_NOOVL(bte, lng);
     816             :                         break;
     817             :                 case TYPE_sht:
     818             :                         AGGR_SUM_NOOVL(sht, lng);
     819             :                         break;
     820             :                 case TYPE_int:
     821             :                         AGGR_SUM_NOOVL(int, lng);
     822             :                         break;
     823             : #else
     824           0 :                 case TYPE_bte:
     825           0 :                         if (ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(bte)) << 3)))
     826           0 :                                 AGGR_SUM_NOOVL(bte, lng);
     827             :                         else
     828           0 :                                 AGGR_SUM(bte, lng);
     829             :                         break;
     830           0 :                 case TYPE_sht:
     831           0 :                         if (ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(sht)) << 3)))
     832           0 :                                 AGGR_SUM_NOOVL(sht, lng);
     833             :                         else
     834           0 :                                 AGGR_SUM(sht, lng);
     835             :                         break;
     836           0 :                 case TYPE_int:
     837           0 :                         if (ncand < ((BUN) 1 << ((sizeof(lng) - sizeof(int)) << 3)))
     838           0 :                                 AGGR_SUM_NOOVL(int, lng);
     839             :                         else
     840           0 :                                 AGGR_SUM(int, lng);
     841             :                         break;
     842             : #endif
     843        7002 :                 case TYPE_lng:
     844     7390672 :                         AGGR_SUM(lng, lng);
     845             :                         break;
     846           0 :                 default:
     847           0 :                         goto unsupported;
     848             :                 }
     849             :                 break;
     850             :         }
     851             : #ifdef HAVE_HGE
     852        2933 :         case TYPE_hge: {
     853             :                 hge *sums = (hge *) results;
     854             :                 switch (tp1) {
     855         236 :                 case TYPE_bte:
     856     3365754 :                         AGGR_SUM_NOOVL(bte, hge);
     857             :                         break;
     858           7 :                 case TYPE_sht:
     859          45 :                         AGGR_SUM_NOOVL(sht, hge);
     860             :                         break;
     861        1904 :                 case TYPE_int:
     862    56479220 :                         AGGR_SUM_NOOVL(int, hge);
     863             :                         break;
     864          21 :                 case TYPE_lng:
     865      655828 :                         AGGR_SUM_NOOVL(lng, hge);
     866             :                         break;
     867         765 :                 case TYPE_hge:
     868    33093149 :                         AGGR_SUM(hge, hge);
     869             :                         break;
     870           0 :                 default:
     871           0 :                         goto unsupported;
     872             :                 }
     873             :                 break;
     874             :         }
     875             : #endif
     876           0 :         default:
     877           0 :                 goto unsupported;
     878             :         }
     879             : 
     880        9945 :         if (nils == 0 && nil_if_empty) {
     881             :                 /* figure out whether there were any empty groups
     882             :                  * (that result in a nil value) */
     883        9946 :                 if (ngrp & 0x1F) {
     884             :                         /* fill last slot */
     885        9911 :                         seen[ngrp >> 5] |= ~0U << (ngrp & 0x1F);
     886             :                 }
     887      235499 :                 for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
     888      226336 :                         if (seen[i] != ~0U) {
     889             :                                 nils = 1;
     890             :                                 break;
     891             :                         }
     892             :                 }
     893             :         }
     894        9945 :         GDKfree(seen);
     895             : 
     896        9944 :         return nils;
     897             : 
     898           0 :   unsupported:
     899           0 :         GDKfree(seen);
     900           0 :         GDKerror("%s: type combination (sum(%s)->%s) not supported.\n",
     901             :                  func, ATOMname(tp1), ATOMname(tp2));
     902           0 :         return BUN_NONE;
     903             : 
     904           0 :   overflow:
     905           0 :         GDKfree(seen);
     906           0 :         GDKerror("22003!overflow in sum aggregate.\n");
     907           0 :         return BUN_NONE;
     908             : 
     909           1 :   bailout:
     910           1 :         GDKfree(seen);
     911           1 :         return BUN_NONE;
     912             : }
     913             : 
     914             : /* calculate group sums with optional candidates list */
     915             : BAT *
     916        7776 : BATgroupsum(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
     917             : {
     918             :         const oid *restrict gids;
     919             :         oid min, max;
     920             :         BUN ngrp;
     921             :         BUN nils;
     922             :         BAT *bn;
     923             :         struct canditer ci;
     924             :         BUN ncand;
     925             :         const char *err;
     926        7776 :         const char *algo = NULL;
     927             :         lng t0 = 0;
     928             : 
     929        7776 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
     930             : 
     931        7776 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
     932           0 :                 GDKerror("%s\n", err);
     933           0 :                 return NULL;
     934             :         }
     935        7773 :         if (g == NULL) {
     936           0 :                 GDKerror("b and g must be aligned\n");
     937           0 :                 return NULL;
     938             :         }
     939             : 
     940        7773 :         if (ncand == 0 || ngrp == 0) {
     941             :                 /* trivial: no sums, so return bat aligned with g with
     942             :                  * nil in the tail */
     943        2673 :                 return BATconstant(ngrp == 0 ? 0 : min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
     944             :         }
     945             : 
     946        5100 :         if ((e == NULL ||
     947        5100 :              (BATcount(e) == ncand && e->hseqbase == ci.hseq)) &&
     948        1744 :             (BATtdense(g) || (g->tkey && g->tnonil))) {
     949             :                 /* trivial: singleton groups, so all results are equal
     950             :                  * to the inputs (but possibly a different type) */
     951        1744 :                 return BATconvert(b, s, tp, abort_on_error, 0, 0, 0);
     952             :         }
     953             : 
     954        3356 :         bn = BATconstant(min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
     955        3356 :         if (bn == NULL) {
     956             :                 return NULL;
     957             :         }
     958             : 
     959        3356 :         if (BATtdense(g))
     960             :                 gids = NULL;
     961             :         else
     962        3174 :                 gids = (const oid *) Tloc(g, 0);
     963             : 
     964        3356 :         BATiter bi = bat_iterator(b);
     965        3356 :         nils = dosum(bi.base, b->tnonil, b->hseqbase, &ci, ncand,
     966        3356 :                      Tloc(bn, 0), ngrp, b->ttype, tp, gids, min, max,
     967             :                      skip_nils, abort_on_error, true, __func__, &algo);
     968        3356 :         bat_iterator_end(&bi);
     969             : 
     970        3354 :         if (nils < BUN_NONE) {
     971        3354 :                 BATsetcount(bn, ngrp);
     972        3355 :                 bn->tkey = BATcount(bn) <= 1;
     973        3355 :                 bn->tsorted = BATcount(bn) <= 1;
     974        3355 :                 bn->trevsorted = BATcount(bn) <= 1;
     975        3355 :                 bn->tnil = nils != 0;
     976        3355 :                 bn->tnonil = nils == 0;
     977             :         } else {
     978           0 :                 BBPunfix(bn->batCacheid);
     979             :                 bn = NULL;
     980             :         }
     981             : 
     982        3355 :         if (algo)
     983        3355 :                 MT_thread_setalgorithm(algo);
     984        3356 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
     985             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
     986             :                   "; start " OIDFMT ", count " BUNFMT " (%s -- " LLFMT " usec)\n",
     987             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
     988             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
     989             :                   ci.seq, ncand, algo ? algo : "", GDKusec() - t0);
     990             :         return bn;
     991             : }
     992             : 
     993             : static BUN
     994           4 : mskCountOnes(BAT *b, struct canditer *ci)
     995             : {
     996           4 :         BUN cnt = 0, ncand = ci->ncand;
     997             : 
     998           4 :         if (ci->s == NULL && mask_cand(b))
     999           4 :                 return BATcount(b);
    1000           0 :         if (ci->tpe == cand_dense && ncand > 0 && !mask_cand(b)) {
    1001           0 :                 BATiter bi = bat_iterator(b);
    1002           0 :                 const uint32_t *restrict src = (const uint32_t *) bi.base + (ci->seq - b->hseqbase) / 32;
    1003           0 :                 int bits = (ci->seq - b->hseqbase) % 32;
    1004           0 :                 if (bits + ncand <= 32) {
    1005           0 :                         if (ncand == 32)
    1006           0 :                                 cnt = candmask_pop(src[0]);
    1007             :                         else
    1008           0 :                                 cnt = candmask_pop(src[0] & (((1U << ncand) - 1) << bits));
    1009           0 :                         bat_iterator_end(&bi);
    1010           0 :                         return cnt;
    1011             :                 }
    1012           0 :                 if (bits != 0) {
    1013           0 :                         cnt = candmask_pop(src[0] & (~0U << bits));
    1014           0 :                         src++;
    1015           0 :                         ncand -= 32 - bits;
    1016             :                 }
    1017           0 :                 while (ncand >= 32) {
    1018           0 :                         cnt += candmask_pop(*src);
    1019           0 :                         src++;
    1020           0 :                         ncand -= 32;
    1021             :                 }
    1022           0 :                 if (ncand > 0)
    1023           0 :                         cnt += candmask_pop(*src & ((1U << ncand) - 1));
    1024           0 :                 bat_iterator_end(&bi);
    1025           0 :                 return cnt;
    1026             :         }
    1027           0 :         for (BUN i = 0; i < ncand; i++) {
    1028           0 :                 BUN x = canditer_next(ci) - b->hseqbase;
    1029           0 :                 cnt += mskGetVal(b, x);
    1030             :         }
    1031             :         return cnt;
    1032             : }
    1033             : 
    1034             : gdk_return
    1035        8043 : BATsum(void *res, int tp, BAT *b, BAT *s, bool skip_nils, bool abort_on_error, bool nil_if_empty)
    1036             : {
    1037             :         oid min, max;
    1038             :         BUN ngrp;
    1039             :         BUN nils;
    1040             :         struct canditer ci;
    1041             :         BUN ncand;
    1042             :         const char *err;
    1043        8043 :         const char *algo = NULL;
    1044             :         lng t0 = 0;
    1045             : 
    1046        8043 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    1047             : 
    1048        8043 :         if ((err = BATgroupaggrinit(b, NULL, NULL, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    1049           0 :                 GDKerror("%s\n", err);
    1050           0 :                 return GDK_FAIL;
    1051             :         }
    1052        8043 :         if (ATOMstorage(b->ttype) == TYPE_msk || mask_cand(b)) {
    1053           4 :                 ncand = mskCountOnes(b, &ci);
    1054           4 :                 switch (tp) {
    1055           0 :                 case TYPE_bte:
    1056           0 :                         if (ncand > GDK_bte_max) {
    1057           0 :                                 GDKerror("22003!overflow in sum aggregate.\n");
    1058           0 :                                 return GDK_FAIL;
    1059             :                         }
    1060           0 :                         * (bte *) res = (bte) ncand;
    1061           0 :                         break;
    1062           0 :                 case TYPE_sht:
    1063           0 :                         if (ncand > GDK_sht_max) {
    1064           0 :                                 GDKerror("22003!overflow in sum aggregate.\n");
    1065           0 :                                 return GDK_FAIL;
    1066             :                         }
    1067           0 :                         * (sht *) res = (sht) ncand;
    1068           0 :                         break;
    1069           0 :                 case TYPE_int:
    1070             : #if SIZEOF_BUN > 4
    1071           0 :                         if (ncand > GDK_int_max) {
    1072           0 :                                 GDKerror("22003!overflow in sum aggregate.\n");
    1073           0 :                                 return GDK_FAIL;
    1074             :                         }
    1075             : #endif
    1076           0 :                         * (int *) res = (int) ncand;
    1077           0 :                         break;
    1078           4 :                 case TYPE_lng:
    1079           4 :                         * (lng *) res = (lng) ncand;
    1080           4 :                         break;
    1081             : #ifdef HAVE_HGE
    1082           0 :                 case TYPE_hge:
    1083           0 :                         * (hge *) res = (hge) ncand;
    1084           0 :                         break;
    1085             : #endif
    1086           0 :                 case TYPE_flt:
    1087           0 :                         * (flt *) res = (flt) ncand;
    1088           0 :                         break;
    1089           0 :                 case TYPE_dbl:
    1090           0 :                         * (dbl *) res = (dbl) ncand;
    1091           0 :                         break;
    1092           0 :                 default:
    1093           0 :                         GDKerror("type combination (sum(%s)->%s) not supported.\n",
    1094             :                                  ATOMname(b->ttype), ATOMname(tp));
    1095           0 :                         return GDK_FAIL;
    1096             :                 }
    1097           4 :                 TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",s=" ALGOOPTBATFMT "; "
    1098             :                           "start " OIDFMT ", count " BUNFMT " (pop count -- " LLFMT " usec)\n",
    1099             :                           ALGOBATPAR(b), ALGOOPTBATPAR(s),
    1100             :                           ci.seq, ci.ncand, GDKusec() - t0);
    1101           4 :                 return GDK_SUCCEED;
    1102             :         }
    1103        8039 :         switch (tp) {
    1104           0 :         case TYPE_bte:
    1105           0 :                 * (bte *) res = nil_if_empty ? bte_nil : 0;
    1106           0 :                 break;
    1107           0 :         case TYPE_sht:
    1108           0 :                 * (sht *) res = nil_if_empty ? sht_nil : 0;
    1109           0 :                 break;
    1110          16 :         case TYPE_int:
    1111          16 :                 * (int *) res = nil_if_empty ? int_nil : 0;
    1112          16 :                 break;
    1113        5959 :         case TYPE_lng:
    1114        5959 :                 * (lng *) res = nil_if_empty ? lng_nil : 0;
    1115        5959 :                 break;
    1116             : #ifdef HAVE_HGE
    1117         940 :         case TYPE_hge:
    1118         940 :                 * (hge *) res = nil_if_empty ? hge_nil : 0;
    1119         940 :                 break;
    1120             : #endif
    1121        1124 :         case TYPE_flt:
    1122             :         case TYPE_dbl:
    1123             :                 switch (b->ttype) {
    1124           1 :                 case TYPE_bte:
    1125             :                 case TYPE_sht:
    1126             :                 case TYPE_int:
    1127             :                 case TYPE_lng:
    1128             : #ifdef HAVE_HGE
    1129             :                 case TYPE_hge:
    1130             : #endif
    1131             :                 {
    1132             :                         /* special case for summing integer types into
    1133             :                          * a floating point: We calculate the average
    1134             :                          * (which is done exactly), and multiply the
    1135             :                          * result by the count to get the sum.  Note
    1136             :                          * that if we just summed into a floating
    1137             :                          * point number, we could loose too much
    1138             :                          * accuracy, and if we summed into lng first,
    1139             :                          * we could get unnecessary overflow. */
    1140             :                         dbl avg;
    1141             :                         BUN cnt;
    1142             : 
    1143           1 :                         if (BATcalcavg(b, s, &avg, &cnt, 0) != GDK_SUCCEED)
    1144             :                                 return GDK_FAIL;
    1145           1 :                         if (cnt == 0) {
    1146           0 :                                 avg = nil_if_empty ? dbl_nil : 0;
    1147             :                         }
    1148           1 :                         if (cnt < ncand && !skip_nils) {
    1149             :                                 /* there were nils */
    1150           0 :                                 avg = dbl_nil;
    1151             :                         }
    1152           1 :                         if (tp == TYPE_flt) {
    1153           0 :                                 if (is_dbl_nil(avg))
    1154           0 :                                         *(flt *) res = flt_nil;
    1155           0 :                                 else if (cnt > 0 &&
    1156           0 :                                          GDK_flt_max / cnt < fabs(avg)) {
    1157           0 :                                         if (abort_on_error) {
    1158           0 :                                                 GDKerror("22003!overflow in sum aggregate.\n");
    1159           0 :                                                 return GDK_FAIL;
    1160             :                                         }
    1161           0 :                                         *(flt *) res = flt_nil;
    1162             :                                 } else {
    1163           0 :                                         *(flt *) res = (flt) avg * cnt;
    1164             :                                 }
    1165             :                         } else {
    1166           1 :                                 if (is_dbl_nil(avg)) {
    1167           0 :                                         *(dbl *) res = dbl_nil;
    1168           1 :                                 } else if (cnt > 0 &&
    1169           1 :                                            GDK_dbl_max / cnt < fabs(avg)) {
    1170           0 :                                         if (abort_on_error) {
    1171           0 :                                                 GDKerror("22003!overflow in sum aggregate.\n");
    1172           0 :                                                 return GDK_FAIL;
    1173             :                                         }
    1174           0 :                                         *(dbl *) res = dbl_nil;
    1175             :                                 } else {
    1176           1 :                                         *(dbl *) res = avg * cnt;
    1177             :                                 }
    1178             :                         }
    1179             :                         return GDK_SUCCEED;
    1180             :                 }
    1181             :                 default:
    1182             :                         break;
    1183             :                 }
    1184        1123 :                 if (b->ttype == TYPE_dbl)
    1185        1110 :                         * (dbl *) res = nil_if_empty ? dbl_nil : 0;
    1186             :                 else
    1187          13 :                         * (flt *) res = nil_if_empty ? flt_nil : 0;
    1188             :                 break;
    1189           0 :         default:
    1190           0 :                 GDKerror("type combination (sum(%s)->%s) not supported.\n",
    1191             :                          ATOMname(b->ttype), ATOMname(tp));
    1192           0 :                 return GDK_FAIL;
    1193             :         }
    1194        8038 :         if (ncand == 0)
    1195             :                 return GDK_SUCCEED;
    1196        7686 :         BATiter bi = bat_iterator(b);
    1197        7686 :         nils = dosum(bi.base, b->tnonil, b->hseqbase, &ci, ncand,
    1198        7686 :                      res, true, b->ttype, tp, &min, min, max,
    1199             :                      skip_nils, abort_on_error, nil_if_empty, __func__, &algo);
    1200        7686 :         bat_iterator_end(&bi);
    1201        7686 :         if (algo)
    1202        7686 :                 MT_thread_setalgorithm(algo);
    1203        7685 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",s=" ALGOOPTBATFMT "; "
    1204             :                   "start " OIDFMT ", count " BUNFMT " (%s -- " LLFMT " usec)\n",
    1205             :                   ALGOBATPAR(b), ALGOOPTBATPAR(s),
    1206             :                   ci.seq, ncand, algo ? algo : "", GDKusec() - t0);
    1207        7685 :         return nils < BUN_NONE ? GDK_SUCCEED : GDK_FAIL;
    1208             : }
    1209             : 
    1210             : /* ---------------------------------------------------------------------- */
    1211             : /* product */
    1212             : 
    1213             : #define AGGR_PROD(TYPE1, TYPE2, TYPE3)                                  \
    1214             :         do {                                                            \
    1215             :                 const TYPE1 *restrict vals = (const TYPE1 *) values;    \
    1216             :                 gid = 0;        /* doesn't change if gidincr == false */ \
    1217             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    1218             :                         i = canditer_next(ci) - seqb;                   \
    1219             :                         if (gids == NULL || !gidincr ||                 \
    1220             :                             (gids[i] >= min && gids[i] <= max)) { \
    1221             :                                 if (gidincr) {                          \
    1222             :                                         if (gids)                       \
    1223             :                                                 gid = gids[i] - min;    \
    1224             :                                         else                            \
    1225             :                                                 gid = (oid) i;          \
    1226             :                                 }                                       \
    1227             :                                 if (is_##TYPE1##_nil(vals[i])) {        \
    1228             :                                         if (!skip_nils) {               \
    1229             :                                                 prods[gid] = TYPE2##_nil; \
    1230             :                                                 nils++;                 \
    1231             :                                         }                               \
    1232             :                                 } else {                                \
    1233             :                                         if (nil_if_empty &&             \
    1234             :                                             !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
    1235             :                                                 seen[gid >> 5] |= 1U << (gid & 0x1F); \
    1236             :                                                 prods[gid] = 1;         \
    1237             :                                         }                               \
    1238             :                                         if (!is_##TYPE2##_nil(prods[gid])) { \
    1239             :                                                 MUL4_WITH_CHECK(        \
    1240             :                                                         vals[i],        \
    1241             :                                                         prods[gid],     \
    1242             :                                                         TYPE2, prods[gid], \
    1243             :                                                         GDK_##TYPE2##_max, \
    1244             :                                                         TYPE3,          \
    1245             :                                                         goto overflow); \
    1246             :                                         }                               \
    1247             :                                 }                                       \
    1248             :                         }                                               \
    1249             :                 }                                                       \
    1250             :                 TIMEOUT_CHECK(timeoffset,                               \
    1251             :                               GOTO_LABEL_TIMEOUT_HANDLER(bailout));     \
    1252             :         } while (0)
    1253             : 
    1254             : #ifdef HAVE_HGE
    1255             : #define AGGR_PROD_HGE(TYPE)                                             \
    1256             :         do {                                                            \
    1257             :                 const TYPE *vals = (const TYPE *) values;               \
    1258             :                 gid = 0;        /* doesn't change if gidincr == false */ \
    1259             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    1260             :                         i = canditer_next(ci) - seqb;                   \
    1261             :                         if (gids == NULL || !gidincr ||                 \
    1262             :                             (gids[i] >= min && gids[i] <= max)) { \
    1263             :                                 if (gidincr) {                          \
    1264             :                                         if (gids)                       \
    1265             :                                                 gid = gids[i] - min;    \
    1266             :                                         else                            \
    1267             :                                                 gid = (oid) i;          \
    1268             :                                 }                                       \
    1269             :                                 if (nil_if_empty &&                     \
    1270             :                                     !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
    1271             :                                         seen[gid >> 5] |= 1U << (gid & 0x1F); \
    1272             :                                         prods[gid] = 1;                 \
    1273             :                                 }                                       \
    1274             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    1275             :                                         if (!skip_nils) {               \
    1276             :                                                 prods[gid] = hge_nil;   \
    1277             :                                                 nils++;                 \
    1278             :                                         }                               \
    1279             :                                 } else if (!is_hge_nil(prods[gid])) {   \
    1280             :                                         HGEMUL_CHECK(vals[i],           \
    1281             :                                                      prods[gid],        \
    1282             :                                                      prods[gid],        \
    1283             :                                                      GDK_hge_max,       \
    1284             :                                                      goto overflow);    \
    1285             :                                 }                                       \
    1286             :                         }                                               \
    1287             :                 }                                                       \
    1288             :                 TIMEOUT_CHECK(timeoffset,                               \
    1289             :                               GOTO_LABEL_TIMEOUT_HANDLER(bailout));     \
    1290             :         } while (0)
    1291             : #else
    1292             : #define AGGR_PROD_LNG(TYPE)                                             \
    1293             :         do {                                                            \
    1294             :                 const TYPE *restrict vals = (const TYPE *) values;      \
    1295             :                 gid = 0;        /* doesn't change if gidincr == false */ \
    1296             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    1297             :                         i = canditer_next(ci) - seqb;                   \
    1298             :                         if (gids == NULL || !gidincr ||                 \
    1299             :                             (gids[i] >= min && gids[i] <= max)) { \
    1300             :                                 if (gidincr) {                          \
    1301             :                                         if (gids)                       \
    1302             :                                                 gid = gids[i] - min;    \
    1303             :                                         else                            \
    1304             :                                                 gid = (oid) i;          \
    1305             :                                 }                                       \
    1306             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    1307             :                                         if (!skip_nils) {               \
    1308             :                                                 prods[gid] = lng_nil;   \
    1309             :                                                 nils++;                 \
    1310             :                                         }                               \
    1311             :                                 } else {                                \
    1312             :                                         if (nil_if_empty &&             \
    1313             :                                             !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
    1314             :                                                 seen[gid >> 5] |= 1U << (gid & 0x1F); \
    1315             :                                                 prods[gid] = 1;         \
    1316             :                                         }                               \
    1317             :                                         if (!is_lng_nil(prods[gid])) {  \
    1318             :                                                 LNGMUL_CHECK(           \
    1319             :                                                         vals[i],        \
    1320             :                                                         prods[gid],     \
    1321             :                                                         prods[gid],     \
    1322             :                                                         GDK_lng_max,    \
    1323             :                                                         goto overflow); \
    1324             :                                         }                               \
    1325             :                                 }                                       \
    1326             :                         }                                               \
    1327             :                 }                                                       \
    1328             :                 TIMEOUT_CHECK(timeoffset,                               \
    1329             :                               GOTO_LABEL_TIMEOUT_HANDLER(bailout));     \
    1330             :         } while (0)
    1331             : #endif
    1332             : 
    1333             : #define AGGR_PROD_FLOAT(TYPE1, TYPE2)                                   \
    1334             :         do {                                                            \
    1335             :                 const TYPE1 *restrict vals = (const TYPE1 *) values;    \
    1336             :                 gid = 0;        /* doesn't change if gidincr == false */ \
    1337             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    1338             :                         i = canditer_next(ci) - seqb;                   \
    1339             :                         if (gids == NULL || !gidincr ||                 \
    1340             :                             (gids[i] >= min && gids[i] <= max)) { \
    1341             :                                 if (gidincr) {                          \
    1342             :                                         if (gids)                       \
    1343             :                                                 gid = gids[i] - min;    \
    1344             :                                         else                            \
    1345             :                                                 gid = (oid) i;          \
    1346             :                                 }                                       \
    1347             :                                 if (is_##TYPE1##_nil(vals[i])) {        \
    1348             :                                         if (!skip_nils) {               \
    1349             :                                                 prods[gid] = TYPE2##_nil; \
    1350             :                                                 nils++;                 \
    1351             :                                         }                               \
    1352             :                                 } else {                                \
    1353             :                                         if (nil_if_empty &&             \
    1354             :                                             !(seen[gid >> 5] & (1U << (gid & 0x1F)))) { \
    1355             :                                                 seen[gid >> 5] |= 1U << (gid & 0x1F); \
    1356             :                                                 prods[gid] = 1;         \
    1357             :                                         }                               \
    1358             :                                         if (!is_##TYPE2##_nil(prods[gid])) { \
    1359             :                                                 if (ABSOLUTE(vals[i]) > 1 && \
    1360             :                                                     GDK_##TYPE2##_max / ABSOLUTE(vals[i]) < ABSOLUTE(prods[gid])) { \
    1361             :                                                         if (abort_on_error) \
    1362             :                                                                 goto overflow; \
    1363             :                                                         prods[gid] = TYPE2##_nil; \
    1364             :                                                         nils++;         \
    1365             :                                                 } else {                \
    1366             :                                                         prods[gid] *= vals[i]; \
    1367             :                                                 }                       \
    1368             :                                         }                               \
    1369             :                                 }                                       \
    1370             :                         }                                               \
    1371             :                 }                                                       \
    1372             :                 TIMEOUT_CHECK(timeoffset,                               \
    1373             :                               GOTO_LABEL_TIMEOUT_HANDLER(bailout));     \
    1374             :         } while (0)
    1375             : 
    1376             : static BUN
    1377          32 : doprod(const void *restrict values, oid seqb, struct canditer *restrict ci, BUN ncand,
    1378             :        void *restrict results, BUN ngrp, int tp1, int tp2,
    1379             :        const oid *restrict gids, bool gidincr, oid min, oid max,
    1380             :        bool skip_nils, bool abort_on_error, bool nil_if_empty, const char *func)
    1381             : {
    1382             :         BUN nils = 0;
    1383             :         BUN i;
    1384             :         oid gid;
    1385             :         unsigned int *restrict seen; /* bitmask for groups that we've seen */
    1386             : 
    1387             :         lng timeoffset = 0;
    1388          32 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    1389          32 :         if (qry_ctx != NULL) {
    1390          32 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    1391             :         }
    1392             : 
    1393             :         /* allocate bitmap for seen group ids */
    1394          32 :         seen = GDKzalloc(((ngrp + 31) / 32) * sizeof(int));
    1395          32 :         if (seen == NULL) {
    1396             :                 return BUN_NONE;
    1397             :         }
    1398             : 
    1399          32 :         switch (tp2) {
    1400           0 :         case TYPE_bte: {
    1401             :                 bte *restrict prods = (bte *) results;
    1402             :                 switch (tp1) {
    1403           0 :                 case TYPE_bte:
    1404           0 :                         AGGR_PROD(bte, bte, sht);
    1405             :                         break;
    1406           0 :                 default:
    1407           0 :                         goto unsupported;
    1408             :                 }
    1409             :                 break;
    1410             :         }
    1411           0 :         case TYPE_sht: {
    1412             :                 sht *restrict prods = (sht *) results;
    1413             :                 switch (tp1) {
    1414           0 :                 case TYPE_bte:
    1415           0 :                         AGGR_PROD(bte, sht, int);
    1416             :                         break;
    1417           0 :                 case TYPE_sht:
    1418           0 :                         AGGR_PROD(sht, sht, int);
    1419             :                         break;
    1420           0 :                 default:
    1421           0 :                         goto unsupported;
    1422             :                 }
    1423             :                 break;
    1424             :         }
    1425           0 :         case TYPE_int: {
    1426             :                 int *restrict prods = (int *) results;
    1427             :                 switch (tp1) {
    1428           0 :                 case TYPE_bte:
    1429           0 :                         AGGR_PROD(bte, int, lng);
    1430             :                         break;
    1431           0 :                 case TYPE_sht:
    1432           0 :                         AGGR_PROD(sht, int, lng);
    1433             :                         break;
    1434           0 :                 case TYPE_int:
    1435           0 :                         AGGR_PROD(int, int, lng);
    1436             :                         break;
    1437           0 :                 default:
    1438           0 :                         goto unsupported;
    1439             :                 }
    1440             :                 break;
    1441             :         }
    1442             : #ifdef HAVE_HGE
    1443           0 :         case TYPE_lng: {
    1444             :                 lng *prods = (lng *) results;
    1445             :                 switch (tp1) {
    1446           0 :                 case TYPE_bte:
    1447           0 :                         AGGR_PROD(bte, lng, hge);
    1448             :                         break;
    1449           0 :                 case TYPE_sht:
    1450           0 :                         AGGR_PROD(sht, lng, hge);
    1451             :                         break;
    1452           0 :                 case TYPE_int:
    1453           0 :                         AGGR_PROD(int, lng, hge);
    1454             :                         break;
    1455           0 :                 case TYPE_lng:
    1456           0 :                         AGGR_PROD(lng, lng, hge);
    1457             :                         break;
    1458           0 :                 default:
    1459           0 :                         goto unsupported;
    1460             :                 }
    1461             :                 break;
    1462             :         }
    1463          25 :         case TYPE_hge: {
    1464             :                 hge *prods = (hge *) results;
    1465             :                 switch (tp1) {
    1466           4 :                 case TYPE_bte:
    1467          14 :                         AGGR_PROD_HGE(bte);
    1468             :                         break;
    1469           4 :                 case TYPE_sht:
    1470          14 :                         AGGR_PROD_HGE(sht);
    1471             :                         break;
    1472           4 :                 case TYPE_int:
    1473          14 :                         AGGR_PROD_HGE(int);
    1474             :                         break;
    1475           8 :                 case TYPE_lng:
    1476          28 :                         AGGR_PROD_HGE(lng);
    1477             :                         break;
    1478           5 :                 case TYPE_hge:
    1479          30 :                         AGGR_PROD_HGE(hge);
    1480             :                         break;
    1481           0 :                 default:
    1482           0 :                         goto unsupported;
    1483             :                 }
    1484             :                 break;
    1485             :         }
    1486             : #else
    1487             :         case TYPE_lng: {
    1488             :                 lng *restrict prods = (lng *) results;
    1489             :                 switch (tp1) {
    1490             :                 case TYPE_bte:
    1491             :                         AGGR_PROD_LNG(bte);
    1492             :                         break;
    1493             :                 case TYPE_sht:
    1494             :                         AGGR_PROD_LNG(sht);
    1495             :                         break;
    1496             :                 case TYPE_int:
    1497             :                         AGGR_PROD_LNG(int);
    1498             :                         break;
    1499             :                 case TYPE_lng:
    1500             :                         AGGR_PROD_LNG(lng);
    1501             :                         break;
    1502             :                 default:
    1503             :                         goto unsupported;
    1504             :                 }
    1505             :                 break;
    1506             :         }
    1507             : #endif
    1508           0 :         case TYPE_flt: {
    1509             :                 flt *restrict prods = (flt *) results;
    1510             :                 switch (tp1) {
    1511           0 :                 case TYPE_bte:
    1512           0 :                         AGGR_PROD_FLOAT(bte, flt);
    1513             :                         break;
    1514           0 :                 case TYPE_sht:
    1515           0 :                         AGGR_PROD_FLOAT(sht, flt);
    1516             :                         break;
    1517           0 :                 case TYPE_int:
    1518           0 :                         AGGR_PROD_FLOAT(int, flt);
    1519             :                         break;
    1520           0 :                 case TYPE_lng:
    1521           0 :                         AGGR_PROD_FLOAT(lng, flt);
    1522             :                         break;
    1523             : #ifdef HAVE_HGE
    1524           0 :                 case TYPE_hge:
    1525           0 :                         AGGR_PROD_FLOAT(hge, flt);
    1526             :                         break;
    1527             : #endif
    1528           0 :                 case TYPE_flt:
    1529           0 :                         AGGR_PROD_FLOAT(flt, flt);
    1530             :                         break;
    1531           0 :                 default:
    1532           0 :                         goto unsupported;
    1533             :                 }
    1534             :                 break;
    1535             :         }
    1536           7 :         case TYPE_dbl: {
    1537             :                 dbl *restrict prods = (dbl *) results;
    1538             :                 switch (tp1) {
    1539           0 :                 case TYPE_bte:
    1540           0 :                         AGGR_PROD_FLOAT(bte, dbl);
    1541             :                         break;
    1542           0 :                 case TYPE_sht:
    1543           0 :                         AGGR_PROD_FLOAT(sht, dbl);
    1544             :                         break;
    1545           0 :                 case TYPE_int:
    1546           0 :                         AGGR_PROD_FLOAT(int, dbl);
    1547             :                         break;
    1548           0 :                 case TYPE_lng:
    1549           0 :                         AGGR_PROD_FLOAT(lng, dbl);
    1550             :                         break;
    1551             : #ifdef HAVE_HGE
    1552           0 :                 case TYPE_hge:
    1553           0 :                         AGGR_PROD_FLOAT(hge, dbl);
    1554             :                         break;
    1555             : #endif
    1556           0 :                 case TYPE_flt:
    1557           0 :                         AGGR_PROD_FLOAT(flt, dbl);
    1558             :                         break;
    1559           7 :                 case TYPE_dbl:
    1560          35 :                         AGGR_PROD_FLOAT(dbl, dbl);
    1561             :                         break;
    1562           0 :                 default:
    1563           0 :                         goto unsupported;
    1564             :                 }
    1565             :                 break;
    1566             :         }
    1567           0 :         default:
    1568           0 :                 goto unsupported;
    1569             :         }
    1570             : 
    1571          32 :         if (nils == 0 && nil_if_empty) {
    1572             :                 /* figure out whether there were any empty groups
    1573             :                  * (that result in a nil value) */
    1574          32 :                 if (ngrp & 0x1F) {
    1575             :                         /* fill last slot */
    1576          32 :                         seen[ngrp >> 5] |= ~0U << (ngrp & 0x1F);
    1577             :                 }
    1578          64 :                 for (i = 0, ngrp = (ngrp + 31) / 32; i < ngrp; i++) {
    1579          32 :                         if (seen[i] != ~0U) {
    1580             :                                 nils = 1;
    1581             :                                 break;
    1582             :                         }
    1583             :                 }
    1584             :         }
    1585          32 :         GDKfree(seen);
    1586             : 
    1587          32 :         return nils;
    1588             : 
    1589           0 :   unsupported:
    1590           0 :         GDKfree(seen);
    1591           0 :         GDKerror("%s: type combination (mul(%s)->%s) not supported.\n",
    1592             :                  func, ATOMname(tp1), ATOMname(tp2));
    1593           0 :         return BUN_NONE;
    1594             : 
    1595           0 :   overflow:
    1596           0 :         GDKfree(seen);
    1597           0 :         GDKerror("22003!overflow in product aggregate.\n");
    1598           0 :         return BUN_NONE;
    1599             : 
    1600           0 :   bailout:
    1601           0 :         GDKfree(seen);
    1602           0 :         return BUN_NONE;
    1603             : }
    1604             : 
    1605             : /* calculate group products with optional candidates list */
    1606             : BAT *
    1607          17 : BATgroupprod(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    1608             : {
    1609             :         const oid *restrict gids;
    1610             :         oid min, max;
    1611             :         BUN ngrp;
    1612             :         BUN nils;
    1613             :         BAT *bn;
    1614             :         struct canditer ci;
    1615             :         BUN ncand;
    1616             :         const char *err;
    1617             :         lng t0 = 0;
    1618             : 
    1619          17 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    1620             : 
    1621          17 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    1622           0 :                 GDKerror("%s\n", err);
    1623           0 :                 return NULL;
    1624             :         }
    1625          17 :         if (g == NULL) {
    1626           0 :                 GDKerror("b and g must be aligned\n");
    1627           0 :                 return NULL;
    1628             :         }
    1629             : 
    1630          17 :         if (ncand == 0 || ngrp == 0) {
    1631             :                 /* trivial: no products, so return bat aligned with g
    1632             :                  * with nil in the tail */
    1633           9 :                 return BATconstant(ngrp == 0 ? 0 : min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
    1634             :         }
    1635             : 
    1636           8 :         if ((e == NULL ||
    1637           8 :              (BATcount(e) == ncand && e->hseqbase == ci.hseq)) &&
    1638           6 :             (BATtdense(g) || (g->tkey && g->tnonil))) {
    1639             :                 /* trivial: singleton groups, so all results are equal
    1640             :                  * to the inputs (but possibly a different type) */
    1641           6 :                 return BATconvert(b, s, tp, abort_on_error, 0, 0, 0);
    1642             :         }
    1643             : 
    1644           2 :         bn = BATconstant(min, tp, ATOMnilptr(tp), ngrp, TRANSIENT);
    1645           2 :         if (bn == NULL) {
    1646             :                 return NULL;
    1647             :         }
    1648             : 
    1649           2 :         if (BATtdense(g))
    1650             :                 gids = NULL;
    1651             :         else
    1652           2 :                 gids = (const oid *) Tloc(g, 0);
    1653             : 
    1654           2 :         BATiter bi = bat_iterator(b);
    1655           2 :         nils = doprod(bi.base, b->hseqbase, &ci, ncand, Tloc(bn, 0), ngrp,
    1656           2 :                       b->ttype, tp, gids, true, min, max, skip_nils,
    1657             :                       abort_on_error, true, __func__);
    1658           2 :         bat_iterator_end(&bi);
    1659             : 
    1660           2 :         if (nils < BUN_NONE) {
    1661           2 :                 BATsetcount(bn, ngrp);
    1662           2 :                 bn->tkey = BATcount(bn) <= 1;
    1663           2 :                 bn->tsorted = BATcount(bn) <= 1;
    1664           2 :                 bn->trevsorted = BATcount(bn) <= 1;
    1665           2 :                 bn->tnil = nils != 0;
    1666           2 :                 bn->tnonil = nils == 0;
    1667             :         } else {
    1668           0 :                 BBPunfix(bn->batCacheid);
    1669             :                 bn = NULL;
    1670             :         }
    1671             : 
    1672           2 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    1673             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    1674             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    1675             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    1676             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    1677             :                   ci.seq, ncand, GDKusec() - t0);
    1678             : 
    1679             :         return bn;
    1680             : }
    1681             : 
    1682             : gdk_return
    1683          30 : BATprod(void *res, int tp, BAT *b, BAT *s, bool skip_nils, bool abort_on_error, bool nil_if_empty)
    1684             : {
    1685             :         oid min, max;
    1686             :         BUN ngrp;
    1687             :         BUN nils;
    1688             :         struct canditer ci;
    1689             :         BUN ncand;
    1690             :         const char *err;
    1691             :         lng t0 = 0;
    1692             : 
    1693          30 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    1694             : 
    1695          30 :         if ((err = BATgroupaggrinit(b, NULL, NULL, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    1696           0 :                 GDKerror("%s\n", err);
    1697           0 :                 return GDK_FAIL;
    1698             :         }
    1699          30 :         switch (tp) {
    1700           0 :         case TYPE_bte:
    1701           0 :                 * (bte *) res = nil_if_empty ? bte_nil : (bte) 1;
    1702           0 :                 break;
    1703           0 :         case TYPE_sht:
    1704           0 :                 * (sht *) res = nil_if_empty ? sht_nil : (sht) 1;
    1705           0 :                 break;
    1706           0 :         case TYPE_int:
    1707           0 :                 * (int *) res = nil_if_empty ? int_nil : (int) 1;
    1708           0 :                 break;
    1709           0 :         case TYPE_lng:
    1710           0 :                 * (lng *) res = nil_if_empty ? lng_nil : (lng) 1;
    1711           0 :                 break;
    1712             : #ifdef HAVE_HGE
    1713          25 :         case TYPE_hge:
    1714          25 :                 * (hge *) res = nil_if_empty ? hge_nil : (hge) 1;
    1715          25 :                 break;
    1716             : #endif
    1717           0 :         case TYPE_flt:
    1718           0 :                 * (flt *) res = nil_if_empty ? flt_nil : (flt) 1;
    1719           0 :                 break;
    1720           5 :         case TYPE_dbl:
    1721           5 :                 * (dbl *) res = nil_if_empty ? dbl_nil : (dbl) 1;
    1722           5 :                 break;
    1723           0 :         default:
    1724           0 :                 GDKerror("type combination (prod(%s)->%s) not supported.\n",
    1725             :                          ATOMname(b->ttype), ATOMname(tp));
    1726           0 :                 return GDK_FAIL;
    1727             :         }
    1728          30 :         if (ncand == 0)
    1729             :                 return GDK_SUCCEED;
    1730          30 :         BATiter bi = bat_iterator(b);
    1731          30 :         nils = doprod(bi.base, b->hseqbase, &ci, ncand, res, true,
    1732          30 :                       b->ttype, tp, &min, false, min, max,
    1733             :                       skip_nils, abort_on_error, nil_if_empty, __func__);
    1734          30 :         bat_iterator_end(&bi);
    1735          30 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",s=" ALGOOPTBATFMT "; "
    1736             :                   "start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    1737             :                   ALGOBATPAR(b), ALGOOPTBATPAR(s),
    1738             :                   ci.seq, ncand, GDKusec() - t0);
    1739          30 :         return nils < BUN_NONE ? GDK_SUCCEED : GDK_FAIL;
    1740             : }
    1741             : 
    1742             : /* ---------------------------------------------------------------------- */
    1743             : /* average */
    1744             : 
    1745             : #define AGGR_AVG(TYPE)                                                  \
    1746             :         do {                                                            \
    1747             :                 const TYPE *restrict vals = (const TYPE *) bi.base;     \
    1748             :                 TYPE *restrict avgs = GDKzalloc(ngrp * sizeof(TYPE));   \
    1749             :                 if (avgs == NULL) {                                     \
    1750             :                         bat_iterator_end(&bi);                              \
    1751             :                         goto alloc_fail;                                \
    1752             :                 }                                                       \
    1753             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    1754             :                         i = canditer_next(&ci) - b->hseqbase;            \
    1755             :                         if (gids == NULL ||                             \
    1756             :                             (gids[i] >= min && gids[i] <= max)) { \
    1757             :                                 if (gids)                               \
    1758             :                                         gid = gids[i] - min;            \
    1759             :                                 else                                    \
    1760             :                                         gid = (oid) i;                  \
    1761             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    1762             :                                         if (!skip_nils)                 \
    1763             :                                                 cnts[gid] = lng_nil;    \
    1764             :                                 } else if (!is_lng_nil(cnts[gid])) {    \
    1765             :                                         AVERAGE_ITER(TYPE, vals[i],     \
    1766             :                                                      avgs[gid],         \
    1767             :                                                      rems[gid],         \
    1768             :                                                      cnts[gid]);        \
    1769             :                                 }                                       \
    1770             :                         }                                               \
    1771             :                 }                                                       \
    1772             :                 TIMEOUT_CHECK(timeoffset,                               \
    1773             :                               GOTO_LABEL_TIMEOUT_HANDLER(alloc_fail));  \
    1774             :                 for (i = 0; i < ngrp; i++) {                         \
    1775             :                         if (cnts[i] == 0 || is_lng_nil(cnts[i])) {      \
    1776             :                                 dbls[i] = dbl_nil;                      \
    1777             :                                 cnts[i] = 0;                            \
    1778             :                                 nils++;                                 \
    1779             :                         } else {                                        \
    1780             :                                 dbls[i] = avgs[i] + (dbl) rems[i] / cnts[i]; \
    1781             :                         }                                               \
    1782             :                 }                                                       \
    1783             :                 GDKfree(avgs);                                          \
    1784             :         } while (0)
    1785             : 
    1786             : #define AGGR_AVG_FLOAT(TYPE)                                            \
    1787             :         do {                                                            \
    1788             :                 const TYPE *restrict vals = (const TYPE *) bi.base;     \
    1789             :                 for (i = 0; i < ngrp; i++)                           \
    1790             :                         dbls[i] = 0;                                    \
    1791             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    1792             :                         i = canditer_next(&ci) - b->hseqbase;            \
    1793             :                         if (gids == NULL ||                             \
    1794             :                             (gids[i] >= min && gids[i] <= max)) { \
    1795             :                                 if (gids)                               \
    1796             :                                         gid = gids[i] - min;            \
    1797             :                                 else                                    \
    1798             :                                         gid = (oid) i;                  \
    1799             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    1800             :                                         if (!skip_nils)                 \
    1801             :                                                 cnts[gid] = lng_nil;    \
    1802             :                                 } else if (!is_lng_nil(cnts[gid])) {    \
    1803             :                                         AVERAGE_ITER_FLOAT(TYPE, vals[i], \
    1804             :                                                            dbls[gid],   \
    1805             :                                                            cnts[gid]);  \
    1806             :                                 }                                       \
    1807             :                         }                                               \
    1808             :                 }                                                       \
    1809             :                 TIMEOUT_CHECK(timeoffset,                               \
    1810             :                               GOTO_LABEL_TIMEOUT_HANDLER(alloc_fail));  \
    1811             :                 for (i = 0; i < ngrp; i++) {                         \
    1812             :                         if (cnts[i] == 0 || is_lng_nil(cnts[i])) {      \
    1813             :                                 dbls[i] = dbl_nil;                      \
    1814             :                                 cnts[i] = 0;                            \
    1815             :                                 nils++;                                 \
    1816             :                         }                                               \
    1817             :                 }                                                       \
    1818             :         } while (0)
    1819             : 
    1820             : /* calculate group averages with optional candidates list */
    1821             : gdk_return
    1822         250 : BATgroupavg(BAT **bnp, BAT **cntsp, BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error, int scale)
    1823             : {
    1824             :         const oid *restrict gids;
    1825             :         oid gid;
    1826             :         oid min, max;
    1827             :         BUN i, ngrp;
    1828             :         BUN nils = 0;
    1829             :         lng *restrict rems = NULL;
    1830             :         lng *restrict cnts = NULL;
    1831             :         dbl *restrict dbls;
    1832             :         BAT *bn = NULL, *cn = NULL;
    1833             :         struct canditer ci;
    1834             :         BUN ncand;
    1835             :         const char *err;
    1836             :         lng t0 = 0;
    1837             : 
    1838             :         lng timeoffset = 0;
    1839         250 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    1840         250 :         if (qry_ctx != NULL) {
    1841         250 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    1842             :         }
    1843             : 
    1844         250 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    1845             : 
    1846         250 :         assert(tp == TYPE_dbl);
    1847             :         (void) tp;              /* compatibility (with other BATgroup*
    1848             :                                  * functions) argument */
    1849             : 
    1850         250 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    1851           0 :                 GDKerror("%s\n", err);
    1852           0 :                 return GDK_FAIL;
    1853             :         }
    1854         250 :         if (g == NULL) {
    1855           0 :                 GDKerror("b and g must be aligned\n");
    1856           0 :                 return GDK_FAIL;
    1857             :         }
    1858             : 
    1859         250 :         if (ncand == 0 || ngrp == 0) {
    1860             :                 /* trivial: no averages, so return bat aligned with g
    1861             :                  * with nil in the tail */
    1862          16 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
    1863          16 :                 if (bn == NULL) {
    1864             :                         return GDK_FAIL;
    1865             :                 }
    1866          16 :                 if (cntsp) {
    1867           9 :                         lng zero = 0;
    1868           9 :                         if ((cn = BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT)) == NULL) {
    1869           0 :                                 BBPreclaim(bn);
    1870           0 :                                 return GDK_FAIL;
    1871             :                         }
    1872           9 :                         *cntsp = cn;
    1873             :                 }
    1874          16 :                 *bnp = bn;
    1875          16 :                 return GDK_SUCCEED;
    1876             :         }
    1877             : 
    1878         234 :         if ((!skip_nils || cntsp == NULL || b->tnonil) &&
    1879         146 :             (e == NULL ||
    1880         146 :              (BATcount(e) == ncand && e->hseqbase == b->hseqbase)) &&
    1881         110 :             (BATtdense(g) || (g->tkey && g->tnonil))) {
    1882             :                 /* trivial: singleton groups, so all results are equal
    1883             :                  * to the inputs (but possibly a different type) */
    1884         110 :                 if ((bn = BATconvert(b, s, TYPE_dbl, abort_on_error, 0, 0, 0)) == NULL)
    1885             :                         return GDK_FAIL;
    1886         110 :                 if (cntsp) {
    1887           6 :                         lng one = 1;
    1888           6 :                         if ((cn = BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &one, ngrp, TRANSIENT)) == NULL) {
    1889           0 :                                 BBPreclaim(bn);
    1890           0 :                                 return GDK_FAIL;
    1891             :                         }
    1892           6 :                         *cntsp = cn;
    1893             :                 }
    1894         110 :                 *bnp = bn;
    1895         110 :                 return GDK_SUCCEED;
    1896             :         }
    1897             : 
    1898             :         /* allocate temporary space to do per group calculations */
    1899         124 :         switch (b->ttype) {
    1900         109 :         case TYPE_bte:
    1901             :         case TYPE_sht:
    1902             :         case TYPE_int:
    1903             :         case TYPE_lng:
    1904             : #ifdef HAVE_HGE
    1905             :         case TYPE_hge:
    1906             : #endif
    1907         109 :                 rems = GDKzalloc(ngrp * sizeof(lng));
    1908         109 :                 if (rems == NULL)
    1909           0 :                         goto alloc_fail;
    1910             :                 break;
    1911             :         default:
    1912             :                 break;
    1913             :         }
    1914         124 :         if (cntsp) {
    1915          94 :                 if ((cn = COLnew(min, TYPE_lng, ngrp, TRANSIENT)) == NULL)
    1916           0 :                         goto alloc_fail;
    1917          94 :                 cnts = (lng *) Tloc(cn, 0);
    1918          94 :                 memset(cnts, 0, ngrp * sizeof(lng));
    1919             :         } else {
    1920          30 :                 cnts = GDKzalloc(ngrp * sizeof(lng));
    1921          30 :                 if (cnts == NULL)
    1922           0 :                         goto alloc_fail;
    1923             :         }
    1924             : 
    1925         124 :         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    1926         124 :         if (bn == NULL)
    1927           0 :                 goto alloc_fail;
    1928         124 :         dbls = (dbl *) Tloc(bn, 0);
    1929             : 
    1930         124 :         if (BATtdense(g))
    1931             :                 gids = NULL;
    1932             :         else
    1933          92 :                 gids = (const oid *) Tloc(g, 0);
    1934             : 
    1935             :         BATiter bi;
    1936         124 :         bi = bat_iterator(b);
    1937         124 :         switch (b->ttype) {
    1938           0 :         case TYPE_bte:
    1939           0 :                 AGGR_AVG(bte);
    1940           0 :                 break;
    1941           4 :         case TYPE_sht:
    1942          17 :                 AGGR_AVG(sht);
    1943           4 :                 break;
    1944          91 :         case TYPE_int:
    1945    22861504 :                 AGGR_AVG(int);
    1946          91 :                 break;
    1947          14 :         case TYPE_lng:
    1948         107 :                 AGGR_AVG(lng);
    1949          14 :                 break;
    1950             : #ifdef HAVE_HGE
    1951           0 :         case TYPE_hge:
    1952           0 :                 AGGR_AVG(hge);
    1953           0 :                 break;
    1954             : #endif
    1955           6 :         case TYPE_flt:
    1956        2965 :                 AGGR_AVG_FLOAT(flt);
    1957             :                 break;
    1958           9 :         case TYPE_dbl:
    1959         298 :                 AGGR_AVG_FLOAT(dbl);
    1960             :                 break;
    1961           0 :         default:
    1962           0 :                 bat_iterator_end(&bi);
    1963           0 :                 GDKfree(rems);
    1964           0 :                 if (cn)
    1965           0 :                         BBPreclaim(cn);
    1966             :                 else
    1967           0 :                         GDKfree(cnts);
    1968           0 :                 BBPunfix(bn->batCacheid);
    1969           0 :                 GDKerror("type (%s) not supported.\n", ATOMname(b->ttype));
    1970           0 :                 return GDK_FAIL;
    1971             :         }
    1972         124 :         bat_iterator_end(&bi);
    1973         124 :         GDKfree(rems);
    1974         124 :         if (cn == NULL)
    1975          30 :                 GDKfree(cnts);
    1976             :         else {
    1977          94 :                 BATsetcount(cn, ngrp);
    1978          94 :                 cn->tkey = BATcount(cn) <= 1;
    1979          94 :                 cn->tsorted = BATcount(cn) <= 1;
    1980          94 :                 cn->trevsorted = BATcount(cn) <= 1;
    1981          94 :                 cn->tnil = false;
    1982          94 :                 cn->tnonil = true;
    1983          94 :                 *cntsp = cn;
    1984             :         }
    1985         124 :         if (scale != 0) {
    1986           0 :                 dbl fac = pow(10.0, (double) scale);
    1987           0 :                 for (i = 0; i < ngrp; i++) {
    1988           0 :                         if (!is_dbl_nil(dbls[i]))
    1989           0 :                                 dbls[i] *= fac;
    1990             :                 }
    1991             :         }
    1992         124 :         BATsetcount(bn, ngrp);
    1993         124 :         bn->tkey = BATcount(bn) <= 1;
    1994         124 :         bn->tsorted = BATcount(bn) <= 1;
    1995         124 :         bn->trevsorted = BATcount(bn) <= 1;
    1996         124 :         bn->tnil = nils != 0;
    1997         124 :         bn->tnonil = nils == 0;
    1998         124 :         *bnp = bn;
    1999         124 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    2000             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    2001             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    2002             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    2003             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    2004             :                   ci.seq, ncand, GDKusec() - t0);
    2005             :         return GDK_SUCCEED;
    2006             : 
    2007           0 :   alloc_fail:
    2008           0 :         if (bn)
    2009           0 :                 BBPunfix(bn->batCacheid);
    2010           0 :         GDKfree(rems);
    2011           0 :         if (cntsp) {
    2012           0 :                 BBPreclaim(*cntsp);
    2013           0 :                 *cntsp = NULL;
    2014           0 :         } else if (cnts) {
    2015           0 :                 GDKfree(cnts);
    2016             :         }
    2017             :         return GDK_FAIL;
    2018             : }
    2019             : 
    2020             : /* An exact numeric average of a bunch of values consists of three parts: the
    2021             :  * average rounded down (towards minus infinity), the number of values that
    2022             :  * participated in the calculation, and the remainder.  The remainder is in the
    2023             :  * range 0 (inclusive) to count (not inclusive).  BATgroupavg3 calculates these
    2024             :  * values for each given group.  The function below, BATgroupavg3combine,
    2025             :  * combines averages calculated this way to correct, rounded or truncated
    2026             :  * towards zero (depending on the symbol TRUNCATE_NUMBERS) averages. */
    2027             : gdk_return
    2028         693 : BATgroupavg3(BAT **avgp, BAT **remp, BAT **cntp, BAT *b, BAT *g, BAT *e, BAT *s, bool skip_nils)
    2029             : {
    2030             :         const char *err;
    2031             :         oid min, max;
    2032             :         BUN ngrp;
    2033             :         struct canditer ci;
    2034             :         BUN ncand;
    2035             :         BAT *bn, *rn, *cn;
    2036             :         BUN i;
    2037             :         oid o;
    2038             : 
    2039         693 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    2040           0 :                 GDKerror("%s\n", err);
    2041           0 :                 return GDK_FAIL;
    2042             :         }
    2043         693 :         if (ncand == 0 || ngrp == 0) {
    2044         100 :                 if (ngrp == 0)
    2045          76 :                         min = 0;
    2046         100 :                 bn = BATconstant(min, b->ttype, ATOMnilptr(b->ttype),
    2047             :                                  ngrp, TRANSIENT);
    2048         100 :                 rn = BATconstant(min, TYPE_lng, &lng_nil, ngrp, TRANSIENT);
    2049         100 :                 cn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
    2050         100 :                 if (bn == NULL || rn == NULL || cn == NULL) {
    2051           0 :                         BBPreclaim(bn);
    2052           0 :                         BBPreclaim(rn);
    2053           0 :                         BBPreclaim(cn);
    2054           0 :                         return GDK_FAIL;
    2055             :                 }
    2056         100 :                 *avgp = bn;
    2057         100 :                 *remp = rn;
    2058         100 :                 *cntp = cn;
    2059         100 :                 return GDK_SUCCEED;
    2060             :         }
    2061             :         ValRecord zero;
    2062         593 :         (void) VALinit(&zero, TYPE_bte, &(bte){0});
    2063         593 :         bn = BATconstant(min, b->ttype, VALconvert(b->ttype, &zero),
    2064             :                          ngrp, TRANSIENT);
    2065         593 :         rn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
    2066         593 :         cn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
    2067         593 :         if (bn == NULL || rn == NULL || cn == NULL) {
    2068           0 :                 BBPreclaim(bn);
    2069           0 :                 BBPreclaim(rn);
    2070           0 :                 BBPreclaim(cn);
    2071           0 :                 return GDK_FAIL;
    2072             :         }
    2073         593 :         lng *rems = Tloc(rn, 0);
    2074         593 :         lng *cnts = Tloc(cn, 0);
    2075         593 :         const oid *gids = g && !BATtdense(g) ? Tloc(g, 0) : NULL;
    2076         593 :         oid gid = ngrp == 1 && gids ? gids[0] - min : 0;
    2077             : 
    2078         593 :         BATiter bi = bat_iterator(b);
    2079             : 
    2080        1186 :         switch (ATOMbasetype(b->ttype)) {
    2081           8 :         case TYPE_bte: {
    2082           8 :                 const bte *vals = (const bte *) bi.base;
    2083           8 :                 bte *avgs = Tloc(bn, 0);
    2084         260 :                 for (i = 0; i < ncand; i++) {
    2085         252 :                         o = canditer_next(&ci) - b->hseqbase;
    2086         252 :                         if (ngrp > 1)
    2087           0 :                                 gid = gids ? gids[o] - min : o;
    2088         252 :                         if (is_bte_nil(vals[o])) {
    2089           0 :                                 if (!skip_nils) {
    2090           0 :                                         avgs[gid] = bte_nil;
    2091           0 :                                         rems[gid] = lng_nil;
    2092           0 :                                         cnts[gid] = lng_nil;
    2093           0 :                                         bn->tnil = true;
    2094           0 :                                         rn->tnil = true;
    2095           0 :                                         cn->tnil = true;
    2096             :                                 }
    2097         252 :                         } else if (!is_lng_nil(cnts[gid])) {
    2098         252 :                                 AVERAGE_ITER(bte, vals[o], avgs[gid], rems[gid], cnts[gid]);
    2099             :                         }
    2100             :                 }
    2101          16 :                 for (i = 0; i < ngrp; i++) {
    2102           8 :                         if (cnts[i] == 0) {
    2103           0 :                                 avgs[i] = bte_nil;
    2104           0 :                                 bn->tnil = true;
    2105             :                         } else
    2106             : #ifdef TRUNCATE_NUMBERS
    2107             :                         if (!is_bte_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
    2108             :                                 avgs[i]++;
    2109             :                                 rems[i] -= cnts[i];
    2110             :                         }
    2111             : #else
    2112           8 :                         if (!is_bte_nil(avgs[i]) && rems[i] > 0) {
    2113           6 :                                 if (avgs[i] < 0) {
    2114           0 :                                         if (2*rems[i] > cnts[i]) {
    2115           0 :                                                 avgs[i]++;
    2116           0 :                                                 rems[i] -= cnts[i];
    2117             :                                         }
    2118             :                                 } else {
    2119           6 :                                         if (2*rems[i] >= cnts[i]) {
    2120           2 :                                                 avgs[i]++;
    2121           2 :                                                 rems[i] -= cnts[i];
    2122             :                                         }
    2123             :                                 }
    2124             :                         }
    2125             : #endif
    2126             :                 }
    2127             :                 break;
    2128             :         }
    2129           1 :         case TYPE_sht: {
    2130           1 :                 const sht *vals = (const sht *) bi.base;
    2131           1 :                 sht *avgs = Tloc(bn, 0);
    2132          20 :                 for (i = 0; i < ncand; i++) {
    2133          19 :                         o = canditer_next(&ci) - b->hseqbase;
    2134          19 :                         if (ngrp > 1)
    2135           0 :                                 gid = gids ? gids[o] - min : o;
    2136          19 :                         if (is_sht_nil(vals[o])) {
    2137           0 :                                 if (!skip_nils) {
    2138           0 :                                         avgs[gid] = sht_nil;
    2139           0 :                                         rems[gid] = lng_nil;
    2140           0 :                                         cnts[gid] = lng_nil;
    2141           0 :                                         bn->tnil = true;
    2142           0 :                                         rn->tnil = true;
    2143           0 :                                         cn->tnil = true;
    2144             :                                 }
    2145          19 :                         } else if (!is_lng_nil(cnts[gid])) {
    2146          19 :                                 AVERAGE_ITER(sht, vals[o], avgs[gid], rems[gid], cnts[gid]);
    2147             :                         }
    2148             :                 }
    2149           2 :                 for (i = 0; i < ngrp; i++) {
    2150           1 :                         if (cnts[i] == 0) {
    2151           0 :                                 avgs[i] = sht_nil;
    2152           0 :                                 bn->tnil = true;
    2153             :                         } else
    2154             : #ifdef TRUNCATE_NUMBERS
    2155             :                         if (!is_sht_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
    2156             :                                 avgs[i]++;
    2157             :                                 rems[i] -= cnts[i];
    2158             :                         }
    2159             : #else
    2160           1 :                         if (!is_sht_nil(avgs[i]) && rems[i] > 0) {
    2161           1 :                                 if (avgs[i] < 0) {
    2162           0 :                                         if (2*rems[i] > cnts[i]) {
    2163           0 :                                                 avgs[i]++;
    2164           0 :                                                 rems[i] -= cnts[i];
    2165             :                                         }
    2166             :                                 } else {
    2167           1 :                                         if (2*rems[i] >= cnts[i]) {
    2168           0 :                                                 avgs[i]++;
    2169           0 :                                                 rems[i] -= cnts[i];
    2170             :                                         }
    2171             :                                 }
    2172             :                         }
    2173             : #endif
    2174             :                 }
    2175             :                 break;
    2176             :         }
    2177         481 :         case TYPE_int: {
    2178         481 :                 const int *vals = (const int *) bi.base;
    2179         481 :                 int *avgs = Tloc(bn, 0);
    2180     4992221 :                 for (i = 0; i < ncand; i++) {
    2181     4991742 :                         o = canditer_next(&ci) - b->hseqbase;
    2182     4991740 :                         if (ngrp > 1)
    2183     1004109 :                                 gid = gids ? gids[o] - min : o;
    2184     4991740 :                         if (is_int_nil(vals[o])) {
    2185      157482 :                                 if (!skip_nils) {
    2186           0 :                                         avgs[gid] = int_nil;
    2187           0 :                                         rems[gid] = lng_nil;
    2188           0 :                                         cnts[gid] = lng_nil;
    2189           0 :                                         bn->tnil = true;
    2190           0 :                                         rn->tnil = true;
    2191           0 :                                         cn->tnil = true;
    2192             :                                 }
    2193     4834258 :                         } else if (!is_lng_nil(cnts[gid])) {
    2194     4863414 :                                 AVERAGE_ITER(int, vals[o], avgs[gid], rems[gid], cnts[gid]);
    2195             :                         }
    2196             :                 }
    2197      511370 :                 for (i = 0; i < ngrp; i++) {
    2198      510891 :                         if (cnts[i] == 0) {
    2199        4776 :                                 avgs[i] = int_nil;
    2200        4776 :                                 bn->tnil = true;
    2201             :                         } else
    2202             : #ifdef TRUNCATE_NUMBERS
    2203             :                         if (!is_int_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
    2204             :                                 avgs[i]++;
    2205             :                                 rems[i] -= cnts[i];
    2206             :                         }
    2207             : #else
    2208      506115 :                         if (!is_int_nil(avgs[i]) && rems[i] > 0) {
    2209      147837 :                                 if (avgs[i] < 0) {
    2210      115183 :                                         if (2*rems[i] > cnts[i]) {
    2211       33468 :                                                 avgs[i]++;
    2212       33468 :                                                 rems[i] -= cnts[i];
    2213             :                                         }
    2214             :                                 } else {
    2215       32654 :                                         if (2*rems[i] >= cnts[i]) {
    2216       24991 :                                                 avgs[i]++;
    2217       24991 :                                                 rems[i] -= cnts[i];
    2218             :                                         }
    2219             :                                 }
    2220             :                         }
    2221             : #endif
    2222             :                 }
    2223             :                 break;
    2224             :         }
    2225          92 :         case TYPE_lng: {
    2226          92 :                 const lng *vals = (const lng *) bi.base;
    2227          92 :                 lng *avgs = Tloc(bn, 0);
    2228      210072 :                 for (i = 0; i < ncand; i++) {
    2229      209980 :                         o = canditer_next(&ci) - b->hseqbase;
    2230      209980 :                         if (ngrp > 1)
    2231      188326 :                                 gid = gids ? gids[o] - min : o;
    2232      209980 :                         if (is_lng_nil(vals[o])) {
    2233         210 :                                 if (!skip_nils) {
    2234           0 :                                         avgs[gid] = lng_nil;
    2235           0 :                                         rems[gid] = lng_nil;
    2236           0 :                                         cnts[gid] = lng_nil;
    2237           0 :                                         bn->tnil = true;
    2238           0 :                                         rn->tnil = true;
    2239           0 :                                         cn->tnil = true;
    2240             :                                 }
    2241      209770 :                         } else if (!is_lng_nil(cnts[gid])) {
    2242      210300 :                                 AVERAGE_ITER(lng, vals[o], avgs[gid], rems[gid], cnts[gid]);
    2243             :                         }
    2244             :                 }
    2245       69285 :                 for (i = 0; i < ngrp; i++) {
    2246       69193 :                         if (cnts[i] == 0) {
    2247         168 :                                 avgs[i] = lng_nil;
    2248         168 :                                 bn->tnil = true;
    2249             :                         } else
    2250             : #ifdef TRUNCATE_NUMBERS
    2251             :                         if (!is_lng_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
    2252             :                                 avgs[i]++;
    2253             :                                 rems[i] -= cnts[i];
    2254             :                         }
    2255             : #else
    2256       69025 :                         if (!is_lng_nil(avgs[i]) && rems[i] > 0) {
    2257         922 :                                 if (avgs[i] < 0) {
    2258         129 :                                         if (2*rems[i] > cnts[i]) {
    2259           0 :                                                 avgs[i]++;
    2260           0 :                                                 rems[i] -= cnts[i];
    2261             :                                         }
    2262             :                                 } else {
    2263         793 :                                         if (2*rems[i] >= cnts[i]) {
    2264         694 :                                                 avgs[i]++;
    2265         694 :                                                 rems[i] -= cnts[i];
    2266             :                                         }
    2267             :                                 }
    2268             :                         }
    2269             : #endif
    2270             :                 }
    2271             :                 break;
    2272             :         }
    2273             : #ifdef HAVE_HGE
    2274          11 :         case TYPE_hge: {
    2275          11 :                 const hge *vals = (const hge *) bi.base;
    2276          11 :                 hge *avgs = Tloc(bn, 0);
    2277    11896076 :                 for (i = 0; i < ncand; i++) {
    2278    11896065 :                         o = canditer_next(&ci) - b->hseqbase;
    2279    11896065 :                         if (ngrp > 1)
    2280      150707 :                                 gid = gids ? gids[o] - min : o;
    2281    11896065 :                         if (is_hge_nil(vals[o])) {
    2282      248594 :                                 if (!skip_nils) {
    2283           0 :                                         avgs[gid] = hge_nil;
    2284           0 :                                         rems[gid] = lng_nil;
    2285           0 :                                         cnts[gid] = lng_nil;
    2286           0 :                                         bn->tnil = true;
    2287           0 :                                         rn->tnil = true;
    2288           0 :                                         cn->tnil = true;
    2289             :                                 }
    2290    11647471 :                         } else if (!is_lng_nil(cnts[gid])) {
    2291    11647471 :                                 AVERAGE_ITER(hge, vals[o], avgs[gid], rems[gid], cnts[gid]);
    2292             :                         }
    2293             :                 }
    2294         135 :                 for (i = 0; i < ngrp; i++) {
    2295         124 :                         if (cnts[i] == 0) {
    2296           0 :                                 avgs[i] = hge_nil;
    2297           0 :                                 bn->tnil = true;
    2298             :                         } else
    2299             : #ifdef TRUNCATE_NUMBERS
    2300             :                         if (!is_hge_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0) {
    2301             :                                 avgs[i]++;
    2302             :                                 rems[i] -= cnts[i];
    2303             :                         }
    2304             : #else
    2305         124 :                         if (!is_hge_nil(avgs[i]) && rems[i] > 0) {
    2306         119 :                                 if (avgs[i] < 0) {
    2307           0 :                                         if (2*rems[i] > cnts[i]) {
    2308           0 :                                                 avgs[i]++;
    2309           0 :                                                 rems[i] -= cnts[i];
    2310             :                                         }
    2311             :                                 } else {
    2312         119 :                                         if (2*rems[i] >= cnts[i]) {
    2313          52 :                                                 avgs[i]++;
    2314          52 :                                                 rems[i] -= cnts[i];
    2315             :                                         }
    2316             :                                 }
    2317             :                         }
    2318             : #endif
    2319             :                 }
    2320             :                 break;
    2321             :         }
    2322             : #endif
    2323             :         }
    2324         591 :         bat_iterator_end(&bi);
    2325         593 :         BATsetcount(bn, ngrp);
    2326         593 :         BATsetcount(rn, ngrp);
    2327         593 :         BATsetcount(cn, ngrp);
    2328         593 :         bn->tnonil = !bn->tnil;
    2329         593 :         rn->tnonil = !rn->tnil;
    2330         593 :         cn->tnonil = !cn->tnil;
    2331         593 :         bn->tkey = rn->tkey = cn->tkey = ngrp == 1;
    2332         593 :         bn->tsorted = rn->tsorted = cn->tsorted = ngrp == 1;
    2333         593 :         bn->trevsorted = rn->trevsorted = cn->trevsorted = ngrp == 1;
    2334         593 :         *avgp = bn;
    2335         593 :         *remp = rn;
    2336         593 :         *cntp = cn;
    2337         593 :         return GDK_SUCCEED;
    2338             : }
    2339             : 
    2340             : #ifdef HAVE_HGE
    2341             : #define BIGINT uhge
    2342             : #else
    2343             : #define BIGINT uint64_t
    2344             : #endif
    2345             : /* calculate the values (a*b)/c and (a*b)%c but without causing overflow
    2346             :  *
    2347             :  * this function is from https://stackoverflow.com/a/8757419, only with the
    2348             :  * slight addition of returning the remainder and changing the first argument
    2349             :  * and return value to BIGINT */
    2350             : static inline BIGINT
    2351           0 : multdiv(BIGINT a, uint64_t b, uint64_t c, uint64_t *rem)
    2352             : {
    2353             :         static BIGINT const base = ((BIGINT)1)<<(sizeof(BIGINT)*4);
    2354             : //      static BIGINT const maxdiv = (base-1)*base + (base-1);
    2355             :         static BIGINT const maxdiv = ~(BIGINT)0;
    2356             : 
    2357             :         // First get the easy thing
    2358           0 :         BIGINT res = (a/c) * b + (a%c) * (b/c);
    2359             :         a %= c;
    2360           0 :         b %= c;
    2361             :         // Are we done?
    2362           0 :         if (a == 0 || b == 0) {
    2363           0 :                 *rem = 0;
    2364           0 :                 return res;
    2365             :         }
    2366             :         // Is it easy to compute what remain to be added?
    2367             :         if (c < base) {
    2368           0 :                 *rem = (uint64_t) ((a*b)%c);
    2369           0 :                 return res + (a*b/c);
    2370             :         }
    2371             :         // Now 0 < a < c, 0 < b < c, c >= 1ULL
    2372             :         // Normalize
    2373             :         BIGINT norm = maxdiv/c;
    2374             :         BIGINT cbig = c * norm; // orig: c *= norm; and below cbig was plain c
    2375             :         a *= norm;
    2376             :         // split into 2 digits
    2377             :         BIGINT ah = a / base, al = a % base;
    2378             :         BIGINT bh = b / base, bl = b % base;
    2379             :         BIGINT ch = cbig / base, cl = cbig % base;
    2380             :         // compute the product
    2381             :         BIGINT p0 = al*bl;
    2382             :         BIGINT p1 = p0 / base + al*bh;
    2383             :         p0 %= base;
    2384             :         BIGINT p2 = p1 / base + ah*bh;
    2385             :         p1 = (p1 % base) + ah * bl;
    2386             :         p2 += p1 / base;
    2387             :         p1 %= base;
    2388             :         // p2 holds 2 digits, p1 and p0 one
    2389             : 
    2390             :         // first digit is easy, not null only in case of overflow
    2391             :         //BIGINT q2 = p2 / cbig;
    2392             :         p2 = p2 % cbig;
    2393             : 
    2394             :         // second digit, estimate
    2395             :         BIGINT q1 = p2 / ch;
    2396             :         // and now adjust
    2397             :         BIGINT rhat = p2 % ch;
    2398             :         // the loop can be unrolled, it will be executed at most twice for
    2399             :         // even bases -- three times for odd one -- due to the normalisation above
    2400             :         while (q1 >= base || (rhat < base && q1*cl > rhat*base+p1)) {
    2401             :                 q1--;
    2402             :                 rhat += ch;
    2403             :         }
    2404             :         // subtract
    2405             :         p1 = ((p2 % base) * base + p1) - q1 * cl;
    2406             :         p2 = (p2 / base * base + p1 / base) - q1 * ch;
    2407             :         p1 = p1 % base + (p2 % base) * base;
    2408             : 
    2409             :         // now p1 hold 2 digits, p0 one and p2 is to be ignored
    2410             :         BIGINT q0 = p1 / ch;
    2411             :         rhat = p1 % ch;
    2412             :         while (q0 >= base || (rhat < base && q0*cl > rhat*base+p0)) {
    2413             :                 q0--;
    2414             :                 rhat += ch;
    2415             :         }
    2416             :         // subtract
    2417             :         p0 = ((p1 % base) * base + p0) - q0 * cl;
    2418             :         p1 = (p1 / base * base + p0 / base) - q0 * ch;
    2419             :         p0 = p0 % base + (p1 % base) * base;
    2420             : 
    2421             :         *rem = p0 / norm;
    2422             :         return res + q0 + q1 * base; // + q2 *base*base
    2423             : }
    2424             : 
    2425             : static inline void
    2426           8 : combine_averages_bte(bte *avgp, lng *remp, lng *cntp,
    2427             :                      bte avg1, lng rem1, lng cnt1,
    2428             :                      bte avg2, lng rem2, lng cnt2)
    2429             : {
    2430           8 :         lng cnt = cnt1 + cnt2;
    2431             : 
    2432           8 :         if (rem2 < 0) {
    2433           2 :                 avg2--;
    2434           2 :                 rem2 += cnt2;
    2435             :         }
    2436           8 :         *cntp = cnt;
    2437           8 :         lng v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
    2438           8 :         bte a = (bte) (v / cnt);
    2439           8 :         v %= cnt;
    2440           8 :         if (v < 0) {
    2441           0 :                 a--;
    2442           0 :                 v += cnt;
    2443             :         }
    2444           8 :         *avgp = a;
    2445           8 :         *remp = v;
    2446           8 : }
    2447             : 
    2448             : static inline void
    2449           1 : combine_averages_sht(sht *avgp, lng *remp, lng *cntp,
    2450             :                      sht avg1, lng rem1, lng cnt1,
    2451             :                      sht avg2, lng rem2, lng cnt2)
    2452             : {
    2453           1 :         lng cnt = cnt1 + cnt2;
    2454             : 
    2455           1 :         if (rem2 < 0) {
    2456           0 :                 avg2--;
    2457           0 :                 rem2 += cnt2;
    2458             :         }
    2459           1 :         *cntp = cnt;
    2460           1 :         lng v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
    2461           1 :         sht a = (sht) (v / cnt);
    2462           1 :         v %= cnt;
    2463           1 :         if (v < 0) {
    2464           0 :                 a--;
    2465           0 :                 v += cnt;
    2466             :         }
    2467           1 :         *avgp = a;
    2468           1 :         *remp = v;
    2469           1 : }
    2470             : 
    2471             : 
    2472             : static inline void
    2473      245886 : combine_averages_int(int *avgp, lng *remp, lng *cntp,
    2474             :                      int avg1, lng rem1, lng cnt1,
    2475             :                      int avg2, lng rem2, lng cnt2)
    2476             : {
    2477      245886 :         lng cnt = cnt1 + cnt2;
    2478             : 
    2479      245886 :         if (rem2 < 0) {
    2480       25904 :                 avg2--;
    2481       25904 :                 rem2 += cnt2;
    2482             :         }
    2483      245886 :         *cntp = cnt;
    2484             : #ifdef HAVE_HGE
    2485      245886 :         hge v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
    2486      245886 :         int a = (int) (v / cnt);
    2487      245886 :         v %= cnt;
    2488      245886 :         if (v < 0) {
    2489      168565 :                 a--;
    2490      168565 :                 v += cnt;
    2491             :         }
    2492      245886 :         *avgp = a;
    2493      245886 :         *remp = (lng) v;
    2494             : #else
    2495             :         if (cnt1 == 0) {
    2496             :                 avg1 = 0;
    2497             :                 rem1 = 0;
    2498             :         }
    2499             :         lng rem = rem1 + rem2;
    2500             :         lng v;
    2501             :         uint64_t r;
    2502             :         if (avg1 < 0) {
    2503             :                 avg1 = (int) multdiv((uint64_t) -avg1, cnt1, cnt, &r);
    2504             :                 if (r > 0) {
    2505             :                         avg1 = -avg1 - 1;
    2506             :                         r = cnt - r;
    2507             :                 } else {
    2508             :                         avg1 = -avg1;
    2509             :                 }
    2510             :         } else {
    2511             :                 avg1 = (int) multdiv((uint64_t) avg1, cnt1, cnt, &r);
    2512             :         }
    2513             :         v = avg1;
    2514             :         rem += r;
    2515             :         if (avg2 < 0) {
    2516             :                 avg2 = (int) multdiv((uint64_t) -avg2, cnt2, cnt, &r);
    2517             :                 if (r > 0) {
    2518             :                         avg2 = -avg2 - 1;
    2519             :                         r = cnt - r;
    2520             :                 } else {
    2521             :                         avg2 = -avg2;
    2522             :                 }
    2523             :         } else {
    2524             :                 avg2 = (int) multdiv((uint64_t) avg2, cnt2, cnt, &r);
    2525             :         }
    2526             :         v += avg2;
    2527             :         rem += r;
    2528             :         while (rem >= cnt) { /* max twice */
    2529             :                 v++;
    2530             :                 rem -= cnt;
    2531             :         }
    2532             :         *avgp = (int) v;
    2533             :         *remp = rem;
    2534             : #endif
    2535      245886 : }
    2536             : 
    2537             : static inline void
    2538          52 : combine_averages_lng(lng *avgp, lng *remp, lng *cntp,
    2539             :                      lng avg1, lng rem1, lng cnt1,
    2540             :                      lng avg2, lng rem2, lng cnt2)
    2541             : {
    2542          52 :         lng cnt = cnt1 + cnt2;
    2543             : 
    2544          52 :         if (rem2 < 0) {
    2545          24 :                 avg2--;
    2546          24 :                 rem2 += cnt2;
    2547             :         }
    2548          52 :         *cntp = cnt;
    2549             : #ifdef HAVE_HGE
    2550          52 :         hge v = avg1 * cnt1 + rem1 + avg2 * cnt2 + rem2;
    2551          52 :         lng a = (lng) (v / cnt);
    2552          52 :         v %= cnt;
    2553          52 :         if (v < 0) {
    2554           0 :                 a--;
    2555           0 :                 v += cnt;
    2556             :         }
    2557          52 :         *avgp = a;
    2558          52 :         *remp = (lng) v;
    2559             : #else
    2560             :         if (cnt1 == 0) {
    2561             :                 avg1 = 0;
    2562             :                 rem1 = 0;
    2563             :         }
    2564             :         lng rem = rem1 + rem2;
    2565             :         lng v;
    2566             :         uint64_t r;
    2567             :         if (avg1 < 0) {
    2568             :                 avg1 = (lng) multdiv((uint64_t) -avg1, cnt1, cnt, &r);
    2569             :                 if (r > 0) {
    2570             :                         avg1 = -avg1 - 1;
    2571             :                         r = cnt - r;
    2572             :                 } else {
    2573             :                         avg1 = -avg1;
    2574             :                 }
    2575             :         } else {
    2576             :                 avg1 = (lng) multdiv((uint64_t) avg1, cnt1, cnt, &r);
    2577             :         }
    2578             :         v = avg1;
    2579             :         rem += r;
    2580             :         if (avg2 < 0) {
    2581             :                 avg2 = (lng) multdiv((uint64_t) -avg2, cnt2, cnt, &r);
    2582             :                 if (r > 0) {
    2583             :                         avg2 = -avg2 - 1;
    2584             :                         r = cnt - r;
    2585             :                 } else {
    2586             :                         avg2 = -avg2;
    2587             :                 }
    2588             :         } else {
    2589             :                 avg2 = (lng) multdiv((uint64_t) avg2, cnt2, cnt, &r);
    2590             :         }
    2591             :         v += avg2;
    2592             :         rem += r;
    2593             :         while (rem >= cnt) { /* max twice */
    2594             :                 v++;
    2595             :                 rem -= cnt;
    2596             :         }
    2597             :         *avgp = v;
    2598             :         *remp = rem;
    2599             : #endif
    2600          52 : }
    2601             : 
    2602             : #ifdef HAVE_HGE
    2603             : static inline void
    2604           0 : combine_averages_hge(hge *avgp, lng *remp, lng *cntp,
    2605             :                      hge avg1, lng rem1, lng cnt1,
    2606             :                      hge avg2, lng rem2, lng cnt2)
    2607             : {
    2608           0 :         if (cnt1 == 0) {
    2609             :                 avg1 = 0;
    2610             :                 rem1 = 0;
    2611             :         }
    2612           0 :         if (rem2 < 0) {
    2613           0 :                 avg2--;
    2614           0 :                 rem2 += cnt2;
    2615             :         }
    2616           0 :         lng cnt = cnt1 + cnt2;
    2617           0 :         lng rem = rem1 + rem2;
    2618             :         hge v;
    2619             :         uint64_t r;
    2620             : 
    2621           0 :         *cntp = cnt;
    2622           0 :         if (avg1 < 0) {
    2623           0 :                 avg1 = (hge) multdiv((uhge) -avg1, cnt1, cnt, &r);
    2624           0 :                 if (r > 0) {
    2625           0 :                         avg1 = -avg1 - 1;
    2626           0 :                         r = cnt - r;
    2627             :                 } else {
    2628           0 :                         avg1 = -avg1;
    2629             :                 }
    2630             :         } else {
    2631           0 :                 avg1 = (hge) multdiv((uhge) avg1, cnt1, cnt, &r);
    2632             :         }
    2633             :         v = avg1;
    2634           0 :         rem += r;
    2635           0 :         if (avg2 < 0) {
    2636           0 :                 avg2 = (hge) multdiv((uhge) -avg2, cnt2, cnt, &r);
    2637           0 :                 if (r > 0) {
    2638           0 :                         avg2 = -avg2 - 1;
    2639           0 :                         r = cnt - r;
    2640             :                 } else {
    2641           0 :                         avg2 = -avg2;
    2642             :                 }
    2643             :         } else {
    2644           0 :                 avg2 = (hge) multdiv((uhge) avg2, cnt2, cnt, &r);
    2645             :         }
    2646           0 :         v += avg2;
    2647           0 :         rem += r;
    2648           0 :         while (rem >= cnt) { /* max twice */
    2649           0 :                 v++;
    2650           0 :                 rem -= cnt;
    2651             :         }
    2652           0 :         *avgp = v;
    2653           0 :         *remp = rem;
    2654           0 : }
    2655             : #endif
    2656             : 
    2657             : BAT *
    2658          42 : BATgroupavg3combine(BAT *avg, BAT *rem, BAT *cnt, BAT *g, BAT *e, bool skip_nils)
    2659             : {
    2660             :         const char *err;
    2661             :         oid min, max;
    2662             :         BUN ngrp;
    2663             :         struct canditer ci;
    2664             :         BUN i, ncand;
    2665             :         BAT *bn, *rn, *cn;
    2666             : 
    2667          42 :         if ((err = BATgroupaggrinit(avg, g, e, NULL, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    2668           0 :                 GDKerror("%s\n", err);
    2669           0 :                 return NULL;
    2670             :         }
    2671          42 :         assert(ci.tpe == cand_dense);
    2672          42 :         if (ncand != BATcount(rem) || ncand != BATcount(cnt)) {
    2673           0 :                 GDKerror("input bats not aligned");
    2674           0 :                 return NULL;
    2675             :         }
    2676          42 :         if (ncand == 0 || ngrp == 0) {
    2677           1 :                 return BATconstant(ngrp == 0 ? 0 : min, avg->ttype,
    2678           1 :                                    ATOMnilptr(avg->ttype), ngrp, TRANSIENT);
    2679             :         }
    2680             :         ValRecord zero;
    2681          41 :         (void) VALinit(&zero, TYPE_bte, &(bte){0});
    2682          41 :         bn = BATconstant(min, avg->ttype, VALconvert(avg->ttype, &zero),
    2683             :                          ngrp, TRANSIENT);
    2684             :         /* rn and cn are temporary storage of intermediates */
    2685          41 :         rn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
    2686          41 :         cn = BATconstant(min, TYPE_lng, &(lng){0}, ngrp, TRANSIENT);
    2687          41 :         if (bn == NULL || rn == NULL || cn == NULL) {
    2688           0 :                 BBPreclaim(bn);
    2689           0 :                 BBPreclaim(rn);
    2690           0 :                 BBPreclaim(cn);
    2691           0 :                 return NULL;
    2692             :         }
    2693          41 :         lng *rems = Tloc(rn, 0);
    2694          41 :         lng *cnts = Tloc(cn, 0);
    2695          41 :         const lng *orems = Tloc(rem, 0);
    2696          41 :         const lng *ocnts = Tloc(cnt, 0);
    2697          41 :         const oid *gids = g && !BATtdense(g) ? Tloc(g, 0) : NULL;
    2698          41 :         oid gid = ngrp == 1 && gids ? gids[0] - min : 0;
    2699             : 
    2700          41 :         BATiter bi = bat_iterator(avg);
    2701             : 
    2702          82 :         switch (ATOMbasetype(avg->ttype)) {
    2703           2 :         case TYPE_bte: {
    2704           2 :                 const bte *vals = (const bte *) bi.base;
    2705           2 :                 bte *avgs = Tloc(bn, 0);
    2706          10 :                 for (i = 0; i < ncand; i++) {
    2707           8 :                         if (ngrp > 1)
    2708           0 :                                 gid = gids ? gids[i] - min : i;
    2709           8 :                         if (is_bte_nil(vals[i])) {
    2710           0 :                                 if (!skip_nils) {
    2711           0 :                                         avgs[gid] = bte_nil;
    2712           0 :                                         bn->tnil = true;
    2713             :                                 }
    2714           8 :                         } else if (!is_bte_nil(avgs[gid])) {
    2715           8 :                                 combine_averages_bte(&avgs[gid], &rems[gid],
    2716             :                                                      &cnts[gid], avgs[gid],
    2717           8 :                                                      rems[gid], cnts[gid],
    2718           8 :                                                      vals[i], orems[i],
    2719           8 :                                                      ocnts[i]);
    2720             :                         }
    2721             :                 }
    2722           4 :                 for (i = 0; i < ngrp; i++) {
    2723           2 :                         if (cnts[i] == 0) {
    2724           0 :                                 avgs[i] = bte_nil;
    2725           0 :                                 bn->tnil = true;
    2726             :                         } else
    2727             : #ifdef TRUNCATE_NUMBERS
    2728             :                         if (!is_bte_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
    2729             :                                 avgs[i]++;
    2730             : #else
    2731           2 :                         if (!is_bte_nil(avgs[i]) && rems[i] > 0) {
    2732           2 :                                 if (avgs[i] < 0) {
    2733           0 :                                         if (2*rems[i] > cnts[i])
    2734           0 :                                                 avgs[i]++;
    2735             :                                 } else {
    2736           2 :                                         if (2*rems[i] >= cnts[i])
    2737           0 :                                                 avgs[i]++;
    2738             :                                 }
    2739             :                         }
    2740             : #endif
    2741             :                 }
    2742             :                 break;
    2743             :         }
    2744           1 :         case TYPE_sht: {
    2745           1 :                 const sht *vals = (const sht *) bi.base;
    2746           1 :                 sht *avgs = Tloc(bn, 0);
    2747           2 :                 for (i = 0; i < ncand; i++) {
    2748           1 :                         if (ngrp > 1)
    2749           0 :                                 gid = gids ? gids[i] - min : i;
    2750           1 :                         if (is_sht_nil(vals[i])) {
    2751           0 :                                 if (!skip_nils) {
    2752           0 :                                         avgs[gid] = sht_nil;
    2753           0 :                                         bn->tnil = true;
    2754             :                                 }
    2755           1 :                         } else if (!is_sht_nil(avgs[gid])) {
    2756           1 :                                 combine_averages_sht(&avgs[gid], &rems[gid],
    2757             :                                                      &cnts[gid], avgs[gid],
    2758           1 :                                                      rems[gid], cnts[gid],
    2759           1 :                                                      vals[i], orems[i],
    2760           1 :                                                      ocnts[i]);
    2761             :                         }
    2762             :                 }
    2763           2 :                 for (i = 0; i < ngrp; i++) {
    2764           1 :                         if (cnts[i] == 0) {
    2765           0 :                                 avgs[i] = sht_nil;
    2766           0 :                                 bn->tnil = true;
    2767             :                         } else
    2768             : #ifdef TRUNCATE_NUMBERS
    2769             :                         if (!is_sht_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
    2770             :                                 avgs[i]++;
    2771             : #else
    2772           1 :                         if (!is_sht_nil(avgs[i]) && rems[i] > 0) {
    2773           1 :                                 if (avgs[i] < 0) {
    2774           0 :                                         if (2*rems[i] > cnts[i])
    2775           0 :                                                 avgs[i]++;
    2776             :                                 } else {
    2777           1 :                                         if (2*rems[i] >= cnts[i])
    2778           0 :                                                 avgs[i]++;
    2779             :                                 }
    2780             :                         }
    2781             : #endif
    2782             :                 }
    2783             :                 break;
    2784             :         }
    2785          34 :         case TYPE_int: {
    2786          34 :                 const int *vals = (const int *) bi.base;
    2787          34 :                 int *avgs = Tloc(bn, 0);
    2788      248283 :                 for (i = 0; i < ncand; i++) {
    2789      248249 :                         if (ngrp > 1)
    2790      247922 :                                 gid = gids ? gids[i] - min : i;
    2791      248249 :                         if (is_int_nil(vals[i])) {
    2792        2363 :                                 if (!skip_nils) {
    2793           0 :                                         avgs[gid] = int_nil;
    2794           0 :                                         bn->tnil = true;
    2795             :                                 }
    2796      245886 :                         } else if (!is_int_nil(avgs[gid])) {
    2797      245886 :                                 combine_averages_int(&avgs[gid], &rems[gid],
    2798             :                                                      &cnts[gid], avgs[gid],
    2799      245886 :                                                      rems[gid], cnts[gid],
    2800      245886 :                                                      vals[i], orems[i],
    2801      245886 :                                                      ocnts[i]);
    2802             :                         }
    2803             :                 }
    2804       61147 :                 for (i = 0; i < ngrp; i++) {
    2805       61113 :                         if (cnts[i] == 0) {
    2806          67 :                                 avgs[i] = int_nil;
    2807          67 :                                 bn->tnil = true;
    2808             :                         } else
    2809             : #ifdef TRUNCATE_NUMBERS
    2810             :                         if (!is_int_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
    2811             :                                 avgs[i]++;
    2812             : #else
    2813       61046 :                         if (!is_int_nil(avgs[i]) && rems[i] > 0) {
    2814       37417 :                                 if (avgs[i] < 0) {
    2815       33434 :                                         if (2*rems[i] > cnts[i])
    2816       16188 :                                                 avgs[i]++;
    2817             :                                 } else {
    2818        3983 :                                         if (2*rems[i] >= cnts[i])
    2819        3202 :                                                 avgs[i]++;
    2820             :                                 }
    2821             :                         }
    2822             : #endif
    2823             :                 }
    2824             :                 break;
    2825             :         }
    2826           4 :         case TYPE_lng: {
    2827           4 :                 const lng *vals = (const lng *) bi.base;
    2828           4 :                 lng *avgs = Tloc(bn, 0);
    2829          56 :                 for (i = 0; i < ncand; i++) {
    2830          52 :                         if (ngrp > 1)
    2831          48 :                                 gid = gids ? gids[i] - min : i;
    2832          52 :                         if (is_lng_nil(vals[i])) {
    2833           0 :                                 if (!skip_nils) {
    2834           0 :                                         avgs[gid] = lng_nil;
    2835           0 :                                         bn->tnil = true;
    2836             :                                 }
    2837          52 :                         } else if (!is_lng_nil(avgs[gid])) {
    2838          52 :                                 combine_averages_lng(&avgs[gid], &rems[gid],
    2839             :                                                      &cnts[gid], avgs[gid],
    2840          52 :                                                      rems[gid], cnts[gid],
    2841          52 :                                                      vals[i], orems[i],
    2842          52 :                                                      ocnts[i]);
    2843             :                         }
    2844             :                 }
    2845          17 :                 for (i = 0; i < ngrp; i++) {
    2846          13 :                         if (cnts[i] == 0) {
    2847           0 :                                 avgs[i] = lng_nil;
    2848           0 :                                 bn->tnil = true;
    2849             :                         } else
    2850             : #ifdef TRUNCATE_NUMBERS
    2851             :                         if (!is_lng_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
    2852             :                                 avgs[i]++;
    2853             : #else
    2854          13 :                         if (!is_lng_nil(avgs[i]) && rems[i] > 0) {
    2855          12 :                                 if (avgs[i] < 0) {
    2856           0 :                                         if (2*rems[i] > cnts[i])
    2857           0 :                                                 avgs[i]++;
    2858             :                                 } else {
    2859          12 :                                         if (2*rems[i] >= cnts[i])
    2860          10 :                                                 avgs[i]++;
    2861             :                                 }
    2862             :                         }
    2863             : #endif
    2864             :                 }
    2865             :                 break;
    2866             :         }
    2867             : #ifdef HAVE_HGE
    2868           0 :         case TYPE_hge: {
    2869           0 :                 const hge *vals = (const hge *) bi.base;
    2870           0 :                 hge *avgs = Tloc(bn, 0);
    2871           0 :                 for (i = 0; i < ncand; i++) {
    2872           0 :                         if (ngrp > 1)
    2873           0 :                                 gid = gids ? gids[i] - min : i;
    2874           0 :                         if (is_hge_nil(vals[i])) {
    2875           0 :                                 if (!skip_nils) {
    2876           0 :                                         avgs[gid] = hge_nil;
    2877           0 :                                         bn->tnil = true;
    2878             :                                 }
    2879           0 :                         } else if (!is_hge_nil(avgs[gid])) {
    2880           0 :                                 combine_averages_hge(&avgs[gid], &rems[gid],
    2881             :                                                      &cnts[gid], avgs[gid],
    2882           0 :                                                      rems[gid], cnts[gid],
    2883           0 :                                                      vals[i], orems[i],
    2884           0 :                                                      ocnts[i]);
    2885             :                         }
    2886             :                 }
    2887           0 :                 for (i = 0; i < ngrp; i++) {
    2888           0 :                         if (cnts[i] == 0) {
    2889           0 :                                 avgs[i] = hge_nil;
    2890           0 :                                 bn->tnil = true;
    2891             :                         } else
    2892             : #ifdef TRUNCATE_NUMBERS
    2893             :                         if (!is_hge_nil(avgs[i]) && rems[i] > 0 && avgs[i] < 0)
    2894             :                                 avgs[i]++;
    2895             : #else
    2896           0 :                         if (!is_hge_nil(avgs[i]) && rems[i] > 0) {
    2897           0 :                                 if (avgs[i] < 0) {
    2898           0 :                                         if (2*rems[i] > cnts[i])
    2899           0 :                                                 avgs[i]++;
    2900             :                                 } else {
    2901           0 :                                         if (2*rems[i] >= cnts[i])
    2902           0 :                                                 avgs[i]++;
    2903             :                                 }
    2904             :                         }
    2905             : #endif
    2906             :                 }
    2907             :                 break;
    2908             :         }
    2909             : #endif
    2910             :         }
    2911          41 :         bat_iterator_end(&bi);
    2912          41 :         BBPreclaim(rn);
    2913          41 :         BBPreclaim(cn);
    2914          41 :         BATsetcount(bn, ngrp);
    2915          41 :         bn->tnonil = !bn->tnil;
    2916          41 :         bn->tkey = ngrp == 1;
    2917          41 :         bn->tsorted = ngrp == 1;
    2918          41 :         bn->trevsorted = ngrp == 1;
    2919          41 :         return bn;
    2920             : }
    2921             : 
    2922             : #define AVERAGE_TYPE_LNG_HGE(TYPE,lng_hge)                              \
    2923             :         do {                                                            \
    2924             :                 TYPE x, a;                                              \
    2925             :                                                                         \
    2926             :                 /* first try to calculate the sum of all values into a */ \
    2927             :                 /* lng/hge */                                           \
    2928             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    2929             :                         i = canditer_next(&ci) - b->hseqbase;            \
    2930             :                         x = ((const TYPE *) src)[i];                    \
    2931             :                         if (is_##TYPE##_nil(x))                         \
    2932             :                                 continue;                               \
    2933             :                         ADD_WITH_CHECK(x, sum,                          \
    2934             :                                        lng_hge, sum,                    \
    2935             :                                        GDK_##lng_hge##_max,             \
    2936             :                                        goto overflow##TYPE);            \
    2937             :                         /* don't count value until after overflow check */ \
    2938             :                         n++;                                            \
    2939             :                 }                                                       \
    2940             :                 TIMEOUT_CHECK(timeoffset,                               \
    2941             :                               TIMEOUT_HANDLER(GDK_FAIL));               \
    2942             :                 /* the sum fit, so now we can calculate the average */  \
    2943             :                 *avg = n > 0 ? (dbl) sum / n : dbl_nil;                      \
    2944             :                 if (0) {                                                \
    2945             :                   overflow##TYPE:                                       \
    2946             :                         /* we get here if sum(x[0],...,x[i]) doesn't */ \
    2947             :                         /* fit in a lng/hge but sum(x[0],...,x[i-1]) did */ \
    2948             :                         /* the variable sum contains that sum */        \
    2949             :                         /* the rest of the calculation is done */       \
    2950             :                         /* according to the loop invariant described */ \
    2951             :                         /* in the below loop */                         \
    2952             :                         /* note that n necessarily is > 0 (else no */        \
    2953             :                         /* overflow possible) */                        \
    2954             :                         assert(n > 0);                                       \
    2955             :                         if (sum >= 0) {                                      \
    2956             :                                 a = (TYPE) (sum / n); /* this fits */   \
    2957             :                                 r = (lng) (sum % n);                    \
    2958             :                         } else {                                        \
    2959             :                                 sum = -sum;                             \
    2960             :                                 a = - (TYPE) (sum / n); /* this fits */ \
    2961             :                                 r = (lng) (sum % n);                    \
    2962             :                                 if (r) {                                \
    2963             :                                         a--;                            \
    2964             :                                         r = n - r;                      \
    2965             :                                 }                                       \
    2966             :                         }                                               \
    2967             :                         while (ncand > 0) {                          \
    2968             :                                 /* loop invariant: */                   \
    2969             :                                 /* a + r/n == average(x[0],...,x[n]); */ \
    2970             :                                 /* 0 <= r < n */                  \
    2971             :                                 ncand--;                                \
    2972             :                                 i = canditer_next(&ci) - b->hseqbase;    \
    2973             :                                 x = ((const TYPE *) src)[i];            \
    2974             :                                 if (is_##TYPE##_nil(x))                 \
    2975             :                                         continue;                       \
    2976             :                                 AVERAGE_ITER(TYPE, x, a, r, n);         \
    2977             :                         }                                               \
    2978             :                         *avg = a + (dbl) r / n;                         \
    2979             :                 }                                                       \
    2980             :         } while (0)
    2981             : 
    2982             : #ifdef HAVE_HGE
    2983             : #define AVERAGE_TYPE(TYPE) AVERAGE_TYPE_LNG_HGE(TYPE,hge)
    2984             : #else
    2985             : #define AVERAGE_TYPE(TYPE) AVERAGE_TYPE_LNG_HGE(TYPE,lng)
    2986             : #endif
    2987             : 
    2988             : #define AVERAGE_FLOATTYPE(TYPE)                                 \
    2989             :         do {                                                    \
    2990             :                 double a = 0;                                   \
    2991             :                 TYPE x;                                         \
    2992             :                 TIMEOUT_LOOP(ncand, timeoffset) {               \
    2993             :                         i = canditer_next(&ci) - b->hseqbase;    \
    2994             :                         x = ((const TYPE *) src)[i];            \
    2995             :                         if (is_##TYPE##_nil(x))                 \
    2996             :                                 continue;                       \
    2997             :                         AVERAGE_ITER_FLOAT(TYPE, x, a, n);      \
    2998             :                 }                                               \
    2999             :                 TIMEOUT_CHECK(timeoffset,                       \
    3000             :                               TIMEOUT_HANDLER(GDK_FAIL));       \
    3001             :                 *avg = n > 0 ? a : dbl_nil;                  \
    3002             :         } while (0)
    3003             : 
    3004             : gdk_return
    3005        4047 : BATcalcavg(BAT *b, BAT *s, dbl *avg, BUN *vals, int scale)
    3006             : {
    3007             :         lng n = 0, r = 0;
    3008             :         BUN i = 0;
    3009             : #ifdef HAVE_HGE
    3010             :         hge sum = 0;
    3011             : #else
    3012             :         lng sum = 0;
    3013             : #endif
    3014             :         struct canditer ci;
    3015             :         BUN ncand;
    3016             :         const void *restrict src;
    3017             :         /* these two needed for ADD_WITH_CHECK macro */
    3018             :         bool abort_on_error = true;
    3019             :         BUN nils = 0;
    3020             : 
    3021             :         lng timeoffset = 0;
    3022        4047 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    3023        4047 :         if (qry_ctx != NULL) {
    3024        4047 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    3025             :         }
    3026             : 
    3027        4047 :         ncand = canditer_init(&ci, b, s);
    3028             : 
    3029        4047 :         BATiter bi = bat_iterator(b);
    3030        4047 :         src = bi.base;
    3031             : 
    3032        4047 :         switch (b->ttype) {
    3033           6 :         case TYPE_bte:
    3034          18 :                 AVERAGE_TYPE(bte);
    3035             :                 break;
    3036           0 :         case TYPE_sht:
    3037           0 :                 AVERAGE_TYPE(sht);
    3038             :                 break;
    3039        3989 :         case TYPE_int:
    3040     4742278 :                 AVERAGE_TYPE(int);
    3041             :                 break;
    3042           9 :         case TYPE_lng:
    3043    10000635 :                 AVERAGE_TYPE(lng);
    3044             :                 break;
    3045             : #ifdef HAVE_HGE
    3046           0 :         case TYPE_hge:
    3047           0 :                 AVERAGE_TYPE(hge);
    3048             :                 break;
    3049             : #endif
    3050           1 :         case TYPE_flt:
    3051   100006105 :                 AVERAGE_FLOATTYPE(flt);
    3052           1 :                 break;
    3053          42 :         case TYPE_dbl:
    3054         155 :                 AVERAGE_FLOATTYPE(dbl);
    3055          42 :                 break;
    3056           0 :         default:
    3057           0 :                 bat_iterator_end(&bi);
    3058           0 :                 GDKerror("average of type %s unsupported.\n",
    3059             :                          ATOMname(b->ttype));
    3060           0 :                 return GDK_FAIL;
    3061             :         }
    3062        4047 :         bat_iterator_end(&bi);
    3063        4047 :         if (scale != 0 && !is_dbl_nil(*avg))
    3064           0 :                 *avg *= pow(10.0, (double) scale);
    3065        4047 :         if (vals)
    3066        4047 :                 *vals = (BUN) n;
    3067             :         return GDK_SUCCEED;
    3068             : }
    3069             : 
    3070             : /* ---------------------------------------------------------------------- */
    3071             : /* count */
    3072             : 
    3073             : #define AGGR_COUNT(TYPE)                                                \
    3074             :         do {                                                            \
    3075             :                 const TYPE *restrict vals = (const TYPE *) bi.base;     \
    3076             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    3077             :                         i = canditer_next(&ci) - b->hseqbase;            \
    3078             :                         if (gids == NULL ||                             \
    3079             :                             (gids[i] >= min && gids[i] <= max)) { \
    3080             :                                 if (gids)                               \
    3081             :                                         gid = gids[i] - min;            \
    3082             :                                 else                                    \
    3083             :                                         gid = (oid) i;                  \
    3084             :                                 if (!is_##TYPE##_nil(vals[i])) {        \
    3085             :                                         cnts[gid]++;                    \
    3086             :                                 }                                       \
    3087             :                         }                                               \
    3088             :                 }                                                       \
    3089             :                 TIMEOUT_CHECK(timeoffset,                               \
    3090             :                               GOTO_LABEL_TIMEOUT_HANDLER(bailout));     \
    3091             :         } while (0)
    3092             : 
    3093             : /* calculate group counts with optional candidates list */
    3094             : BAT *
    3095       13404 : BATgroupcount(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    3096             : {
    3097             :         const oid *restrict gids;
    3098             :         oid gid;
    3099             :         oid min, max;
    3100             :         BUN i, ngrp;
    3101             :         lng *restrict cnts;
    3102             :         BAT *bn = NULL;
    3103             :         int t;
    3104             :         const void *nil;
    3105             :         int (*atomcmp)(const void *, const void *);
    3106             :         struct canditer ci;
    3107             :         BUN ncand;
    3108             :         const char *err;
    3109             :         lng t0 = 0;
    3110             : 
    3111             :         lng timeoffset = 0;
    3112       13404 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    3113       13404 :         if (qry_ctx != NULL) {
    3114       13404 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    3115             :         }
    3116             : 
    3117       13404 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3118             : 
    3119       13404 :         assert(tp == TYPE_lng);
    3120             :         (void) tp;              /* compatibility (with other BATgroup* */
    3121             :         (void) abort_on_error;  /* functions) argument */
    3122             : 
    3123       13404 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    3124           0 :                 GDKerror("%s\n", err);
    3125           0 :                 return NULL;
    3126             :         }
    3127       13404 :         if (g == NULL) {
    3128           0 :                 GDKerror("b and g must be aligned\n");
    3129           0 :                 return NULL;
    3130             :         }
    3131             : 
    3132       13404 :         if (ncand == 0 || ngrp == 0) {
    3133             :                 /* trivial: no counts, so return bat aligned with g
    3134             :                  * with zero in the tail */
    3135        5840 :                 lng zero = 0;
    3136        5840 :                 return BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT);
    3137             :         }
    3138             : 
    3139        7564 :         bn = COLnew(min, TYPE_lng, ngrp, TRANSIENT);
    3140        7564 :         if (bn == NULL)
    3141             :                 return NULL;
    3142        7564 :         cnts = (lng *) Tloc(bn, 0);
    3143        7564 :         memset(cnts, 0, ngrp * sizeof(lng));
    3144             : 
    3145        7564 :         if (BATtdense(g))
    3146             :                 gids = NULL;
    3147             :         else
    3148        5625 :                 gids = (const oid *) Tloc(g, 0);
    3149             : 
    3150        7564 :         if (!skip_nils || b->tnonil) {
    3151             :                 /* if nils are nothing special, or if there are no
    3152             :                  * nils, we don't need to look at the values at all */
    3153        7464 :                 if (gids) {
    3154    16854940 :                         while (ncand > 0) {
    3155    16849361 :                                 ncand--;
    3156    16849361 :                                 i = canditer_next(&ci) - b->hseqbase;
    3157    16849361 :                                 if (gids[i] >= min && gids[i] <= max)
    3158    16857734 :                                         cnts[gids[i] - min]++;
    3159             :                         }
    3160             :                 } else {
    3161      101527 :                         while (ncand > 0) {
    3162       99642 :                                 ncand--;
    3163       99642 :                                 i = canditer_next(&ci) - b->hseqbase;
    3164       99642 :                                 cnts[i] = 1;
    3165             :                         }
    3166             :                 }
    3167             :         } else {
    3168         100 :                 t = b->ttype;
    3169         100 :                 nil = ATOMnilptr(t);
    3170         100 :                 atomcmp = ATOMcompare(t);
    3171         100 :                 t = ATOMbasetype(t);
    3172             : 
    3173         100 :                 BATiter bi = bat_iterator(b);
    3174             : 
    3175         100 :                 switch (t) {
    3176           6 :                 case TYPE_bte:
    3177          20 :                         AGGR_COUNT(bte);
    3178             :                         break;
    3179           0 :                 case TYPE_sht:
    3180           0 :                         AGGR_COUNT(sht);
    3181             :                         break;
    3182          74 :                 case TYPE_int:
    3183       15771 :                         AGGR_COUNT(int);
    3184             :                         break;
    3185           4 :                 case TYPE_lng:
    3186          12 :                         AGGR_COUNT(lng);
    3187             :                         break;
    3188             : #ifdef HAVE_HGE
    3189           0 :                 case TYPE_hge:
    3190           0 :                         AGGR_COUNT(hge);
    3191             :                         break;
    3192             : #endif
    3193           0 :                 case TYPE_flt:
    3194           0 :                         AGGR_COUNT(flt);
    3195             :                         break;
    3196           5 :                 case TYPE_dbl:
    3197         155 :                         AGGR_COUNT(dbl);
    3198             :                         break;
    3199             :                 default:
    3200         175 :                         while (ncand > 0) {
    3201         164 :                                 ncand--;
    3202         164 :                                 i = canditer_next(&ci) - b->hseqbase;
    3203         164 :                                 if (gids == NULL ||
    3204         162 :                                     (gids[i] >= min && gids[i] <= max)) {
    3205         164 :                                         if (gids)
    3206         162 :                                                 gid = gids[i] - min;
    3207             :                                         else
    3208             :                                                 gid = (oid) i;
    3209         164 :                                         if ((*atomcmp)(BUNtail(bi, i), nil) != 0) {
    3210          88 :                                                 cnts[gid]++;
    3211             :                                         }
    3212             :                                 }
    3213             :                         }
    3214             :                         break;
    3215             :                 }
    3216         100 :                 bat_iterator_end(&bi);
    3217             :         }
    3218        7564 :         BATsetcount(bn, ngrp);
    3219        7564 :         bn->tkey = BATcount(bn) <= 1;
    3220        7564 :         bn->tsorted = BATcount(bn) <= 1;
    3221        7564 :         bn->trevsorted = BATcount(bn) <= 1;
    3222        7564 :         bn->tnil = false;
    3223        7564 :         bn->tnonil = true;
    3224        7564 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    3225             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    3226             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    3227             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    3228             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    3229             :                   ci.seq, ncand, GDKusec() - t0);
    3230             :         return bn;
    3231             : 
    3232             :   bailout:
    3233           0 :         BBPreclaim(bn);
    3234           0 :         return NULL;
    3235             : }
    3236             : 
    3237             : /* calculate group sizes (number of TRUE values) with optional
    3238             :  * candidates list */
    3239             : BAT *
    3240           0 : BATgroupsize(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    3241             : {
    3242             :         const oid *restrict gids;
    3243             :         oid min, max;
    3244             :         BUN i, ngrp;
    3245             :         const bit *restrict bits;
    3246             :         lng *restrict cnts;
    3247             :         BAT *bn = NULL;
    3248             :         struct canditer ci;
    3249             :         BUN ncand;
    3250             :         const char *err;
    3251             :         lng t0 = 0;
    3252             : 
    3253           0 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3254             : 
    3255           0 :         assert(tp == TYPE_lng);
    3256           0 :         assert(b->ttype == TYPE_bit);
    3257             :         /* compatibility arguments */
    3258             :         (void) tp;
    3259             :         (void) abort_on_error;
    3260             :         (void) skip_nils;
    3261             : 
    3262           0 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    3263           0 :                 GDKerror("%s\n", err);
    3264           0 :                 return NULL;
    3265             :         }
    3266           0 :         if (g == NULL) {
    3267           0 :                 GDKerror("b and g must be aligned\n");
    3268           0 :                 return NULL;
    3269             :         }
    3270             : 
    3271           0 :         if (ncand == 0 || ngrp == 0) {
    3272             :                 /* trivial: no products, so return bat aligned with g
    3273             :                  * with zero in the tail */
    3274           0 :                 lng zero = 0;
    3275           0 :                 return BATconstant(ngrp == 0 ? 0 : min, TYPE_lng, &zero, ngrp, TRANSIENT);
    3276             :         }
    3277             : 
    3278           0 :         bn = COLnew(min, TYPE_lng, ngrp, TRANSIENT);
    3279           0 :         if (bn == NULL)
    3280             :                 return NULL;
    3281           0 :         cnts = (lng *) Tloc(bn, 0);
    3282           0 :         memset(cnts, 0, ngrp * sizeof(lng));
    3283             : 
    3284           0 :         if (BATtdense(g))
    3285             :                 gids = NULL;
    3286             :         else
    3287           0 :                 gids = (const oid *) Tloc(g, 0);
    3288             : 
    3289           0 :         BATiter bi = bat_iterator(b);
    3290             : 
    3291           0 :         bits = (const bit *) bi.base;
    3292             : 
    3293           0 :         while (ncand > 0) {
    3294           0 :                 ncand--;
    3295           0 :                 i = canditer_next(&ci) - b->hseqbase;
    3296           0 :                 if (bits[i] == 1 &&
    3297           0 :                     (gids == NULL || (gids[i] >= min && gids[i] <= max))) {
    3298           0 :                         cnts[gids ? gids[i] - min : (oid) i]++;
    3299             :                 }
    3300             :         }
    3301           0 :         bat_iterator_end(&bi);
    3302           0 :         BATsetcount(bn, ngrp);
    3303           0 :         bn->tkey = BATcount(bn) <= 1;
    3304           0 :         bn->tsorted = BATcount(bn) <= 1;
    3305           0 :         bn->trevsorted = BATcount(bn) <= 1;
    3306           0 :         bn->tnil = false;
    3307           0 :         bn->tnonil = true;
    3308           0 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    3309             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    3310             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    3311             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    3312             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    3313             :                   ci.seq, ncand, GDKusec() - t0);
    3314             :         return bn;
    3315             : }
    3316             : 
    3317             : /* ---------------------------------------------------------------------- */
    3318             : /* min and max */
    3319             : 
    3320             : #define AGGR_CMP(TYPE, OP)                                              \
    3321             :         do {                                                            \
    3322             :                 const TYPE *restrict vals = (const TYPE *) bi->base; \
    3323             :                 if (ngrp == ncand) {                                    \
    3324             :                         /* single element groups */                     \
    3325             :                         TIMEOUT_LOOP(ncand, timeoffset) {               \
    3326             :                                 i = canditer_next(ci) - hseq;           \
    3327             :                                 if (!skip_nils ||                       \
    3328             :                                     !is_##TYPE##_nil(vals[i])) {        \
    3329             :                                         oids[i] = i + hseq;             \
    3330             :                                         nils--;                         \
    3331             :                                 }                                       \
    3332             :                         }                                               \
    3333             :                         TIMEOUT_CHECK(timeoffset,                       \
    3334             :                                       TIMEOUT_HANDLER(BUN_NONE));       \
    3335             :                 } else {                                                \
    3336             :                         gid = 0; /* in case gids == NULL */             \
    3337             :                         TIMEOUT_LOOP(ncand, timeoffset) {               \
    3338             :                                 i = canditer_next(ci) - hseq;           \
    3339             :                                 if (gids == NULL ||                     \
    3340             :                                     (gids[i] >= min && gids[i] <= max)) { \
    3341             :                                         if (gids)                       \
    3342             :                                                 gid = gids[i] - min;    \
    3343             :                                         if (!skip_nils || !is_##TYPE##_nil(vals[i])) { \
    3344             :                                                 if (is_oid_nil(oids[gid])) { \
    3345             :                                                         oids[gid] = i + hseq; \
    3346             :                                                         nils--;         \
    3347             :                                                 } else if (!is_##TYPE##_nil(vals[oids[gid] - hseq]) && \
    3348             :                                                            (is_##TYPE##_nil(vals[i]) || \
    3349             :                                                             OP(vals[i], vals[oids[gid] - hseq]))) \
    3350             :                                                         oids[gid] = i + hseq; \
    3351             :                                         }                               \
    3352             :                                 }                                       \
    3353             :                         }                                               \
    3354             :                         TIMEOUT_CHECK(timeoffset,                       \
    3355             :                                       TIMEOUT_HANDLER(BUN_NONE));       \
    3356             :                 }                                                       \
    3357             :         } while (0)
    3358             : 
    3359             : /* calculate group minimums with optional candidates list
    3360             :  *
    3361             :  * note that this functions returns *positions* of where the minimum
    3362             :  * values occur */
    3363             : static BUN
    3364         683 : do_groupmin(oid *restrict oids, BATiter *bi, const oid *restrict gids, BUN ngrp,
    3365             :             oid min, oid max, struct canditer *restrict ci, BUN ncand,
    3366             :             bool skip_nils, bool gdense)
    3367             : {
    3368             :         oid gid;
    3369             :         BUN i, nils;
    3370             :         int t;
    3371             :         const void *nil;
    3372             :         int (*atomcmp)(const void *, const void *);
    3373             : 
    3374             :         lng timeoffset = 0;
    3375         683 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    3376         683 :         if (qry_ctx != NULL) {
    3377         683 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    3378             :         }
    3379             : 
    3380             :         nils = ngrp;
    3381       31319 :         for (i = 0; i < ngrp; i++)
    3382       30636 :                 oids[i] = oid_nil;
    3383         683 :         if (ncand == 0)
    3384             :                 return nils;
    3385             : 
    3386         683 :         t = bi->b->ttype;
    3387         683 :         nil = ATOMnilptr(t);
    3388         683 :         atomcmp = ATOMcompare(t);
    3389         683 :         t = ATOMbasetype(t);
    3390         683 :         oid hseq = bi->b->hseqbase;
    3391             : 
    3392         683 :         switch (t) {
    3393          34 :         case TYPE_bte:
    3394        4281 :                 AGGR_CMP(bte, LT);
    3395             :                 break;
    3396           1 :         case TYPE_sht:
    3397        1322 :                 AGGR_CMP(sht, LT);
    3398             :                 break;
    3399         470 :         case TYPE_int:
    3400      219764 :                 AGGR_CMP(int, LT);
    3401             :                 break;
    3402          59 :         case TYPE_lng:
    3403       50152 :                 AGGR_CMP(lng, LT);
    3404             :                 break;
    3405             : #ifdef HAVE_HGE
    3406           4 :         case TYPE_hge:
    3407          18 :                 AGGR_CMP(hge, LT);
    3408             :                 break;
    3409             : #endif
    3410           0 :         case TYPE_flt:
    3411           0 :                 AGGR_CMP(flt, LT);
    3412             :                 break;
    3413          37 :         case TYPE_dbl:
    3414         363 :                 AGGR_CMP(dbl, LT);
    3415             :                 break;
    3416           3 :         case TYPE_void:
    3417           3 :                 if (!skip_nils || !is_oid_nil(bi->tseq)) {
    3418           0 :                         if (!gdense && gids == NULL) {
    3419           0 :                                 oids[0] = canditer_next(ci);
    3420           0 :                                 nils--;
    3421           0 :                         } else if (gdense) {
    3422             :                                 /* single element groups */
    3423           0 :                                 while (ncand > 0) {
    3424           0 :                                         ncand--;
    3425           0 :                                         i = canditer_next(ci);
    3426           0 :                                         oids[i - hseq] = i;
    3427           0 :                                         nils--;
    3428             :                                 }
    3429             :                         } else {
    3430           0 :                                 while (ncand > 0) {
    3431           0 :                                         ncand--;
    3432           0 :                                         i = canditer_next(ci);
    3433           0 :                                         if (is_oid_nil(oids[i - hseq])) {
    3434           0 :                                                 oids[i - hseq] = i;
    3435           0 :                                                 nils--;
    3436             :                                         }
    3437             :                                 }
    3438             :                         }
    3439             :                 }
    3440             :                 break;
    3441          75 :         default:
    3442          75 :                 assert(bi->b->ttype != TYPE_oid);
    3443             : 
    3444          75 :                 if (gdense) {
    3445             :                         /* single element groups */
    3446           3 :                         TIMEOUT_LOOP(ncand, timeoffset) {
    3447           1 :                                 i = canditer_next(ci) - hseq;
    3448           2 :                                 if (!skip_nils ||
    3449           1 :                                     (*atomcmp)(BUNtail(*bi, i), nil) != 0) {
    3450           0 :                                         oids[i] = i + hseq;
    3451           0 :                                         nils--;
    3452             :                                 }
    3453             :                         }
    3454           1 :                         TIMEOUT_CHECK(timeoffset,
    3455             :                                       TIMEOUT_HANDLER(BUN_NONE));
    3456             :                 } else {
    3457             :                         gid = 0; /* in case gids == NULL */
    3458      393176 :                         TIMEOUT_LOOP(ncand, timeoffset) {
    3459      393004 :                                 i = canditer_next(ci) - hseq;
    3460      395051 :                                 if (gids == NULL ||
    3461          51 :                                     (gids[i] >= min && gids[i] <= max)) {
    3462      395051 :                                         const void *v = BUNtail(*bi, i);
    3463      395909 :                                         if (gids)
    3464         669 :                                                 gid = gids[i] - min;
    3465      831894 :                                         if (!skip_nils ||
    3466      395315 :                                             (*atomcmp)(v, nil) != 0) {
    3467      433558 :                                                 if (is_oid_nil(oids[gid])) {
    3468          85 :                                                         oids[gid] = i + hseq;
    3469          85 :                                                         nils--;
    3470      433473 :                                                 } else if (t != TYPE_void) {
    3471      437505 :                                                         const void *g = BUNtail(*bi, (BUN) (oids[gid] - hseq));
    3472      791027 :                                                         if ((*atomcmp)(g, nil) != 0 &&
    3473      777233 :                                                             ((*atomcmp)(v, nil) == 0 ||
    3474      353522 :                                                              LT((*atomcmp)(v, g), 0)))
    3475          75 :                                                                 oids[gid] = i + hseq;
    3476             :                                                 }
    3477             :                                         }
    3478             :                                 }
    3479             :                         }
    3480          74 :                         TIMEOUT_CHECK(timeoffset,
    3481             :                                       TIMEOUT_HANDLER(BUN_NONE));
    3482             :                 }
    3483             :                 break;
    3484             :         }
    3485             : 
    3486             :         return nils;
    3487             : }
    3488             : 
    3489             : /* calculate group maximums with optional candidates list
    3490             :  *
    3491             :  * note that this functions returns *positions* of where the maximum
    3492             :  * values occur */
    3493             : static BUN
    3494         639 : do_groupmax(oid *restrict oids, BATiter *bi, const oid *restrict gids, BUN ngrp,
    3495             :             oid min, oid max, struct canditer *restrict ci, BUN ncand,
    3496             :             bool skip_nils, bool gdense)
    3497             : {
    3498             :         oid gid;
    3499             :         BUN i, nils;
    3500             :         int t;
    3501             :         const void *nil;
    3502             :         int (*atomcmp)(const void *, const void *);
    3503             : 
    3504             :         lng timeoffset = 0;
    3505         639 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    3506         639 :         if (qry_ctx != NULL) {
    3507         639 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    3508             :         }
    3509             : 
    3510             :         nils = ngrp;
    3511      677601 :         for (i = 0; i < ngrp; i++)
    3512      676962 :                 oids[i] = oid_nil;
    3513         639 :         if (ncand == 0)
    3514             :                 return nils;
    3515             : 
    3516         639 :         t = bi->b->ttype;
    3517         639 :         nil = ATOMnilptr(t);
    3518         639 :         atomcmp = ATOMcompare(t);
    3519         639 :         t = ATOMbasetype(t);
    3520         639 :         oid hseq = bi->b->hseqbase;
    3521             : 
    3522         639 :         switch (t) {
    3523          45 :         case TYPE_bte:
    3524        4347 :                 AGGR_CMP(bte, GT);
    3525             :                 break;
    3526           1 :         case TYPE_sht:
    3527         130 :                 AGGR_CMP(sht, GT);
    3528             :                 break;
    3529         407 :         case TYPE_int:
    3530      333536 :                 AGGR_CMP(int, GT);
    3531             :                 break;
    3532          74 :         case TYPE_lng:
    3533       57790 :                 AGGR_CMP(lng, GT);
    3534             :                 break;
    3535             : #ifdef HAVE_HGE
    3536          10 :         case TYPE_hge:
    3537    34144874 :                 AGGR_CMP(hge, GT);
    3538             :                 break;
    3539             : #endif
    3540           0 :         case TYPE_flt:
    3541           0 :                 AGGR_CMP(flt, GT);
    3542             :                 break;
    3543          29 :         case TYPE_dbl:
    3544         283 :                 AGGR_CMP(dbl, GT);
    3545             :                 break;
    3546           0 :         case TYPE_void:
    3547           0 :                 if (!skip_nils || !is_oid_nil(bi->tseq)) {
    3548           0 :                         if (!gdense && gids == NULL) {
    3549           0 :                                 oids[0] = canditer_last(ci);
    3550           0 :                                 nils--;
    3551           0 :                         } else if (gdense) {
    3552             :                                 /* single element groups */
    3553           0 :                                 while (ncand > 0) {
    3554           0 :                                         ncand--;
    3555           0 :                                         i = canditer_next(ci);
    3556           0 :                                         oids[i - hseq] = i;
    3557           0 :                                         nils--;
    3558             :                                 }
    3559             :                         } else {
    3560           0 :                                 while (ncand > 0) {
    3561           0 :                                         ncand--;
    3562           0 :                                         i = canditer_next(ci);
    3563           0 :                                         if (is_oid_nil(oids[i - hseq]))
    3564           0 :                                                 nils--;
    3565           0 :                                         oids[i - hseq] = i;
    3566             :                                 }
    3567             :                         }
    3568             :                 }
    3569             :                 break;
    3570          73 :         default:
    3571          73 :                 assert(bi->b->ttype != TYPE_oid);
    3572             : 
    3573          73 :                 if (gdense) {
    3574             :                         /* single element groups */
    3575           9 :                         TIMEOUT_LOOP(ncand, timeoffset) {
    3576           3 :                                 i = canditer_next(ci) - hseq;
    3577           6 :                                 if (!skip_nils ||
    3578           3 :                                     (*atomcmp)(BUNtail(*bi, i), nil) != 0) {
    3579           3 :                                         oids[i] = i + hseq;
    3580           3 :                                         nils--;
    3581             :                                 }
    3582             :                         }
    3583           3 :                         TIMEOUT_CHECK(timeoffset,
    3584             :                                       TIMEOUT_HANDLER(BUN_NONE));
    3585             :                 } else {
    3586             :                         gid = 0; /* in case gids == NULL */
    3587      403749 :                         TIMEOUT_LOOP(ncand, timeoffset) {
    3588      403585 :                                 i = canditer_next(ci) - hseq;
    3589      403871 :                                 if (gids == NULL ||
    3590          26 :                                     (gids[i] >= min && gids[i] <= max)) {
    3591      403871 :                                         const void *v = BUNtail(*bi, i);
    3592      439652 :                                         if (gids)
    3593       35508 :                                                 gid = gids[i] - min;
    3594      891033 :                                         if (!skip_nils ||
    3595      435092 :                                             (*atomcmp)(v, nil) != 0) {
    3596      453193 :                                                 if (is_oid_nil(oids[gid])) {
    3597          69 :                                                         oids[gid] = i + hseq;
    3598          69 :                                                         nils--;
    3599             :                                                 } else {
    3600      453124 :                                                         const void *g = BUNtail(*bi, (BUN) (oids[gid] - hseq));
    3601      843779 :                                                         if (t == TYPE_void ||
    3602      807454 :                                                             ((*atomcmp)(g, nil) != 0 &&
    3603      793437 :                                                              ((*atomcmp)(v, nil) == 0 ||
    3604      354592 :                                                               GT((*atomcmp)(v, g), 0))))
    3605           0 :                                                                 oids[gid] = i + hseq;
    3606             :                                                 }
    3607             :                                         }
    3608             :                                 }
    3609             :                         }
    3610          70 :                         TIMEOUT_CHECK(timeoffset,
    3611             :                                       TIMEOUT_HANDLER(BUN_NONE));
    3612             :                 }
    3613             :                 break;
    3614             :         }
    3615             : 
    3616             :         return nils;
    3617             : }
    3618             : 
    3619             : static BAT *
    3620         851 : BATgroupminmax(BAT *b, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils,
    3621             :                bool abort_on_error,
    3622             :                BUN (*minmax)(oid *restrict, BATiter *, const oid *restrict, BUN,
    3623             :                              oid, oid, struct canditer *restrict, BUN,
    3624             :                              bool, bool),
    3625             :                const char *name)
    3626             : {
    3627             :         const oid *restrict gids;
    3628             :         oid min, max;
    3629             :         BUN ngrp;
    3630             :         oid *restrict oids;
    3631             :         BAT *bn = NULL;
    3632             :         BUN nils;
    3633             :         struct canditer ci;
    3634             :         BUN ncand;
    3635             :         const char *err;
    3636             :         lng t0 = 0;
    3637             : 
    3638         851 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3639             : 
    3640         851 :         assert(tp == TYPE_oid);
    3641             :         (void) tp;              /* compatibility (with other BATgroup* */
    3642             :         (void) abort_on_error;  /* functions) argument */
    3643             : 
    3644         851 :         if (!ATOMlinear(b->ttype)) {
    3645           0 :                 GDKerror("%s: cannot determine minimum on "
    3646             :                          "non-linear type %s\n", name, ATOMname(b->ttype));
    3647           0 :                 return NULL;
    3648             :         }
    3649             : 
    3650         851 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    3651           0 :                 GDKerror("%s: %s\n", name, err);
    3652           0 :                 return NULL;
    3653             :         }
    3654             : 
    3655         851 :         if (ncand == 0 || ngrp == 0) {
    3656             :                 /* trivial: no minimums, so return bat aligned with g
    3657             :                  * with nil in the tail */
    3658         136 :                 return BATconstant(ngrp == 0 ? 0 : min, TYPE_oid, &oid_nil, ngrp, TRANSIENT);
    3659             :         }
    3660             : 
    3661         715 :         bn = COLnew(min, TYPE_oid, ngrp, TRANSIENT);
    3662         714 :         if (bn == NULL)
    3663             :                 return NULL;
    3664         714 :         oids = (oid *) Tloc(bn, 0);
    3665             : 
    3666         714 :         if (g == NULL || BATtdense(g))
    3667             :                 gids = NULL;
    3668             :         else
    3669         375 :                 gids = (const oid *) Tloc(g, 0);
    3670             : 
    3671         714 :         BATiter bi = bat_iterator(b);
    3672         715 :         nils = (*minmax)(oids, &bi, gids, ngrp, min, max, &ci, ncand,
    3673         715 :                          skip_nils, g && BATtdense(g));
    3674         715 :         bat_iterator_end(&bi);
    3675         715 :         if (nils == BUN_NONE) {
    3676           0 :                 BBPreclaim(bn);
    3677           0 :                 return NULL;
    3678             :         }
    3679             : 
    3680         715 :         BATsetcount(bn, ngrp);
    3681             : 
    3682         715 :         bn->tkey = BATcount(bn) <= 1;
    3683         715 :         bn->tsorted = BATcount(bn) <= 1;
    3684         715 :         bn->trevsorted = BATcount(bn) <= 1;
    3685         715 :         bn->tnil = nils != 0;
    3686         715 :         bn->tnonil = nils == 0;
    3687         715 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    3688             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT " -> " ALGOOPTBATFMT
    3689             :                   "; start " OIDFMT ", count " BUNFMT " (%s -- " LLFMT " usec)\n",
    3690             :                   ALGOBATPAR(b), ALGOOPTBATPAR(g), ALGOOPTBATPAR(e),
    3691             :                   ALGOOPTBATPAR(s), ALGOOPTBATPAR(bn),
    3692             :                   ci.seq, ncand, name, GDKusec() - t0);
    3693             :         return bn;
    3694             : }
    3695             : 
    3696             : BAT *
    3697         392 : BATgroupmin(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    3698             :             bool skip_nils, bool abort_on_error)
    3699             : {
    3700         392 :         return BATgroupminmax(b, g, e, s, tp, skip_nils, abort_on_error,
    3701             :                               do_groupmin, __func__);
    3702             : }
    3703             : 
    3704             : /* return pointer to smallest non-nil value in b, or pointer to nil if
    3705             :  * there is no such value (no values at all, or only nil) */
    3706             : void *
    3707         960 : BATmin_skipnil(BAT *b, void *aggr, bit skipnil)
    3708             : {
    3709             :         const void *res = NULL;
    3710             :         size_t s;
    3711             :         lng t0 = 0;
    3712             : 
    3713         960 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3714             : 
    3715         960 :         if (!ATOMlinear(b->ttype)) {
    3716             :                 /* there is no such thing as a smallest value if you
    3717             :                  * can't compare values */
    3718           0 :                 GDKerror("non-linear type");
    3719           0 :                 return NULL;
    3720             :         }
    3721         960 :         BATiter bi = bat_iterator(b);
    3722         960 :         if (bi.count == 0) {
    3723         223 :                 res = ATOMnilptr(b->ttype);
    3724         737 :         } else if (bi.minpos != BUN_NONE) {
    3725         409 :                 res = BUNtail(bi, bi.minpos);
    3726             :         } else {
    3727             :                 oid pos;
    3728             :                 BAT *pb = NULL;
    3729             :                 Heap *oidxh = NULL;
    3730             : 
    3731         328 :                 if (BATcheckorderidx(b)) {
    3732           0 :                         MT_lock_set(&b->batIdxLock);
    3733           0 :                         oidxh = b->torderidx;
    3734           0 :                         if (oidxh != NULL)
    3735           0 :                                 HEAPincref(oidxh);
    3736           0 :                         MT_lock_unset(&b->batIdxLock);
    3737             :                 }
    3738         204 :                 if (oidxh == NULL &&
    3739         328 :                     VIEWtparent(b) &&
    3740         408 :                     (pb = BBP_cache(VIEWtparent(b))) != NULL &&
    3741         204 :                     BATcheckorderidx(pb)) {
    3742             :                         /* no lock on b needed since it's a view */
    3743           0 :                         MT_lock_set(&pb->batIdxLock);
    3744           0 :                         MT_lock_set(&pb->theaplock);
    3745           0 :                         if (pb->tbaseoff == b->tbaseoff &&
    3746           0 :                             BATcount(pb) == bi.count &&
    3747           0 :                             pb->hseqbase == b->hseqbase &&
    3748           0 :                             (oidxh = pb->torderidx) != NULL) {
    3749           0 :                                 HEAPincref(oidxh);
    3750             :                         }
    3751           0 :                         MT_lock_unset(&pb->batIdxLock);
    3752           0 :                         MT_lock_unset(&pb->theaplock);
    3753             :                 }
    3754         328 :                 if (oidxh != NULL) {
    3755           0 :                         const oid *ords = (const oid *) oidxh->base + ORDERIDXOFF;
    3756             :                         BUN r;
    3757           0 :                         if (!b->tnonil) {
    3758           0 :                                 MT_thread_setalgorithm(pb ? "binsearch on parent oidx" : "binsearch on oidx");
    3759           0 :                                 r = binsearch(ords, 0, b->ttype, bi.base,
    3760           0 :                                               bi.vh ? bi.vh->base : NULL,
    3761           0 :                                               bi.width, 0, bi.count,
    3762           0 :                                               ATOMnilptr(b->ttype), 1, 1);
    3763           0 :                                 if (r == 0) {
    3764           0 :                                         b->tnonil = true;
    3765           0 :                                         b->batDirtydesc = true;
    3766             :                                 }
    3767             :                         } else {
    3768             :                                 r = 0;
    3769             :                         }
    3770           0 :                         if (r == bi.count) {
    3771             :                                 /* no non-nil values */
    3772           0 :                                 pos = oid_nil;
    3773             :                         } else {
    3774           0 :                                 MT_thread_setalgorithm(pb ? "using parent oidx" : "using oidx");
    3775           0 :                                 pos = ords[r];
    3776             :                         }
    3777           0 :                         HEAPdecref(oidxh, false);
    3778             :                 } else {
    3779             :                         Imprints *imprints = NULL;
    3780         328 :                         if ((VIEWtparent(b) == 0 ||
    3781         413 :                              bi.count == BATcount(BBP_cache(VIEWtparent(b)))) &&
    3782         209 :                             BATcheckimprints(b)) {
    3783           0 :                                 if (VIEWtparent(b)) {
    3784           0 :                                         BAT *pb = BBP_cache(VIEWtparent(b));
    3785           0 :                                         MT_lock_set(&pb->batIdxLock);
    3786           0 :                                         imprints = pb->timprints;
    3787           0 :                                         if (imprints != NULL)
    3788           0 :                                                 IMPSincref(imprints);
    3789             :                                         else
    3790             :                                                 imprints = NULL;
    3791           0 :                                         MT_lock_unset(&pb->batIdxLock);
    3792             :                                 } else {
    3793           0 :                                         MT_lock_set(&b->batIdxLock);
    3794           0 :                                         imprints = b->timprints;
    3795           0 :                                         if (imprints != NULL)
    3796           0 :                                                 IMPSincref(imprints);
    3797             :                                         else
    3798             :                                                 imprints = NULL;
    3799           0 :                                         MT_lock_unset(&b->batIdxLock);
    3800             :                                 }
    3801             :                         }
    3802           0 :                         if (imprints) {
    3803             :                                 int i;
    3804             : 
    3805           0 :                                 MT_thread_setalgorithm(VIEWtparent(b) ? "using parent imprints" : "using imprints");
    3806           0 :                                 pos = oid_nil;
    3807             :                                 /* find first non-empty bin */
    3808           0 :                                 for (i = 0; i < imprints->bits; i++) {
    3809           0 :                                         if (imprints->stats[i + 128]) {
    3810           0 :                                                 pos = imprints->stats[i] + b->hseqbase;
    3811           0 :                                                 break;
    3812             :                                         }
    3813             :                                 }
    3814           0 :                                 IMPSdecref(imprints, false);
    3815             :                         } else {
    3816             :                                 struct canditer ci;
    3817         328 :                                 BUN ncand = canditer_init(&ci, b, NULL);
    3818         328 :                                 (void) do_groupmin(&pos, &bi, NULL, 1, 0, 0, &ci, ncand,
    3819             :                                                    skipnil, false);
    3820             :                         }
    3821             :                 }
    3822         328 :                 if (is_oid_nil(pos)) {
    3823          39 :                         res = ATOMnilptr(b->ttype);
    3824             :                 } else {
    3825         289 :                         res = BUNtail(bi, pos - b->hseqbase);
    3826             :                 }
    3827             :         }
    3828         960 :         if (aggr == NULL) {
    3829         534 :                 s = ATOMlen(b->ttype, res);
    3830         534 :                 aggr = GDKmalloc(s);
    3831             :         } else {
    3832         852 :                 s = ATOMsize(ATOMtype(b->ttype));
    3833             :         }
    3834         960 :         if (aggr != NULL)       /* else: malloc error */
    3835         960 :                 memcpy(aggr, res, s);
    3836         960 :         bat_iterator_end(&bi);
    3837         960 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
    3838             :                   ALGOBATPAR(b), skipnil, GDKusec() - t0);
    3839             :         return aggr;
    3840             : }
    3841             : 
    3842             : void *
    3843         497 : BATmin(BAT *b, void *aggr)
    3844             : {
    3845         497 :         return BATmin_skipnil(b, aggr, 1);
    3846             : }
    3847             : 
    3848             : BAT *
    3849         459 : BATgroupmax(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    3850             :             bool skip_nils, bool abort_on_error)
    3851             : {
    3852         459 :         return BATgroupminmax(b, g, e, s, tp, skip_nils, abort_on_error,
    3853             :                               do_groupmax, __func__);
    3854             : }
    3855             : 
    3856             : void *
    3857        1114 : BATmax_skipnil(BAT *b, void *aggr, bit skipnil)
    3858             : {
    3859             :         const void *res = NULL;
    3860             :         size_t s;
    3861             :         BATiter bi;
    3862             :         lng t0 = 0;
    3863             : 
    3864        1114 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    3865             : 
    3866        1114 :         if (!ATOMlinear(b->ttype)) {
    3867           0 :                 GDKerror("non-linear type");
    3868           0 :                 return NULL;
    3869             :         }
    3870        1114 :         bi = bat_iterator(b);
    3871        1114 :         if (bi.count == 0) {
    3872         193 :                 res = ATOMnilptr(b->ttype);
    3873         921 :         } else if (bi.maxpos != BUN_NONE) {
    3874         642 :                 res = BUNtail(bi, bi.maxpos);
    3875             :         } else {
    3876             :                 oid pos;
    3877             :                 BAT *pb = NULL;
    3878             :                 Heap *oidxh = NULL;
    3879             : 
    3880         279 :                 if (BATcheckorderidx(b)) {
    3881           0 :                         MT_lock_set(&b->batIdxLock);
    3882           0 :                         oidxh = b->torderidx;
    3883           0 :                         if (oidxh != NULL)
    3884           0 :                                 HEAPincref(oidxh);
    3885           0 :                         MT_lock_unset(&b->batIdxLock);
    3886             :                 }
    3887         237 :                 if (oidxh == NULL &&
    3888         279 :                     VIEWtparent(b) &&
    3889         474 :                     (pb = BBP_cache(VIEWtparent(b))) != NULL &&
    3890         237 :                     BATcheckorderidx(pb)) {
    3891             :                         /* no lock on b needed since it's a view */
    3892           0 :                         MT_lock_set(&pb->batIdxLock);
    3893           0 :                         MT_lock_set(&pb->theaplock);
    3894           0 :                         if (pb->tbaseoff == b->tbaseoff &&
    3895           0 :                             BATcount(pb) == bi.count &&
    3896           0 :                             pb->hseqbase == b->hseqbase &&
    3897           0 :                             (oidxh = pb->torderidx) != NULL) {
    3898           0 :                                 HEAPincref(oidxh);
    3899             :                         }
    3900           0 :                         MT_lock_unset(&pb->batIdxLock);
    3901           0 :                         MT_lock_unset(&pb->theaplock);
    3902             :                 }
    3903         279 :                 if (oidxh != NULL) {
    3904           0 :                         const oid *ords = (const oid *) oidxh->base + ORDERIDXOFF;
    3905             : 
    3906           0 :                         MT_thread_setalgorithm(pb ? "using parent oidx" : "using oids");
    3907           0 :                         pos = ords[bi.count - 1];
    3908             :                         /* nils are first, ie !skipnil, check for nils */
    3909           0 :                         if (!skipnil) {
    3910           0 :                                 BUN z = ords[0];
    3911             : 
    3912           0 :                                 res = BUNtail(bi, z - b->hseqbase);
    3913             : 
    3914           0 :                                 if (ATOMcmp(b->ttype, res, ATOMnilptr(b->ttype)) == 0)
    3915           0 :                                         pos = z;
    3916             :                         }
    3917           0 :                         HEAPdecref(oidxh, false);
    3918             :                 } else {
    3919             :                         Imprints *imprints = NULL;
    3920         279 :                         if ((VIEWtparent(b) == 0 ||
    3921         367 :                              BATcount(b) == BATcount(BBP_cache(VIEWtparent(b)))) &&
    3922         130 :                             BATcheckimprints(b)) {
    3923           0 :                                 if (VIEWtparent(b)) {
    3924           0 :                                         BAT *pb = BBP_cache(VIEWtparent(b));
    3925           0 :                                         MT_lock_set(&pb->batIdxLock);
    3926           0 :                                         imprints = pb->timprints;
    3927           0 :                                         if (imprints != NULL)
    3928           0 :                                                 IMPSincref(imprints);
    3929             :                                         else
    3930             :                                                 imprints = NULL;
    3931           0 :                                         MT_lock_unset(&pb->batIdxLock);
    3932             :                                 } else {
    3933           0 :                                         MT_lock_set(&b->batIdxLock);
    3934           0 :                                         imprints = b->timprints;
    3935           0 :                                         if (imprints != NULL)
    3936           0 :                                                 IMPSincref(imprints);
    3937             :                                         else
    3938             :                                                 imprints = NULL;
    3939           0 :                                         MT_lock_unset(&b->batIdxLock);
    3940             :                                 }
    3941             :                         }
    3942           0 :                         if (imprints) {
    3943             :                                 int i;
    3944             : 
    3945           0 :                                 MT_thread_setalgorithm(VIEWtparent(b) ? "using parent imprints" : "using imprints");
    3946           0 :                                 pos = oid_nil;
    3947             :                                 /* find last non-empty bin */
    3948           0 :                                 for (i = imprints->bits - 1; i >= 0; i--) {
    3949           0 :                                         if (imprints->stats[i + 128]) {
    3950           0 :                                                 pos = imprints->stats[i + 64] + b->hseqbase;
    3951           0 :                                                 break;
    3952             :                                         }
    3953             :                                 }
    3954           0 :                                 IMPSdecref(imprints, false);
    3955             :                         } else {
    3956             :                                 struct canditer ci;
    3957         279 :                                 BUN ncand = canditer_init(&ci, b, NULL);
    3958         279 :                                 (void) do_groupmax(&pos, &bi, NULL, 1, 0, 0, &ci, ncand,
    3959             :                                                    skipnil, false);
    3960             :                         }
    3961             :                 }
    3962         279 :                 if (is_oid_nil(pos)) {
    3963          34 :                         res = ATOMnilptr(b->ttype);
    3964             :                 } else {
    3965         245 :                         res = BUNtail(bi, pos - b->hseqbase);
    3966             :                 }
    3967             :         }
    3968        1114 :         if (aggr == NULL) {
    3969         508 :                 s = ATOMlen(b->ttype, res);
    3970         508 :                 aggr = GDKmalloc(s);
    3971             :         } else {
    3972        1212 :                 s = ATOMsize(ATOMtype(b->ttype));
    3973             :         }
    3974        1114 :         if (aggr != NULL)       /* else: malloc error */
    3975        1114 :                 memcpy(aggr, res, s);
    3976        1114 :         bat_iterator_end(&bi);
    3977        1114 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",skipnil=%d; (" LLFMT " usec)\n",
    3978             :                   ALGOBATPAR(b), skipnil, GDKusec() - t0);
    3979             :         return aggr;
    3980             : }
    3981             : 
    3982             : void *
    3983         811 : BATmax(BAT *b, void *aggr)
    3984             : {
    3985         811 :         return BATmax_skipnil(b, aggr, 1);
    3986             : }
    3987             : 
    3988             : 
    3989             : /* ---------------------------------------------------------------------- */
    3990             : /* quantiles/median */
    3991             : 
    3992             : #if SIZEOF_OID == SIZEOF_INT
    3993             : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_int(indir, offset, (const int *) vals, lo, hi, (int) (v), ordering, last)
    3994             : #endif
    3995             : #if SIZEOF_OID == SIZEOF_LNG
    3996             : #define binsearch_oid(indir, offset, vals, lo, hi, v, ordering, last) binsearch_lng(indir, offset, (const lng *) vals, lo, hi, (lng) (v), ordering, last)
    3997             : #endif
    3998             : 
    3999             : #define DO_QUANTILE_AVG(TPE)                                            \
    4000             :         do {                                                            \
    4001             :                 TPE low = *(TPE*) BUNtail(bi, r + (BUN) hi);            \
    4002             :                 TPE high = *(TPE*) BUNtail(bi, r + (BUN) lo);           \
    4003             :                 if (is_##TPE##_nil(low) || is_##TPE##_nil(high)) {      \
    4004             :                         val = dbl_nil;                                  \
    4005             :                         nils++;                                         \
    4006             :                 } else {                                                \
    4007             :                         val = (f - lo) * low + (lo + 1 - f) * high;     \
    4008             :                 }                                                       \
    4009             :         } while (0)
    4010             : 
    4011             : static BAT *
    4012          63 : doBATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
    4013             :                    bool skip_nils, bool abort_on_error, bool average)
    4014             : {
    4015             :         BAT *origb = b;
    4016             :         BAT *origg = g;
    4017             :         oid min, max;
    4018             :         BUN ngrp;
    4019             :         BUN nils = 0;
    4020             :         BAT *bn = NULL;
    4021             :         struct canditer ci;
    4022             :         BUN ncand;
    4023             :         BAT *t1, *t2;
    4024             :         BATiter bi;
    4025             :         const void *v;
    4026          63 :         const void *nil = ATOMnilptr(tp);
    4027             :         const void *dnil = nil;
    4028             :         dbl val;                /* only used for average */
    4029          63 :         int (*atomcmp)(const void *, const void *) = ATOMcompare(tp);
    4030             :         const char *err;
    4031             :         (void) abort_on_error;
    4032             :         lng t0 = 0;
    4033             : 
    4034             :         size_t counter = 0;
    4035             :         lng timeoffset = 0;
    4036          63 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    4037          63 :         if (qry_ctx != NULL) {
    4038          63 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    4039             :         }
    4040             : 
    4041          63 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    4042             : 
    4043          63 :         if (average) {
    4044          12 :                 switch (ATOMbasetype(b->ttype)) {
    4045             :                 case TYPE_bte:
    4046             :                 case TYPE_sht:
    4047             :                 case TYPE_int:
    4048             :                 case TYPE_lng:
    4049             : #ifdef HAVE_HGE
    4050             :                 case TYPE_hge:
    4051             : #endif
    4052             :                 case TYPE_flt:
    4053             :                 case TYPE_dbl:
    4054             :                         break;
    4055           0 :                 default:
    4056           0 :                         GDKerror("incompatible type\n");
    4057           0 :                         return NULL;
    4058             :                 }
    4059             :                 dnil = &dbl_nil;
    4060             :         }
    4061          63 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    4062           0 :                 GDKerror("%s\n", err);
    4063           0 :                 return NULL;
    4064             :         }
    4065          63 :         assert(tp == b->ttype);
    4066          63 :         if (!ATOMlinear(tp)) {
    4067           0 :                 GDKerror("cannot determine quantile on "
    4068             :                          "non-linear type %s\n", ATOMname(tp));
    4069           0 :                 return NULL;
    4070             :         }
    4071          63 :         if (quantile < 0 || quantile > 1) {
    4072           0 :                 GDKerror("cannot determine quantile for "
    4073             :                          "p=%f (p has to be in [0,1])\n", quantile);
    4074           0 :                 return NULL;
    4075             :         }
    4076             : 
    4077          63 :         if (BATcount(b) == 0 || ngrp == 0 || is_dbl_nil(quantile)) {
    4078             :                 /* trivial: no values, thus also no quantiles,
    4079             :                  * so return bat aligned with e with nil in the tail
    4080             :                  * The same happens for a NULL quantile */
    4081          24 :                 return BATconstant(ngrp == 0 ? 0 : min, average ? TYPE_dbl : tp, dnil, ngrp, TRANSIENT);
    4082             :         }
    4083             : 
    4084          51 :         if (s) {
    4085             :                 /* there is a candidate list, replace b (and g, if
    4086             :                  * given) with just the values we're interested in */
    4087           0 :                 b = BATproject(s, b);
    4088           0 :                 if (b == NULL)
    4089             :                         return NULL;
    4090           0 :                 if (g) {
    4091           0 :                         g = BATproject(s, g);
    4092           0 :                         if (g == NULL)
    4093           0 :                                 goto bunins_failed;
    4094             :                 }
    4095             :         }
    4096             : 
    4097             :         /* we want to sort b so that we can figure out the quantile, but
    4098             :          * if g is given, sort g and subsort b so that we can get the
    4099             :          * quantile for each group */
    4100          51 :         if (g) {
    4101             :                 const oid *restrict grps;
    4102             :                 oid prev;
    4103             :                 BUN p, q, r;
    4104             : 
    4105          19 :                 if (BATtdense(g)) {
    4106             :                         /* singleton groups, so calculating quantile is
    4107             :                          * easy */
    4108           0 :                         if (average)
    4109           0 :                                 bn = BATconvert(b, NULL, TYPE_dbl, abort_on_error, 0, 0, 0);
    4110             :                         else
    4111           0 :                                 bn = COLcopy(b, tp, false, TRANSIENT);
    4112           0 :                         BAThseqbase(bn, g->tseqbase); /* deals with NULL */
    4113           0 :                         if (b != origb)
    4114           0 :                                 BBPunfix(b->batCacheid);
    4115           0 :                         if (g != origg)
    4116           0 :                                 BBPunfix(g->batCacheid);
    4117           0 :                         return bn;
    4118             :                 }
    4119          19 :                 if (BATsort(&t1, &t2, NULL, g, NULL, NULL, false, false, false) != GDK_SUCCEED)
    4120           0 :                         goto bunins_failed;
    4121          19 :                 if (g != origg)
    4122           0 :                         BBPunfix(g->batCacheid);
    4123          19 :                 g = t1;
    4124             : 
    4125          19 :                 if (BATsort(&t1, NULL, NULL, b, t2, g, false, false, false) != GDK_SUCCEED) {
    4126           0 :                         BBPunfix(t2->batCacheid);
    4127           0 :                         goto bunins_failed;
    4128             :                 }
    4129          19 :                 if (b != origb)
    4130           0 :                         BBPunfix(b->batCacheid);
    4131          19 :                 b = t1;
    4132          19 :                 BBPunfix(t2->batCacheid);
    4133             : 
    4134          19 :                 if (average)
    4135           2 :                         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    4136             :                 else
    4137          17 :                         bn = COLnew(min, tp, ngrp, TRANSIENT);
    4138          19 :                 if (bn == NULL)
    4139           0 :                         goto bunins_failed;
    4140             : 
    4141          19 :                 bi = bat_iterator(b);
    4142             : 
    4143          19 :                 grps = (const oid *) Tloc(g, 0);
    4144             :                 /* for each group (r and p are the beginning and end
    4145             :                  * of the current group, respectively) */
    4146          78 :                 for (r = 0, q = BATcount(g); r < q; r = p) {
    4147          59 :                         GDK_CHECK_TIMEOUT(timeoffset, counter,
    4148             :                                         GOTO_LABEL_TIMEOUT_HANDLER(bunins_failed));
    4149             :                         BUN qindex;
    4150          59 :                         prev = grps[r];
    4151             :                         /* search for end of current group (grps is
    4152             :                          * sorted so we can use binary search) */
    4153          59 :                         p = binsearch_oid(NULL, 0, grps, r, q - 1, prev, 1, 1);
    4154          59 :                         if (skip_nils && !b->tnonil) {
    4155             :                                 /* within group, locate start of non-nils */
    4156          32 :                                 r = binsearch(NULL, 0, tp, bi.base,
    4157          32 :                                               bi.vh ? bi.vh->base : NULL,
    4158          32 :                                               bi.width, r, p, nil,
    4159             :                                               1, 1);
    4160             :                         }
    4161          59 :                         if (r == p) {
    4162             :                                 /* no non-nils found */
    4163             :                                 v = dnil;
    4164          12 :                                 nils++;
    4165          47 :                         } else if (average) {
    4166           4 :                                 double f = (p - r - 1) * quantile;
    4167           4 :                                 double lo = floor(f);
    4168           4 :                                 double hi = ceil(f);
    4169           4 :                                 switch (ATOMbasetype(tp)) {
    4170           4 :                                 case TYPE_bte:
    4171           4 :                                         DO_QUANTILE_AVG(bte);
    4172             :                                         break;
    4173           0 :                                 case TYPE_sht:
    4174           0 :                                         DO_QUANTILE_AVG(sht);
    4175             :                                         break;
    4176           0 :                                 case TYPE_int:
    4177           0 :                                         DO_QUANTILE_AVG(int);
    4178             :                                         break;
    4179           0 :                                 case TYPE_lng:
    4180           0 :                                         DO_QUANTILE_AVG(lng);
    4181             :                                         break;
    4182             : #ifdef HAVE_HGE
    4183           0 :                                 case TYPE_hge:
    4184           0 :                                         DO_QUANTILE_AVG(hge);
    4185             :                                         break;
    4186             : #endif
    4187           0 :                                 case TYPE_flt:
    4188           0 :                                         DO_QUANTILE_AVG(flt);
    4189             :                                         break;
    4190           0 :                                 case TYPE_dbl:
    4191           0 :                                         DO_QUANTILE_AVG(dbl);
    4192             :                                         break;
    4193             :                                 }
    4194             :                                 v = &val;
    4195             :                         } else {
    4196             :                                 /* round *down* to nearest integer */
    4197          43 :                                 double f = (p - r - 1) * quantile;
    4198          43 :                                 qindex = r + p - (BUN) (p + 0.5 - f);
    4199             :                                 /* be a little paranoid about the index */
    4200          43 :                                 assert(qindex >= r && qindex <  p);
    4201          43 :                                 v = BUNtail(bi, qindex);
    4202          43 :                                 if (!skip_nils && !b->tnonil)
    4203           0 :                                         nils += (*atomcmp)(v, dnil) == 0;
    4204             :                         }
    4205          59 :                         if (bunfastapp_nocheck(bn, v) != GDK_SUCCEED) {
    4206           0 :                                 bat_iterator_end(&bi);
    4207           0 :                                 goto bunins_failed;
    4208             :                         }
    4209             :                 }
    4210          19 :                 bat_iterator_end(&bi);
    4211          19 :                 nils += ngrp - BATcount(bn);
    4212          19 :                 while (BATcount(bn) < ngrp) {
    4213           0 :                         if (bunfastapp_nocheck(bn, dnil) != GDK_SUCCEED)
    4214           0 :                                 goto bunins_failed;
    4215             :                 }
    4216          19 :                 bn->theap->dirty = true;
    4217          19 :                 BBPunfix(g->batCacheid);
    4218             :         } else {
    4219          32 :                 BUN index, r, p = BATcount(b);
    4220             :                 BAT *pb = NULL;
    4221             :                 const oid *ords;
    4222             :                 Heap *oidxh = NULL;
    4223             : 
    4224          60 :                 bn = COLnew(0, average ? TYPE_dbl : tp, 1, TRANSIENT);
    4225          32 :                 if (bn == NULL)
    4226           0 :                         goto bunins_failed;
    4227             : 
    4228          32 :                 t1 = NULL;
    4229             : 
    4230          32 :                 if (BATcheckorderidx(b)) {
    4231           0 :                         MT_lock_set(&b->batIdxLock);
    4232           0 :                         oidxh = b->torderidx;
    4233           0 :                         if (oidxh != NULL)
    4234           0 :                                 HEAPincref(oidxh);
    4235           0 :                         MT_lock_unset(&b->batIdxLock);
    4236             :                 }
    4237           4 :                 if (oidxh == NULL &&
    4238          32 :                     VIEWtparent(b) &&
    4239           8 :                     (pb = BBP_cache(VIEWtparent(b))) != NULL &&
    4240           4 :                     BATcheckorderidx(pb)) {
    4241             :                         /* no lock on b needed since it's a view */
    4242           0 :                         MT_lock_set(&pb->batIdxLock);
    4243           0 :                         MT_lock_set(&pb->theaplock);
    4244           0 :                         if (pb->tbaseoff == b->tbaseoff &&
    4245           0 :                             BATcount(pb) == BATcount(b) &&
    4246           0 :                             pb->hseqbase == b->hseqbase &&
    4247           0 :                             (oidxh = pb->torderidx) != NULL) {
    4248           0 :                                 HEAPincref(oidxh);
    4249             :                         }
    4250           0 :                         MT_lock_unset(&pb->batIdxLock);
    4251           0 :                         MT_lock_unset(&pb->theaplock);
    4252             :                 }
    4253          32 :                 if (oidxh != NULL) {
    4254           0 :                         MT_thread_setalgorithm(pb ? "using parent oidx" : "using oids");
    4255           0 :                         ords = (const oid *) oidxh->base + ORDERIDXOFF;
    4256             :                 } else {
    4257          32 :                         if (BATsort(NULL, &t1, NULL, b, NULL, g, false, false, false) != GDK_SUCCEED)
    4258           0 :                                 goto bunins_failed;
    4259          32 :                         if (BATtdense(t1))
    4260             :                                 ords = NULL;
    4261             :                         else
    4262          22 :                                 ords = (const oid *) Tloc(t1, 0);
    4263             :                 }
    4264             : 
    4265          32 :                 bi = bat_iterator(b);
    4266             : 
    4267          32 :                 if (skip_nils && !b->tnonil)
    4268          15 :                         r = binsearch(ords, 0, tp, bi.base,
    4269          15 :                                       bi.vh ? bi.vh->base : NULL,
    4270          15 :                                       bi.width, 0, p,
    4271             :                                       nil, 1, 1);
    4272             :                 else
    4273             :                         r = 0;
    4274             : 
    4275          32 :                 if (r == p) {
    4276             :                         /* no non-nil values, so quantile is nil */
    4277             :                         v = dnil;
    4278             :                         nils++;
    4279          29 :                 } else if (average) {
    4280           4 :                         double f = (p - r - 1) * quantile;
    4281           4 :                         double lo = floor(f);
    4282           4 :                         double hi = ceil(f);
    4283           4 :                         switch (ATOMbasetype(tp)) {
    4284           4 :                         case TYPE_bte:
    4285           4 :                                 DO_QUANTILE_AVG(bte);
    4286             :                                 break;
    4287           0 :                         case TYPE_sht:
    4288           0 :                                 DO_QUANTILE_AVG(sht);
    4289             :                                 break;
    4290           0 :                         case TYPE_int:
    4291           0 :                                 DO_QUANTILE_AVG(int);
    4292             :                                 break;
    4293           0 :                         case TYPE_lng:
    4294           0 :                                 DO_QUANTILE_AVG(lng);
    4295             :                                 break;
    4296             : #ifdef HAVE_HGE
    4297           0 :                         case TYPE_hge:
    4298           0 :                                 DO_QUANTILE_AVG(hge);
    4299             :                                 break;
    4300             : #endif
    4301           0 :                         case TYPE_flt:
    4302           0 :                                 DO_QUANTILE_AVG(flt);
    4303             :                                 break;
    4304           0 :                         case TYPE_dbl:
    4305           0 :                                 DO_QUANTILE_AVG(dbl);
    4306             :                                 break;
    4307             :                         }
    4308             :                         v = &val;
    4309             :                 } else {
    4310             :                         double f;
    4311             :                         /* round (p-r-1)*quantile *down* to nearest
    4312             :                          * integer (i.e., 1.49 and 1.5 are rounded to
    4313             :                          * 1, 1.51 is rounded to 2) */
    4314          25 :                         f = (p - r - 1) * quantile;
    4315          25 :                         index = r + p - (BUN) (p + 0.5 - f);
    4316          25 :                         if (ords)
    4317          20 :                                 index = ords[index] - b->hseqbase;
    4318             :                         else
    4319           5 :                                 index = index + t1->tseqbase;
    4320          25 :                         v = BUNtail(bi, index);
    4321          25 :                         nils += (*atomcmp)(v, dnil) == 0;
    4322             :                 }
    4323          32 :                 if (oidxh != NULL)
    4324           0 :                         HEAPdecref(oidxh, false);
    4325          32 :                 if (t1)
    4326          32 :                         BBPunfix(t1->batCacheid);
    4327          32 :                 gdk_return rc = BUNappend(bn, v, false);
    4328          32 :                 bat_iterator_end(&bi);
    4329          32 :                 if (rc != GDK_SUCCEED)
    4330           0 :                         goto bunins_failed;
    4331             :         }
    4332             : 
    4333          51 :         if (b != origb)
    4334          19 :                 BBPunfix(b->batCacheid);
    4335             : 
    4336          51 :         bn->tkey = BATcount(bn) <= 1;
    4337          51 :         bn->tsorted = BATcount(bn) <= 1;
    4338          51 :         bn->trevsorted = BATcount(bn) <= 1;
    4339          51 :         bn->tnil = nils != 0;
    4340          51 :         bn->tnonil = nils == 0;
    4341          51 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOOPTBATFMT ","
    4342             :                   "e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
    4343             :                   ",quantile=%g,average=%s -> " ALGOOPTBATFMT
    4344             :                   "; start " OIDFMT ", count " BUNFMT " (" LLFMT " usec)\n",
    4345             :                   ALGOBATPAR(origb), ALGOOPTBATPAR(origg), ALGOOPTBATPAR(e),
    4346             :                   ALGOOPTBATPAR(s), quantile, average ? "true" : "false",
    4347             :                   ALGOOPTBATPAR(bn), ci.seq, ncand, GDKusec() - t0);
    4348             :         return bn;
    4349             : 
    4350           0 :   bunins_failed:
    4351           0 :         if (b && b != origb)
    4352           0 :                 BBPunfix(b->batCacheid);
    4353           0 :         if (g && g != origg)
    4354           0 :                 BBPunfix(g->batCacheid);
    4355           0 :         if (bn)
    4356           0 :                 BBPunfix(bn->batCacheid);
    4357             :         return NULL;
    4358             : }
    4359             : 
    4360             : BAT *
    4361          28 : BATgroupmedian(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    4362             :                bool skip_nils, bool abort_on_error)
    4363             : {
    4364          28 :         return doBATgroupquantile(b, g, e, s, tp, 0.5,
    4365             :                                   skip_nils, abort_on_error, false);
    4366             : }
    4367             : 
    4368             : BAT *
    4369          29 : BATgroupquantile(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
    4370             :                  bool skip_nils, bool abort_on_error)
    4371             : {
    4372          29 :         return doBATgroupquantile(b, g, e, s, tp, quantile,
    4373             :                                   skip_nils, abort_on_error, false);
    4374             : }
    4375             : 
    4376             : BAT *
    4377           3 : BATgroupmedian_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    4378             :                    bool skip_nils, bool abort_on_error)
    4379             : {
    4380           3 :         return doBATgroupquantile(b, g, e, s, tp, 0.5,
    4381             :                                   skip_nils, abort_on_error, true);
    4382             : }
    4383             : 
    4384             : BAT *
    4385           3 : BATgroupquantile_avg(BAT *b, BAT *g, BAT *e, BAT *s, int tp, double quantile,
    4386             :                      bool skip_nils, bool abort_on_error)
    4387             : {
    4388           3 :         return doBATgroupquantile(b, g, e, s, tp, quantile,
    4389             :                                   skip_nils, abort_on_error, true);
    4390             : }
    4391             : 
    4392             : /* ---------------------------------------------------------------------- */
    4393             : /* standard deviation (both biased and non-biased) */
    4394             : 
    4395             : #define AGGR_STDEV_SINGLE(TYPE)                                 \
    4396             :         do {                                                    \
    4397             :                 TYPE x;                                         \
    4398             :                 TIMEOUT_LOOP_IDX(i, cnt, timeoffset) {          \
    4399             :                         x = ((const TYPE *) values)[i];         \
    4400             :                         if (is_##TYPE##_nil(x))                 \
    4401             :                                 continue;                       \
    4402             :                         n++;                                    \
    4403             :                         delta = (dbl) x - mean;                 \
    4404             :                         mean += delta / n;                      \
    4405             :                         m2 += delta * ((dbl) x - mean);         \
    4406             :                         if (isinf(m2))                          \
    4407             :                                 goto overflow;                  \
    4408             :                 }                                               \
    4409             :                 TIMEOUT_CHECK(timeoffset,                       \
    4410             :                               TIMEOUT_HANDLER(dbl_nil));        \
    4411             :         } while (0)
    4412             : 
    4413             : static dbl
    4414          26 : calcvariance(dbl *restrict avgp, const void *restrict values, BUN cnt, int tp, bool issample)
    4415             : {
    4416             :         BUN n = 0, i;
    4417             :         dbl mean = 0;
    4418             :         dbl m2 = 0;
    4419             :         dbl delta;
    4420             : 
    4421             :         lng timeoffset = 0;
    4422          26 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    4423          26 :         if (qry_ctx != NULL) {
    4424          26 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    4425             :         }
    4426             : 
    4427          26 :         switch (tp) {
    4428             :         case TYPE_bte:
    4429          18 :                 AGGR_STDEV_SINGLE(bte);
    4430             :                 break;
    4431             :         case TYPE_sht:
    4432          18 :                 AGGR_STDEV_SINGLE(sht);
    4433             :                 break;
    4434             :         case TYPE_int:
    4435          47 :                 AGGR_STDEV_SINGLE(int);
    4436             :                 break;
    4437             :         case TYPE_lng:
    4438          36 :                 AGGR_STDEV_SINGLE(lng);
    4439             :                 break;
    4440             : #ifdef HAVE_HGE
    4441             :         case TYPE_hge:
    4442           0 :                 AGGR_STDEV_SINGLE(hge);
    4443             :                 break;
    4444             : #endif
    4445             :         case TYPE_flt:
    4446           0 :                 AGGR_STDEV_SINGLE(flt);
    4447             :                 break;
    4448             :         case TYPE_dbl:
    4449          39 :                 AGGR_STDEV_SINGLE(dbl);
    4450             :                 break;
    4451           0 :         default:
    4452           0 :                 GDKerror("type (%s) not supported.\n", ATOMname(tp));
    4453           0 :                 return dbl_nil;
    4454             :         }
    4455          24 :         if (n <= (BUN) issample) {
    4456           7 :                 if (avgp)
    4457           0 :                         *avgp = dbl_nil;
    4458           7 :                 return dbl_nil;
    4459             :         }
    4460          17 :         if (avgp)
    4461           0 :                 *avgp = mean;
    4462          17 :         return m2 / (n - issample);
    4463           2 :   overflow:
    4464           2 :         GDKerror("22003!overflow in calculation.\n");
    4465           2 :         return dbl_nil;
    4466             : }
    4467             : 
    4468             : dbl
    4469          13 : BATcalcstdev_population(dbl *avgp, BAT *b)
    4470             : {
    4471             :         lng t0 = 0;
    4472             : 
    4473          13 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    4474          13 :         BATiter bi = bat_iterator(b);
    4475          13 :         dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, false);
    4476          13 :         bat_iterator_end(&bi);
    4477          13 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
    4478             :                   ALGOBATPAR(b), GDKusec() - t0);
    4479          13 :         return is_dbl_nil(v) ? dbl_nil : sqrt(v);
    4480             : }
    4481             : 
    4482             : dbl
    4483           7 : BATcalcstdev_sample(dbl *avgp, BAT *b)
    4484             : {
    4485             :         lng t0 = 0;
    4486             : 
    4487           7 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    4488           7 :         BATiter bi = bat_iterator(b);
    4489           7 :         dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, true);
    4490           7 :         bat_iterator_end(&bi);
    4491           7 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
    4492             :                   ALGOBATPAR(b), GDKusec() - t0);
    4493           7 :         return is_dbl_nil(v) ? dbl_nil : sqrt(v);
    4494             : }
    4495             : 
    4496             : dbl
    4497           4 : BATcalcvariance_population(dbl *avgp, BAT *b)
    4498             : {
    4499             :         lng t0 = 0;
    4500             : 
    4501           4 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    4502           4 :         BATiter bi = bat_iterator(b);
    4503           4 :         dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, false);
    4504           4 :         bat_iterator_end(&bi);
    4505           4 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
    4506             :                   ALGOBATPAR(b), GDKusec() - t0);
    4507           4 :         return v;
    4508             : }
    4509             : 
    4510             : dbl
    4511           2 : BATcalcvariance_sample(dbl *avgp, BAT *b)
    4512             : {
    4513             :         lng t0 = 0;
    4514             : 
    4515           2 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    4516           2 :         BATiter bi = bat_iterator(b);
    4517           2 :         dbl v = calcvariance(avgp, bi.base, bi.count, bi.type, true);
    4518           2 :         bat_iterator_end(&bi);
    4519           2 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT " (" LLFMT " usec)\n",
    4520             :                   ALGOBATPAR(b), GDKusec() - t0);
    4521           2 :         return v;
    4522             : }
    4523             : 
    4524             : #define AGGR_COVARIANCE_SINGLE(TYPE)                                    \
    4525             :         do {                                                            \
    4526             :                 TYPE x, y;                                              \
    4527             :                 TIMEOUT_LOOP_IDX(i, cnt, timeoffset) {                  \
    4528             :                         x = ((const TYPE *) v1)[i];                     \
    4529             :                         y = ((const TYPE *) v2)[i];                     \
    4530             :                         if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y))   \
    4531             :                                 continue;                               \
    4532             :                         n++;                                            \
    4533             :                         delta1 = (dbl) x - mean1;                       \
    4534             :                         mean1 += delta1 / n;                            \
    4535             :                         delta2 = (dbl) y - mean2;                       \
    4536             :                         mean2 += delta2 / n;                            \
    4537             :                         m2 += delta1 * ((dbl) y - mean2);               \
    4538             :                         if (isinf(m2))                                  \
    4539             :                                 goto overflow;                          \
    4540             :                 }                                                       \
    4541             :                 TIMEOUT_CHECK(timeoffset,                               \
    4542             :                               TIMEOUT_HANDLER(dbl_nil));                \
    4543             :         } while (0)
    4544             : 
    4545             : static dbl
    4546          12 : calccovariance(const void *v1, const void *v2, BUN cnt, int tp, bool issample)
    4547             : {
    4548             :         BUN n = 0, i;
    4549             :         dbl mean1 = 0, mean2 = 0, m2 = 0, delta1, delta2;
    4550             : 
    4551             :         lng timeoffset = 0;
    4552          12 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    4553          13 :         if (qry_ctx != NULL) {
    4554          13 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    4555             :         }
    4556             : 
    4557             : 
    4558          13 :         switch (tp) {
    4559             :         case TYPE_bte:
    4560           0 :                 AGGR_COVARIANCE_SINGLE(bte);
    4561             :                 break;
    4562             :         case TYPE_sht:
    4563           0 :                 AGGR_COVARIANCE_SINGLE(sht);
    4564             :                 break;
    4565             :         case TYPE_int:
    4566         103 :                 AGGR_COVARIANCE_SINGLE(int);
    4567             :                 break;
    4568             :         case TYPE_lng:
    4569           0 :                 AGGR_COVARIANCE_SINGLE(lng);
    4570             :                 break;
    4571             : #ifdef HAVE_HGE
    4572             :         case TYPE_hge:
    4573           0 :                 AGGR_COVARIANCE_SINGLE(hge);
    4574             :                 break;
    4575             : #endif
    4576             :         case TYPE_flt:
    4577           0 :                 AGGR_COVARIANCE_SINGLE(flt);
    4578             :                 break;
    4579             :         case TYPE_dbl:
    4580          14 :                 AGGR_COVARIANCE_SINGLE(dbl);
    4581             :                 break;
    4582           0 :         default:
    4583           0 :                 GDKerror("type (%s) not supported.\n", ATOMname(tp));
    4584           0 :                 return dbl_nil;
    4585             :         }
    4586          12 :         if (n <= (BUN) issample)
    4587           1 :                 return dbl_nil;
    4588          11 :         return m2 / (n - issample);
    4589           1 :   overflow:
    4590           1 :         GDKerror("22003!overflow in calculation.\n");
    4591           1 :         return dbl_nil;
    4592             : }
    4593             : 
    4594             : dbl
    4595           7 : BATcalccovariance_population(BAT *b1, BAT *b2)
    4596             : {
    4597             :         lng t0 = 0;
    4598             : 
    4599           7 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    4600           7 :         BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
    4601           7 :         dbl v = calccovariance(b1i.base, b2i.base, b1i.count, b1i.type, false);
    4602           7 :         bat_iterator_end(&b1i);
    4603           8 :         bat_iterator_end(&b2i);
    4604           8 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
    4605             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
    4606           8 :         return v;
    4607             : }
    4608             : 
    4609             : dbl
    4610           5 : BATcalccovariance_sample(BAT *b1, BAT *b2)
    4611             : {
    4612             :         lng t0 = 0;
    4613             : 
    4614           5 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    4615           5 :         BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
    4616           5 :         dbl v = calccovariance(b1i.base, b2i.base, b1i.count, b1i.type, true);
    4617           5 :         bat_iterator_end(&b1i);
    4618           5 :         bat_iterator_end(&b2i);
    4619           5 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
    4620             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
    4621           5 :         return v;
    4622             : }
    4623             : 
    4624             : #define AGGR_CORRELATION_SINGLE(TYPE)                                   \
    4625             :         do {                                                            \
    4626             :                 TYPE x, y;                                              \
    4627             :                 TIMEOUT_LOOP_IDX(i, cnt, timeoffset) {                  \
    4628             :                         x = ((const TYPE *) v1)[i];                     \
    4629             :                         y = ((const TYPE *) v2)[i];                     \
    4630             :                         if (is_##TYPE##_nil(x) || is_##TYPE##_nil(y))   \
    4631             :                                 continue;                               \
    4632             :                         n++;                                            \
    4633             :                         delta1 = (dbl) x - mean1;                       \
    4634             :                         mean1 += delta1 / n;                            \
    4635             :                         delta2 = (dbl) y - mean2;                       \
    4636             :                         mean2 += delta2 / n;                            \
    4637             :                         aux = (dbl) y - mean2;                          \
    4638             :                         up += delta1 * aux;                             \
    4639             :                         down1 += delta1 * ((dbl) x - mean1);            \
    4640             :                         down2 += delta2 * aux;                          \
    4641             :                         if (isinf(up) || isinf(down1) || isinf(down2))  \
    4642             :                                 goto overflow;                          \
    4643             :                 }                                                       \
    4644             :                 TIMEOUT_CHECK(timeoffset,                               \
    4645             :                               TIMEOUT_HANDLER(dbl_nil));                \
    4646             :         } while (0)
    4647             : 
    4648             : dbl
    4649          19 : BATcalccorrelation(BAT *b1, BAT *b2)
    4650             : {
    4651          19 :         BUN n = 0, i, cnt = BATcount(b1);
    4652             :         dbl mean1 = 0, mean2 = 0, up = 0, down1 = 0, down2 = 0, delta1, delta2, aux;
    4653          19 :         BATiter b1i = bat_iterator(b1), b2i = bat_iterator(b2);
    4654          19 :         const void *v1 = b1i.base, *v2 = b2i.base;
    4655          19 :         int tp = b1i.type;
    4656             :         lng t0 = 0;
    4657             : 
    4658             :         lng timeoffset = 0;
    4659          19 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    4660          19 :         if (qry_ctx != NULL) {
    4661          19 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    4662             :         }
    4663             : 
    4664          19 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    4665             : 
    4666          19 :         switch (tp) {
    4667             :         case TYPE_bte:
    4668          11 :                 AGGR_CORRELATION_SINGLE(bte);
    4669             :                 break;
    4670             :         case TYPE_sht:
    4671           0 :                 AGGR_CORRELATION_SINGLE(sht);
    4672             :                 break;
    4673             :         case TYPE_int:
    4674          95 :                 AGGR_CORRELATION_SINGLE(int);
    4675             :                 break;
    4676             :         case TYPE_lng:
    4677           0 :                 AGGR_CORRELATION_SINGLE(lng);
    4678             :                 break;
    4679             : #ifdef HAVE_HGE
    4680             :         case TYPE_hge:
    4681           0 :                 AGGR_CORRELATION_SINGLE(hge);
    4682             :                 break;
    4683             : #endif
    4684             :         case TYPE_flt:
    4685           0 :                 AGGR_CORRELATION_SINGLE(flt);
    4686             :                 break;
    4687             :         case TYPE_dbl:
    4688          17 :                 AGGR_CORRELATION_SINGLE(dbl);
    4689             :                 break;
    4690           0 :         default:
    4691           0 :                 bat_iterator_end(&b1i);
    4692           0 :                 bat_iterator_end(&b2i);
    4693           0 :                 GDKerror("type (%s) not supported.\n", ATOMname(tp));
    4694           0 :                 return dbl_nil;
    4695             :         }
    4696          18 :         bat_iterator_end(&b1i);
    4697          18 :         bat_iterator_end(&b2i);
    4698          18 :         if (n != 0 && down1 != 0 && down2 != 0)
    4699           9 :                 aux = (up / n) / (sqrt(down1 / n) * sqrt(down2 / n));
    4700             :         else
    4701           9 :                 aux = dbl_nil;
    4702          18 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT " (" LLFMT " usec)\n",
    4703             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), GDKusec() - t0);
    4704             :         return aux;
    4705           1 :   overflow:
    4706           1 :         bat_iterator_end(&b1i);
    4707           1 :         bat_iterator_end(&b2i);
    4708           1 :         GDKerror("22003!overflow in calculation.\n");
    4709           1 :         return dbl_nil;
    4710             : }
    4711             : 
    4712             : #define AGGR_STDEV(TYPE)                                                \
    4713             :         do {                                                            \
    4714             :                 const TYPE *restrict vals = (const TYPE *) bi.base;     \
    4715             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    4716             :                         i = canditer_next(&ci) - b->hseqbase;            \
    4717             :                         if (gids == NULL ||                             \
    4718             :                             (gids[i] >= min && gids[i] <= max)) { \
    4719             :                                 if (gids)                               \
    4720             :                                         gid = gids[i] - min;            \
    4721             :                                 else                                    \
    4722             :                                         gid = (oid) i;                  \
    4723             :                                 if (is_##TYPE##_nil(vals[i])) {         \
    4724             :                                         if (!skip_nils)                 \
    4725             :                                                 cnts[gid] = BUN_NONE;   \
    4726             :                                 } else if (cnts[gid] != BUN_NONE) {     \
    4727             :                                         cnts[gid]++;                    \
    4728             :                                         delta[gid] = (dbl) vals[i] - mean[gid]; \
    4729             :                                         mean[gid] += delta[gid] / cnts[gid]; \
    4730             :                                         m2[gid] += delta[gid] * ((dbl) vals[i] - mean[gid]); \
    4731             :                                 }                                       \
    4732             :                         }                                               \
    4733             :                 }                                                       \
    4734             :                 TIMEOUT_CHECK(timeoffset,                               \
    4735             :                               GOTO_LABEL_TIMEOUT_HANDLER(alloc_fail));  \
    4736             :                 for (i = 0; i < ngrp; i++) {                         \
    4737             :                         if (cnts[i] == 0 || cnts[i] == BUN_NONE) {      \
    4738             :                                 dbls[i] = dbl_nil;                      \
    4739             :                                 mean[i] = dbl_nil;                      \
    4740             :                                 nils++;                                 \
    4741             :                         } else if (cnts[i] == 1) {                      \
    4742             :                                 dbls[i] = issample ? dbl_nil : 0;       \
    4743             :                                 nils2++;                                \
    4744             :                         } else if (isinf(m2[i])) {                      \
    4745             :                                 goto overflow;                          \
    4746             :                         } else if (variance) {                          \
    4747             :                                 dbls[i] = m2[i] / (cnts[i] - issample); \
    4748             :                         } else {                                        \
    4749             :                                 dbls[i] = sqrt(m2[i] / (cnts[i] - issample)); \
    4750             :                         }                                               \
    4751             :                 }                                                       \
    4752             :         } while (0)
    4753             : 
    4754             : /* Calculate group standard deviation (population (i.e. biased) or
    4755             :  * sample (i.e. non-biased)) with optional candidates list.
    4756             :  *
    4757             :  * Note that this helper function is prepared to return two BATs: one
    4758             :  * (as return value) with the standard deviation per group, and one
    4759             :  * (as return argument) with the average per group.  This isn't
    4760             :  * currently used since it doesn't fit into the mold of grouped
    4761             :  * aggregates. */
    4762             : static BAT *
    4763          18 : dogroupstdev(BAT **avgb, BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    4764             :              bool skip_nils, bool issample, bool variance, const char *func)
    4765             : {
    4766             :         const oid *restrict gids;
    4767             :         oid gid;
    4768             :         oid min, max;
    4769             :         BUN i, ngrp;
    4770             :         BUN nils = 0, nils2 = 0;
    4771             :         BUN *restrict cnts = NULL;
    4772             :         dbl *restrict dbls, *restrict mean, *restrict delta, *restrict m2;
    4773             :         BAT *bn = NULL, *an = NULL;
    4774             :         struct canditer ci;
    4775             :         BUN ncand;
    4776             :         const char *err;
    4777             :         lng t0 = 0;
    4778             : 
    4779             :         lng timeoffset = 0;
    4780          18 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    4781          18 :         if (qry_ctx != NULL) {
    4782          18 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    4783             :         }
    4784             : 
    4785          18 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    4786             : 
    4787          18 :         assert(tp == TYPE_dbl);
    4788             :         (void) tp;              /* compatibility (with other BATgroup*
    4789             :                                  * functions) argument */
    4790             : 
    4791          18 :         if ((err = BATgroupaggrinit(b, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    4792           0 :                 GDKerror("%s: %s\n", func, err);
    4793           0 :                 return NULL;
    4794             :         }
    4795          18 :         if (g == NULL) {
    4796           0 :                 GDKerror("%s: b and g must be aligned\n", func);
    4797           0 :                 return NULL;
    4798             :         }
    4799             : 
    4800          18 :         if (BATcount(b) == 0 || ngrp == 0) {
    4801             :                 /* trivial: no products, so return bat aligned with g
    4802             :                  * with nil in the tail */
    4803           3 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
    4804           3 :                 goto doreturn;
    4805             :         }
    4806             : 
    4807          15 :         if ((e == NULL ||
    4808          15 :              (BATcount(e) == ncand && e->hseqbase == ci.hseq)) &&
    4809           7 :             (BATtdense(g) || (g->tkey && g->tnonil)) &&
    4810           4 :             (issample || b->tnonil)) {
    4811             :                 /* trivial: singleton groups, so all results are equal
    4812             :                  * to zero (population) or nil (sample) */
    4813           4 :                 dbl v = issample ? dbl_nil : 0;
    4814           4 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
    4815           4 :                 goto doreturn;
    4816             :         }
    4817             : 
    4818          11 :         delta = GDKmalloc(ngrp * sizeof(dbl));
    4819          11 :         m2 = GDKmalloc(ngrp * sizeof(dbl));
    4820          11 :         cnts = GDKzalloc(ngrp * sizeof(BUN));
    4821          11 :         if (avgb) {
    4822           0 :                 an = COLnew(0, TYPE_dbl, ngrp, TRANSIENT);
    4823           0 :                 *avgb = an;
    4824           0 :                 if (an == NULL) {
    4825             :                         mean = NULL;
    4826           0 :                         goto alloc_fail;
    4827             :                 }
    4828           0 :                 mean = (dbl *) Tloc(an, 0);
    4829             :         } else {
    4830          11 :                 mean = GDKmalloc(ngrp * sizeof(dbl));
    4831             :         }
    4832          11 :         if (mean == NULL || delta == NULL || m2 == NULL || cnts == NULL)
    4833           0 :                 goto alloc_fail;
    4834             : 
    4835          11 :         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    4836          11 :         if (bn == NULL)
    4837           0 :                 goto alloc_fail;
    4838          11 :         dbls = (dbl *) Tloc(bn, 0);
    4839             : 
    4840     1080049 :         for (i = 0; i < ngrp; i++) {
    4841     1080038 :                 mean[i] = 0;
    4842     1080038 :                 delta[i] = 0;
    4843     1080038 :                 m2[i] = 0;
    4844             :         }
    4845             : 
    4846             :         BATiter bi;
    4847          11 :         bi = bat_iterator(b);
    4848          11 :         if (BATtdense(g))
    4849             :                 gids = NULL;
    4850             :         else
    4851           9 :                 gids = (const oid *) Tloc(g, 0);
    4852             : 
    4853          11 :         switch (b->ttype) {
    4854           0 :         case TYPE_bte:
    4855           0 :                 AGGR_STDEV(bte);
    4856             :                 break;
    4857           0 :         case TYPE_sht:
    4858           0 :                 AGGR_STDEV(sht);
    4859             :                 break;
    4860           4 :         case TYPE_int:
    4861     5760304 :                 AGGR_STDEV(int);
    4862             :                 break;
    4863           0 :         case TYPE_lng:
    4864           0 :                 AGGR_STDEV(lng);
    4865             :                 break;
    4866             : #ifdef HAVE_HGE
    4867           0 :         case TYPE_hge:
    4868           0 :                 AGGR_STDEV(hge);
    4869             :                 break;
    4870             : #endif
    4871           0 :         case TYPE_flt:
    4872           0 :                 AGGR_STDEV(flt);
    4873             :                 break;
    4874           7 :         case TYPE_dbl:
    4875          94 :                 AGGR_STDEV(dbl);
    4876             :                 break;
    4877           0 :         default:
    4878           0 :                 bat_iterator_end(&bi);
    4879           0 :                 if (an)
    4880           0 :                         BBPreclaim(an);
    4881             :                 else
    4882           0 :                         GDKfree(mean);
    4883           0 :                 GDKfree(delta);
    4884           0 :                 GDKfree(m2);
    4885           0 :                 GDKfree(cnts);
    4886           0 :                 BBPunfix(bn->batCacheid);
    4887           0 :                 GDKerror("%s: type (%s) not supported.\n",
    4888             :                          func, ATOMname(b->ttype));
    4889           0 :                 return NULL;
    4890             :         }
    4891           9 :         bat_iterator_end(&bi);
    4892           9 :         if (an) {
    4893           0 :                 BATsetcount(an, ngrp);
    4894           0 :                 an->tkey = ngrp <= 1;
    4895           0 :                 an->tsorted = ngrp <= 1;
    4896           0 :                 an->trevsorted = ngrp <= 1;
    4897           0 :                 an->tnil = nils != 0;
    4898           0 :                 an->tnonil = nils == 0;
    4899             :         } else {
    4900           9 :                 GDKfree(mean);
    4901             :         }
    4902           9 :         if (issample)
    4903           3 :                 nils += nils2;
    4904           9 :         GDKfree(delta);
    4905           9 :         GDKfree(m2);
    4906           9 :         GDKfree(cnts);
    4907           9 :         BATsetcount(bn, ngrp);
    4908           9 :         bn->tkey = ngrp <= 1;
    4909           9 :         bn->tsorted = ngrp <= 1;
    4910           9 :         bn->trevsorted = ngrp <= 1;
    4911           9 :         bn->tnil = nils != 0;
    4912           9 :         bn->tnonil = nils == 0;
    4913          16 :   doreturn:
    4914          16 :         TRC_DEBUG(ALGO, "b=" ALGOBATFMT ",g=" ALGOBATFMT ",e=" ALGOOPTBATFMT
    4915             :                   ",s=" ALGOOPTBATFMT
    4916             :                   ",skip_nils=%s,issample=%s,variance=%s -> " ALGOOPTBATFMT
    4917             :                   ",avgb=" ALGOOPTBATFMT " (%s -- " LLFMT " usec)\n",
    4918             :                   ALGOBATPAR(b), ALGOBATPAR(g), ALGOOPTBATPAR(e),
    4919             :                   ALGOOPTBATPAR(s),
    4920             :                   skip_nils ? "true" : "false",
    4921             :                   issample ? "true" : "false",
    4922             :                   variance ? "true" : "false",
    4923             :                   ALGOOPTBATPAR(bn), ALGOOPTBATPAR(an),
    4924             :                   func, GDKusec() - t0);
    4925             :         return bn;
    4926           2 :   overflow:
    4927           2 :         bat_iterator_end(&bi);
    4928           2 :         GDKerror("22003!overflow in calculation.\n");
    4929           2 :   alloc_fail:
    4930           2 :         if (an)
    4931           0 :                 BBPreclaim(an);
    4932             :         else
    4933           2 :                 GDKfree(mean);
    4934           2 :         BBPreclaim(bn);
    4935           2 :         GDKfree(delta);
    4936           2 :         GDKfree(m2);
    4937           2 :         GDKfree(cnts);
    4938           2 :         return NULL;
    4939             : }
    4940             : 
    4941             : BAT *
    4942           7 : BATgroupstdev_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    4943             :                      bool skip_nils, bool abort_on_error)
    4944             : {
    4945             :         (void) abort_on_error;
    4946           7 :         return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, false,
    4947             :                             __func__);
    4948             : }
    4949             : 
    4950             : BAT *
    4951           5 : BATgroupstdev_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    4952             :                          bool skip_nils, bool abort_on_error)
    4953             : {
    4954             :         (void) abort_on_error;
    4955           5 :         return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, false,
    4956             :                             __func__);
    4957             : }
    4958             : 
    4959             : BAT *
    4960           1 : BATgroupvariance_sample(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    4961             :                         bool skip_nils, bool abort_on_error)
    4962             : {
    4963             :         (void) abort_on_error;
    4964           1 :         return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, true, true,
    4965             :                             __func__);
    4966             : }
    4967             : 
    4968             : BAT *
    4969           5 : BATgroupvariance_population(BAT *b, BAT *g, BAT *e, BAT *s, int tp,
    4970             :                             bool skip_nils, bool abort_on_error)
    4971             : {
    4972             :         (void) abort_on_error;
    4973           5 :         return dogroupstdev(NULL, b, g, e, s, tp, skip_nils, false, true,
    4974             :                             __func__);
    4975             : }
    4976             : 
    4977             : #define AGGR_COVARIANCE(TYPE)                                           \
    4978             :         do {                                                            \
    4979             :                 const TYPE *vals1 = (const TYPE *) b1i.base;            \
    4980             :                 const TYPE *vals2 = (const TYPE *) b2i.base;            \
    4981             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    4982             :                         i = canditer_next(&ci) - b1->hseqbase;           \
    4983             :                         if (gids == NULL ||                             \
    4984             :                             (gids[i] >= min && gids[i] <= max)) { \
    4985             :                                 if (gids)                               \
    4986             :                                         gid = gids[i] - min;            \
    4987             :                                 else                                    \
    4988             :                                         gid = (oid) i;                  \
    4989             :                                 if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
    4990             :                                         if (!skip_nils)                 \
    4991             :                                                 cnts[gid] = BUN_NONE;   \
    4992             :                                 } else if (cnts[gid] != BUN_NONE) {     \
    4993             :                                         cnts[gid]++;                    \
    4994             :                                         delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
    4995             :                                         mean1[gid] += delta1[gid] / cnts[gid]; \
    4996             :                                         delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
    4997             :                                         mean2[gid] += delta2[gid] / cnts[gid]; \
    4998             :                                         m2[gid] += delta1[gid] * ((dbl) vals2[i] - mean2[gid]); \
    4999             :                                 }                                       \
    5000             :                         }                                               \
    5001             :                 }                                                       \
    5002             :                 TIMEOUT_CHECK(timeoffset,                               \
    5003             :                               GOTO_LABEL_TIMEOUT_HANDLER(alloc_fail));  \
    5004             :                 for (i = 0; i < ngrp; i++) {                         \
    5005             :                         if (cnts[i] == 0 || cnts[i] == BUN_NONE) {      \
    5006             :                                 dbls[i] = dbl_nil;                      \
    5007             :                                 nils++;                                 \
    5008             :                         } else if (cnts[i] == 1) {                      \
    5009             :                                 dbls[i] = issample ? dbl_nil : 0;       \
    5010             :                                 nils2++;                                \
    5011             :                         } else if (isinf(m2[i])) {                      \
    5012             :                                 goto overflow;                          \
    5013             :                         } else {                                        \
    5014             :                                 dbls[i] = m2[i] / (cnts[i] - issample); \
    5015             :                         }                                               \
    5016             :                 }                                                       \
    5017             :         } while (0)
    5018             : 
    5019             : static BAT *
    5020          60 : dogroupcovariance(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp,
    5021             :                   bool skip_nils, bool issample, const char *func)
    5022             : {
    5023             :         const oid *restrict gids;
    5024             :         oid gid, min, max;
    5025             :         BUN i, ngrp, nils = 0, nils2 = 0, ncand;
    5026             :         BUN *restrict cnts = NULL;
    5027             :         dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict m2;
    5028             :         BAT *bn = NULL;
    5029             :         struct canditer ci;
    5030             :         const char *err;
    5031             :         lng t0 = 0;
    5032             : 
    5033             :         lng timeoffset = 0;
    5034          60 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    5035          60 :         if (qry_ctx != NULL) {
    5036          60 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    5037             :         }
    5038             : 
    5039             : 
    5040          60 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    5041             : 
    5042         180 :         assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
    5043             :         (void) tp;
    5044             : 
    5045          60 :         if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    5046           0 :                 GDKerror("%s: %s\n", func, err);
    5047           0 :                 return NULL;
    5048             :         }
    5049          60 :         if (g == NULL) {
    5050           0 :                 GDKerror("%s: b1, b2 and g must be aligned\n", func);
    5051           0 :                 return NULL;
    5052             :         }
    5053             : 
    5054          60 :         if (BATcount(b1) == 0 || ngrp == 0) {
    5055           0 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
    5056           0 :                 goto doreturn;
    5057             :         }
    5058             : 
    5059          60 :         if ((e == NULL ||
    5060          60 :              (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
    5061          20 :             (BATtdense(g) || (g->tkey && g->tnonil)) &&
    5062          10 :             (issample || (b1->tnonil && b2->tnonil))) {
    5063             :                 /* trivial: singleton groups, so all results are equal
    5064             :                  * to zero (population) or nil (sample) */
    5065          11 :                 dbl v = issample ? dbl_nil : 0;
    5066          11 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &v, ngrp, TRANSIENT);
    5067          11 :                 goto doreturn;
    5068             :         }
    5069             : 
    5070          49 :         delta1 = GDKmalloc(ngrp * sizeof(dbl));
    5071          49 :         delta2 = GDKmalloc(ngrp * sizeof(dbl));
    5072          49 :         m2 = GDKmalloc(ngrp * sizeof(dbl));
    5073          49 :         cnts = GDKzalloc(ngrp * sizeof(BUN));
    5074          49 :         mean1 = GDKmalloc(ngrp * sizeof(dbl));
    5075          49 :         mean2 = GDKmalloc(ngrp * sizeof(dbl));
    5076             : 
    5077          49 :         if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || m2 == NULL || cnts == NULL)
    5078           0 :                 goto alloc_fail;
    5079             : 
    5080         379 :         for (i = 0; i < ngrp; i++) {
    5081         330 :                 m2[i] = 0;
    5082         330 :                 mean1[i] = 0;
    5083         330 :                 mean2[i] = 0;
    5084             :         }
    5085             : 
    5086          49 :         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    5087          49 :         if (bn == NULL)
    5088           0 :                 goto alloc_fail;
    5089          49 :         dbls = (dbl *) Tloc(bn, 0);
    5090             : 
    5091          49 :         if (!g || BATtdense(g))
    5092             :                 gids = NULL;
    5093             :         else
    5094          49 :                 gids = (const oid *) Tloc(g, 0);
    5095             : 
    5096             :         BATiter b1i, b2i;
    5097          49 :         b1i = bat_iterator(b1);
    5098          49 :         b2i = bat_iterator(b2);
    5099          49 :         switch (b1->ttype) {
    5100          14 :         case TYPE_bte:
    5101         260 :                 AGGR_COVARIANCE(bte);
    5102             :                 break;
    5103           0 :         case TYPE_sht:
    5104           0 :                 AGGR_COVARIANCE(sht);
    5105             :                 break;
    5106          35 :         case TYPE_int:
    5107         654 :                 AGGR_COVARIANCE(int);
    5108             :                 break;
    5109           0 :         case TYPE_lng:
    5110           0 :                 AGGR_COVARIANCE(lng);
    5111             :                 break;
    5112             : #ifdef HAVE_HGE
    5113           0 :         case TYPE_hge:
    5114           0 :                 AGGR_COVARIANCE(hge);
    5115             :                 break;
    5116             : #endif
    5117           0 :         case TYPE_flt:
    5118           0 :                 AGGR_COVARIANCE(flt);
    5119             :                 break;
    5120           0 :         case TYPE_dbl:
    5121           0 :                 AGGR_COVARIANCE(dbl);
    5122             :                 break;
    5123           0 :         default:
    5124           0 :                 bat_iterator_end(&b1i);
    5125           0 :                 bat_iterator_end(&b2i);
    5126           0 :                 GDKerror("%s: type (%s) not supported.\n", func, ATOMname(b1->ttype));
    5127           0 :                 goto alloc_fail;
    5128             :         }
    5129          49 :         bat_iterator_end(&b1i);
    5130          49 :         bat_iterator_end(&b2i);
    5131          49 :         GDKfree(mean1);
    5132          49 :         GDKfree(mean2);
    5133             : 
    5134          49 :         if (issample)
    5135          20 :                 nils += nils2;
    5136          49 :         GDKfree(delta1);
    5137          49 :         GDKfree(delta2);
    5138          49 :         GDKfree(m2);
    5139          49 :         GDKfree(cnts);
    5140          49 :         BATsetcount(bn, ngrp);
    5141          49 :         bn->tkey = ngrp <= 1;
    5142          49 :         bn->tsorted = ngrp <= 1;
    5143          49 :         bn->trevsorted = ngrp <= 1;
    5144          49 :         bn->tnil = nils != 0;
    5145          49 :         bn->tnonil = nils == 0;
    5146          60 :   doreturn:
    5147          60 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
    5148             :                   ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
    5149             :                   ",skip_nils=%s,issample=%s -> " ALGOOPTBATFMT
    5150             :                   " (%s -- " LLFMT " usec)\n",
    5151             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
    5152             :                   ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
    5153             :                   skip_nils ? "true" : "false",
    5154             :                   issample ? "true" : "false",
    5155             :                   ALGOOPTBATPAR(bn),
    5156             :                   func, GDKusec() - t0);
    5157             :         return bn;
    5158           0 :   overflow:
    5159           0 :         bat_iterator_end(&b1i);
    5160           0 :         bat_iterator_end(&b2i);
    5161           0 :         GDKerror("22003!overflow in calculation.\n");
    5162           0 :   alloc_fail:
    5163           0 :         BBPreclaim(bn);
    5164           0 :         GDKfree(mean1);
    5165           0 :         GDKfree(mean2);
    5166           0 :         GDKfree(delta1);
    5167           0 :         GDKfree(delta2);
    5168           0 :         GDKfree(m2);
    5169           0 :         GDKfree(cnts);
    5170           0 :         return NULL;
    5171             : }
    5172             : 
    5173             : BAT *
    5174          30 : BATgroupcovariance_sample(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    5175             : {
    5176             :         (void) abort_on_error;
    5177          30 :         return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, true,
    5178             :                                  __func__);
    5179             : }
    5180             : 
    5181             : BAT *
    5182          30 : BATgroupcovariance_population(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    5183             : {
    5184             :         (void) abort_on_error;
    5185          30 :         return dogroupcovariance(b1, b2, g, e, s, tp, skip_nils, false,
    5186             :                                  __func__);
    5187             : }
    5188             : 
    5189             : #define AGGR_CORRELATION(TYPE)                                          \
    5190             :         do {                                                            \
    5191             :                 const TYPE *vals1 = (const TYPE *) b1i.base;            \
    5192             :                 const TYPE *vals2 = (const TYPE *) b2i.base;            \
    5193             :                 TIMEOUT_LOOP(ncand, timeoffset) {                       \
    5194             :                         i = canditer_next(&ci) - b1->hseqbase;           \
    5195             :                         if (gids == NULL ||                             \
    5196             :                             (gids[i] >= min && gids[i] <= max)) { \
    5197             :                                 if (gids)                               \
    5198             :                                         gid = gids[i] - min;            \
    5199             :                                 else                                    \
    5200             :                                         gid = (oid) i;                  \
    5201             :                                 if (is_##TYPE##_nil(vals1[i]) || is_##TYPE##_nil(vals2[i])) { \
    5202             :                                         if (!skip_nils)                 \
    5203             :                                                 cnts[gid] = BUN_NONE;   \
    5204             :                                 } else if (cnts[gid] != BUN_NONE) {     \
    5205             :                                         cnts[gid]++;                    \
    5206             :                                         delta1[gid] = (dbl) vals1[i] - mean1[gid]; \
    5207             :                                         mean1[gid] += delta1[gid] / cnts[gid]; \
    5208             :                                         delta2[gid] = (dbl) vals2[i] - mean2[gid]; \
    5209             :                                         mean2[gid] += delta2[gid] / cnts[gid]; \
    5210             :                                         aux = (dbl) vals2[i] - mean2[gid]; \
    5211             :                                         up[gid] += delta1[gid] * aux;   \
    5212             :                                         down1[gid] += delta1[gid] * ((dbl) vals1[i] - mean1[gid]); \
    5213             :                                         down2[gid] += delta2[gid] * aux; \
    5214             :                                 }                                       \
    5215             :                         }                                               \
    5216             :                 }                                                       \
    5217             :                 TIMEOUT_CHECK(timeoffset,                               \
    5218             :                               GOTO_LABEL_TIMEOUT_HANDLER(alloc_fail));  \
    5219             :                 for (i = 0; i < ngrp; i++) {                         \
    5220             :                         if (cnts[i] <= 1 || cnts[i] == BUN_NONE || down1[i] == 0 || down2[i] == 0) { \
    5221             :                                 dbls[i] = dbl_nil;                      \
    5222             :                                 nils++;                                 \
    5223             :                         } else if (isinf(up[i]) || isinf(down1[i]) || isinf(down2[i])) { \
    5224             :                                 goto overflow;                          \
    5225             :                         } else {                                        \
    5226             :                                 dbls[i] = (up[i] / cnts[i]) / (sqrt(down1[i] / cnts[i]) * sqrt(down2[i] / cnts[i])); \
    5227             :                                 assert(!is_dbl_nil(dbls[i]));           \
    5228             :                         }                                               \
    5229             :                 }                                                       \
    5230             :         } while (0)
    5231             : 
    5232             : BAT *
    5233          46 : BATgroupcorrelation(BAT *b1, BAT *b2, BAT *g, BAT *e, BAT *s, int tp, bool skip_nils, bool abort_on_error)
    5234             : {
    5235             :         const oid *restrict gids;
    5236             :         oid gid, min, max;
    5237             :         BUN i, ngrp, nils = 0, ncand;
    5238             :         BUN *restrict cnts = NULL;
    5239             :         dbl *restrict dbls, *restrict mean1, *restrict mean2, *restrict delta1, *restrict delta2, *restrict up, *restrict down1, *restrict down2, aux;
    5240             :         BAT *bn = NULL;
    5241             :         struct canditer ci;
    5242             :         const char *err;
    5243             :         lng t0 = 0;
    5244             : 
    5245             :         lng timeoffset = 0;
    5246          46 :         QryCtx *qry_ctx = MT_thread_get_qry_ctx();
    5247          46 :         if (qry_ctx != NULL) {
    5248          46 :                 timeoffset = (qry_ctx->starttime && qry_ctx->querytimeout) ? (qry_ctx->starttime + qry_ctx->querytimeout) : 0;
    5249             :         }
    5250             : 
    5251          46 :         TRC_DEBUG_IF(ALGO) t0 = GDKusec();
    5252             : 
    5253         138 :         assert(tp == TYPE_dbl && BATcount(b1) == BATcount(b2) && b1->ttype == b2->ttype && BATtdense(b1) == BATtdense(b2));
    5254             :         (void) tp;
    5255             :         (void) abort_on_error;
    5256             : 
    5257          46 :         if ((err = BATgroupaggrinit(b1, g, e, s, &min, &max, &ngrp, &ci, &ncand)) != NULL) {
    5258           0 :                 GDKerror("%s\n", err);
    5259           0 :                 return NULL;
    5260             :         }
    5261          46 :         if (g == NULL) {
    5262           0 :                 GDKerror("b1, b2 and g must be aligned\n");
    5263           0 :                 return NULL;
    5264             :         }
    5265             : 
    5266          46 :         if (BATcount(b1) == 0 || ngrp == 0) {
    5267           1 :                 bn = BATconstant(ngrp == 0 ? 0 : min, TYPE_dbl, &dbl_nil, ngrp, TRANSIENT);
    5268           1 :                 goto doreturn;
    5269             :         }
    5270             : 
    5271          45 :         if ((e == NULL ||
    5272          45 :              (BATcount(e) == BATcount(b1) && (e->hseqbase == b1->hseqbase || e->hseqbase == b2->hseqbase))) &&
    5273          11 :             (BATtdense(g) || (g->tkey && g->tnonil))) {
    5274          11 :                 dbl v = dbl_nil;
    5275          11 :                 bn = BATconstant(min, TYPE_dbl, &v, ngrp, TRANSIENT);
    5276          11 :                 goto doreturn;
    5277             :         }
    5278             : 
    5279          34 :         delta1 = GDKmalloc(ngrp * sizeof(dbl));
    5280          34 :         delta2 = GDKmalloc(ngrp * sizeof(dbl));
    5281          34 :         up = GDKmalloc(ngrp * sizeof(dbl));
    5282          34 :         down1 = GDKmalloc(ngrp * sizeof(dbl));
    5283          34 :         down2 = GDKmalloc(ngrp * sizeof(dbl));
    5284          34 :         cnts = GDKzalloc(ngrp * sizeof(BUN));
    5285          34 :         mean1 = GDKmalloc(ngrp * sizeof(dbl));
    5286          34 :         mean2 = GDKmalloc(ngrp * sizeof(dbl));
    5287             : 
    5288          34 :         if (mean1 == NULL || mean2 == NULL || delta1 == NULL || delta2 == NULL || up == NULL || down1 == NULL || down2 == NULL || cnts == NULL)
    5289           0 :                 goto alloc_fail;
    5290             : 
    5291         239 :         for (i = 0; i < ngrp; i++) {
    5292         205 :                 up[i] = 0;
    5293         205 :                 down1[i] = 0;
    5294         205 :                 down2[i] = 0;
    5295         205 :                 mean1[i] = 0;
    5296         205 :                 mean2[i] = 0;
    5297             :         }
    5298             : 
    5299          34 :         bn = COLnew(min, TYPE_dbl, ngrp, TRANSIENT);
    5300          34 :         if (bn == NULL)
    5301           0 :                 goto alloc_fail;
    5302          34 :         dbls = (dbl *) Tloc(bn, 0);
    5303             : 
    5304          34 :         if (!g || BATtdense(g))
    5305             :                 gids = NULL;
    5306             :         else
    5307          34 :                 gids = (const oid *) Tloc(g, 0);
    5308             : 
    5309             :         BATiter b1i, b2i;
    5310          34 :         b1i = bat_iterator(b1);
    5311          34 :         b2i = bat_iterator(b2);
    5312          34 :         switch (b1->ttype) {
    5313           6 :         case TYPE_bte:
    5314         108 :                 AGGR_CORRELATION(bte);
    5315             :                 break;
    5316           0 :         case TYPE_sht:
    5317           0 :                 AGGR_CORRELATION(sht);
    5318             :                 break;
    5319          23 :         case TYPE_int:
    5320         825 :                 AGGR_CORRELATION(int);
    5321             :                 break;
    5322           4 :         case TYPE_lng:
    5323         316 :                 AGGR_CORRELATION(lng);
    5324             :                 break;
    5325             : #ifdef HAVE_HGE
    5326           0 :         case TYPE_hge:
    5327           0 :                 AGGR_CORRELATION(hge);
    5328             :                 break;
    5329             : #endif
    5330           0 :         case TYPE_flt:
    5331           0 :                 AGGR_CORRELATION(flt);
    5332             :                 break;
    5333           1 :         case TYPE_dbl:
    5334           4 :                 AGGR_CORRELATION(dbl);
    5335             :                 break;
    5336           0 :         default:
    5337           0 :                 bat_iterator_end(&b1i);
    5338           0 :                 bat_iterator_end(&b2i);
    5339           0 :                 GDKerror("type (%s) not supported.\n", ATOMname(b1->ttype));
    5340           0 :                 goto alloc_fail;
    5341             :         }
    5342          33 :         bat_iterator_end(&b1i);
    5343          33 :         bat_iterator_end(&b2i);
    5344          33 :         GDKfree(mean1);
    5345          33 :         GDKfree(mean2);
    5346          33 :         GDKfree(delta1);
    5347          33 :         GDKfree(delta2);
    5348          33 :         GDKfree(up);
    5349          33 :         GDKfree(down1);
    5350          33 :         GDKfree(down2);
    5351          33 :         GDKfree(cnts);
    5352          33 :         BATsetcount(bn, ngrp);
    5353          33 :         bn->tkey = ngrp <= 1;
    5354          33 :         bn->tsorted = ngrp <= 1;
    5355          33 :         bn->trevsorted = ngrp <= 1;
    5356          33 :         bn->tnil = nils != 0;
    5357          33 :         bn->tnonil = nils == 0;
    5358          45 :   doreturn:
    5359          45 :         TRC_DEBUG(ALGO, "b1=" ALGOBATFMT ",b2=" ALGOBATFMT ",g=" ALGOBATFMT
    5360             :                   ",e=" ALGOOPTBATFMT ",s=" ALGOOPTBATFMT
    5361             :                   ",skip_nils=%s -> " ALGOOPTBATFMT
    5362             :                   " (" LLFMT " usec)\n",
    5363             :                   ALGOBATPAR(b1), ALGOBATPAR(b2), ALGOBATPAR(g),
    5364             :                   ALGOOPTBATPAR(e), ALGOOPTBATPAR(s),
    5365             :                   skip_nils ? "true" : "false",
    5366             :                   ALGOOPTBATPAR(bn),
    5367             :                   GDKusec() - t0);
    5368             :         return bn;
    5369           1 :   overflow:
    5370           1 :         bat_iterator_end(&b1i);
    5371           1 :         bat_iterator_end(&b2i);
    5372           1 :         GDKerror("22003!overflow in calculation.\n");
    5373           1 :   alloc_fail:
    5374           1 :         BBPreclaim(bn);
    5375           1 :         GDKfree(mean1);
    5376           1 :         GDKfree(mean2);
    5377           1 :         GDKfree(delta1);
    5378           1 :         GDKfree(delta2);
    5379           1 :         GDKfree(up);
    5380           1 :         GDKfree(down1);
    5381           1 :         GDKfree(down2);
    5382           1 :         GDKfree(cnts);
    5383           1 :         return NULL;
    5384             : }

Generated by: LCOV version 1.14