/*

  Teem: Tools to process and visualize scientific data and images             .

  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago

  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann

  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah

  This library is free software; you can redistribute it and/or

  modify it under the terms of the GNU Lesser General Public License

  (LGPL) as published by the Free Software Foundation; either

  version 2.1 of the License, or (at your option) any later version.

  The terms of redistributing and/or modifying this software also

  include exceptions to the LGPL that facilitate static linking.

  This library is distributed in the hope that it will be useful,

  but WITHOUT ANY WARRANTY; without even the implied warranty of

  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public License

  along with this library; if not, write to Free Software Foundation, Inc.,

  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

*/

#include "air.h"

/*

** by the way, the organization of functions into files is a little

** arbitrary around here

*/

/*

** from: Gavin C Cawley, "On a Fast, Compact Approximation of the

** Exponential Function" Neural Computation, 2000, 12(9), 2009-2012.

**

** which in turn is based on: N N Schraudolph, "A fast, compact approximation

** of the exponential function." Neural Computation, 1999, 11(4), 853-862.

** http://cnl.salk.edu/~schraudo/pubs/Schraudolph99.pdf

**

** some possible code which does not depend on endian

** http://www.machinedlearnings.com/2011/06/fast-approximate-logarithm-exponential.html

*/

typedef union {

  double dd;

  int nn[2];

} eco_t;

#define EXPA (1048576/0.69314718055994530942)

#define EXPC 60801

double

airFastExp(double val) {

  double ret;

  int tmpI, EXPI;

  /* HEY: COPY AND PASTE from airMyEndian */

  tmpI = 1;

  EXPI = *(AIR_CAST(char*, &tmpI));

         ? eco.dd

         /* seems that only times this happens is when the real exp()

            returns either 0 or +inf */

}

#undef EXPA

#undef EXPC

/*

** based on MiniMaxApproximation, but this has failed in its goal

** of being faster than libc's exp(), so should be considered

** a work-in-progress

*/

double

airExp(double x) {

  double num, den, ee;

  unsigned int pp;

    pp = 0;

    }

    }

    pp = 0;

    }

    }

  } else {

    /* x is probably not finite */

  }

}

/*

******** airNormalRand

**

** generates two numbers with normal distribution (mean 0, stdv 1)

** using the Box-Muller transform.

**

** on (seemingly sound) advice of

** <http://www.taygeta.com/random/gaussian.html>,

** I'm using the polar form of the Box-Muller transform, instead of the

** Cartesian one described at

** <http://mathworld.wolfram.com/Box-MullerTransformation.html>

**

** this is careful to not write into given NULL pointers

*/

void

airNormalRand(double *z1, double *z2) {

  double w, x1, x2;

  return;

void

airNormalRand_r(double *z1, double *z2, airRandMTState *state) {

  double w, x1, x2;

  return;

/*

******** airShuffle()

**

** generates a random permutation of integers [0..N-1] if perm is non-zero,

** otherwise, just fills buff with [0..N-1] in order

**

** see http://en.wikipedia.org/wiki/Knuth_shuffle

*/

void

airShuffle(unsigned int *buff, unsigned int N, int perm) {

  unsigned int i;

  }

  }

    unsigned int swp, tmp;

    }

  }

}

void

airShuffle_r(airRandMTState *state,

             unsigned int *buff, unsigned int N, int perm) {

  unsigned int i;

  /* HEY !!! COPY AND PASTE !!!! */

  }

  }

    unsigned int swp, tmp;

    }

  }

  /* HEY !!! COPY AND PASTE !!!! */

}

double

airSgnPow(double v, double p) {

          ? v

}

double

airFlippedSgnPow(double vv, double pp) {

  double sn=1.0;

  }

    sn = -1.0;

  }

}

double

airIntPow(double v, int p) {

  double sq, ret;

    sq = v;

      /* while the low bit is zero */

    }

    /* must terminate because we know p != 0, and when

       it terminates we know that the low bit is 1 */

    ret = sq;

      /* while there are any non-zero bits in p left */

      }

    }

    ret = 1.0;

  }

}

/*

******** airLog2()

**

** silly little function which returns log_2(n) if n is exactly a power of 2,

** and -1 otherwise

*/

int

