CbmRoot
PolynomRealRoots.h
Go to the documentation of this file.
1 #include <cmath>
2 
3 //*************************************************************************
4 float polinom(int n, float x, float* k) {
5  float s = 1;
6  for (int i = n - 1; i >= 0; i--)
7  s = s * x + k[i];
8  return s;
9 } //polinom
10 
11 float dihot(int degree, float edgeNegativ, float edgePositiv, float* kf) {
12  for (;;) {
13  float x = 0.5 * (edgeNegativ + edgePositiv);
14  if (std::abs(x - edgeNegativ) < 1e-3 || std::abs(x - edgePositiv) < 1e-3)
15  return x;
16  if (polinom(degree, x, kf) < 0)
17  edgeNegativ = x;
18  else
19  edgePositiv = x;
20  }
21 } //dihot
22 
23 void stepUp(int level, float** A, float** B, int* currentRootsCount) {
24 
25  float major = 0;
26  for (int i = 0; i < level; i++) {
27  float s = fabs(A[level][i]);
28  if (s > major) major = s;
29  }
30  major += 1.0;
31 
32  currentRootsCount[level] = 0;
33 
34  for (int i = 0; i <= currentRootsCount[level - 1]; i++) {
35  int signLeft, signRight;
36  float edgeLeft, edgeRight;
37  float edgeNegativ, edgePositiv;
38 
39  if (i == 0)
40  edgeLeft = -major;
41  else
42  edgeLeft = B[level - 1][i - 1];
43 
44  float rb = polinom(level, edgeLeft, A[level]);
45 
46  if (rb == 0) {
47  B[level][currentRootsCount[level]] = edgeLeft;
48  currentRootsCount[level]++;
49  continue;
50  }
51 
52  if (rb > 0)
53  signLeft = 1;
54  else
55  signLeft = -1;
56 
57  if (i == currentRootsCount[level - 1])
58  edgeRight = major;
59  else
60  edgeRight = B[level - 1][i];
61 
62  rb = polinom(level, edgeRight, A[level]);
63 
64  if (rb == 0) {
65  B[level][currentRootsCount[level]] = edgeRight;
66  currentRootsCount[level]++;
67  continue;
68  }
69 
70  if (rb > 0)
71  signRight = 1;
72  else
73  signRight = -1;
74  if (signLeft == signRight) continue;
75 
76  if (signLeft < 0) {
77  edgeNegativ = edgeLeft;
78  edgePositiv = edgeRight;
79  } else {
80  edgeNegativ = edgeRight;
81  edgePositiv = edgeLeft;
82  }
83 
84  B[level][currentRootsCount[level]] =
85  dihot(level, edgeNegativ, edgePositiv, A[level]);
86  currentRootsCount[level]++;
87  }
88  return;
89 } //stepUp
90 
91 void polynomRealRoots(float* rootsArray, int n, float* kf_, int& rootsCount) {
92 
93  float* kf = new float[n + 1];
94 
95  float** A = new float*[n + 1];
96  float** B = new float*[n + 1];
97 
98  int* currentRootsCount = new int[n + 1];
99 
100  for (int i = 1; i <= n; i++) {
101  A[i] = new float[i];
102  B[i] = new float[i];
103  }
104 
105  for (int i = 0; i <= n; i++)
106  kf[i] = kf_[n - i];
107 
108  for (int i = 0; i < n; i++)
109  A[n][i] = kf[i] / kf[n];
110 
111  for (int i1 = n, i = n - 1; i > 0; i1 = i, i--) {
112  for (int j1 = i, j = i - 1; j >= 0; j1 = j, j--) {
113  A[i][j] = A[i1][j1] * j1 / i1;
114  }
115  }
116 
117  currentRootsCount[1] = 1;
118  B[1][0] = -A[1][0];
119 
120  for (int i = 2; i <= n; i++)
121  stepUp(i, A, B, currentRootsCount);
122 
123  rootsCount = currentRootsCount[n];
124  for (int i = 0; i < rootsCount; i++)
125  rootsArray[i] = B[n][i];
126 
127  for (int i = 1; i <= n; i++) {
128  delete[] A[i];
129  delete[] B[i];
130  }
131  delete[] A;
132  delete[] B;
133  delete[] currentRootsCount;
134  delete[] kf;
135 } //polynomRealRoots
136 
137 //*************************************************************************
polynomRealRoots
void polynomRealRoots(float *rootsArray, int n, float *kf_, int &rootsCount)
Definition: PolynomRealRoots.h:91
polinom
float polinom(int n, float x, float *k)
Definition: PolynomRealRoots.h:4
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
stepUp
void stepUp(int level, float **A, float **B, int *currentRootsCount)
Definition: PolynomRealRoots.h:23
dihot
float dihot(int degree, float edgeNegativ, float edgePositiv, float *kf)
Definition: PolynomRealRoots.h:11
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60