source: git/coeffs/mpr_numeric.cc @ 634d93b

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