text: /*-------------------------------------------------------------------------
   *
   * hyperloglog.c
   *	  HyperLogLog cardinality estimator
   *
   * Portions Copyright (c) 2014-2020, PostgreSQL Global Development Group
   *
   * Based on Hideaki Ohno's C++ implementation.  This is probably not ideally
   * suited to estimating the cardinality of very large sets;  in particular, we
   * have not attempted to further optimize the implementation as described in
   * the Heule, Nunkesser and Hall paper "HyperLogLog in Practice: Algorithmic
   * Engineering of a State of The Art Cardinality Estimation Algorithm".
   *
   * A sparse representation of HyperLogLog state is used, with fixed space
   * overhead.
   *
   * The copyright terms of Ohno's original version (the MIT license) follow.
   *
   * IDENTIFICATION
   *	  src/backend/lib/hyperloglog.c
   *
   *-------------------------------------------------------------------------
   */
  
  /*
   * Copyright (c) 2013 Hideaki Ohno <hide.o.j55{at}gmail.com>
   *
   * Permission is hereby granted, free of charge, to any person obtaining a copy
   * of this software and associated documentation files (the 'Software'), to
   * deal in the Software without restriction, including without limitation the
   * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
   * sell copies of the Software, and to permit persons to whom the Software is
   * furnished to do so, subject to the following conditions:
   *
   * The above copyright notice and this permission notice shall be included in
   * all copies or substantial portions of the Software.
   *
   * THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
   * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
   * IN THE SOFTWARE.
   */
  
  #include "postgres.h"
  
  #include <math.h>
  
  #include "lib/hyperloglog.h"
  #include "port/pg_bitutils.h"
  
  #define POW_2_32			(4294967296.0)
  #define NEG_POW_2_32		(-4294967296.0)
  
  static inline uint8 rho(uint32 x, uint8 b);
  
  /*
   * Initialize HyperLogLog track state, by bit width
   *
   * bwidth is bit width (so register size will be 2 to the power of bwidth).
   * Must be between 4 and 16 inclusive.
   */
  void
  initHyperLogLog(hyperLogLogState *cState, uint8 bwidth)
  {
  	double		alpha;
  
  	if (bwidth < 4 || bwidth > 16)
  		elog(ERROR, "bit width must be between 4 and 16 inclusive");
  
  	cState->registerWidth = bwidth;
  	cState->nRegisters = (Size) 1 << bwidth;
  	cState->arrSize = sizeof(uint8) * cState->nRegisters + 1;
  
  	/*
  	 * Initialize hashes array to zero, not negative infinity, per discussion
  	 * of the coupon collector problem in the HyperLogLog paper
  	 */
  	cState->hashesArr = palloc0(cState->arrSize);
  
  	/*
  	 * "alpha" is a value that for each possible number of registers (m) is
  	 * used to correct a systematic multiplicative bias present in m ^ 2 Z (Z
  	 * is "the indicator function" through which we finally compute E,
  	 * estimated cardinality).
  	 */
  	switch (cState->nRegisters)
  	{
  		case 16:
  			alpha = 0.673;
  			break;
  		case 32:
  			alpha = 0.697;
  			break;
  		case 64:
  			alpha = 0.709;
  			break;
  		default:
  			alpha = 0.7213 / (1.0 + 1.079 / cState->nRegisters);
  	}
  
  	/*
  	 * Precalculate alpha m ^ 2, later used to generate "raw" HyperLogLog
  	 * estimate E
  	 */
  	cState->alphaMM = alpha * cState->nRegisters * cState->nRegisters;
  }
  
  /*
   * Initialize HyperLogLog track state, by error rate
   *
   * Instead of specifying bwidth (number of bits used for addressing the
   * register), this method allows sizing the counter for particular error
   * rate using a simple formula from the paper:
   *
   *	 e = 1.04 / sqrt(m)
   *
   * where 'm' is the number of registers, i.e. (2^bwidth). The method
   * finds the lowest bwidth with 'e' below the requested error rate, and
   * then uses it to initialize the counter.
   *
   * As bwidth has to be between 4 and 16, the worst possible error rate
   * is between ~25% (bwidth=4) and 0.4% (bwidth=16).
   */
  void
  initHyperLogLogError(hyperLogLogState *cState, double error)
  {
  	uint8		bwidth = 4;
  
  	while (bwidth < 16)
  	{
  		double		m = (Size) 1 << bwidth;
  
  		if (1.04 / sqrt(m) < error)
  			break;
  		bwidth++;
  	}
  
  	initHyperLogLog(cState, bwidth);
  }
  
  /*
   * Free HyperLogLog track state
   *
   * Releases allocated resources, but not the state itself (in case it's not
   * allocated by palloc).
   */
  void
  freeHyperLogLog(hyperLogLogState *cState)
  {
  	Assert(cState->hashesArr != NULL);
  	pfree(cState->hashesArr);
  }
  
  /*
   * Adds element to the estimator, from caller-supplied hash.
   *
   * It is critical that the hash value passed be an actual hash value, typically
   * generated using hash_any().  The algorithm relies on a specific bit-pattern
   * observable in conjunction with stochastic averaging.  There must be a
   * uniform distribution of bits in hash values for each distinct original value
   * observed.
   */
  void
  addHyperLogLog(hyperLogLogState *cState, uint32 hash)
  {
  	uint8		count;
  	uint32		index;
  
  	/* Use the first "k" (registerWidth) bits as a zero based index */
  	index = hash >> (BITS_PER_BYTE * sizeof(uint32) - cState->registerWidth);
  
  	/* Compute the rank of the remaining 32 - "k" (registerWidth) bits */
  	count = rho(hash << cState->registerWidth,
  				BITS_PER_BYTE * sizeof(uint32) - cState->registerWidth);
  
  	cState->hashesArr[index] = Max(count, cState->hashesArr[index]);
  }
  
  /*
   * Estimates cardinality, based on elements added so far
   */
  double
  estimateHyperLogLog(hyperLogLogState *cState)
  {
  	double		result;
  	double		sum = 0.0;
  	int			i;
  
  	for (i = 0; i < cState->nRegisters; i++)
  	{
  		sum += 1.0 / pow(2.0, cState->hashesArr[i]);
  	}
  
  	/* result set to "raw" HyperLogLog estimate (E in the HyperLogLog paper) */
  	result = cState->alphaMM / sum;
  
  	if (result <= (5.0 / 2.0) * cState->nRegisters)
  	{
  		/* Small range correction */
  		int			zero_count = 0;
  
  		for (i = 0; i < cState->nRegisters; i++)
  		{
  			if (cState->hashesArr[i] == 0)
  				zero_count++;
  		}
  
  		if (zero_count != 0)
  			result = cState->nRegisters * log((double) cState->nRegisters /
  											  zero_count);
  	}
  	else if (result > (1.0 / 30.0) * POW_2_32)
  	{
  		/* Large range correction */
  		result = NEG_POW_2_32 * log(1.0 - (result / POW_2_32));
  	}
  
  	return result;
  }
  
  /*
   * Worker for addHyperLogLog().
   *
   * Calculates the position of the first set bit in first b bits of x argument
   * starting from the first, reading from most significant to least significant
   * bits.
   *
   * Example (when considering fist 10 bits of x):
   *
   * rho(x = 0b1000000000)   returns 1
   * rho(x = 0b0010000000)   returns 3
   * rho(x = 0b0000000000)   returns b + 1
   *
   * "The binary address determined by the first b bits of x"
   *
   * Return value "j" used to index bit pattern to watch.
   */
  static inline uint8
  rho(uint32 x, uint8 b)
  {
  	uint8		j = 1;
  
  	if (x == 0)
  		return b + 1;
  
  	j = 32 - pg_leftmost_one_pos32(x);
  
  	if (j > b)
  		return b + 1;
  
  	return j;
  }
  