← Back to branch summary

~ubuntu-branches/ubuntu/saucy/gfan/saucy-proposed

« back to all changes in this revision

Viewing changes to nbody.cpp

  • Committer: Package Import Robot
  • Author(s): Cédric Boutillier
  • Date: 2013-07-09 10:44:01 UTC
  • mfrom: (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130709104401-5q66ozz5j5af0dak
Tags: 0.5+dfsg-3
* Upload to unstable.
* modify remove_failing_tests_on_32bits.patch to replace command of
  0009RenderStairCase test with an empty one instead of deleting it.
* remove lintian override about spelling error

Show diffs side-by-side

added added

removed removed

Lines of Context:
nbody.cpp
 
1
#include "polynomial.h"
 
2
#include "field_rationals.h"
 
3
#include "printer.h"
 
4
#include "log.h"
 
5
 
 
6
#include <sstream>
 
7
#include <iostream>
 
8
 
 
9
static int rIndex(int i, int j, int N, bool withMasses)
 
10
{
 
11
  if(j<i){i^=j;j^=i;i^=j;}
 
12
  int r=withMasses*N;
 
13
  for(int I=0;I<N;I++)
 
14
    for(int J=I+1;J<N;J++)
 
15
      {
 
16
        if((I==i)&&(J==j))return r;
 
17
        r++;
 
18
      }
 
19
  assert(0);
 
20
}
 
21
 
 
22
static int sIndex(int i, int j, int N, bool withMasses)
 
23
{
 
24
  return rIndex(i,j,N,withMasses)+N*(N-1)/2;
 
25
}
 
26
 
 
27
static Polynomial S(PolynomialRing const &r, int N, bool withMasses, int i, int j, bool withSVariables)
 
28
{
 
29
  Polynomial ret(r);
 
30
  int n=r.getNumberOfVariables();
 
31
 
 
32
  if(withSVariables)
 
33
    {
 
34
      ret+=Term(r.getField().zHomomorphism(1),Monomial(r,IntegerVector::standardVector(n,sIndex(i,j,N,withMasses))));
 
35
    }
 
36
  else
 
37
    {
 
38
      ret+=Term(r.getField().zHomomorphism(-1),Monomial(r,IntegerVector(n)));
 
39
      ret+=Term(r.getField().zHomomorphism(1),Monomial(r,-3*IntegerVector::standardVector(n,rIndex(i,j,N,withMasses))));
 
40
    }
 
41
 
 
42
  return ret;
 
43
}
 
44
 
 
45
 
 
46
static Polynomial rPolynomial(PolynomialRing const &r, int N, bool withMasses, int i, int j, int exponent=2)
 
47
{
 
48
  if(i==j)return Polynomial(r);
 
49
  return Term(r.getField().zHomomorphism(1),Monomial(r,exponent*IntegerVector::standardVector(r.getNumberOfVariables(),rIndex(i,j,N,withMasses))));
 
50
}
 
51
 
 
52
 
 
53
static vector<string> variableNames(int N, bool withMasses, bool withSVariables)
 
54
{
 
55
  vector<string> ret;
 
56
 
 
57
  if(withMasses)
 
58
    {
 
59
      for(int i=0;i<N;i++)
 
60
        {
 
61
          stringstream s;
 
62
          s<<"m"<<i+1;
 
63
          ret.push_back(s.str());
 
64
        }
 
65
    }
 
66
 
 
67
  for(int i=0;i<N;i++)
 
68
    for(int j=i+1;j<N;j++)
 
69
      {
 
70
        stringstream s;
 
71
        s<<"r"<<i+1<<j+1;
 
72
        ret.push_back(s.str());
 
73
      }
 
74
 
 
75
  if(withSVariables)
 
76
    {
 
77
      for(int i=0;i<N;i++)
 
78
        for(int j=i+1;j<N;j++)
 
79
          {
 
80
            stringstream s;
 
81
            s<<"S"<<i+1<<j+1;
 
82
            ret.push_back(s.str());
 
83
          }
 
84
    }
 
85
 
 
86
  return ret;
 
87
}
 
88
 
 
89
 
 
90
Polynomial AlbouyChencinerEquation(PolynomialRing const &r, int N, bool withMasses, int i, int j, bool symmetric, bool withSVariables, bool saturate)
 
91
{
 
92
  int n=r.getNumberOfVariables();
 
93
  Polynomial ret(r);
 
94
  for(int k=0;k<N;k++)
 
95
    {
 
96
      Term massTerm=Term(r.getField().zHomomorphism(1),Monomial(r,withMasses*IntegerVector::standardVector(n,k)));
 
97
      if(i!=k)
 
98
        ret+=massTerm*
 
99
          (
 
100
           S(r,N,withMasses,i,k,withSVariables)*(
 
101
                       rPolynomial(r,N,withMasses,j,k)
 
102
                       -rPolynomial(r,N,withMasses,i,k)
 
103
                       -rPolynomial(r,N,withMasses,i,j)
 
104
                       )
 
105
           );
 
106
      if(symmetric)
 
107
        if(j!=k)
 
108
          ret+=massTerm*
 
109
            (
 
110
             S(r,N,withMasses,j,k,withSVariables)*(
 
111
                         rPolynomial(r,N,withMasses,i,k)
 
112
                         -rPolynomial(r,N,withMasses,j,k)
 
113
                         -rPolynomial(r,N,withMasses,i,j)
 
114
                         )
 
115
             );
 
116
    }
 
117
  if(saturate)ret.saturate();
 
118
  return ret;
 
119
}
 
120
 
 
121
 
 
122
PolynomialSet AlbouyChencinerEquations(int N, bool withMasses, bool symmetric, bool withSVariables, bool saturate)
 
123
{
 
124
  PolynomialRing r(Q,variableNames(N,withMasses,withSVariables));
 
125
 
 
126
  PolynomialSet ret(r);
 
127
  for(int i=0;i<N;i++)
 
128
    {
 
129
      if(!symmetric)
 
130
        for(int j=0;j<i;j++)
 
131
          ret.push_back(AlbouyChencinerEquation(r,N,withMasses,i,j,symmetric,withSVariables,saturate));
 
132
 
 
133
      for(int j=i+1;j<N;j++)
 
134
        ret.push_back(AlbouyChencinerEquation(r,N,withMasses,i,j,symmetric,withSVariables,saturate));
 
135
    }
 
136
  return ret;
 
137
}
 
138
 
 
139
 
 
140
PolynomialSet DziobekEquations(PolynomialRing const &r, int N, bool withMasses, bool withSVariables, bool saturate)
 
141
{
 
142
  PolynomialSet ret(r);
 
143
 
 
144
  list<int> a;
 
145
  for(int i=4;i<N;i++)
 
146
    a.push_back(0);
 
147
  a.push_back(1);
 
148
  a.push_back(1);
 
149
  a.push_back(1);
 
150
  a.push_back(1);
 
151
 
 
152
  do
 
153
    {
 
154
      list<int>::const_iterator i=a.begin();
 
155
      int I=0;
 
156
      while((*i)==0){i++;I++;}
 
157
      int first=I;
 
158
      do{i++;I++;}while((*i)==0);
 
159
      int second=I;
 
160
      do{i++;I++;}while((*i)==0);
 
161
      int third=I;
 
162
      do{i++;I++;}while((*i)==0);
 
163
      int fourth=I;
 
164
 
 
165
      //      cerr<<first<<second<<third<<fourth;
 
166
      Polynomial f1=S(r,N,withMasses,first,second,withSVariables)*S(r,N,withMasses,third,fourth,withSVariables)-S(r,N,withMasses,first,third,withSVariables)*S(r,N,withMasses,second,fourth,withSVariables);
 
167
      if(saturate)f1.saturate();
 
168
      Polynomial f2=S(r,N,withMasses,first,third,withSVariables)*S(r,N,withMasses,second,fourth,withSVariables)-S(r,N,withMasses,first,fourth,withSVariables)*S(r,N,withMasses,second,third,withSVariables);
 
169
      if(saturate)f2.saturate();
 
170
      Polynomial f3=S(r,N,withMasses,first,second,withSVariables)*S(r,N,withMasses,third,fourth,withSVariables)-S(r,N,withMasses,first,fourth,withSVariables)*S(r,N,withMasses,second,third,withSVariables);
 
171
      if(saturate)f3.saturate();
 
172
      ret.push_back(f1);
 
173
      ret.push_back(f2);
 
174
      ret.push_back(f3);
 
175
    }
 
176
  while(next_permutation(a.begin(),a.end()));
 
177
 
 
178
  return ret;
 
179
}
 
180
 
 
181
 
 
182
static Polynomial mlookup(PolynomialRing const &r, int N, bool withMasses,int i,int j)
 
183
{
 
184
  if(i==j)return r.zero();
 
185
  if(i==0 || j==0)return r.one();
 
186
  return Polynomial(Term(r.getField().zHomomorphism(1),Monomial(r,2*IntegerVector::standardVector(r.getNumberOfVariables(),rIndex(i-1,j-1,N,withMasses)))));
 
187
}
 
188
 
 
189
 
 
190
PolynomialSet nbodyDeterminants(PolynomialRing const &r, int N, bool withMasses, int determinantSize)
 
191
{
 
192
  vector<int> l;
 
193
 
 
194
  for(int i=0;i<N+1-determinantSize;i++)
 
195
    l.push_back(1);
 
196
  for(int i=0;i<determinantSize-1;i++)
 
197
    l.push_back(0);
 
198
 
 
199
  PolynomialSet ret(r);
 
200
 
 
201
  if(determinantSize==N+2)return ret;
 
202
 
 
203
  do
 
204
    {
 
205
      vector<int> indexList;
 
206
      indexList.push_back(0);
 
207
      for(int i=0;i<l.size();i++)if(l[i]==0)indexList.push_back(i+1);
 
208
 
 
209
      log1
 
210
      {
 
211
          for(int A=0;A<determinantSize;A++)
 
212
          {
 
213
                  for(int B=0;B<determinantSize;B++)
 
214
                          AsciiPrinter(Stderr)<<mlookup(r,N,withMasses,indexList[A],indexList[B])<<";";
 
215
                                  cerr<<endl;
 
216
          }
 
217
      }
 
218
 
 
219
      vector<int> perm;
 
220
      for(int i=0;i<indexList.size();i++)
 
221
        perm.push_back(i);
 
222
 
 
223
      Polynomial p(r);
 
224
 
 
225
      do
 
226
        {
 
227
          Polynomial prod=r.one();
 
228
          for(int j=0;j<perm.size();j++)
 
229
            {
 
230
              prod*=mlookup(r,N,withMasses,indexList[j],indexList[perm[j]]);
 
231
            }
 
232
          int s=1;
 
233
          for(int x=0;x<perm.size();x++)
 
234
            for(int y=0;y<x;y++)
 
235
              if(perm[y]>perm[x])s*=-1;
 
236
          if(s==1)
 
237
            p+=prod;
 
238
          else
 
239
            p-=prod;
 
240
        }
 
241
      while(next_permutation(perm.begin(),perm.end()));
 
242
      ret.push_back(p);
 
243
    }
 
244
  while(prev_permutation(l.begin(),l.end()));
 
245
 
 
246
  return ret;
 
247
}
 
248
 
 
249
 
 
250
Polynomial massEquation(PolynomialRing const &r, int N, bool withMasses, bool saturate)
 
251
{
 
252
  Polynomial ret(r);
 
253
  for(int i=0;i<N;i++)
 
254
    for(int j=0;j<i;j++)
 
255
      {
 
256
        Polynomial mm=r.one();
 
257
        if(withMasses)mm=Polynomial(Term(r.getField().zHomomorphism(1),Monomial(r,IntegerVector::standardVector(r.getNumberOfVariables(),i)+IntegerVector::standardVector(r.getNumberOfVariables(),j))));
 
258
        ret+=(rPolynomial(r,N,withMasses,i,j,2)-rPolynomial(r,N,withMasses,i,j,-1))*mm;
 
259
      }
 
260
  if(saturate)ret.saturate();
 
261
  return ret;
 
262
}
 
263
 
 
264
 
 
265
Polynomial SEquation(PolynomialRing const &r, int N, bool withMasses, int i, int j, bool withSVariables, bool saturate=true)
 
266
{
 
267
  Polynomial ret(r);
 
268
  int n=r.getNumberOfVariables();
 
269
 
 
270
  ret-=Term(r.getField().zHomomorphism(1),Monomial(r,IntegerVector::standardVector(n,sIndex(i,j,N,withMasses))));
 
271
 
 
272
  ret+=Term(r.getField().zHomomorphism(-1),Monomial(r,IntegerVector(n)));
 
273
  ret+=Term(r.getField().zHomomorphism(1),Monomial(r,-3*IntegerVector::standardVector(n,rIndex(i,j,N,withMasses))));
 
274
 
 
275
  if(saturate)ret.saturate();
 
276
 
 
277
  return ret;
 
278
}
 
279
 
 
280
 
 
281
PolynomialSet SEquations(PolynomialRing const &r, int N, bool withMasses, bool saturate)
 
282
{
 
283
  PolynomialSet ret(r);
 
284
  for(int i=0;i<N;i++)
 
285
    for(int j=0;j<i;j++)
 
286
      {
 
287
        ret.push_back(SEquation(r,N,withMasses,i,j,true,saturate));
 
288
      }
 
289
  return ret;
 
290
}