00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef __PJ_MATH_H__
00022 #define __PJ_MATH_H__
00023
00029 #include <pj/string.h>
00030 #include <pj/compat/high_precision.h>
00031
00032 PJ_BEGIN_DECL
00033
00048 #define PJ_PI 3.14159265358979323846
00049 #define PJ_1_PI 0.318309886183790671538
00050
00054 #define PJ_ABS(x) ((x) > 0 ? (x) : -(x))
00055 #define PJ_MAX(x, y) ((x) > (y)? (x) : (y))
00056 #define PJ_MIN(x, y) ((x) < (y)? (x) : (y))
00057
00061 typedef struct pj_math_stat
00062 {
00063 int n;
00064 int max;
00065 int min;
00066 int last;
00067 int mean;
00068
00069
00070 #if PJ_HAS_FLOATING_POINT
00071 float fmean_;
00072 #else
00073 int mean_res_;
00074 #endif
00075 pj_highprec_t m2_;
00076 } pj_math_stat;
00077
00085 PJ_INLINE(unsigned) pj_isqrt(unsigned i)
00086 {
00087 unsigned res = 1, prev;
00088
00089
00090 prev = i >> 2;
00091 while (prev) {
00092 prev >>= 2;
00093 res <<= 1;
00094 }
00095
00096
00097 do {
00098 prev = res;
00099 res = (prev + i/prev) >> 1;
00100 } while ((prev+res)>>1 != res);
00101
00102 return res;
00103 }
00104
00110 PJ_INLINE(void) pj_math_stat_init(pj_math_stat *stat)
00111 {
00112 pj_bzero(stat, sizeof(pj_math_stat));
00113 }
00114
00121 PJ_INLINE(void) pj_math_stat_update(pj_math_stat *stat, int val)
00122 {
00123 #if PJ_HAS_FLOATING_POINT
00124 float delta;
00125 #else
00126 int delta;
00127 #endif
00128
00129 stat->last = val;
00130
00131 if (stat->n++) {
00132 if (stat->min > val)
00133 stat->min = val;
00134 if (stat->max < val)
00135 stat->max = val;
00136 } else {
00137 stat->min = stat->max = val;
00138 }
00139
00140 #if PJ_HAS_FLOATING_POINT
00141 delta = val - stat->fmean_;
00142 stat->fmean_ += delta/stat->n;
00143
00144
00145 stat->mean = (int) (stat->fmean_ + 0.5);
00146
00147 stat->m2_ += (int)(delta * (val-stat->fmean_));
00148 #else
00149 delta = val - stat->mean;
00150 stat->mean += delta/stat->n;
00151 stat->mean_res_ += delta % stat->n;
00152 if (stat->mean_res_ >= stat->n) {
00153 ++stat->mean;
00154 stat->mean_res_ -= stat->n;
00155 } else if (stat->mean_res_ <= -stat->n) {
00156 --stat->mean;
00157 stat->mean_res_ += stat->n;
00158 }
00159
00160 stat->m2_ += delta * (val-stat->mean);
00161 #endif
00162 }
00163
00171 PJ_INLINE(unsigned) pj_math_stat_get_stddev(const pj_math_stat *stat)
00172 {
00173 if (stat->n == 0) return 0;
00174 return (pj_isqrt((unsigned)(stat->m2_/stat->n)));
00175 }
00176
00186 PJ_INLINE(void) pj_math_stat_set_stddev(pj_math_stat *stat, unsigned dev)
00187 {
00188 if (stat->n == 0)
00189 stat->n = 1;
00190 stat->m2_ = dev*dev*stat->n;
00191 }
00192
00195 PJ_END_DECL
00196
00197 #endif