source: git/Singular/mpr_numeric.cc @ 50cbdc

fieker-DuValspielwiese
Last change on this file since 50cbdc was 50cbdc, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes: merge-2-0-2 git-svn-id: file:///usr/local/Singular/svn/trunk@5619 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 28.4 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5/* $Id: mpr_numeric.cc,v 1.16 2001-08-27 14:47:15 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#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      20
268#define MAXIT   (MT*MR)   // max number of iterations in laguer root finder
269
270// set these values according to gmp_default_prec_bits and gmp_equalupto_bits!
271#define EPS     2.0e-34   // used by rootContainer::laguer_driver(), imag() == 0.0 ???
272//<-
273
274//-> rootContainer::rootContainer()
275rootContainer::rootContainer()
276{
277  rt=none;
278
279  coeffs= NULL;
280  ievpoint= NULL;
281  theroots= NULL;
282
283  found_roots= false;
284}
285//<-
286
287//-> rootContainer::~rootContainer()
288rootContainer::~rootContainer()
289{
290  int i;
291  int n= pVariables;
292
293  // free coeffs, ievpoint
294  if ( ievpoint != NULL )
295  {
296    for ( i=0; i < anz+2; i++ ) nDelete( ievpoint + i );
297    omFreeSize( (ADDRESS)ievpoint, (anz+2) * sizeof( number ) );
298  }
299
300  for ( i=0; i <= tdg; i++ ) nDelete( coeffs + i );
301  omFreeSize( (ADDRESS)coeffs, (tdg+1) * sizeof( number ) );
302
303  // theroots löschen
304  for ( i=0; i < tdg; i++ ) delete theroots[i];
305  omFreeSize( (ADDRESS) theroots, (tdg)*sizeof(gmp_complex*) );
306
307  mprPROTnl("~rootContainer()");
308}
309//<-
310
311//-> void rootContainer::fillContainer( ... )
312void rootContainer::fillContainer( number *_coeffs, number *_ievpoint,
313                                   const int _var, const int _tdg,
314                                   const rootType  _rt, const int _anz )
315{
316  int i;
317  number nn= nInit(0);
318  var=_var;
319  tdg=_tdg;
320  coeffs=_coeffs;
321  rt=_rt;
322  anz=_anz;
323
324  for ( i=0; i <= tdg; i++ )
325  {
326    if ( nEqual(coeffs[i],nn) )
327    {
328      nDelete( &coeffs[i] );
329      coeffs[i]=NULL;
330    }
331  }
332  nDelete( &nn );
333
334  if ( rt == cspecialmu && _ievpoint ) { // copy ievpoint
335    ievpoint= (number *)omAlloc( (anz+2) * sizeof( number ) );
336    for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] );
337  }
338
339  theroots= NULL;
340  found_roots= false;
341}
342//<-
343
344//-> poly rootContainer::getPoly()
345poly rootContainer::getPoly()
346{
347  int i;
348
349  poly result= NULL;
350  poly ppos;
351
352  if ( (rt == cspecial) || ( rt == cspecialmu ) )
353  {
354    for ( i= tdg; i >= 0; i-- )
355    {
356      if ( coeffs[i] )
357      {
358        poly p= pOne();
359        //pSetExp( p, var+1, i);
360        pSetExp( p, 1, i);
361        pSetCoeff( p, nCopy( coeffs[i] ) );
362        pSetm( p );
363        if (result)
364        {
365          ppos->next=p;
366          ppos=ppos->next;
367        }
368        else
369        {
370          result=p;
371          ppos=p;
372        }
373
374      }
375    }
376    pSetm( result );
377  }
378
379  return result;
380}
381//<-
382
383//-> const gmp_complex & rootContainer::opterator[] ( const int i )
384// this is now inline!
385//  gmp_complex & rootContainer::operator[] ( const int i )
386//  {
387//    if ( found_roots && ( i >= 0) && ( i < tdg ) )
388//    {
389//      return *theroots[i];
390//    }
391//    // warning
392//    Warn("rootContainer::getRoot: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
393//    gmp_complex *tmp= new gmp_complex();
394//    return *tmp;
395//  }
396//<-
397
398//-> gmp_complex & rootContainer::evPointCoord( int i )
399gmp_complex & rootContainer::evPointCoord( const int i )
400{
401  if (! ((i >= 0) && (i < anz+2) ) )
402    WarnS("rootContainer::evPointCoord: index out of range");
403  if (ievpoint == NULL)
404    WarnS("rootContainer::evPointCoord: ievpoint == NULL");
405
406  if ( (rt == cspecialmu) && found_roots ) {  // FIX ME
407    if ( ievpoint[i] != NULL )
408    {
409      gmp_complex *tmp= new gmp_complex();
410      *tmp= numberToComplex(ievpoint[i]);
411      return *tmp;
412    }
413    else
414    {
415      Warn("rootContainer::evPointCoord: NULL index %d",i);
416    }
417  }
418
419  // warning
420  Warn("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
421  gmp_complex *tmp= new gmp_complex();
422  return *tmp;
423}
424//<-
425
426//-> bool rootContainer::changeRoots( int from, int to )
427bool rootContainer::swapRoots( const int from, const int to )
428{
429  if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) )
430  {
431    if ( to != from )
432    {
433      gmp_complex tmp( *theroots[from] );
434      *theroots[from]= *theroots[to];
435      *theroots[to]= tmp;
436    }
437    return true;
438  }
439
440  // warning
441  Warn(" rootContainer::changeRoots: Wrong index %d, %d",from,to);
442  return false;
443}
444//<-
445
446//-> void rootContainer::solver()
447bool rootContainer::solver( const int polishmode )
448{
449  int i;
450
451  // there are maximal tdg roots, so *roots ranges form 0 to tdg-1.
452  theroots= (gmp_complex**)omAlloc( tdg*sizeof(gmp_complex*) );
453  for ( i=0; i < tdg; i++ ) theroots[i]= new gmp_complex();
454
455  // copy the coefficients of type number to type gmp_complex
456  gmp_complex **ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
457  for ( i=0; i <= tdg; i++ )
458  {
459    ad[i]= new gmp_complex();
460    if ( coeffs[i] ) *ad[i] = numberToComplex( coeffs[i] );
461  }
462
463  // now solve
464  switch (polishmode)
465  {
466  case PM_NONE:
467  case PM_POLISH:
468    found_roots= laguer_driver( ad, theroots, polishmode == PM_POLISH );
469    if (!found_roots)
470    {
471      WarnS("rootContainer::solver: No roots found!");
472      goto solverend;
473    }
474    break;
475  case PM_CORRUPT:
476    found_roots= laguer_driver( ad, theroots, false );
477    // corrupt the roots
478    for ( i= 0; i < tdg; i++ )
479      *theroots[i]= *theroots[i] * (gmp_float)(1.0+0.01*(mprfloat)i);
480    // and interpolate again
481    found_roots= laguer_driver( ad, theroots, true );
482    if (!found_roots)
483    {
484      WarnS("rootContainer::solver: No roots found!");
485      goto solverend;
486    }
487    break;
488  default:
489    Warn("rootContainer::solver: Unsupported polish mode %d! Valid are [0,1,2].",polishmode);
490    found_roots= false;
491  } // switch
492
493 solverend:
494  for ( i=0; i <= tdg; i++ ) delete ad[i];
495  omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof(gmp_complex*) );
496
497  return found_roots;
498}
499//<-
500
501//-> gmp_complex* rootContainer::laguer_driver( bool polish )
502bool rootContainer::laguer_driver(gmp_complex ** a, gmp_complex ** roots, bool polish )
503{
504  int i,j,jj;
505  int its;
506  gmp_complex x,b,c;
507  bool ret= true;
508
509  gmp_complex ** ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
510  for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] );
511
512  for ( j= tdg; j >= 1; j-- )
513  {
514    x= gmp_complex();
515
516    // run laguer alg
517    laguer(ad, j, &x, &its);
518
519    mprSTICKYPROT(ST_ROOTS_LGSTEP);
520    if ( its > MAXIT ) {  // error
521      WarnS("rootContainer::laguer_driver: To many iterations!");
522      ret= false;
523      goto theend;
524    }
525    if ( abs(x.imag()) <= (gmp_float)(2.0*EPS)*abs(x.real()))
526    {
527      x= gmp_complex( x.real() );
528    }
529    *roots[j-1]= x;
530    b= *ad[j];
531    for ( jj= j-1; jj >= 0; jj-- )
532    {
533      c= *ad[jj];
534      *ad[jj]= b;
535      b= ( x * b ) + c;
536    }
537  }
538
539  if ( polish )
540  {
541    mprSTICKYPROT(ST_ROOTS_LGPOLISH);
542    for ( i=0; i <= tdg; i++ ) *ad[i]=*a[i];
543
544    for ( j= 1; j <= tdg; j++ )
545    {
546      // run laguer alg with corrupted roots
547      laguer( ad, tdg, roots[j-1], &its );
548
549      mprSTICKYPROT(ST_ROOTS_LGSTEP);
550      if ( its > MAXIT )
551      {  // error
552        WarnS("rootContainer::laguer_driver: To many iterations!");
553        ret= false;
554        goto theend;
555      }
556    }
557    for ( j= 2; j <= tdg; j++ )
558    {
559      // sort root by their absolute real parts by straight insertion
560      x= *roots[j-1];
561      for ( i= j-1; i >= 1; i-- )
562      {
563        if ( abs(roots[i-1]->real()) <= abs(x.real()) ) break;
564        *roots[i]= *roots[i-1];
565      }
566      *roots[i]= x;
567    }
568  }
569
570 theend:
571  mprSTICKYPROT("\n");
572  for ( i=0; i <= tdg; i++ ) delete ad[i];
573  omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof( gmp_complex* ));
574
575  return ret;
576}
577//<-
578
579//-> void rootContainer::laguer(...)
580void rootContainer::laguer(gmp_complex ** a, int m, gmp_complex *x, int *its)
581{
582  int iter,j;
583  gmp_float abx_g, abp_g, abm_g, err_g;
584  gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2;
585  static gmp_float frac_g[MR+1] = { 0.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0 };
586
587  gmp_float epss(1.0/pow(10.0,(int)(gmp_output_digits+gmp_output_digits/4)));
588
589  for ( iter= 1; iter <= MAXIT; iter++ )
590  {
591    mprSTICKYPROT(ST_ROOTS_LG);
592
593    *its=iter;
594
595    b= *a[m];
596    err_g= abs(b);
597    d= gmp_complex();
598    f= gmp_complex();
599    abx_g= abs(*x);
600
601// gcc 2.95.2 on the dec alpha chokes on this
602#if defined(__GNUC__)
603#if ! (defined(__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 95)
604    for ( j= m-1; j >= 0; j-- )
605    {
606      f= ( *x * f ) + d;
607      d= ( *x * d ) + b;
608      b= ( *x * b ) + *a[j];
609      err_g= abs( b ) + ( abx_g * err_g );
610    }
611
612    err_g *= epss; // EPSS;
613
614    if ( abs(b) <= err_g ) return;
615
616    g= d / b;
617    g2 = g * g;
618    h= g2 - ( ( f / b ) * (gmp_float)2.0 );
619    sq= sqrt(( ( h * (gmp_float)m ) - g2 ) * (gmp_float)(m - 1));
620    gp= g + sq;
621    gm= g - sq;
622    abp_g= abs( gp );
623    abm_g= abs( gm );
624
625    if ( abp_g < abm_g ) gp= gm;
626
627    dx = ( (max(abp_g,abm_g) > (gmp_float)0.0)
628           ? ( gmp_complex( (mprfloat)m ) / gp )
629           : ( gmp_complex( cos((mprfloat)iter),sin((mprfloat)iter))
630               * exp(log((gmp_float)1.0+abx_g))) );
631
632    x1= *x - dx;
633
634    if ( (*x == x1) ) return;
635
636    if ( iter % MT ) *x= x1;
637    else *x -= ( dx * frac_g[ iter / MT ] );
638#endif
639#endif
640  }
641
642  *its= MAXIT+1;
643  return;
644}
645//<-
646
647//-----------------------------------------------------------------------------
648//-------------- rootArranger -------------------------------------------------
649//-----------------------------------------------------------------------------
650
651//-> rootArranger::rootArranger(...)
652rootArranger::rootArranger( rootContainer ** _roots,
653                            rootContainer ** _mu,
654                            const int _howclean )
655  : roots(_roots), mu(_mu), howclean(_howclean)
656{
657  found_roots=false;
658}
659//<-
660
661//-> void rootArranger::solve_all()
662void rootArranger::solve_all()
663{
664  int i;
665  found_roots= true;
666
667  // find roots of polys given by coeffs in roots
668  rc= roots[0]->getAnzElems();
669  for ( i= 0; i < rc; i++ )
670    if ( !roots[i]->solver( howclean ) )
671    {
672      found_roots= false;
673      return;
674    }
675  // find roots of polys given by coeffs in mu
676  mc= mu[0]->getAnzElems();
677  for ( i= 0; i < mc; i++ )
678    if ( ! mu[i]->solver( howclean ) )
679    {
680      found_roots= false;
681      return;
682    }
683}
684//<-
685
686//-> void rootArranger::arrange()
687void rootArranger::arrange()
688{
689  gmp_complex tmp,zwerg;
690  int anzm= mu[0]->getAnzElems();
691  int anzr= roots[0]->getAnzRoots();
692  int xkoord, r, rtest, xk, mtest;
693  bool found;
694  //gmp_complex mprec(1.0/pow(10,gmp_output_digits-5),1.0/pow(10,gmp_output_digits-5));
695  gmp_float mprec(1.0/pow(10.0,(int)(gmp_output_digits/3)));
696
697  for ( xkoord= 0; xkoord < anzm; xkoord++ ) {    // für x1,x2, x1,x2,x3, x1,x2,...,xn
698    for ( r= 0; r < anzr; r++ ) {                 // für jede Nullstelle
699      // (x1-koordinate) * evp[1] + (x2-koordinate) * evp[2] +
700      //                                  ... + (xkoord-koordinate) * evp[xkoord]
701      tmp= gmp_complex();
702      for ( xk =0; xk <= xkoord; xk++ )
703      {
704        tmp -= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); //xk+1
705      }
706      found= false;
707      for ( rtest= r; rtest < anzr; rtest++ ) {   // für jede Nullstelle
708        zwerg = tmp - (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xk+1, xkoord+2
709        for ( mtest= 0; mtest < anzr; mtest++ )
710        {
711          //          if ( tmp == (*mu[xkoord])[mtest] )
712          //          {
713          if ( ((zwerg.real() <= (*mu[xkoord])[mtest].real() + mprec) &&
714                (zwerg.real() >= (*mu[xkoord])[mtest].real() - mprec)) &&
715               ((zwerg.imag() <= (*mu[xkoord])[mtest].imag() + mprec) &&
716                (zwerg.imag() >= (*mu[xkoord])[mtest].imag() - mprec)) )
717           {
718             roots[xk]->swapRoots( r, rtest );
719             found= true;
720             break;
721           }
722        }
723      } // rtest
724      if ( !found )
725      {
726        Warn("rootArranger::arrange: No match? coord %d, root %d.",xkoord,r);
727#ifdef mprDEBUG_PROT
728        WarnS("One of these ...");
729        for ( rtest= r; rtest < anzr; rtest++ )
730        {
731          tmp= gmp_complex();
732          for ( xk =0; xk <= xkoord; xk++ )
733          {
734            tmp-= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1);
735          }
736          tmp-= (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xkoord+2
737          Warn("  %s",complexToStr(tmp,gmp_output_digits+1),rtest);
738        }
739        WarnS(" ... must match to one of these:");
740        for ( mtest= 0; mtest < anzr; mtest++ )
741        {
742          Warn("                  %s",complexToStr((*mu[xkoord])[mtest],gmp_output_digits+1));
743        }
744#endif
745      }
746    } // r
747  } // xkoord
748}
749//<-
750
751//-> lists rootArranger::listOfRoots( int oprec )
752lists rootArranger::listOfRoots( const unsigned int oprec )
753{
754  int i,j,tr;
755  int count= roots[0]->getAnzRoots(); // number of roots
756  int elem= roots[0]->getAnzElems();  // number of koordinates per root
757
758  lists listofroots= (lists)omAlloc( sizeof(slists) ); // must be done this way!
759
760  if ( found_roots )
761  {
762    listofroots->Init( count );
763
764    for (i=0; i < count; i++)
765    {
766      lists onepoint= (lists)omAlloc(sizeof(slists)); // must be done this way!
767      onepoint->Init(elem);
768      for ( j= 0; j < elem; j++ )
769      {
770        if ( !rField_is_long_C() )
771        {
772          onepoint->m[j].rtyp=STRING_CMD;
773          onepoint->m[j].data=(void *)complexToStr((*roots[j])[i],oprec);
774        }
775        else
776        {
777          onepoint->m[j].rtyp=NUMBER_CMD;
778          onepoint->m[j].data=(void *)nCopy((number)(roots[j]->getRoot(i)));
779        }
780        onepoint->m[j].next= NULL;
781        onepoint->m[j].name= NULL;
782      }
783      listofroots->m[i].rtyp=LIST_CMD;
784      listofroots->m[i].data=(void *)onepoint;
785      listofroots->m[j].next= NULL;
786      listofroots->m[j].name= NULL;
787    }
788
789  }
790  else
791  {
792    listofroots->Init( 0 );
793  }
794
795  return listofroots;
796}
797//<-
798
799//-----------------------------------------------------------------------------
800//-------------- simplex ----- ------------------------------------------------
801//-----------------------------------------------------------------------------
802
803//  #ifdef mprDEBUG_PROT
804//  #define error(a) a
805//  #else
806//  #define error(a)
807//  #endif
808
809#define error(a) a
810
811#define MAXPOINTS      1000
812
813//-> simplex::*
814//
815simplex::simplex( int rows, int cols )
816   : LiPM_cols(cols), LiPM_rows(rows)
817{
818  int i;
819
820  LiPM_rows=LiPM_rows+3;
821  LiPM_cols=LiPM_cols+2;
822
823  LiPM = (mprfloat **)omAlloc( LiPM_rows * sizeof(mprfloat *) );  // LP matrix
824  for( i= 0; i < LiPM_rows; i++ )
825  {
826    // Mem must be allocated aligned, also for type double!
827    LiPM[i] = (mprfloat *)omAlloc0Aligned( LiPM_cols * sizeof(mprfloat) );
828  }
829
830  iposv = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
831  izrov = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
832
833  m=n=m1=m2=m3=icase=0;
834
835#ifdef mprDEBUG_ALL
836  Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols);
837#endif
838}
839
840simplex::~simplex()
841{
842  // clean up
843  int i;
844  for( i= 0; i < LiPM_rows; i++ )
845  {
846    omFreeSize( (ADDRESS) LiPM[i], LiPM_cols * sizeof(mprfloat) );
847  }
848  omFreeSize( (ADDRESS) LiPM, LiPM_rows * sizeof(mprfloat *) );
849
850  omFreeSize( (ADDRESS) iposv, 2*LiPM_rows*sizeof(int) );
851  omFreeSize( (ADDRESS) izrov, 2*LiPM_rows*sizeof(int) );
852}
853
854BOOLEAN simplex::mapFromMatrix( matrix m )
855{
856  int i,j;
857//    if ( MATROWS( m ) > LiPM_rows ||  MATCOLS( m ) > LiPM_cols ) {
858//      WarnS("");
859//      return FALSE;
860//    }
861
862  number coef;
863  for ( i= 1; i <= MATROWS( m ); i++ )
864  {
865     for ( j= 1; j <= MATCOLS( m ); j++ )
866     {
867        if ( MATELEM(m,i,j) != NULL )
868        {
869           coef= pGetCoeff( MATELEM(m,i,j) );
870           if ( coef != NULL && !nIsZero(coef) )
871              LiPM[i][j]= (double)(*(gmp_float*)coef);
872           //#ifdef mpr_DEBUG_PROT
873           //Print("%f ",LiPM[i][j]);
874           //#endif
875        }
876     }
877     //     PrintLn();
878  }
879
880  return TRUE;
881}
882
883matrix simplex::mapToMatrix( matrix m )
884{
885  int i,j;
886//    if ( MATROWS( m ) < LiPM_rows-3 ||  MATCOLS( m ) < LiPM_cols-2 ) {
887//      WarnS("");
888//      return NULL;
889//    }
890
891//Print(" %d x %d\n",MATROWS( m ),MATCOLS( m ));
892
893  number coef;
894  gmp_float * bla;
895  for ( i= 1; i <= MATROWS( m ); i++ )
896  {
897    for ( j= 1; j <= MATCOLS( m ); j++ )
898    {
899       pDelete( &(MATELEM(m,i,j)) );
900       MATELEM(m,i,j)= NULL;
901//Print(" %3.0f ",LiPM[i][j]);
902       if ( LiPM[i][j] != 0.0 )
903       {
904          bla= new gmp_float(LiPM[i][j]);
905          coef= (number)bla;
906          MATELEM(m,i,j)= pOne();
907          pSetCoeff( MATELEM(m,i,j), coef );
908       }
909    }
910//PrintLn();
911  }
912
913  return m;
914}
915
916intvec * simplex::posvToIV()
917{
918   int i;
919   intvec * iv = new intvec( m );
920   for ( i= 1; i <= m; i++ )
921   {
922      IMATELEM(*iv,i,1)= iposv[i];
923   }
924   return iv;
925}
926
927intvec * simplex::zrovToIV()
928{
929   int i;
930   intvec * iv = new intvec( n );
931   for ( i= 1; i <= n; i++ )
932   {
933      IMATELEM(*iv,i,1)= izrov[i];
934   }
935   return iv;
936}
937
938void simplex::compute()
939{
940  int i,ip,ir,is,k,kh,kp,m12,nl1,nl2;
941  int *l1,*l2,*l3;
942  mprfloat q1, bmax;
943
944  if ( m != (m1+m2+m3) )
945  {
946    // error: bad input
947    error(WarnS("simplex::compute: Bad input constraint counts!");)
948    icase=-2;
949    return;
950  }
951
952  l1= (int *) omAlloc0( (n+1) * sizeof(int) );
953  l2= (int *) omAlloc0( (m+1) * sizeof(int) );
954  l3= (int *) omAlloc0( (m+1) * sizeof(int) );
955
956  nl1= n;
957  for ( k=1; k<=n; k++ ) l1[k]=izrov[k]=k;
958  nl2=m;
959  for ( i=1; i<=m; i++ )
960  {
961    if ( LiPM[i+1][1] < 0.0 )
962    {
963      // error: bad input
964      error(WarnS("simplex::compute: Bad input tableau!");)
965      error(Warn("simplex::compute: in input Matrix row %d, column 1, value %f",i+1,LiPM[i+1][1]);)
966      icase=-2;
967      // free mem l1,l2,l3;
968      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
969      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
970      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
971      return;
972    }
973    l2[i]= i;
974    iposv[i]= n+i;
975  }
976  for ( i=1; i<=m2; i++) l3[i]= 1;
977  ir= 0;
978  if (m2+m3)
979  {
980    ir=1;
981    for ( k=1; k <= (n+1); k++ )
982    {
983      q1=0.0;
984      for ( i=m1+1; i <= m; i++ ) q1+= LiPM[i+1][k];
985      LiPM[m+2][k]= -q1;
986    }
987
988    do
989    {
990      simp1(LiPM,m+1,l1,nl1,0,&kp,&bmax);
991      if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] < -SIMPLEX_EPS )
992      {
993        icase= -1; // no solution found
994        // free mem l1,l2,l3;
995        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
996        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
997        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
998        return;
999      }
1000      else if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] <= SIMPLEX_EPS )
1001      {
1002        m12= m1+m2+1;
1003        if ( m12 <= m )
1004        {
1005          for ( ip= m12; ip <= m; ip++ )
1006          {
1007            if ( iposv[ip] == (ip+n) )
1008            {
1009              simp1(LiPM,ip,l1,nl1,1,&kp,&bmax);
1010              if ( fabs(bmax) >= SIMPLEX_EPS)
1011                goto one;
1012            }
1013          }
1014        }
1015        ir= 0;
1016        --m12;
1017        if ( m1+1 <= m12 )
1018          for ( i=m1+1; i <= m12; i++ )
1019            if ( l3[i-m1] == 1 )
1020              for ( k=1; k <= n+1; k++ )
1021                LiPM[i+1][k] = -(LiPM[i+1][k]);
1022        break;
1023      }
1024      //#if DEBUG
1025      //print_bmat( a, m+2, n+3);
1026      //#endif
1027      simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1028      if ( ip == 0 )
1029      {
1030        icase = -1; // no solution found
1031        // free mem l1,l2,l3;
1032        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1033        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1034        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1035        return;
1036      }
1037    one: simp3(LiPM,m+1,n,ip,kp);
1038      // #if DEBUG
1039      // print_bmat(a,m+2,n+3);
1040      // #endif
1041      if ( iposv[ip] >= (n+m1+m2+1))
1042      {
1043        for ( k= 1; k <= nl1; k++ )
1044          if ( l1[k] == kp ) break;
1045        --nl1;
1046        for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1];
1047        ++(LiPM[m+2][kp+1]);
1048        for ( i= 1; i <= m+2; i++ ) LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1049      }
1050      else
1051      {
1052        if ( iposv[ip] >= (n+m1+1) )
1053        {
1054          kh= iposv[ip]-m1-n;
1055          if ( l3[kh] )
1056          {
1057            l3[kh]= 0;
1058            ++(LiPM[m+2][kp+1]);
1059            for ( i=1; i<= m+2; i++ )
1060              LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1061          }
1062        }
1063      }
1064      is= izrov[kp];
1065      izrov[kp]= iposv[ip];
1066      iposv[ip]= is;
1067    } while (ir);
1068  }
1069  /* end of phase 1, have feasible sol, now optimize it */
1070  loop
1071  {
1072    // #if DEBUG
1073    // print_bmat( a, m+1, n+5);
1074    // #endif
1075    simp1(LiPM,0,l1,nl1,0,&kp,&bmax);
1076    if (bmax <= /*SIMPLEX_EPS*/0.0)
1077    {
1078      icase=0; // finite solution found
1079      // free mem l1,l2,l3
1080      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1081      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1082      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1083      return;
1084    }
1085    simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1086    if (ip == 0)
1087    {
1088      //printf("Unbounded:");
1089      // #if DEBUG
1090      //       print_bmat( a, m+1, n+1);
1091      // #endif
1092      icase=1;                /* unbounded */
1093      // free mem
1094      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1095      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1096      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1097      return;
1098    }
1099    simp3(LiPM,m,n,ip,kp);
1100    is= izrov[kp];
1101    izrov[kp]= iposv[ip];
1102    iposv[ip]= is;
1103  }/*for ;;*/
1104}
1105
1106void simplex::simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax )
1107{
1108  int k;
1109  mprfloat test;
1110
1111  if( nll <= 0)
1112  {                        /* init'tion: fixed */
1113    *bmax = 0.0;
1114    return;
1115  }
1116  *kp=ll[1];
1117  *bmax=a[mm+1][*kp+1];
1118  for (k=2;k<=nll;k++)
1119  {
1120    if (iabf == 0)
1121    {
1122      test=a[mm+1][ll[k]+1]-(*bmax);
1123      if (test > 0.0)
1124      {
1125        *bmax=a[mm+1][ll[k]+1];
1126        *kp=ll[k];
1127      }
1128    }
1129    else
1130    {                        /* abs values: have fixed it */
1131      test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
1132      if (test > 0.0)
1133      {
1134        *bmax=a[mm+1][ll[k]+1];
1135        *kp=ll[k];
1136      }
1137    }
1138  }
1139}
1140
1141void simplex::simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 )
1142{
1143  int k,ii,i;
1144  mprfloat qp,q0,q;
1145
1146  *ip= 0;
1147  for ( i=1; i <= nl2; i++ )
1148  {
1149    if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS )
1150    {
1151      *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1];
1152      *ip= l2[i];
1153      for ( i= i+1; i <= nl2; i++ )
1154      {
1155        ii= l2[i];
1156        if (a[ii+1][kp+1] < -SIMPLEX_EPS)
1157        {
1158          q= -a[ii+1][1] / a[ii+1][kp+1];
1159          if (q - *q1 < -SIMPLEX_EPS)
1160          {
1161            *ip=ii;
1162            *q1=q;
1163          }
1164          else if (q - *q1 < SIMPLEX_EPS)
1165          {
1166            for ( k=1; k<= n; k++ )
1167            {
1168              qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
1169              q0= -a[ii+1][k+1]/a[ii+1][kp+1];
1170              if ( q0 != qp ) break;
1171            }
1172            if ( q0 < qp ) *ip= ii;
1173          }
1174        }
1175      }
1176    }
1177  }
1178}
1179
1180void simplex::simp3( mprfloat **a, int i1, int k1, int ip, int kp )
1181{
1182  int kk,ii;
1183  mprfloat piv;
1184
1185  piv= 1.0 / a[ip+1][kp+1];
1186  for ( ii=1; ii <= i1+1; ii++ )
1187  {
1188    if ( ii -1 != ip )
1189    {
1190      a[ii][kp+1] *= piv;
1191      for ( kk=1; kk <= k1+1; kk++ )
1192        if ( kk-1 != kp )
1193          a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
1194    }
1195  }
1196  for ( kk=1; kk<= k1+1; kk++ )
1197    if ( kk-1 != kp ) a[ip+1][kk] *= -piv;
1198  a[ip+1][kp+1]= piv;
1199}
1200//<-
1201
1202//-----------------------------------------------------------------------------
1203
1204//#endif // HAVE_MPR
1205
1206// local Variables: ***
1207// folded-file: t ***
1208// compile-command-1: "make installg" ***
1209// compile-command-2: "make install" ***
1210// End: ***
Note: See TracBrowser for help on using the repository browser.