source: git/kernel/mpr_numeric.cc @ c30c46

spielwiese
Last change on this file since c30c46 was c30c46, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: use !isZero() instead of !=zero: PORTING stuff git-svn-id: file:///usr/local/Singular/svn/trunk@8028 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 30.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5/* $Id: mpr_numeric.cc,v 1.3 2005-05-04 07:23:17 Singular Exp $ */
6
7/*
8* ABSTRACT - multipolynomial resultants - numeric stuff
9*            ( root finder, vandermonde system solver, simplex )
10*/
11
12#include "mod2.h"
13//#ifdef HAVE_MPR
14
15//#define mprDEBUG_ALL
16
17//-> includes
18#include "structs.h"
19#include "febase.h"
20#include "omalloc.h"
21#include "numbers.h"
22#include "polys.h"
23#include "ideals.h"
24#include "intvec.h"
25#include "matpol.h"
26#include "ring.h"
27//#include "longrat.h"
28
29#include <math.h>
30#include "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    ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) );
333    for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] );
334  }
335
336  theroots= NULL;
337  found_roots= false;
338}
339//<-
340
341//-> poly rootContainer::getPoly()
342poly rootContainer::getPoly()
343{
344  int i;
345
346  poly result= NULL;
347  poly ppos;
348
349  if ( (rt == cspecial) || ( rt == cspecialmu ) )
350  {
351    for ( i= tdg; i >= 0; i-- )
352    {
353      if ( coeffs[i] )
354      {
355        poly p= pOne();
356        //pSetExp( p, var+1, i);
357        pSetExp( p, 1, i);
358        pSetCoeff( p, nCopy( coeffs[i] ) );
359        pSetm( p );
360        if (result)
361        {
362          ppos->next=p;
363          ppos=ppos->next;
364        }
365        else
366        {
367          result=p;
368          ppos=p;
369        }
370
371      }
372    }
373    if (result!=NULL) pSetm( result );
374  }
375
376  return result;
377}
378//<-
379
380//-> const gmp_complex & rootContainer::opterator[] ( const int i )
381// this is now inline!
382//  gmp_complex & rootContainer::operator[] ( const int i )
383//  {
384//    if ( found_roots && ( i >= 0) && ( i < tdg ) )
385//    {
386//      return *theroots[i];
387//    }
388//    // warning
389//    Warn("rootContainer::getRoot: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
390//    gmp_complex *tmp= new gmp_complex();
391//    return *tmp;
392//  }
393//<-
394
395//-> gmp_complex & rootContainer::evPointCoord( int i )
396gmp_complex & rootContainer::evPointCoord( const int i )
397{
398  if (! ((i >= 0) && (i < anz+2) ) )
399    WarnS("rootContainer::evPointCoord: index out of range");
400  if (ievpoint == NULL)
401    WarnS("rootContainer::evPointCoord: ievpoint == NULL");
402
403  if ( (rt == cspecialmu) && found_roots ) {  // FIX ME
404    if ( ievpoint[i] != NULL )
405    {
406      gmp_complex *tmp= new gmp_complex();
407      *tmp= numberToComplex(ievpoint[i]);
408      return *tmp;
409    }
410    else
411    {
412      Warn("rootContainer::evPointCoord: NULL index %d",i);
413    }
414  }
415
416  // warning
417  Warn("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
418  gmp_complex *tmp= new gmp_complex();
419  return *tmp;
420}
421//<-
422
423//-> bool rootContainer::changeRoots( int from, int to )
424bool rootContainer::swapRoots( const int from, const int to )
425{
426  if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) )
427  {
428    if ( to != from )
429    {
430      gmp_complex tmp( *theroots[from] );
431      *theroots[from]= *theroots[to];
432      *theroots[to]= tmp;
433    }
434    return true;
435  }
436
437  // warning
438  Warn(" rootContainer::changeRoots: Wrong index %d, %d",from,to);
439  return false;
440}
441//<-
442
443//-> void rootContainer::solver()
444bool rootContainer::solver( const int polishmode )
445{
446  int i;
447
448  // there are maximal tdg roots, so *roots ranges form 0 to tdg-1.
449  theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) );
450  for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex();
451
452  // copy the coefficients of type number to type gmp_complex
453  gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
454  for ( i=0; i <= tdg; i++ )
455  {
456    ad[i]= new gmp_complex();
457    if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i] );
458  }
459
460  // now solve
461  found_roots= laguer_driver( ad, theroots, polishmode != 0 );
462  if (!found_roots)
463    WarnS("rootContainer::solver: No roots found!");
464 
465  // free memory
466  for ( i=0; i <= tdg; i++ ) delete ad[i];
467  omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof(gmp_complex*) );
468
469  return found_roots;
470}
471//<-
472
473//-> gmp_complex* rootContainer::laguer_driver( bool polish )
474bool rootContainer::laguer_driver(gmp_complex ** a, gmp_complex ** roots, bool polish )
475{
476  int i,j,k,its;
477  gmp_float zero(0.0);
478  gmp_complex x(0.0),o(1.0);
479  bool ret= true, isf=isfloat(a), type=true;
480
481  gmp_complex ** ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
482  for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] );
483
484  k = 0;
485  i = tdg;
486  j = i-1;
487  while (i>2)
488  {
489    // run laguer alg
490    x = zero;
491    laguer(ad, i, &x, &its, type);
492    if ( its > MAXIT )
493    {
494      type = !type;
495      x = zero;
496      laguer(ad, i, &x, &its, type);
497    }
498
499    mprSTICKYPROT(ST_ROOTS_LGSTEP);
500    if ( its > MAXIT )
501    {  // error
502      WarnS("rootContainer::laguer_driver: To many iterations!");
503      ret= false;
504      goto theend;
505    }
506    if ( polish )
507    {
508      laguer( a, tdg, &x, &its , type);
509      if ( its > MAXIT )
510      {  // error
511        WarnS("rootContainer::laguer_driver: To many iterations in polish!");
512        ret= false;
513        goto theend;
514      }
515    }
516    if ((!type)&&(!((x.real()==zero)&&(x.imag()==zero)))) x = o/x;
517    if (x.imag() == zero)
518    {
519      *roots[k] = x;
520      k++;
521      divlin(ad,x,i);
522      i--;
523    }
524    else
525    {
526      if(isf)
527      {
528        *roots[j] = x;
529        *roots[j-1]= gmp_complex(x.real(),-x.imag());
530        j -= 2;
531        divquad(ad,x,i);
532        i -= 2;
533      }
534      else
535      {
536        *roots[j] = x;
537        j--;
538        divlin(ad,x,i);
539        i--;
540      }
541    }
542    type = !type;
543  }
544  solvequad(ad,roots,k,j);
545  sortroots(roots,k,j,isf);
546
547theend:
548  mprSTICKYPROT("\n");
549  for ( i=0; i <= tdg; i++ ) delete ad[i];
550  omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof( gmp_complex* ));
551
552  return ret;
553}
554//<-
555
556//-> void rootContainer::laguer(...)
557void rootContainer::laguer(gmp_complex ** a, int m, gmp_complex *x, int *its, bool type)
558{
559  int iter,j;
560  gmp_float zero(0.0),one(1.0),deg(m);
561  gmp_float abx_g, err_g, fabs;
562  gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2;
563  gmp_float frac_g[MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0 };
564
565  gmp_float epss(0.1);
566  mpf_pow_ui(*epss._mpfp(),*epss.mpfp(),getGMPFloatDigits());
567
568  for ( iter= 1; iter <= MAXIT; iter++ )
569  {
570    mprSTICKYPROT(ST_ROOTS_LG);
571    *its=iter;
572    if (type)
573      computefx(a,*x,m,b,d,f,abx_g,err_g);
574    else
575      computegx(a,*x,m,b,d,f,abx_g,err_g);
576    err_g *= epss; // EPSS;
577
578    fabs = abs(b);
579    if (fabs <= err_g)
580    {
581      if ((fabs==zero) || (abs(d)==zero)) return;
582      *x -= (b/d); // a last newton-step
583      goto ende;
584    }
585
586    g= d / b;
587    g2 = g * g;
588    h= g2 - (((f+f) / b ));
589    sq= sqrt(( ( h * deg ) - g2 ) * (deg - one));
590    gp= g + sq;
591    gm= g - sq;
592    if (abs(gp)<abs(gm))
593    {
594      dx = deg/gm;
595    }
596    else
597    {
598      if((gp.real()==zero)&&(gp.imag()==zero))
599      {
600        dx.real(cos((mprfloat)iter));
601        dx.imag(sin((mprfloat)iter));
602        dx = dx*(one+abx_g);
603      }
604      else
605      {
606        dx = deg/gp;
607      }
608    }
609    x1= *x - dx;
610
611    if (*x == x1) goto ende;
612
613    j = iter%MMOD;
614    if (j==0) j=MT;
615    if ( j % MT ) *x= x1;
616    else *x -= ( dx * frac_g[ j / MT ] );
617  }
618
619  *its= MAXIT+1;
620ende:
621  checkimag(x,epss);
622}
623
624void rootContainer::checkimag(gmp_complex *x, gmp_float &e)
625{
626  if(abs(x->imag())<abs(x->real())*e)
627  {
628    x->imag(0.0);
629  }
630}
631
632bool rootContainer::isfloat(gmp_complex **a)
633{
634  gmp_float z(0.0);
635  gmp_complex *b;
636  for (int i=tdg; i >= 0; i-- )
637  {
638    b = &(*a[i]);
639    if (!(b->imag()==z))
640      return false;
641  }
642  return true;
643}
644
645void rootContainer::divlin(gmp_complex **a, gmp_complex x, int j)
646{
647  int i;
648  gmp_float o(1.0);
649
650  if (abs(x)<o)
651  {
652    for (i= j-1; i > 0; i-- )
653      *a[i] += (*a[i+1]*x);
654    for (i= 0; i < j; i++ )
655      *a[i] = *a[i+1];
656  }
657  else
658  {
659    gmp_complex y(o/x);
660    for (i= 1; i < j; i++)
661      *a[i] += (*a[i-1]*y);
662  }
663}
664
665void rootContainer::divquad(gmp_complex **a, gmp_complex x, int j)
666{
667  int i;
668  gmp_float o(1.0),p(x.real()+x.real()),
669            q((x.real()*x.real())+(x.imag()*x.imag()));
670
671  if (abs(x)<o)
672  {
673    *a[j-1] += (*a[j]*p);
674    for (i= j-2; i > 1; i-- )
675      *a[i] += ((*a[i+1]*p)-(*a[i+2]*q));
676    for (i= 0; i < j-1; i++ )
677      *a[i] = *a[i+2];
678  }
679  else
680  {
681    p = p/q;
682    q = o/q;
683    *a[1] += (*a[0]*p);
684    for (i= 2; i < j-1; i++)
685      *a[i] += ((*a[i-1]*p)-(*a[i-2]*q));
686  }
687}
688
689void rootContainer::solvequad(gmp_complex **a, gmp_complex **r, int &k, int &j)
690{
691  gmp_float zero(0.0);
692
693  if ((j>k)
694  &&((!(*a[2]).real().isZero())||(!(*a[2]).imag().isZero())))
695  {
696    gmp_complex sq(zero);
697    gmp_complex h1(*a[1]/(*a[2] + *a[2])), h2(*a[0] / *a[2]);
698    gmp_complex disk((h1 * h1) - h2);
699    if (disk.imag()==zero)
700    {
701      if (disk.real()<zero)
702      {
703        sq.real(zero);
704        sq.imag(sqrt(-disk.real()));
705      }
706      else
707        sq = (gmp_complex)sqrt(disk.real());
708    }
709    else
710      sq = sqrt(disk);
711    *r[k+1] = sq - h1;
712    sq += h1;
713    *r[k] = (gmp_complex)0.0-sq;
714    if(sq.imag()==zero)
715    {
716      k = j;
717      j++;
718    }
719    else
720    {
721      j = k;
722      k--;
723    }
724  }
725  else
726  {
727    if (((*a[1]).real()==zero) && ((*a[1]).imag()==zero))
728    {
729      WerrorS("precision lost, try again with higher precision");
730    }
731    else
732    {
733      *r[k]= (gmp_complex)0.0-(*a[0] / *a[1]);
734      if(r[k]->imag()==zero)
735        j++;
736      else
737        k--;
738    }
739  }
740}
741
742void rootContainer::sortroots(gmp_complex **ro, int r, int c, bool isf)
743{
744  int j;
745
746  for (j=0; j<r; j++) // the real roots
747    sortre(ro, j, r, 1);
748  if (c>=tdg) return;
749  if (isf)
750  {
751    for (j=c; j+2<tdg; j+=2) // the complex roots for a real poly
752      sortre(ro, j, tdg-1, 2);
753  }
754  else
755  {
756    for (j=c; j+1<tdg; j++) // the complex roots for a general poly
757      sortre(ro, j, tdg-1, 1);
758  }
759}
760
761void rootContainer::sortre(gmp_complex **r, int l, int u, int inc)
762{
763  int pos,i;
764  gmp_complex *x,*y;
765
766  pos = l;
767  x = r[pos];
768  for (i=l+inc; i<=u; i+=inc)
769  {
770    if (r[i]->real()<x->real())
771    {
772      pos = i;
773      x = r[pos];
774    }
775  }
776  if (pos>l)
777  {
778    if (inc==1)
779    {
780      for (i=pos; i>l; i--)
781        r[i] = r[i-1];
782      r[l] = x;
783    }
784    else
785    {
786      y = r[pos+1];
787      for (i=pos+1; i+1>l; i--)
788        r[i] = r[i-2];
789      if (x->imag()>y->imag())
790      {
791        r[l] = x;
792        r[l+1] = y;
793      }
794      else
795      {
796        r[l] = y;
797        r[l+1] = x;
798      }
799    }
800  }
801  else if ((inc==2)&&(x->imag()<r[l+1]->imag()))
802  {
803    r[l] = r[l+1];
804    r[l+1] = x;
805  }
806}
807
808void rootContainer::computefx(gmp_complex **a, gmp_complex x, int m,
809                  gmp_complex &f0, gmp_complex &f1, gmp_complex &f2,
810                  gmp_float &ex, gmp_float &ef)
811{
812  int k;
813
814  f0= *a[m];
815  ef= abs(f0);
816  f1= gmp_complex(0.0);
817  f2= f1;
818  ex= abs(x);
819
820  for ( k= m-1; k >= 0; k-- )
821  {
822    f2 = ( x * f2 ) + f1;
823    f1 = ( x * f1 ) + f0;
824    f0 = ( x * f0 ) + *a[k];
825    ef = abs( f0 ) + ( ex * ef );
826  }
827}
828
829void rootContainer::computegx(gmp_complex **a, gmp_complex x, int m,
830                  gmp_complex &f0, gmp_complex &f1, gmp_complex &f2,
831                  gmp_float &ex, gmp_float &ef)
832{
833  int k;
834
835  f0= *a[0];
836  ef= abs(f0);
837  f1= gmp_complex(0.0);
838  f2= f1;
839  ex= abs(x);
840
841  for ( k= 1; k <= m; k++ )
842  {
843    f2 = ( x * f2 ) + f1;
844    f1 = ( x * f1 ) + f0;
845    f0 = ( x * f0 ) + *a[k];
846    ef = abs( f0 ) + ( ex * ef );
847  }
848}
849
850//-----------------------------------------------------------------------------
851//-------------- rootArranger -------------------------------------------------
852//-----------------------------------------------------------------------------
853
854//-> rootArranger::rootArranger(...)
855rootArranger::rootArranger( rootContainer ** _roots,
856                            rootContainer ** _mu,
857                            const int _howclean )
858  : roots(_roots), mu(_mu), howclean(_howclean)
859{
860  found_roots=false;
861}
862//<-
863
864//-> void rootArranger::solve_all()
865void rootArranger::solve_all()
866{
867  int i;
868  found_roots= true;
869
870  // find roots of polys given by coeffs in roots
871  rc= roots[0]->getAnzElems();
872  for ( i= 0; i < rc; i++ )
873    if ( !roots[i]->solver( howclean ) )
874    {
875      found_roots= false;
876      return;
877    }
878  // find roots of polys given by coeffs in mu
879  mc= mu[0]->getAnzElems();
880  for ( i= 0; i < mc; i++ )
881    if ( ! mu[i]->solver( howclean ) )
882    {
883      found_roots= false;
884      return;
885    }
886}
887//<-
888
889//-> void rootArranger::arrange()
890void rootArranger::arrange()
891{
892  gmp_complex tmp,zwerg;
893  int anzm= mu[0]->getAnzElems();
894  int anzr= roots[0]->getAnzRoots();
895  int xkoord, r, rtest, xk, mtest;
896  bool found;
897  //gmp_complex mprec(1.0/pow(10,gmp_output_digits-5),1.0/pow(10,gmp_output_digits-5));
898  gmp_float mprec(1.0/pow(10.0,(int)(gmp_output_digits/3)));
899
900  for ( xkoord= 0; xkoord < anzm; xkoord++ ) {    // für x1,x2, x1,x2,x3, x1,x2,...,xn
901    for ( r= 0; r < anzr; r++ ) {                 // für jede Nullstelle
902      // (x1-koordinate) * evp[1] + (x2-koordinate) * evp[2] +
903      //                                  ... + (xkoord-koordinate) * evp[xkoord]
904      tmp= gmp_complex();
905      for ( xk =0; xk <= xkoord; xk++ )
906      {
907        tmp -= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); //xk+1
908      }
909      found= false;
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        Warn("rootArranger::arrange: No match? coord %d, root %d.",xkoord,r);
930#ifdef mprDEBUG_PROT
931        WarnS("One of these ...");
932        for ( rtest= r; rtest < anzr; rtest++ )
933        {
934          tmp= gmp_complex();
935          for ( xk =0; xk <= xkoord; xk++ )
936          {
937            tmp-= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1);
938          }
939          tmp-= (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xkoord+2
940          Warn("  %s",complexToStr(tmp,gmp_output_digits+1),rtest);
941        }
942        WarnS(" ... must match to one of these:");
943        for ( mtest= 0; mtest < anzr; mtest++ )
944        {
945          Warn("                  %s",complexToStr((*mu[xkoord])[mtest],gmp_output_digits+1));
946        }
947#endif
948      }
949    } // r
950  } // xkoord
951}
952//<-
953
954//-----------------------------------------------------------------------------
955//-------------- simplex ----- ------------------------------------------------
956//-----------------------------------------------------------------------------
957
958//  #ifdef mprDEBUG_PROT
959//  #define error(a) a
960//  #else
961//  #define error(a)
962//  #endif
963
964#define error(a) a
965
966#define MAXPOINTS      1000
967
968//-> simplex::*
969//
970simplex::simplex( int rows, int cols )
971   : LiPM_cols(cols), LiPM_rows(rows)
972{
973  int i;
974
975  LiPM_rows=LiPM_rows+3;
976  LiPM_cols=LiPM_cols+2;
977
978  LiPM = (mprfloat **)omAlloc( LiPM_rows * sizeof(mprfloat *) );  // LP matrix
979  for( i= 0; i < LiPM_rows; i++ )
980  {
981    // Mem must be allocated aligned, also for type double!
982    LiPM[i] = (mprfloat *)omAlloc0Aligned( LiPM_cols * sizeof(mprfloat) );
983  }
984
985  iposv = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
986  izrov = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
987
988  m=n=m1=m2=m3=icase=0;
989
990#ifdef mprDEBUG_ALL
991  Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols);
992#endif
993}
994
995simplex::~simplex()
996{
997  // clean up
998  int i;
999  for( i= 0; i < LiPM_rows; i++ )
1000  {
1001    omFreeSize( (ADDRESS) LiPM[i], LiPM_cols * sizeof(mprfloat) );
1002  }
1003  omFreeSize( (ADDRESS) LiPM, LiPM_rows * sizeof(mprfloat *) );
1004
1005  omFreeSize( (ADDRESS) iposv, 2*LiPM_rows*sizeof(int) );
1006  omFreeSize( (ADDRESS) izrov, 2*LiPM_rows*sizeof(int) );
1007}
1008
1009BOOLEAN simplex::mapFromMatrix( matrix m )
1010{
1011  int i,j;
1012//    if ( MATROWS( m ) > LiPM_rows ||  MATCOLS( m ) > LiPM_cols ) {
1013//      WarnS("");
1014//      return FALSE;
1015//    }
1016
1017  number coef;
1018  for ( i= 1; i <= MATROWS( m ); i++ )
1019  {
1020     for ( j= 1; j <= MATCOLS( m ); j++ )
1021     {
1022        if ( MATELEM(m,i,j) != NULL )
1023        {
1024           coef= pGetCoeff( MATELEM(m,i,j) );
1025           if ( coef != NULL && !nIsZero(coef) )
1026              LiPM[i][j]= (double)(*(gmp_float*)coef);
1027           //#ifdef mpr_DEBUG_PROT
1028           //Print("%f ",LiPM[i][j]);
1029           //#endif
1030        }
1031     }
1032     //     PrintLn();
1033  }
1034
1035  return TRUE;
1036}
1037
1038matrix simplex::mapToMatrix( matrix m )
1039{
1040  int i,j;
1041//    if ( MATROWS( m ) < LiPM_rows-3 ||  MATCOLS( m ) < LiPM_cols-2 ) {
1042//      WarnS("");
1043//      return NULL;
1044//    }
1045
1046//Print(" %d x %d\n",MATROWS( m ),MATCOLS( m ));
1047
1048  number coef;
1049  gmp_float * bla;
1050  for ( i= 1; i <= MATROWS( m ); i++ )
1051  {
1052    for ( j= 1; j <= MATCOLS( m ); j++ )
1053    {
1054       pDelete( &(MATELEM(m,i,j)) );
1055       MATELEM(m,i,j)= NULL;
1056//Print(" %3.0f ",LiPM[i][j]);
1057       if ( LiPM[i][j] != 0.0 )
1058       {
1059          bla= new gmp_float(LiPM[i][j]);
1060          coef= (number)bla;
1061          MATELEM(m,i,j)= pOne();
1062          pSetCoeff( MATELEM(m,i,j), coef );
1063       }
1064    }
1065//PrintLn();
1066  }
1067
1068  return m;
1069}
1070
1071intvec * simplex::posvToIV()
1072{
1073   int i;
1074   intvec * iv = new intvec( m );
1075   for ( i= 1; i <= m; i++ )
1076   {
1077      IMATELEM(*iv,i,1)= iposv[i];
1078   }
1079   return iv;
1080}
1081
1082intvec * simplex::zrovToIV()
1083{
1084   int i;
1085   intvec * iv = new intvec( n );
1086   for ( i= 1; i <= n; i++ )
1087   {
1088      IMATELEM(*iv,i,1)= izrov[i];
1089   }
1090   return iv;
1091}
1092
1093void simplex::compute()
1094{
1095  int i,ip,ir,is,k,kh,kp,m12,nl1,nl2;
1096  int *l1,*l2,*l3;
1097  mprfloat q1, bmax;
1098
1099  if ( m != (m1+m2+m3) )
1100  {
1101    // error: bad input
1102    error(WarnS("simplex::compute: Bad input constraint counts!");)
1103    icase=-2;
1104    return;
1105  }
1106
1107  l1= (int *) omAlloc0( (n+1) * sizeof(int) );
1108  l2= (int *) omAlloc0( (m+1) * sizeof(int) );
1109  l3= (int *) omAlloc0( (m+1) * sizeof(int) );
1110
1111  nl1= n;
1112  for ( k=1; k<=n; k++ ) l1[k]=izrov[k]=k;
1113  nl2=m;
1114  for ( i=1; i<=m; i++ )
1115  {
1116    if ( LiPM[i+1][1] < 0.0 )
1117    {
1118      // error: bad input
1119      error(WarnS("simplex::compute: Bad input tableau!");)
1120      error(Warn("simplex::compute: in input Matrix row %d, column 1, value %f",i+1,LiPM[i+1][1]);)
1121      icase=-2;
1122      // free mem l1,l2,l3;
1123      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1124      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1125      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1126      return;
1127    }
1128    l2[i]= i;
1129    iposv[i]= n+i;
1130  }
1131  for ( i=1; i<=m2; i++) l3[i]= 1;
1132  ir= 0;
1133  if (m2+m3)
1134  {
1135    ir=1;
1136    for ( k=1; k <= (n+1); k++ )
1137    {
1138      q1=0.0;
1139      for ( i=m1+1; i <= m; i++ ) q1+= LiPM[i+1][k];
1140      LiPM[m+2][k]= -q1;
1141    }
1142
1143    do
1144    {
1145      simp1(LiPM,m+1,l1,nl1,0,&kp,&bmax);
1146      if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] < -SIMPLEX_EPS )
1147      {
1148        icase= -1; // no solution found
1149        // free mem l1,l2,l3;
1150        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1151        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1152        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1153        return;
1154      }
1155      else if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] <= SIMPLEX_EPS )
1156      {
1157        m12= m1+m2+1;
1158        if ( m12 <= m )
1159        {
1160          for ( ip= m12; ip <= m; ip++ )
1161          {
1162            if ( iposv[ip] == (ip+n) )
1163            {
1164              simp1(LiPM,ip,l1,nl1,1,&kp,&bmax);
1165              if ( fabs(bmax) >= SIMPLEX_EPS)
1166                goto one;
1167            }
1168          }
1169        }
1170        ir= 0;
1171        --m12;
1172        if ( m1+1 <= m12 )
1173          for ( i=m1+1; i <= m12; i++ )
1174            if ( l3[i-m1] == 1 )
1175              for ( k=1; k <= n+1; k++ )
1176                LiPM[i+1][k] = -(LiPM[i+1][k]);
1177        break;
1178      }
1179      //#if DEBUG
1180      //print_bmat( a, m+2, n+3);
1181      //#endif
1182      simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1183      if ( ip == 0 )
1184      {
1185        icase = -1; // no solution found
1186        // free mem l1,l2,l3;
1187        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1188        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1189        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1190        return;
1191      }
1192    one: simp3(LiPM,m+1,n,ip,kp);
1193      // #if DEBUG
1194      // print_bmat(a,m+2,n+3);
1195      // #endif
1196      if ( iposv[ip] >= (n+m1+m2+1))
1197      {
1198        for ( k= 1; k <= nl1; k++ )
1199          if ( l1[k] == kp ) break;
1200        --nl1;
1201        for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1];
1202        ++(LiPM[m+2][kp+1]);
1203        for ( i= 1; i <= m+2; i++ ) LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1204      }
1205      else
1206      {
1207        if ( iposv[ip] >= (n+m1+1) )
1208        {
1209          kh= iposv[ip]-m1-n;
1210          if ( l3[kh] )
1211          {
1212            l3[kh]= 0;
1213            ++(LiPM[m+2][kp+1]);
1214            for ( i=1; i<= m+2; i++ )
1215              LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1216          }
1217        }
1218      }
1219      is= izrov[kp];
1220      izrov[kp]= iposv[ip];
1221      iposv[ip]= is;
1222    } while (ir);
1223  }
1224  /* end of phase 1, have feasible sol, now optimize it */
1225  loop
1226  {
1227    // #if DEBUG
1228    // print_bmat( a, m+1, n+5);
1229    // #endif
1230    simp1(LiPM,0,l1,nl1,0,&kp,&bmax);
1231    if (bmax <= /*SIMPLEX_EPS*/0.0)
1232    {
1233      icase=0; // finite solution found
1234      // free mem l1,l2,l3
1235      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1236      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1237      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1238      return;
1239    }
1240    simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1241    if (ip == 0)
1242    {
1243      //printf("Unbounded:");
1244      // #if DEBUG
1245      //       print_bmat( a, m+1, n+1);
1246      // #endif
1247      icase=1;                /* unbounded */
1248      // free mem
1249      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1250      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1251      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1252      return;
1253    }
1254    simp3(LiPM,m,n,ip,kp);
1255    is= izrov[kp];
1256    izrov[kp]= iposv[ip];
1257    iposv[ip]= is;
1258  }/*for ;;*/
1259}
1260
1261void simplex::simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax )
1262{
1263  int k;
1264  mprfloat test;
1265
1266  if( nll <= 0)
1267  {                        /* init'tion: fixed */
1268    *bmax = 0.0;
1269    return;
1270  }
1271  *kp=ll[1];
1272  *bmax=a[mm+1][*kp+1];
1273  for (k=2;k<=nll;k++)
1274  {
1275    if (iabf == 0)
1276    {
1277      test=a[mm+1][ll[k]+1]-(*bmax);
1278      if (test > 0.0)
1279      {
1280        *bmax=a[mm+1][ll[k]+1];
1281        *kp=ll[k];
1282      }
1283    }
1284    else
1285    {                        /* abs values: have fixed it */
1286      test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
1287      if (test > 0.0)
1288      {
1289        *bmax=a[mm+1][ll[k]+1];
1290        *kp=ll[k];
1291      }
1292    }
1293  }
1294}
1295
1296void simplex::simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 )
1297{
1298  int k,ii,i;
1299  mprfloat qp,q0,q;
1300
1301  *ip= 0;
1302  for ( i=1; i <= nl2; i++ )
1303  {
1304    if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS )
1305    {
1306      *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1];
1307      *ip= l2[i];
1308      for ( i= i+1; i <= nl2; i++ )
1309      {
1310        ii= l2[i];
1311        if (a[ii+1][kp+1] < -SIMPLEX_EPS)
1312        {
1313          q= -a[ii+1][1] / a[ii+1][kp+1];
1314          if (q - *q1 < -SIMPLEX_EPS)
1315          {
1316            *ip=ii;
1317            *q1=q;
1318          }
1319          else if (q - *q1 < SIMPLEX_EPS)
1320          {
1321            for ( k=1; k<= n; k++ )
1322            {
1323              qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
1324              q0= -a[ii+1][k+1]/a[ii+1][kp+1];
1325              if ( q0 != qp ) break;
1326            }
1327            if ( q0 < qp ) *ip= ii;
1328          }
1329        }
1330      }
1331    }
1332  }
1333}
1334
1335void simplex::simp3( mprfloat **a, int i1, int k1, int ip, int kp )
1336{
1337  int kk,ii;
1338  mprfloat piv;
1339
1340  piv= 1.0 / a[ip+1][kp+1];
1341  for ( ii=1; ii <= i1+1; ii++ )
1342  {
1343    if ( ii -1 != ip )
1344    {
1345      a[ii][kp+1] *= piv;
1346      for ( kk=1; kk <= k1+1; kk++ )
1347        if ( kk-1 != kp )
1348          a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
1349    }
1350  }
1351  for ( kk=1; kk<= k1+1; kk++ )
1352    if ( kk-1 != kp ) a[ip+1][kk] *= -piv;
1353  a[ip+1][kp+1]= piv;
1354}
1355//<-
1356
1357//-----------------------------------------------------------------------------
1358
1359//#endif // HAVE_MPR
1360
1361// local Variables: ***
1362// folded-file: t ***
1363// compile-command-1: "make installg" ***
1364// compile-command-2: "make install" ***
1365// End: ***
Note: See TracBrowser for help on using the repository browser.