HepLib
Rational.cpp
Go to the documentation of this file.
1 // Rational Polynomial ReConstruction
2 #include "BASIC.h"
3 #include <vector>
4 
5 namespace HepLib {
6 
7  // FROM FIRE6.m
8  numeric RationalReconstruct(numeric a, numeric p) {
9  numeric g = a, s = 1, t = 0, g1 = p, s1 = 0, t1 = 1;
10  while(g*g>p) {
11  numeric o[6];
12  int i = 0;
13  o[i++] = g1;
14  o[i++] = s1;
15  o[i++] = t1;
16  o[i++] = g - iquo(g,g1)*g1;
17  o[i++] = s - iquo(g,g1)*s1;
18  o[i++] = t - iquo(g,g1)*t1;
19  i = 0;
20  g = o[i++];
21  s = o[i++];
22  t = o[i++];
23  g1 = o[i++];
24  s1 = o[i++];
25  t1 = o[i++];
26  }
27  return g/s;
28  }
29 
30  // FROM FIRE6.m
31  numeric mulInv(numeric a, numeric b) {
32  numeric b0 = b;
33  numeric x0 = 0;
34  numeric x1 = 1;
35 
36  if (b == 1) return 1;
37 
38  while (a > 1) {
39  numeric q = iquo(a, b);
40  numeric amb = mod(a, b);
41  a = b;
42  b = amb;
43 
44  numeric xqx = x1 - q * x0;
45  x1 = x0;
46  x0 = xqx;
47  }
48  if (x1 < 0) x1 += b0;
49 
50  return x1;
51  }
52 
53  numeric ChineseRemainder(std::vector<numeric> a, std::vector<numeric> n) {
54  numeric prod = 1;
55  for (int i = 0; i < n.size(); i++) prod *= n[i];
56 
57  numeric sm = 0;
58  for (int i = 0; i < n.size(); i++) {
59  numeric p = iquo(prod, n[i]);
60  sm += a[i] * mulInv(p, n[i]) * p;
61  }
62 
63  return mod(sm, prod);
64  }
65 
66  numeric RationalReconstruct(vector<numeric> aa, vector<numeric> pp) {
67  numeric prod = 1, res = 0, resQ = 1, i, temp, steps = 0, resQprev = 0;
68  for(int i=0; i<aa.size(); i++) {
69  res = ChineseRemainder({res, aa[i]}, {prod, pp[i]});
70  prod = prod*pp[i];
71  resQ = RationalReconstruct(res, prod);
72  if(resQ == resQprev) return resQ;
73  resQprev = resQ;
74  steps++;
75  }
76  throw Error("RationalReconstruct unstable!");
77  return resQ;
78  }
79 
80  // FROM FIRE6.m
81  ex Thiele(const exvector & keys, const exvector & values, const ex & d) {
82  int n = keys.size();
83  exvector coeffs(n);
84  for(int i=0; i<n; i++) coeffs[i] = 0;
85  for(int i=0; i<n; i++) {
86  auto temp = values[i];
87  for(int j=0; j<i; j++) {
88  if(temp.is_equal(coeffs[j])) {
89  n = i-1;
90  goto out;
91  }
92  temp = (keys[i] - keys[j])/(temp - coeffs[j]);
93  }
94  coeffs[i] = temp;
95  }
96  throw Error("Thiele unstable!");
97  out: ;
98  ex result = coeffs[n];
99  for(int i=n-1; i>=0; i--) result = coeffs[i] + (d - keys[i])/result;
100  result = normal(result);
101  //result = exnormal(result);
102  return result;
103  }
104 
105  ex Newton(const exvector & keys, const exvector & values, const ex & d, const ex factor) {
106  int n = keys.size();
107  exvector coeffs(n);
108  for(int i=0; i<n; i++) coeffs[i] = 0;
109  for(int i=0; i<n; i++) {
110  auto temp = values[i] * factor.subs(d==keys[i]);
111  temp = normal(temp);
112  //temp = exnormal(temp);
113  for(int j=0; j<i; j++) {
114  if(temp.is_equal(coeffs[j])) {
115  n = i-1;
116  goto out;
117  }
118  temp = (temp - coeffs[j])/(keys[i] - keys[j]);
119  temp = normal(temp);
120  //temp = exnormal(temp);
121  }
122  coeffs[i] = temp;
123  }
124  throw Error("Newton unstable!");
125  out: ;
126  ex result = coeffs[n];
127  for(int i=n-1; i>=0; i--) result = coeffs[i] + (d - keys[i]) * result;
128  result = normal(result/factor);
129  //result = exnormal(result/factor);
130  return result;
131  }
132 
133 }
int * a
Definition: Functions.cpp:234
Basic header file.
class used to wrap error message
Definition: BASIC.h:242
HepLib namespace.
Definition: BASIC.cpp:17
numeric RationalReconstruct(numeric a, numeric p)
Definition: Rational.cpp:8
ex Thiele(const exvector &keys, const exvector &values, const ex &d)
Definition: Rational.cpp:81
numeric ChineseRemainder(std::vector< numeric > a, std::vector< numeric > n)
Definition: Rational.cpp:53
const Symbol d
numeric mulInv(numeric a, numeric b)
Definition: Rational.cpp:31
ex Newton(const exvector &keys, const exvector &values, const ex &d, const ex factor=1)
Definition: Rational.cpp:105