← Back to branch summary

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

« back to all changes in this revision

Viewing changes to app_realroots.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:
app_realroots.cpp
 
1
#include "dimension.h"
 
2
#include "printer.h"
 
3
#include "parser.h"
 
4
#include "gfanapplication.h"
 
5
#include "division.h"
 
6
#include "field_rationals.h"
 
7
#include "buchberger.h"
 
8
 
 
9
class RealRootsApplication : public GFanApplication
 
10
{
 
11
public:
 
12
PolynomialSet sturmPolynomials(Polynomial  f1)
 
13
{
 
14
    Polynomial f2=f1.derivative();
 
15
 
 
16
    PolynomialRing theRing=f1.getRing();
 
17
    PolynomialSet result(theRing);
 
18
    result.push_back(f1);
 
19
    while(!f2.isZero())
 
20
      {
 
21
        result.push_back(f2);
 
22
        PolynomialSet g(theRing);
 
23
        Polynomial temp=f2;
 
24
        g.push_back(f2);
 
25
        g.markAndScale(LexicographicTermOrder());
 
26
        f2=(f1-f1)-division(f1,g,LexicographicTermOrder());
 
27
        f1=temp;
 
28
      }
 
29
    return result;
 
30
}
 
31
 
 
32
 
 
33
 
 
34
int numberOfSignChangesAtMinusInfinity(PolynomialSet const &l)
 
35
{
 
36
        int ret=0;
 
37
        int sign=0;
 
38
        for(PolynomialSet::const_iterator i=l.begin();i!=l.end();i++)
 
39
        {
 
40
                Polynomial p=*i;
 
41
                p.mark(LexicographicTermOrder());
 
42
                int newSign=p.getMarked().c.sign()*(1-2*(p.getMarked().m.exponent[0]&1));
 
43
                if(newSign && (newSign!=sign))
 
44
                {
 
45
                        ret++;
 
46
                        sign=newSign;
 
47
                }
 
48
        }
 
49
        return ret;
 
50
}
 
51
 
 
52
 
 
53
int numberOfSignChangesAtInfinity(PolynomialSet const &l)
 
54
{
 
55
        int ret=0;
 
56
        int sign=0;
 
57
        for(PolynomialSet::const_iterator i=l.begin();i!=l.end();i++)
 
58
        {
 
59
                Polynomial p=*i;
 
60
                p.mark(LexicographicTermOrder());
 
61
                int newSign=p.getMarked().c.sign();
 
62
                if(newSign && (newSign!=sign))
 
63
                {
 
64
                        ret++;
 
65
                        sign=newSign;
 
66
                }
 
67
        }
 
68
        return ret;
 
69
}
 
70
 
 
71
 
 
72
int numberOfSignChanges(PolynomialSet const &l, FieldElement const &a)
 
73
{
 
74
        int ret=0;
 
75
        int sign=0;
 
76
        for(PolynomialSet::const_iterator i=l.begin();i!=l.end();i++)
 
77
        {
 
78
                FieldElement v=i->evaluate(a);
 
79
                int newSign=v.sign();
 
80
                if(newSign && (newSign!=sign))
 
81
                {
 
82
                        ret++;
 
83
                        sign=newSign;
 
84
                }
 
85
        }
 
86
        return ret;
 
87
}
 
88
 
 
89
/**
 
90
 * Returns a negative number less than all roots of the polynomial whose Sturm sequence is given.
 
91
 */
 
92
FieldElement lowerBoundForRoots(PolynomialSet const &l)
 
93
{
 
94
        FieldElement ret=Q.zHomomorphism(-1);
 
95
        while(numberOfSignChangesAtMinusInfinity(l)!=numberOfSignChanges(l,ret))ret*=Q.zHomomorphism(2);
 
96
        return ret-Q.zHomomorphism(1);
 
97
}
 
98
 
 
99
FieldElement upperBoundForRoots(PolynomialSet const &l)
 
100
{
 
101
        FieldElement ret=Q.zHomomorphism(1);
 
102
        while(numberOfSignChangesAtInfinity(l)!=numberOfSignChanges(l,ret))ret*=Q.zHomomorphism(2);
 
103
        return ret+Q.zHomomorphism(1);
 
104
}
 
105
 
 
106
list<FieldElement> intervals(PolynomialSet const &l)
 
107
{
 
108
        list<FieldElement> ret;
 
109
        FieldElement lo(lowerBoundForRoots(l));
 
110
 
 
111
//      ret.push_back(a);
 
112
        while(1)
 
113
        {
 
114
                FieldElement hi=upperBoundForRoots(l);
 
115
        int cLo=numberOfSignChanges(l,lo);
 
116
        int cHi=numberOfSignChanges(l,hi);
 
117
        if(cLo<=cHi)break;
 
118
        {
 
119
                while(cHi!=cLo-1)
 
120
                {
 
121
                        FieldElement med=(hi+lo)*Q.zHomomorphism(2).inverse();
 
122
                        int cMed=numberOfSignChanges(l,med);
 
123
                        if(cMed==cLo)
 
124
                        {
 
125
                                lo=med;
 
126
                                cLo=cMed;
 
127
                        }
 
128
                        else
 
129
                        {
 
130
                                hi=med;
 
131
                                cHi=cMed;
 
132
                        }
 
133
                }
 
134
                ret.push_back(lo);
 
135
                ret.push_back(hi);
 
136
        }
 
137
        lo=hi;
 
138
        }
 
139
        return ret;
 
140
}
 
141
 
 
142
list<FieldElement> narrow(PolynomialSet const &l, list<FieldElement> intervals, FieldElement epsilon)
 
143
{
 
144
        list<FieldElement> ret;
 
145
 
 
146
        for(list<FieldElement>::const_iterator i=intervals.begin();i!=intervals.end();i++)
 
147
        {
 
148
                FieldElement lo=*i;i++;
 
149
                FieldElement hi=*i;
 
150
 
 
151
                while((hi-lo-epsilon).sign()>0)
 
152
                {
 
153
                        FieldElement med=(hi+lo)*Q.zHomomorphism(2).inverse();
 
154
                        if(numberOfSignChanges(l,med)==numberOfSignChanges(l,hi))
 
155
                                hi=med;
 
156
                        else
 
157
                                lo=med;
 
158
                }
 
159
                ret.push_back(lo);
 
160
                ret.push_back(hi);
 
161
        }
 
162
        return ret;
 
163
}
 
164
 
 
165
//Returns a length that fist in between any two consequtive intervals
 
166
FieldElement smallestDistance(list<FieldElement> intervals)
 
167
{
 
168
        FieldElement ret=Q.zHomomorphism(1000);
 
169
        if(intervals.size()>=4)
 
170
        {
 
171
                list<FieldElement>::const_iterator i=intervals.begin();
 
172
                i++;
 
173
                for(;i!=intervals.end();i++)
 
174
                {
 
175
                        FieldElement a=*i;i++;
 
176
                        FieldElement diff=*i-a;
 
177
                        if((diff-ret).sign()<0)ret=diff;
 
178
                }
 
179
        }
 
180
        return ret;
 
181
}
 
182
/*int numberOfRootsBetweenMinusInfinityAndHere(FieldElement )
 
183
{
 
184
 
 
185
}*/
 
186
 
 
187
bool includeInDefaultInstallation()
 
188
  {
 
189
    return false;
 
190
  }
 
191
  const char *helpText()
 
192
  {
 
193
    return "Not working yet. Given generators for a zero-dimensional ideal this program will compute all real points on the variety.\n";
 
194
  }
 
195
  RealRootsApplication() {
 
196
    registerOptions();
 
197
  }
 
198
  const char *name()
 
199
  {
 
200
    return "_realroots";
 
201
  }
 
202
  int main()
 
203
  {
 
204
    FileParser P(Stdin);
 
205
 
 
206
    PolynomialSet g=P.parsePolynomialSetWithRing();
 
207
    PolynomialRing r=g.getRing();
 
208
    int n=r.getNumberOfVariables();
 
209
//    Polynomial f1=P.parsePolynomialWithRing();
 
210
    WeightReverseLexicographicTermOrder T(IntegerVector::allOnes(n));
 
211
    buchberger(&g,T);
 
212
    debug<<g;
 
213
    int d=krullDimension(g);
 
214
    if(0!=d)
 
215
    {
 
216
                debug<<"Input ideal is not zero-dimensional.\n";
 
217
        assert(0);
 
218
    }
 
219
 
 
220
    PolynomialRing r2(r.getField(),1);
 
221
    PolynomialSet projectionPolys(r2);
 
222
 
 
223
    for(int i=0;i<n;i++)
 
224
    {
 
225
        LexicographicTermOrder T((i+1)%n);
 
226
 
 
227
        buchberger(&g,T);
 
228
 
 
229
        debug<<g;
 
230
 
 
231
        list<int> l;
 
232
        l.push_back(i);
 
233
        PolynomialSet intersection=g.polynomialRingIntersection(r2,&l);
 
234
        assert(intersection.size()==1);
 
235
        projectionPolys.push_back(*intersection.begin());
 
236
    }
 
237
    debug<<projectionPolys;
 
238
 
 
239
    PolynomialSetList sturmPolys;
 
240
    for(PolynomialSet::const_iterator i=projectionPolys.begin();i!=projectionPolys.end();i++)
 
241
    {
 
242
                        sturmPolys.push_back(sturmPolynomials(*i));
 
243
        }
 
244
    debug.printPolynomialSetList(sturmPolys);
 
245
 
 
246
 
 
247
    FieldElement bound=Q.zHomomorphism(1);
 
248
    FieldElement distance=Q.zHomomorphism(1);
 
249
    for(PolynomialSetList::const_iterator i=sturmPolys.begin();i!=sturmPolys.end();i++)
 
250
    {
 
251
        debug<<"Lower "<< lowerBoundForRoots(*i)<<"\n";
 
252
        debug<<"Upper "<< upperBoundForRoots(*i)<<"\n";
 
253
        debug<<"Sign changes "<< numberOfSignChangesAtMinusInfinity(*i)<<"\n";
 
254
        FieldElement s=Q.zHomomorphism(-10000);
 
255
        debug<<"Sign changes "<< numberOfSignChanges(*i,s)<<"\n";
 
256
        FieldElement t=Q.zHomomorphism(100000);
 
257
        debug<<"Sign changes "<< numberOfSignChanges(*i,t)<<"\n";
 
258
        debug<<"Sign changes "<< numberOfSignChangesAtInfinity(*i)<<"\n";
 
259
 
 
260
        list<FieldElement> l=intervals(*i);
 
261
        FieldElement epsilon=(upperBoundForRoots(*i)-lowerBoundForRoots(*i))*Q.zHomomorphism(100).inverse();
 
262
        l=narrow(*i,l,epsilon);
 
263
 
 
264
 
 
265
        for(list<FieldElement>::const_iterator j=l.begin();j!=l.end();j++)
 
266
                debug<<*j<<" "<< fieldElementToFloatingPoint(*j)<<"\n";
 
267
 
 
268
        FieldElement b=upperBoundForRoots(*i)-lowerBoundForRoots(*i);
 
269
//      if(bound<b)bound=b;
 
270
//      FieldElement d=smallestDistance(l);
 
271
//      if(d<distance)distance=d;
 
272
    }
 
273
/*
 
274
    FieldElement epsilon2=distance/bound; //<------------------fix this
 
275
        l=narrow(*i,l,epsilon);
 
276
 
 
277
        FieldElement delta=epsilon/bound;//<--------------------fix this
 
278
 
 
279
        PolynomialRing r3=r.withVariablesAppended("W");
 
280
        Polynomial f=-r3.ithVariable(n);
 
281
        FieldElement multiplier=Q.zHomomorphism(1);
 
282
        for(int i=0;i<n;i++)
 
283
        {
 
284
                f+=multiplier*r3.ithVariable(i);
 
285
                multiplier*=delta;
 
286
        }
 
287
        PolynomialSet g2=g.embeddedInto(r3);
 
288
        g2.push_back(f);
 
289
*/
 
290
 
 
291
 
 
292
/*    FieldRationalFunctions k(r.getField(),"t");
 
293
    PolynomialRing r2(k,n+1);
 
294
    Polynomial generic(r2)=-r2.ithVariable(n);
 
295
    for(int i=0;i<n;i++)
 
296
        generic+=k.exponent(i)*r2.ithVariable(i);
 
297
 
 
298
    PolynomialSet g2(r2);
 
299
    g2.push_back(generic);
 
300
 
 
301
    for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
 
302
    {
 
303
        Polynomial f(r2);
 
304
        for(TermMap::const_iterator j=g.begin();j!=g.end();j++)
 
305
        {
 
306
                f+=Term(k.fromCoefficientField(),Monomial(r2,j->first.m.v));
 
307
        }
 
308
        g2.push_back(f);
 
309
    }
 
310
    pout << g2;
 
311
    LexicographicTermOrder T2;
 
312
    buchberger(g2,T2);
 
313
    Polynomial p(r2);
 
314
    for(PolynomialSet::const_iterator i=g2.begin();i!=g2.end();i++)
 
315
    {
 
316
        if(i->numberOfVariablesInUseConsecutive()==1)p=*i;
 
317
    }
 
318
 
 
319
    PolynomialRing r3(k,1);
 
320
 
 
321
*/
 
322
 
 
323
 
 
324
 
 
325
//AsciiPrinter(Stdout).printPolynomialSet(result);
 
326
 
 
327
/*      if(evaluateOption.getValue())
 
328
      {
 
329
        IntegerVector v=P.parseIntegerVector();
 
330
        for(int i=0;i<v.size();i++)
 
331
          {
 
332
            FieldElement x=Q.zHomomorphism(v[i]);
 
333
            AsciiPrinter(Stdout).printString("Evaluating in ");
 
334
            AsciiPrinter(Stdout).printFieldElement(x);
 
335
            AsciiPrinter(Stdout).printNewLine();
 
336
            for(PolynomialSet::const_iterator j=result.begin();j!=result.end();j++)
 
337
              {
 
338
                AsciiPrinter(Stdout).printFieldElement(j->evaluate(x));
 
339
                AsciiPrinter(Stdout).printNewLine();
 
340
              }
 
341
          }
 
342
      }
 
343
*/
 
344
    return 0;
 
345
  }
 
346
};
 
347
 
 
348
static RealRootsApplication theApplication;