airLog2(size_t _nn) {

  size_t nn;

  int ret;

  nn = _nn;

    /* 0 is not a power of 2 */

    ret = -1;

    int alog = 0;

    /* divide by 2 while its non-zero and the low bit is off */

    }

      ret = alog;

      ret = -1;

    }

  }

}

int

airSgn(double v) {

          ? 1

             ? -1

             : 0));

}

/*

******** airCbrt

**

** cbrt() isn't in C89, so any hacks to create a stand-in for cbrt()

** are done here.

*/

double

airCbrt(double v) {

#if defined(_WIN32) || defined(__STRICT_ANSI__)

  return (v < 0.0 ? -pow(-v,1.0/3.0) : pow(v,1.0/3.0));

#else

#endif

}

/*

** skewness of three numbers, scaled to fit in [-1,+1]

** -1: small, big, big

** +1: small, small, big

*/

double

airMode3_d(const double _v[3]) {

  double num, den, mean, v[3];

}

double

airMode3(double v0, double v1, double v2) {

}

double

airGaussian(double x, double mean, double stdv) {

}

/*

** The function approximations below were done by GLK in Mathematica,

** using its MiniMaxApproximation[] function.  The functional forms

** used for the Bessel functions were copied from Numerical Recipes

** (which were in turn copied from the "Handbook of Mathematical

** Functions with Formulas, Graphs, and Mathematical Tables" by

** Abramowitz and Stegun), but the coefficients here represent

** an increase in accuracy.

**

** The rational functions (crudely copy/paste from Mathematica into

** this file) upon which the approximations are based have a relative

** error of less than 10^-9, at least on the intervals on which they

** were created (in the case of the second branch of the

** approximation, the lower end of the interval was chosen as close to

** zero as possible). The relative error of the total approximation

** may be greater.

*/

double

airErfc(double x) {

  double ax, y, b;

  }

}

double

airErf(double x) {

}

/*

******** airBesselI0

**

** modified Bessel function of the first kind, order 0

*/

double

airBesselI0(double x) {

  double b, ax, y;

  }

}

/*

******** airBesselI0ExpScaled

**

** modified Bessel function of the first kind, order 0,

** scaled by exp(-abs(x)).

*/

double

airBesselI0ExpScaled(double x) {

  double b, ax, y;

  }

}

/*

******** airBesselI1

**

** modified Bessel function of the first kind, order 1

*/

double

airBesselI1(double x) {

  double b, ax, y;

  }

}

/*

******** airBesselI1ExpScaled

**

** modified Bessel function of the first kind, order 1,

** scaled by exp(-abs(x))

*/

double

airBesselI1ExpScaled(double x) {

  double b, ax, y;

  }

}

/*

******** airLogBesselI0

**

** natural logarithm of airBesselI0

*/

double

airLogBesselI0(double x) {

  double b, ax, y;

  }

}

/*

******** airLogRician

**

** natural logarithm of Rician distribution

** mes is measured value

** tru is "true" underlying value

** sig is sigma of 2-D Gaussian

*/

double

airLogRician(double mes, double tru, double sig) {

  double lb, ss;

}

double

airRician(double mes, double tru, double sig) {

}

/*

******** airBesselI1By0

**

** the quotient airBesselI1(x)/airBesselI0(x)

*/

double

airBesselI1By0(double x) {

  double q, ax, y;

  }

}

/*

******** airBesselIn

**

** modified Bessel function of the first kind, order n.

**

*/

double

airBesselIn(int nn, double xx) {

  double tax, bb, bi, bim, bip;

  int ii, an, top;

  }

  }

  bip = bb = 0.0;

  bi = 1.0;

    bip = bi;

    bi = bim;

    }

      bb = bip;

    }

  }

}

/*

******** airBesselInExpScaled

**

** modified Bessel function of the first kind, order n,

** scaled by exp(-abs(x))

**

*/

double

airBesselInExpScaled(int nn, double xx) {

  double tax, bb, bi, bim, bip, eps;

  int top, ii, an;

  }

  }

  bip = bb = 0.0;

  bi = 1.0;

  /* HEY: GLK tried to increase sqrt(40.0*an) to sqrt(100.0*an) to avoid

     jagged discontinuities in (e.g.) airBesselInExpScaled(n, 17*17); the

     problem was detected because of glitches in the highest blurring

     levels for scale-space feature detection; but that didn't quite

     work either; this will have to be debugged further! */

  eps = 1.0e-10;

    bip = bi;

    bi = bim;

      bb = bip;

  }

