source: git/kernel/numeric/mpr_numeric.cc @ 4f8fd1d

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