hgabord/math.c at main · darkon42/hgabord · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
/*..........................................................................*/
/*                                                                          */
/*      ------------------------------------------------------------*/
/*      (C) 1993 Copyright New York University, All Right Reserved.         */
/*                                                                          */
/*..........................................................................*/
/****************************************************************************/
/*                                                                          */
/*  math.c           Some useful mathematical functions                     */
/*                                                                          */
/****************************************************************************/
#include "mpp.h"

//#define M_LN2      0.693147180559945309417
#define M_2_SQRTPI 1.12837916709551257390
//#define M_SQRT1_2  0.707106781186547524401

#ifndef sun4
/* log in base 2 */
double log2(double x)
{
  double y;
  y= ((double)log(x))/ M_LN2;
  return(y);
} /* this function is defined (as double) in math.h on the sun4 */
#endif


/* 2^j  (j > 0) */
int iexp2(int j)
{
  return( 1 << j);
}

/* 2^j  (j can be >= 0 or < 0 ) */
double fexp2(int j)
{
  int k;
  double s;
  s = 1.0;
  if (j >= 0) {
    return( (double)(1 << j));
  }
  else {
    for (k = j; k < 0 ; k++)
      s /= 2.0;
    return(s);
  }
}


/*
 * gauss function g(x) = exp(-(a*(x-b))**2)
  */
double gaussian(double x,double a,double b)
{
   double y;

   y=(double)exp(-(double)(a*a*(x-b)*(x-b)));
   return(y);
}
/*
 *
 * compute the L2 normalized gaussian:
 *	 g(t)=1/sqrt(2*pi*sigma) exp(-t*t/(2*sigma*sigma)
 *
 * M_2_SQRTPI is 2/sqrt(pi) defined in math.h
 * M_SQRT1_2 is 1/sqrt(2) defined in math.h
 *
 *
 */
double gaussianL2(double t,double sigma)
{
    double y, tmp;

    tmp = sqrt(sigma);
    y = exp(-t*t/(2.0*sigma*sigma))*
		(double)(M_2_SQRTPI*M_SQRT1_2)/(2.0*tmp);
/*
    return((double)(exp(-t*t/(2.0*sigma*sigma))*
		(double)(M_2_SQRTPI*M_SQRT1_2)/(2.0*sqrt(sigma))));
*/
    return((double)y);
}
/*
 * singauss function g(x) = [exp(-(a*(x+b))**2) -
                             exp(-(a*(x - b))**2)] / 2.0
 */
double singauss(double x,double a,double b)
{
double y;

	y=(double)exp(-(double)(a*a*(x+b)*(x+b)));
	y -= (double)exp(-(double)(a*a*(x-b)*(x-b)));
	y /= 2.0;
return(y);
}


int find2power(int n)
{
   long m, m2;

   m = 0;
   m2 = 1<<m; /* 2 to the power of m */
   while (m2-n < 0) {
	m++;
	m2 <<= 1; /* m2 = m2*2 */
   }
   return(m);
}
#ifdef hp
int rint(r)
double r;
{
    return((int)r);
}
#endif
#ifndef sun4
int nint(double r)
{
    return((int)r);
}
#endif
#ifdef i386sco
int random()
{
    return((int)Irand48());
}
#endif
/*
 * end of math.c
 */