/*

** based on: T. Lindeberg. "Effective Scale: A Natural Unit For

** Measuring Scale-Space Lifetime" IEEE PAMI, 15:1068-1074 (1993)

**

** which includes tau(tee) as equation (29),

** here taking A'' == 0 and B'' == 1, with

** a0 and a1 as defined by eqs (22) and (23)

**

** Used MiniMaxApproximation[] from Mathematica (see

** ~gk/papers/ssp/nb/effective-scale-TauOfTee.nb) to get functions,

** but the specific functions and their domains could certainly be

** improved upon.  Also, the absence of conversions directly between

** tau and sigma is quite unfortunate: going through tee loses

** precision and takes more time.

**

** ACCURATE: can set this to 0 or 1:

** 0: use a quick-and-dirty approximation to tau(tee), which

** uses a straight line segment for small scales, and then

** has a C^1-continuous transition to the large-scale approximation eq (33)

** 1: careful approximation based on the MiniMaxApproximation[]s

*/

#define ACCURATE 1

double

airTauOfTime(double tee) {

  double tau;

    tau = 0;

#if ACCURATE

    /* mmtau0tweaked */

    /* mmtau1 */

    /* mmtau2 */

    /* mmtau3 */

#else /* quick and dirty approximation */

  } else if (tee < 1.3741310015234351) {

    tau = 0.600069568241882*tee/1.3741310015234351;

#endif

    /* lindtau = eq (33) in paper */

  }

}

double

airTimeOfTau(double tau) {

  double tee;

  /* the number of branches here is not good; needs re-working */

    tee = 0;

#if ACCURATE

    /* mmtee0tweaked */

    /* mmtee1 */

    /* mmtee2 */

    /* mmtee3 */

    /* mmtee4 */

    /* mmtee5 */

    /* mmtee6 */

#else /* quick and dirty approximation */

  } else if (tau < 0.600069568241882) {

    tee = 1.3741310015234351*tau/0.600069568241882;

#endif

    /* lindtee = lindtau^{-1} */

    double ee;

  }

}

double

airSigmaOfTau(double tau) {

}

double

airTauOfSigma(double sig) {

}

/*  http://en.wikipedia.org/wiki/Halton_sequences */

double

airVanDerCorput(unsigned int indx, unsigned int base) {

  double result=0.0, ff;

  unsigned int ii;

  ii = indx;

  }

}

void

airHalton(double *out, unsigned int indx,

          const unsigned int *base, unsigned int num) {

  unsigned int nn, bb;

    double ff, result = 0.0;

    unsigned int ii;

    ii = indx;

    }

  }

  return;

/* via "Table[Prime[n], {n, 1000}]" in Mathematica */

const unsigned int airPrimeList[AIR_PRIME_NUM] = {

  2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,

  67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137,

  139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211,

  223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283,

  293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379,

  383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461,

  463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563,

  569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643,

  647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739,

  743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829,

  839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937,

  941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021,

  1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093,

  1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,

  1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259,

  1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,

  1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433,

  1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493,

  1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579,

  1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,

  1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741,

  1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831,

  1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913,

  1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,

  2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087,

  2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,

  2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269,

  2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347,

  2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417,

  2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,

  2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621,

  2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,

  2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767,

  2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,

  2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953,

  2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041,

  3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163,

  3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251,

  3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329,

  3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,

  3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517,

  3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583,

  3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673,

  3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767,

  3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853,

  3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931,

  3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027,

  4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129,

  4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229,

  4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,

  4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421,

  4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513,

  4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603,

  4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691,

  4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793,

  4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909,

  4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987,

  4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077,

  5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171,

  5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279,

  5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393,

  5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471,

  5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557,

  5563, 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653,

  5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741,

  5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839,

  5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903,

  5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043,

  6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131,

  6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221,

  6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311,

  6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379,

  6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521,

  6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607,

  6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703,

  6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803,

  6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899,

  6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983,

  6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079,

  7103, 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,

  7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307,

  7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433,

  7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523,

  7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589,

  7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687,