#ifndef gauss_table_H
#define gauss_table_H

/* Tools for Gaussian kernel tables */
/* Last edited on 2008-10-05 11:10:32 by stolfi */

#include <bool.h>

#define gauss_table_BIG_ARG (8.6)
  /* Effective support radius of a Gaussian with unit deviation in {double}
    arithmentic. Namely, the value of {z} such that {exp(-z^2/2)}
    is less than {1.0e-16}. */

#define gauss_table_TINY_DEV (1.0/gauss_table_BIG_ARG)
  /* Standard deviation of a Gaussian that is effectively zero at
    argument value 1.0. */

double *gauss_table_make(int n, double avg, double dev, bool_t normSum);
  /* Allocates a vector {w} with {n} elements and sets element {w[i]}
    to {g((i-avg)/dev)}, where {g(z)} is the Gaussian bell {exp(-z*z/2)}.
    
    The bell will be wrapped around the edges of the vector. That is,
    {w[i]} is actually set to {SUM { g((i-avg+k*n)/dev) : k in Z }}.
    The fold-over is not significant if the Gaussian is far from the
    table ends --- that is, if {BIG*dev <= avg <= n - BIG*dev} where
    {BIG = gauss_table_BIG_ARG}.
    
    In particular, if {dev} is very small (specifically, less than
    {gauss_table_TINY_DEV}), entry {w[i]} is set to 1 if {i} is
    congruent to {avg} modulo {n}, ad to 0 otherwise. If {dev} is
    {+INF}, all entries are set to 1.
    
    These raw values are normalized before the procedure returns. If
    {normSum} is TRUE, the elements {w[0..n-1]} are scaled so that
    their sum is 1.0. Otherwise they are scaled so that the weight at
    index {avg} (interpolated if {avg} is fractional, and including
    any wrap-around terms) is 1.0. */

double gauss_table_folded_bell(double z, double dev, int n);
  /* Returns the sum of the Gaussian bell {exp(-(x/dev)^2/2)} evaluated 
    at all arguments {x} that are congruent to {z} modulo {n}.
    
    As special cases, if {dev} is zero, returns 1 when {z} is
    congruent to 0 modulo {n}, and 0 otherwise. If {dev} is {+INF},
    returns 1. If {n} is zero, assumes {n = +oo}, meaning
    no fold-over. */

#endif