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

Generated by: LCOV version 1.14