source: git/Singular/mpr_numeric.cc @ 5ca9807

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