source: git/kernel/mpr_numeric.cc @ 599326

spielwiese
Last change on this file since 599326 was 599326, checked in by Kai Krüger <krueger@…>, 14 years ago
Anne, Kai, Frank: - changes to #include "..." statements to allow cleaner build structure - affected directories: omalloc, kernel, Singular - not yet done: IntergerProgramming git-svn-id: file:///usr/local/Singular/svn/trunk@13032 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 30.4 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5/* $Id$ */
6
7/*
8* ABSTRACT - multipolynomial resultants - numeric stuff
9*            ( root finder, vandermonde system solver, simplex )
10*/
11
12#include <kernel/mod2.h>
13//#ifdef HAVE_MPR
14
15//#define mprDEBUG_ALL
16
17//-> includes
18#include <kernel/options.h>
19#include <kernel/febase.h>
20#include <omalloc.h>
21#include <kernel/numbers.h>
22#include <kernel/polys.h>
23#include <kernel/ideals.h>
24#include <kernel/intvec.h>
25#include <kernel/matpol.h>
26#include <kernel/ring.h>
27//#include "longrat.h"
28
29#include <math.h>
30#include <kernel/mpr_numeric.h>
31
32extern size_t gmp_output_digits;
33//<-
34
35extern void nPrint(number n);  // for debugging output
36
37//-----------------------------------------------------------------------------
38//-------------- vandermonde system solver ------------------------------------
39//-----------------------------------------------------------------------------
40
41//-> vandermonde::*
42vandermonde::vandermonde( const long _cn, const long _n, const long _maxdeg,
43                          number *_p, const bool _homog )
44  : n(_n), cn(_cn), maxdeg(_maxdeg), p(_p), homog(_homog)
45{
46  long j;
47  l= (long)pow((double)maxdeg+1,(int)n);
48  x= (number *)omAlloc( cn * sizeof(number) );
49  for ( j= 0; j < cn; j++ ) x[j]= nInit(1);
50  init();
51}
52
53vandermonde::~vandermonde()
54{
55  int j;
56  for ( j= 0; j < cn; j++ ) nDelete( x+j );
57  omFreeSize( (ADDRESS)x, cn * sizeof( number ) );
58}
59
60void vandermonde::init()
61{
62  int j;
63  long i,c,sum;
64  number tmp,tmp1;
65
66  c=0;
67  sum=0;
68
69  intvec exp( n );
70  for ( j= 0; j < n; j++ ) exp[j]=0;
71
72  for ( i= 0; i < l; i++ )
73  {
74    if ( !homog || (sum == maxdeg) )
75    {
76      for ( j= 0; j < n; j++ )
77      {
78        nPower( p[j], exp[j], &tmp );
79        tmp1 = nMult( tmp, x[c] );
80        x[c]= tmp1;
81        nDelete( &tmp );
82      }
83      c++;
84    }
85    exp[0]++;
86    sum=0;
87    for ( j= 0; j < n - 1; j++ )
88    {
89      if ( exp[j] > maxdeg )
90      {
91        exp[j]= 0;
92        exp[j + 1]++;
93        }
94      sum+= exp[j];
95    }
96    sum+= exp[n - 1];
97  }
98}
99
100poly vandermonde::numvec2poly( const number * q )
101{
102  int j;
103  long i,c,sum;
104  number tmp;
105
106  poly pnew,pit=NULL;
107
108  c=0;
109  sum=0;
110
111  int *exp= (int *) omAlloc( (n+1) * sizeof(int) );
112
113  for ( j= 0; j < n+1; j++ ) exp[j]=0;
114
115  for ( i= 0; i < l; i++ )
116  {
117    if ( (!homog || (sum == maxdeg)) && q[i] && !nIsZero(q[i]) )
118    {
119      pnew= pOne();
120      pSetCoeff(pnew,q[i]);
121      pSetExpV(pnew,exp);
122      if ( pit )
123      {
124        pNext(pnew)= pit;
125        pit= pnew;
126      }
127      else
128      {
129        pit= pnew;
130        pNext(pnew)= NULL;
131      }
132      pSetm(pit);
133    }
134    exp[1]++;
135    sum=0;
136    for ( j= 1; j < n; j++ )
137    {
138      if ( exp[j] > maxdeg )
139      {
140        exp[j]= 0;
141        exp[j + 1]++;
142        }
143      sum+= exp[j];
144    }
145    sum+= exp[n];
146  }
147
148  omFreeSize( (ADDRESS) exp, (n+1) * sizeof(int) );
149
150  pSortAdd(pit);
151  return pit;
152}
153
154number * vandermonde::interpolateDense( const number * q )
155{
156  int i,j,k;
157  number newnum,tmp1;
158  number b,t,xx,s;
159  number *c;
160  number *w;
161
162  b=t=xx=s=tmp1=NULL;
163
164  w= (number *)omAlloc( cn * sizeof(number) );
165  c= (number *)omAlloc( cn * sizeof(number) );
166  for ( j= 0; j < cn; j++ )
167  {
168    w[j]= nInit(0);
169    c[j]= nInit(0);
170  }
171
172  if ( cn == 1 )
173  {
174    nDelete( &w[0] );
175    w[0]= nCopy(q[0]);
176  }
177  else
178  {
179    nDelete( &c[cn-1] );
180    c[cn-1]= nCopy(x[0]);
181    c[cn-1]= nNeg(c[cn-1]);              // c[cn]= -x[1]
182
183    for ( i= 1; i < cn; i++ ) {              // i=2; i <= cn
184      nDelete( &xx );
185      xx= nCopy(x[i]);
186      xx= nNeg(xx);               // xx= -x[i]
187
188      for ( j= (cn-i-1); j <= (cn-2); j++) { // j=(cn+1-i); j <= (cn-1)
189        nDelete( &tmp1 );
190        tmp1= nMult( xx, c[j+1] );           // c[j]= c[j] + (xx * c[j+1])
191        newnum= nAdd( c[j], tmp1 );
192        nDelete( &c[j] );
193        c[j]= newnum;
194      }
195
196      newnum= nAdd( xx, c[cn-1] );           // c[cn-1]= c[cn-1] + xx
197      nDelete( &c[cn-1] );
198      c[cn-1]= newnum;
199    }
200
201    for ( i= 0; i < cn; i++ ) {              // i=1; i <= cn
202      nDelete( &xx );
203      xx= nCopy(x[i]);                     // xx= x[i]
204
205      nDelete( &t );
206      t= nInit( 1 );                         // t= b= 1
207      nDelete( &b );
208      b= nInit( 1 );
209      nDelete( &s );                         // s= q[cn-1]
210      s= nCopy( q[cn-1] );
211
212      for ( k= cn-1; k >= 1; k-- ) {         // k=cn; k >= 2
213        nDelete( &tmp1 );
214        tmp1= nMult( xx, b );                // b= c[k] + (xx * b)
215        nDelete( &b );
216        b= nAdd( c[k], tmp1 );
217
218        nDelete( &tmp1 );
219        tmp1= nMult( q[k-1], b );            // s= s + (q[k-1] * b)
220        newnum= nAdd( s, tmp1 );
221        nDelete( &s );
222        s= newnum;
223
224        nDelete( &tmp1 );
225        tmp1= nMult( xx, t );                // t= (t * xx) + b
226        newnum= nAdd( tmp1, b );
227        nDelete( &t );
228        t= newnum;
229      }
230
231      if (!nIsZero(t))
232      {
233        nDelete( &w[i] );                      // w[i]= s/t
234        w[i]= nDiv( s, t );
235        nNormalize( w[i] );
236      }
237
238      mprSTICKYPROT(ST_VANDER_STEP);
239    }
240  }
241  mprSTICKYPROT("\n");
242
243  // free mem
244  for ( j= 0; j < cn; j++ ) nDelete( c+j );
245  omFreeSize( (ADDRESS)c, cn * sizeof( number ) );
246
247  nDelete( &tmp1 );
248  nDelete( &s );
249  nDelete( &t );
250  nDelete( &b );
251  nDelete( &xx );
252
253  // makes quotiens smaller
254  for ( j= 0; j < cn; j++ ) nNormalize( w[j] );
255
256  return w;
257}
258//<-
259
260//-----------------------------------------------------------------------------
261//-------------- rootContainer ------------------------------------------------
262//-----------------------------------------------------------------------------
263
264//-> definitions
265#define MR       8        // never change this value
266#define MT      5
267#define MMOD    (MT*MR)
268#define MAXIT   (5*MMOD)   // max number of iterations in laguer root finder
269
270
271//-> rootContainer::rootContainer()
272rootContainer::rootContainer()
273{
274  rt=none;
275
276  coeffs= NULL;
277  ievpoint= NULL;
278  theroots= NULL;
279
280  found_roots= false;
281}
282//<-
283
284//-> rootContainer::~rootContainer()
285rootContainer::~rootContainer()
286{
287  int i;
288  int n= pVariables;
289
290  // free coeffs, ievpoint
291  if ( ievpoint != NULL )
292  {
293    for ( i=0; i < anz+2; i++ ) nDelete( ievpoint + i );
294    omFreeSize( (ADDRESS)ievpoint, (anz+2) * sizeof( number ) );
295  }
296
297  for ( i=0; i <= tdg; i++ ) nDelete( coeffs + i );
298  omFreeSize( (ADDRESS)coeffs, (tdg+1) * sizeof( number ) );
299
300  // theroots löschen
301  for ( i=0; i < tdg; i++ ) delete theroots[i];
302  omFreeSize( (ADDRESS) theroots, (tdg)*sizeof(gmp_complex*) );
303
304  //mprPROTnl("~rootContainer()");
305}
306//<-
307
308//-> void rootContainer::fillContainer( ... )
309void rootContainer::fillContainer( number *_coeffs, number *_ievpoint,
310                                   const int _var, const int _tdg,
311                                   const rootType  _rt, const int _anz )
312{
313  int i;
314  number nn= nInit(0);
315  var=_var;
316  tdg=_tdg;
317  coeffs=_coeffs;
318  rt=_rt;
319  anz=_anz;
320
321  for ( i=0; i <= tdg; i++ )
322  {
323    if ( nEqual(coeffs[i],nn) )
324    {
325      nDelete( &coeffs[i] );
326      coeffs[i]=NULL;
327    }
328  }
329  nDelete( &nn );
330
331  if ( rt == cspecialmu && _ievpoint )  // copy ievpoint
332  {
333    ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) );
334    for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] );
335  }
336
337  theroots= NULL;
338  found_roots= false;
339}
340//<-
341
342//-> poly rootContainer::getPoly()
343poly rootContainer::getPoly()
344{
345  int i;
346
347  poly result= NULL;
348  poly ppos;
349
350  if ( (rt == cspecial) || ( rt == cspecialmu ) )
351  {
352    for ( i= tdg; i >= 0; i-- )
353    {
354      if ( coeffs[i] )
355      {
356        poly p= pOne();
357        //pSetExp( p, var+1, i);
358        pSetExp( p, 1, i);
359        pSetCoeff( p, nCopy( coeffs[i] ) );
360        pSetm( p );
361        if (result)
362        {
363          ppos->next=p;
364          ppos=ppos->next;
365        }
366        else
367        {
368          result=p;
369          ppos=p;
370        }
371
372      }
373    }
374    if (result!=NULL) pSetm( result );
375  }
376
377  return result;
378}
379//<-
380
381//-> const gmp_complex & rootContainer::opterator[] ( const int i )
382// this is now inline!
383//  gmp_complex & rootContainer::operator[] ( const int i )
384//  {
385//    if ( found_roots && ( i >= 0) && ( i < tdg ) )
386//    {
387//      return *theroots[i];
388//    }
389//    // warning
390//    Warn("rootContainer::getRoot: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
391//    gmp_complex *tmp= new gmp_complex();
392//    return *tmp;
393//  }
394//<-
395
396//-> gmp_complex & rootContainer::evPointCoord( int i )
397gmp_complex & rootContainer::evPointCoord( const int i )
398{
399  if (! ((i >= 0) && (i < anz+2) ) )
400    WarnS("rootContainer::evPointCoord: index out of range");
401  if (ievpoint == NULL)
402    WarnS("rootContainer::evPointCoord: ievpoint == NULL");
403
404  if ( (rt == cspecialmu) && found_roots ) // FIX ME
405  {
406    if ( ievpoint[i] != NULL )
407    {
408      gmp_complex *tmp= new gmp_complex();
409      *tmp= numberToComplex(ievpoint[i]);
410      return *tmp;
411    }
412    else
413    {
414      Warn("rootContainer::evPointCoord: NULL index %d",i);
415    }
416  }
417
418  // warning
419  Warn("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
420  gmp_complex *tmp= new gmp_complex();
421  return *tmp;
422}
423//<-
424
425//-> bool rootContainer::changeRoots( int from, int to )
426bool rootContainer::swapRoots( const int from, const int to )
427{
428  if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) )
429  {
430    if ( to != from )
431    {
432      gmp_complex tmp( *theroots[from] );
433      *theroots[from]= *theroots[to];
434      *theroots[to]= tmp;
435    }
436    return true;
437  }
438
439  // warning
440  Warn(" rootContainer::changeRoots: Wrong index %d, %d",from,to);
441  return false;
442}
443//<-
444
445//-> void rootContainer::solver()
446bool rootContainer::solver( const int polishmode )
447{
448  int i;
449
450  // there are maximal tdg roots, so *roots ranges form 0 to tdg-1.
451  theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) );
452  for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex();
453
454  // copy the coefficients of type number to type gmp_complex
455  gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
456  for ( i=0; i <= tdg; i++ )
457  {
458    ad[i]= new gmp_complex();
459    if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i] );
460  }
461
462  // now solve
463  found_roots= laguer_driver( ad, theroots, polishmode != 0 );
464  if (!found_roots)
465    WarnS("rootContainer::solver: No roots found!");
466 
467  // free memory
468  for ( i=0; i <= tdg; i++ ) delete ad[i];
469  omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof(gmp_complex*) );
470
471  return found_roots;
472}
473//<-
474
475//-> gmp_complex* rootContainer::laguer_driver( bool polish )
476bool rootContainer::laguer_driver(gmp_complex ** a, gmp_complex ** roots, bool polish )
477{
478  int i,j,k,its;
479  gmp_float zero(0.0);
480  gmp_complex x(0.0),o(1.0);
481  bool ret= true, isf=isfloat(a), type=true;
482
483  gmp_complex ** ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
484  for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] );
485
486  k = 0;
487  i = tdg;
488  j = i-1;
489  while (i>2)
490  {
491    // run laguer alg
492    x = zero;
493    laguer(ad, i, &x, &its, type);
494    if ( its > MAXIT )
495    {
496      type = !type;
497      x = zero;
498      laguer(ad, i, &x, &its, type);
499    }
500
501    mprSTICKYPROT(ST_ROOTS_LGSTEP);
502    if ( its > MAXIT )
503    {  // error
504      WarnS("Laguerre solver: Too many iterations!");
505      ret= false;
506      goto theend;
507    }
508    if ( polish )
509    {
510      laguer( a, tdg, &x, &its , type);
511      if ( its > MAXIT )
512      {  // error
513        WarnS("Laguerre solver: Too many iterations in polish!");
514        ret= false;
515        goto theend;
516      }
517    }
518    if ((!type)&&(!((x.real()==zero)&&(x.imag()==zero)))) x = o/x;
519    if (x.imag() == zero)
520    {
521      *roots[k] = x;
522      k++;
523      divlin(ad,x,i);
524      i--;
525    }
526    else
527    {
528      if(isf)
529      {
530        *roots[j] = x;
531        *roots[j-1]= gmp_complex(x.real(),-x.imag());
532        j -= 2;
533        divquad(ad,x,i);
534        i -= 2;
535      }
536      else
537      {
538        *roots[j] = x;
539        j--;
540        divlin(ad,x,i);
541        i--;
542      }
543    }
544    type = !type;
545  }
546  solvequad(ad,roots,k,j);
547  sortroots(roots,k,j,isf);
548
549theend:
550  mprSTICKYPROT("\n");
551  for ( i=0; i <= tdg; i++ ) delete ad[i];
552  omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof( gmp_complex* ));
553
554  return ret;
555}
556//<-
557
558//-> void rootContainer::laguer(...)
559void rootContainer::laguer(gmp_complex ** a, int m, gmp_complex *x, int *its, bool type)
560{
561  int iter,j;
562  gmp_float zero(0.0),one(1.0),deg(m);
563  gmp_float abx_g, err_g, fabs;
564  gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2;
565  gmp_float frac_g[MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 };
566
567  gmp_float epss(0.1);
568  mpf_pow_ui(*epss._mpfp(),*epss.mpfp(),getGMPFloatDigits());
569
570  for ( iter= 1; iter <= MAXIT; iter++ )
571  {
572    mprSTICKYPROT(ST_ROOTS_LG);
573    *its=iter;
574    if (type)
575      computefx(a,*x,m,b,d,f,abx_g,err_g);
576    else
577      computegx(a,*x,m,b,d,f,abx_g,err_g);
578    err_g *= epss; // EPSS;
579
580    fabs = abs(b);
581    if (fabs <= err_g)
582    {
583      if ((fabs==zero) || (abs(d)==zero)) return;
584      *x -= (b/d); // a last newton-step
585      goto ende;
586    }
587
588    g= d / b;
589    g2 = g * g;
590    h= g2 - (((f+f) / b ));
591    sq= sqrt(( ( h * deg ) - g2 ) * (deg - one));
592    gp= g + sq;
593    gm= g - sq;
594    if (abs(gp)<abs(gm))
595    {
596      dx = deg/gm;
597    }
598    else
599    {
600      if((gp.real()==zero)&&(gp.imag()==zero))
601      {
602        dx.real(cos((mprfloat)iter));
603        dx.imag(sin((mprfloat)iter));
604        dx = dx*(one+abx_g);
605      }
606      else
607      {
608        dx = deg/gp;
609      }
610    }
611    x1= *x - dx;
612
613    if (*x == x1) goto ende;
614
615    j = iter%MMOD;
616    if (j==0) j=MT;
617    if ( j % MT ) *x= x1;
618    else *x -= ( dx * frac_g[ j / MT ] );
619  }
620
621  *its= MAXIT+1;
622ende:
623  checkimag(x,epss);
624}
625
626void rootContainer::checkimag(gmp_complex *x, gmp_float &e)
627{
628  if(abs(x->imag())<abs(x->real())*e)
629  {
630    x->imag(0.0);
631  }
632}
633
634bool rootContainer::isfloat(gmp_complex **a)
635{
636  gmp_float z(0.0);
637  gmp_complex *b;
638  for (int i=tdg; i >= 0; i-- )
639  {
640    b = &(*a[i]);
641    if (!(b->imag()==z))
642      return false;
643  }
644  return true;
645}
646
647void rootContainer::divlin(gmp_complex **a, gmp_complex x, int j)
648{
649  int i;
650  gmp_float o(1.0);
651
652  if (abs(x)<o)
653  {
654    for (i= j-1; i > 0; i-- )
655      *a[i] += (*a[i+1]*x);
656    for (i= 0; i < j; i++ )
657      *a[i] = *a[i+1];
658  }
659  else
660  {
661    gmp_complex y(o/x);
662    for (i= 1; i < j; i++)
663      *a[i] += (*a[i-1]*y);
664  }
665}
666
667void rootContainer::divquad(gmp_complex **a, gmp_complex x, int j)
668{
669  int i;
670  gmp_float o(1.0),p(x.real()+x.real()),
671            q((x.real()*x.real())+(x.imag()*x.imag()));
672
673  if (abs(x)<o)
674  {
675    *a[j-1] += (*a[j]*p);
676    for (i= j-2; i > 1; i-- )
677      *a[i] += ((*a[i+1]*p)-(*a[i+2]*q));
678    for (i= 0; i < j-1; i++ )
679      *a[i] = *a[i+2];
680  }
681  else
682  {
683    p = p/q;
684    q = o/q;
685    *a[1] += (*a[0]*p);
686    for (i= 2; i < j-1; i++)
687      *a[i] += ((*a[i-1]*p)-(*a[i-2]*q));
688  }
689}
690
691void rootContainer::solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
692{
693  gmp_float zero(0.0);
694
695  if ((j>k)
696  &&((!(*a[2]).real().isZero())||(!(*a[2]).imag().isZero())))
697  {
698    gmp_complex sq(zero);
699    gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]);
700    gmp_complex disk((h1 * h1) - h2);
701    if (disk.imag().isZero())
702    {
703      if (disk.real()<zero)
704      {
705        sq.real(zero);
706        sq.imag(sqrt(-disk.real()));
707      }
708      else
709        sq = (gmp_complex)sqrt(disk.real());
710    }
711    else
712      sq = sqrt(disk);
713    *r[k+1] = sq - h1;
714    sq += h1;
715    *r[k] = (gmp_complex)0.0-sq;
716    if(sq.imag().isZero())
717    {
718      k = j;
719      j++;
720    }
721    else
722    {
723      j = k;
724      k--;
725    }
726  }
727  else
728  {
729    if (((*a[1]).real().isZero()) && ((*a[1]).imag().isZero()))
730    {
731      WerrorS("precision lost, try again with higher precision");
732    }
733    else
734    {
735      *r[k]= (gmp_complex)0.0-(*a[0] / *a[1]);
736      if(r[k]->imag().isZero())
737        j++;
738      else
739        k--;
740    }
741  }
742}
743
744void rootContainer::sortroots(gmp_complex **ro, int r, int c, bool isf)
745{
746  int j;
747
748  for (j=0; j<r; j++) // the real roots
749    sortre(ro, j, r, 1);
750  if (c>=tdg) return;
751  if (isf)
752  {
753    for (j=c; j+2<tdg; j+=2) // the complex roots for a real poly
754      sortre(ro, j, tdg-1, 2);
755  }
756  else
757  {
758    for (j=c; j+1<tdg; j++) // the complex roots for a general poly
759      sortre(ro, j, tdg-1, 1);
760  }
761}
762
763void rootContainer::sortre(gmp_complex **r, int l, int u, int inc)
764{
765  int pos,i;
766  gmp_complex *x,*y;
767
768  pos = l;
769  x = r[pos];
770  for (i=l+inc; i<=u; i+=inc)
771  {
772    if (r[i]->real()<x->real())
773    {
774      pos = i;
775      x = r[pos];
776    }
777  }
778  if (pos>l)
779  {
780    if (inc==1)
781    {
782      for (i=pos; i>l; i--)
783        r[i] = r[i-1];
784      r[l] = x;
785    }
786    else
787    {
788      y = r[pos+1];
789      for (i=pos+1; i+1>l; i--)
790        r[i] = r[i-2];
791      if (x->imag()>y->imag())
792      {
793        r[l] = x;
794        r[l+1] = y;
795      }
796      else
797      {
798        r[l] = y;
799        r[l+1] = x;
800      }
801    }
802  }
803  else if ((inc==2)&&(x->imag()<r[l+1]->imag()))
804  {
805    r[l] = r[l+1];
806    r[l+1] = x;
807  }
808}
809
810void rootContainer::computefx(gmp_complex **a, gmp_complex x, int m,
811                  gmp_complex &f0, gmp_complex &f1, gmp_complex &f2,
812                  gmp_float &ex, gmp_float &ef)
813{
814  int k;
815
816  f0= *a[m];
817  ef= abs(f0);
818  f1= gmp_complex(0.0);
819  f2= f1;
820  ex= abs(x);
821
822  for ( k= m-1; k >= 0; k-- )
823  {
824    f2 = ( x * f2 ) + f1;
825    f1 = ( x * f1 ) + f0;
826    f0 = ( x * f0 ) + *a[k];
827    ef = abs( f0 ) + ( ex * ef );
828  }
829}
830
831void rootContainer::computegx(gmp_complex **a, gmp_complex x, int m,
832                  gmp_complex &f0, gmp_complex &f1, gmp_complex &f2,
833                  gmp_float &ex, gmp_float &ef)
834{
835  int k;
836
837  f0= *a[0];
838  ef= abs(f0);
839  f1= gmp_complex(0.0);
840  f2= f1;
841  ex= abs(x);
842
843  for ( k= 1; k <= m; k++ )
844  {
845    f2 = ( x * f2 ) + f1;
846    f1 = ( x * f1 ) + f0;
847    f0 = ( x * f0 ) + *a[k];
848    ef = abs( f0 ) + ( ex * ef );
849  }
850}
851
852//-----------------------------------------------------------------------------
853//-------------- rootArranger -------------------------------------------------
854//-----------------------------------------------------------------------------
855
856//-> rootArranger::rootArranger(...)
857rootArranger::rootArranger( rootContainer ** _roots,
858                            rootContainer ** _mu,
859                            const int _howclean )
860  : roots(_roots), mu(_mu), howclean(_howclean)
861{
862  found_roots=false;
863}
864//<-
865
866//-> void rootArranger::solve_all()
867void rootArranger::solve_all()
868{
869  int i;
870  found_roots= true;
871
872  // find roots of polys given by coeffs in roots
873  rc= roots[0]->getAnzElems();
874  for ( i= 0; i < rc; i++ )
875    if ( !roots[i]->solver( howclean ) )
876    {
877      found_roots= false;
878      return;
879    }
880  // find roots of polys given by coeffs in mu
881  mc= mu[0]->getAnzElems();
882  for ( i= 0; i < mc; i++ )
883    if ( ! mu[i]->solver( howclean ) )
884    {
885      found_roots= false;
886      return;
887    }
888}
889//<-
890
891//-> void rootArranger::arrange()
892void rootArranger::arrange()
893{
894  gmp_complex tmp,zwerg;
895  int anzm= mu[0]->getAnzElems();
896  int anzr= roots[0]->getAnzRoots();
897  int xkoord, r, rtest, xk, mtest;
898  bool found;
899  //gmp_complex mprec(1.0/pow(10,gmp_output_digits-5),1.0/pow(10,gmp_output_digits-5));
900
901  for ( xkoord= 0; xkoord < anzm; xkoord++ ) {    // für x1,x2, x1,x2,x3, x1,x2,...,xn
902    gmp_float mprec(1.0/pow(10.0,(int)(gmp_output_digits/3)));
903    for ( r= 0; r < anzr; r++ ) {                 // für jede Nullstelle
904      // (x1-koordinate) * evp[1] + (x2-koordinate) * evp[2] +
905      //                                  ... + (xkoord-koordinate) * evp[xkoord]
906      tmp= gmp_complex();
907      for ( xk =0; xk <= xkoord; xk++ )
908      {
909        tmp -= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); //xk+1
910      }
911      found= false;
912      do { // while not found
913        for ( rtest= r; rtest < anzr; rtest++ ) {   // für jede Nullstelle
914           zwerg = tmp - (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xk+1, xkoord+2
915          for ( mtest= 0; mtest < anzr; mtest++ )
916          {
917            //          if ( tmp == (*mu[xkoord])[mtest] )
918            //          {
919            if ( ((zwerg.real() <= (*mu[xkoord])[mtest].real() + mprec) &&
920                  (zwerg.real() >= (*mu[xkoord])[mtest].real() - mprec)) &&
921                 ((zwerg.imag() <= (*mu[xkoord])[mtest].imag() + mprec) &&
922                  (zwerg.imag() >= (*mu[xkoord])[mtest].imag() - mprec)) )
923            {
924              roots[xk]->swapRoots( r, rtest );
925              found= true;
926              break;
927            }
928          }
929        } // rtest
930        if (!found) 
931        {
932          WarnS("rootArranger::arrange: precision lost");
933          mprec*=10;
934        }
935      } while(!found);
936#if 0
937      if ( !found )
938      {
939        Warn("rootArranger::arrange: No match? coord %d, root %d.",xkoord,r);
940//#ifdef mprDEBUG_PROT
941        WarnS("One of these ...");
942        for ( rtest= r; rtest < anzr; rtest++ )
943        {
944          tmp= gmp_complex();
945          for ( xk =0; xk <= xkoord; xk++ )
946          {
947            tmp-= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1);
948          }
949          tmp-= (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xkoord+2
950          Warn("  %s",complexToStr(tmp,gmp_output_digits+1),rtest);
951        }
952        WarnS(" ... must match to one of these:");
953        for ( mtest= 0; mtest < anzr; mtest++ )
954        {
955          Warn("                  %s",complexToStr((*mu[xkoord])[mtest],gmp_output_digits+1));
956        }
957//#endif
958      }
959#endif     
960    } // r
961  } // xkoord
962}
963//<-
964
965//-----------------------------------------------------------------------------
966//-------------- simplex ----- ------------------------------------------------
967//-----------------------------------------------------------------------------
968
969//  #ifdef mprDEBUG_PROT
970//  #define error(a) a
971//  #else
972//  #define error(a)
973//  #endif
974
975#define error(a) a
976
977#define MAXPOINTS      1000
978
979//-> simplex::*
980//
981simplex::simplex( int rows, int cols )
982   : LiPM_cols(cols), LiPM_rows(rows)
983{
984  int i;
985
986  LiPM_rows=LiPM_rows+3;
987  LiPM_cols=LiPM_cols+2;
988
989  LiPM = (mprfloat **)omAlloc( LiPM_rows * sizeof(mprfloat *) );  // LP matrix
990  for( i= 0; i < LiPM_rows; i++ )
991  {
992    // Mem must be allocated aligned, also for type double!
993    LiPM[i] = (mprfloat *)omAlloc0Aligned( LiPM_cols * sizeof(mprfloat) );
994  }
995
996  iposv = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
997  izrov = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
998
999  m=n=m1=m2=m3=icase=0;
1000
1001#ifdef mprDEBUG_ALL
1002  Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols);
1003#endif
1004}
1005
1006simplex::~simplex()
1007{
1008  // clean up
1009  int i;
1010  for( i= 0; i < LiPM_rows; i++ )
1011  {
1012    omFreeSize( (ADDRESS) LiPM[i], LiPM_cols * sizeof(mprfloat) );
1013  }
1014  omFreeSize( (ADDRESS) LiPM, LiPM_rows * sizeof(mprfloat *) );
1015
1016  omFreeSize( (ADDRESS) iposv, 2*LiPM_rows*sizeof(int) );
1017  omFreeSize( (ADDRESS) izrov, 2*LiPM_rows*sizeof(int) );
1018}
1019
1020BOOLEAN simplex::mapFromMatrix( matrix m )
1021{
1022  int i,j;
1023//    if ( MATROWS( m ) > LiPM_rows ||  MATCOLS( m ) > LiPM_cols ) {
1024//      WarnS("");
1025//      return FALSE;
1026//    }
1027
1028  number coef;
1029  for ( i= 1; i <= MATROWS( m ); i++ )
1030  {
1031     for ( j= 1; j <= MATCOLS( m ); j++ )
1032     {
1033        if ( MATELEM(m,i,j) != NULL )
1034        {
1035           coef= pGetCoeff( MATELEM(m,i,j) );
1036           if ( coef != NULL && !nIsZero(coef) )
1037              LiPM[i][j]= (double)(*(gmp_float*)coef);
1038           //#ifdef mpr_DEBUG_PROT
1039           //Print("%f ",LiPM[i][j]);
1040           //#endif
1041        }
1042     }
1043     //     PrintLn();
1044  }
1045
1046  return TRUE;
1047}
1048
1049matrix simplex::mapToMatrix( matrix m )
1050{
1051  int i,j;
1052//    if ( MATROWS( m ) < LiPM_rows-3 ||  MATCOLS( m ) < LiPM_cols-2 ) {
1053//      WarnS("");
1054//      return NULL;
1055//    }
1056
1057//Print(" %d x %d\n",MATROWS( m ),MATCOLS( m ));
1058
1059  number coef;
1060  gmp_float * bla;
1061  for ( i= 1; i <= MATROWS( m ); i++ )
1062  {
1063    for ( j= 1; j <= MATCOLS( m ); j++ )
1064    {
1065       pDelete( &(MATELEM(m,i,j)) );
1066       MATELEM(m,i,j)= NULL;
1067//Print(" %3.0f ",LiPM[i][j]);
1068       if ( LiPM[i][j] != 0.0 )
1069       {
1070          bla= new gmp_float(LiPM[i][j]);
1071          coef= (number)bla;
1072          MATELEM(m,i,j)= pOne();
1073          pSetCoeff( MATELEM(m,i,j), coef );
1074       }
1075    }
1076//PrintLn();
1077  }
1078
1079  return m;
1080}
1081
1082intvec * simplex::posvToIV()
1083{
1084   int i;
1085   intvec * iv = new intvec( m );
1086   for ( i= 1; i <= m; i++ )
1087   {
1088      IMATELEM(*iv,i,1)= iposv[i];
1089   }
1090   return iv;
1091}
1092
1093intvec * simplex::zrovToIV()
1094{
1095   int i;
1096   intvec * iv = new intvec( n );
1097   for ( i= 1; i <= n; i++ )
1098   {
1099      IMATELEM(*iv,i,1)= izrov[i];
1100   }
1101   return iv;
1102}
1103
1104void simplex::compute()
1105{
1106  int i,ip,ir,is,k,kh,kp,m12,nl1,nl2;
1107  int *l1,*l2,*l3;
1108  mprfloat q1, bmax;
1109
1110  if ( m != (m1+m2+m3) )
1111  {
1112    // error: bad input
1113    error(WarnS("simplex::compute: Bad input constraint counts!");)
1114    icase=-2;
1115    return;
1116  }
1117
1118  l1= (int *) omAlloc0( (n+1) * sizeof(int) );
1119  l2= (int *) omAlloc0( (m+1) * sizeof(int) );
1120  l3= (int *) omAlloc0( (m+1) * sizeof(int) );
1121
1122  nl1= n;
1123  for ( k=1; k<=n; k++ ) l1[k]=izrov[k]=k;
1124  nl2=m;
1125  for ( i=1; i<=m; i++ )
1126  {
1127    if ( LiPM[i+1][1] < 0.0 )
1128    {
1129      // error: bad input
1130      error(WarnS("simplex::compute: Bad input tableau!");)
1131      error(Warn("simplex::compute: in input Matrix row %d, column 1, value %f",i+1,LiPM[i+1][1]);)
1132      icase=-2;
1133      // free mem l1,l2,l3;
1134      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1135      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1136      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1137      return;
1138    }
1139    l2[i]= i;
1140    iposv[i]= n+i;
1141  }
1142  for ( i=1; i<=m2; i++) l3[i]= 1;
1143  ir= 0;
1144  if (m2+m3)
1145  {
1146    ir=1;
1147    for ( k=1; k <= (n+1); k++ )
1148    {
1149      q1=0.0;
1150      for ( i=m1+1; i <= m; i++ ) q1+= LiPM[i+1][k];
1151      LiPM[m+2][k]= -q1;
1152    }
1153
1154    do
1155    {
1156      simp1(LiPM,m+1,l1,nl1,0,&kp,&bmax);
1157      if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] < -SIMPLEX_EPS )
1158      {
1159        icase= -1; // no solution found
1160        // free mem l1,l2,l3;
1161        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1162        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1163        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1164        return;
1165      }
1166      else if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] <= SIMPLEX_EPS )
1167      {
1168        m12= m1+m2+1;
1169        if ( m12 <= m )
1170        {
1171          for ( ip= m12; ip <= m; ip++ )
1172          {
1173            if ( iposv[ip] == (ip+n) )
1174            {
1175              simp1(LiPM,ip,l1,nl1,1,&kp,&bmax);
1176              if ( fabs(bmax) >= SIMPLEX_EPS)
1177                goto one;
1178            }
1179          }
1180        }
1181        ir= 0;
1182        --m12;
1183        if ( m1+1 <= m12 )
1184          for ( i=m1+1; i <= m12; i++ )
1185            if ( l3[i-m1] == 1 )
1186              for ( k=1; k <= n+1; k++ )
1187                LiPM[i+1][k] = -(LiPM[i+1][k]);
1188        break;
1189      }
1190      //#if DEBUG
1191      //print_bmat( a, m+2, n+3);
1192      //#endif
1193      simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1194      if ( ip == 0 )
1195      {
1196        icase = -1; // no solution found
1197        // free mem l1,l2,l3;
1198        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1199        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1200        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1201        return;
1202      }
1203    one: simp3(LiPM,m+1,n,ip,kp);
1204      // #if DEBUG
1205      // print_bmat(a,m+2,n+3);
1206      // #endif
1207      if ( iposv[ip] >= (n+m1+m2+1))
1208      {
1209        for ( k= 1; k <= nl1; k++ )
1210          if ( l1[k] == kp ) break;
1211        --nl1;
1212        for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1];
1213        ++(LiPM[m+2][kp+1]);
1214        for ( i= 1; i <= m+2; i++ ) LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1215      }
1216      else
1217      {
1218        if ( iposv[ip] >= (n+m1+1) )
1219        {
1220          kh= iposv[ip]-m1-n;
1221          if ( l3[kh] )
1222          {
1223            l3[kh]= 0;
1224            ++(LiPM[m+2][kp+1]);
1225            for ( i=1; i<= m+2; i++ )
1226              LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1227          }
1228        }
1229      }
1230      is= izrov[kp];
1231      izrov[kp]= iposv[ip];
1232      iposv[ip]= is;
1233    } while (ir);
1234  }
1235  /* end of phase 1, have feasible sol, now optimize it */
1236  loop
1237  {
1238    // #if DEBUG
1239    // print_bmat( a, m+1, n+5);
1240    // #endif
1241    simp1(LiPM,0,l1,nl1,0,&kp,&bmax);
1242    if (bmax <= /*SIMPLEX_EPS*/0.0)
1243    {
1244      icase=0; // finite solution found
1245      // free mem l1,l2,l3
1246      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1247      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1248      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1249      return;
1250    }
1251    simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1252    if (ip == 0)
1253    {
1254      //printf("Unbounded:");
1255      // #if DEBUG
1256      //       print_bmat( a, m+1, n+1);
1257      // #endif
1258      icase=1;                /* unbounded */
1259      // free mem
1260      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1261      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1262      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1263      return;
1264    }
1265    simp3(LiPM,m,n,ip,kp);
1266    is= izrov[kp];
1267    izrov[kp]= iposv[ip];
1268    iposv[ip]= is;
1269  }/*for ;;*/
1270}
1271
1272void simplex::simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax )
1273{
1274  int k;
1275  mprfloat test;
1276
1277  if( nll <= 0)
1278  {                        /* init'tion: fixed */
1279    *bmax = 0.0;
1280    return;
1281  }
1282  *kp=ll[1];
1283  *bmax=a[mm+1][*kp+1];
1284  for (k=2;k<=nll;k++)
1285  {
1286    if (iabf == 0)
1287    {
1288      test=a[mm+1][ll[k]+1]-(*bmax);
1289      if (test > 0.0)
1290      {
1291        *bmax=a[mm+1][ll[k]+1];
1292        *kp=ll[k];
1293      }
1294    }
1295    else
1296    {                        /* abs values: have fixed it */
1297      test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
1298      if (test > 0.0)
1299      {
1300        *bmax=a[mm+1][ll[k]+1];
1301        *kp=ll[k];
1302      }
1303    }
1304  }
1305}
1306
1307void simplex::simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 )
1308{
1309  int k,ii,i;
1310  mprfloat qp,q0,q;
1311
1312  *ip= 0;
1313  for ( i=1; i <= nl2; i++ )
1314  {
1315    if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS )
1316    {
1317      *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1];
1318      *ip= l2[i];
1319      for ( i= i+1; i <= nl2; i++ )
1320      {
1321        ii= l2[i];
1322        if (a[ii+1][kp+1] < -SIMPLEX_EPS)
1323        {
1324          q= -a[ii+1][1] / a[ii+1][kp+1];
1325          if (q - *q1 < -SIMPLEX_EPS)
1326          {
1327            *ip=ii;
1328            *q1=q;
1329          }
1330          else if (q - *q1 < SIMPLEX_EPS)
1331          {
1332            for ( k=1; k<= n; k++ )
1333            {
1334              qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
1335              q0= -a[ii+1][k+1]/a[ii+1][kp+1];
1336              if ( q0 != qp ) break;
1337            }
1338            if ( q0 < qp ) *ip= ii;
1339          }
1340        }
1341      }
1342    }
1343  }
1344}
1345
1346void simplex::simp3( mprfloat **a, int i1, int k1, int ip, int kp )
1347{
1348  int kk,ii;
1349  mprfloat piv;
1350
1351  piv= 1.0 / a[ip+1][kp+1];
1352  for ( ii=1; ii <= i1+1; ii++ )
1353  {
1354    if ( ii -1 != ip )
1355    {
1356      a[ii][kp+1] *= piv;
1357      for ( kk=1; kk <= k1+1; kk++ )
1358        if ( kk-1 != kp )
1359          a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
1360    }
1361  }
1362  for ( kk=1; kk<= k1+1; kk++ )
1363    if ( kk-1 != kp ) a[ip+1][kk] *= -piv;
1364  a[ip+1][kp+1]= piv;
1365}
1366//<-
1367
1368//-----------------------------------------------------------------------------
1369
1370//#endif // HAVE_MPR
1371
1372// local Variables: ***
1373// folded-file: t ***
1374// compile-command-1: "make installg" ***
1375// compile-command-2: "make install" ***
1376// End: ***
Note: See TracBrowser for help on using the repository browser.