source: git/Singular/mpr_numeric.cc @ b39d4d

spielwiese
Last change on this file since b39d4d was b39d4d, checked in by Olaf Bachmann <obachman@…>, 23 years ago
Dec alpha port git-svn-id: file:///usr/local/Singular/svn/trunk@4933 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 28.3 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5/* $Id: mpr_numeric.cc,v 1.12 2000-12-18 15:44:42 obachman 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 "longalg.h"
26#include "matpol.h"
27#include "ring.h"
28//#include "longrat.h"
29#include "lists.h"
30
31#include <math.h>
32#include "mpr_numeric.h"
33
34extern size_t gmp_output_digits;
35//<-
36
37extern void nPrint(number n);  // for debugging output
38
39//-----------------------------------------------------------------------------
40//-------------- vandermonde system solver ------------------------------------
41//-----------------------------------------------------------------------------
42
43//-> vandermonde::*
44vandermonde::vandermonde( const long _cn, const long _n, const long _maxdeg,
45                          number *_p, const bool _homog )
46  : n(_n), cn(_cn), maxdeg(_maxdeg), p(_p), homog(_homog)
47{
48  long j;
49  l= (long)pow((double)maxdeg+1,(int)n);
50  x= (number *)omAlloc( cn * sizeof(number) );
51  for ( j= 0; j < cn; j++ ) x[j]= nInit(1);
52  init();
53}
54
55vandermonde::~vandermonde()
56{
57  int j;
58  for ( j= 0; j < cn; j++ ) nDelete( x+j );
59  omFreeSize( (ADDRESS)x, cn * sizeof( number ) );
60}
61
62void vandermonde::init()
63{
64  int j;
65  long i,c,sum;
66  number tmp,tmp1;
67
68  c=0;
69  sum=0;
70
71  intvec exp( n );
72  for ( j= 0; j < n; j++ ) exp[j]=0;
73
74  for ( i= 0; i < l; i++ )
75  {
76    if ( !homog || (sum == maxdeg) )
77    {
78      for ( j= 0; j < n; j++ )
79      {
80        nPower( p[j], exp[j], &tmp );
81        tmp1 = nMult( tmp, x[c] );
82        x[c]= tmp1;
83        nDelete( &tmp );
84      }
85      c++;
86    }
87    exp[0]++;
88    sum=0;
89    for ( j= 0; j < n - 1; j++ )
90    {
91      if ( exp[j] > maxdeg )
92      {
93        exp[j]= 0;
94        exp[j + 1]++;
95        }
96      sum+= exp[j];
97    }
98    sum+= exp[n - 1];
99  }
100}
101
102poly vandermonde::numvec2poly( const number * q )
103{
104  int j;
105  long i,c,sum;
106  number tmp;
107
108  poly pnew,pit=NULL;
109
110  c=0;
111  sum=0;
112
113  Exponent_t *exp= (Exponent_t *) omAlloc( (n+1) * sizeof(Exponent_t) );
114
115  for ( j= 0; j < n+1; j++ ) exp[j]=0;
116
117  for ( i= 0; i < l; i++ )
118  {
119    if ( (!homog || (sum == maxdeg)) && q[i] && !nIsZero(q[i]) )
120    {
121      pnew= pOne();
122      pSetCoeff(pnew,q[i]);
123      pSetExpV(pnew,exp);
124      if ( pit )
125      {
126        pNext(pnew)= pit;
127        pit= pnew;
128      }
129      else
130      {
131        pit= pnew;
132        pNext(pnew)= NULL;
133      }
134      pSetm(pit);
135    }
136    exp[1]++;
137    sum=0;
138    for ( j= 1; j < n; j++ )
139    {
140      if ( exp[j] > maxdeg )
141      {
142        exp[j]= 0;
143        exp[j + 1]++;
144        }
145      sum+= exp[j];
146    }
147    sum+= exp[n];
148  }
149
150  omFreeSize( (ADDRESS) exp, (n+1) * sizeof(Exponent_t) );
151
152  pOrdPolyMerge(pit);
153  return pit;
154}
155
156number * vandermonde::interpolateDense( const number * q )
157{
158  int i,j,k;
159  number newnum,tmp1;
160  number b,t,xx,s;
161  number *c;
162  number *w;
163
164  b=t=xx=s=tmp1=NULL;
165
166  w= (number *)omAlloc( cn * sizeof(number) );
167  c= (number *)omAlloc( cn * sizeof(number) );
168  for ( j= 0; j < cn; j++ )
169  {
170    w[j]= nInit(0);
171    c[j]= nInit(0);
172  }
173
174  if ( cn == 1 )
175  {
176    nDelete( &w[0] );
177    w[0]= nCopy(q[0]);
178  }
179  else
180  {
181    nDelete( &c[cn-1] );
182    c[cn-1]= nCopy(x[0]);
183    c[cn-1]= nNeg(c[cn-1]);              // c[cn]= -x[1]
184
185    for ( i= 1; i < cn; i++ ) {              // i=2; i <= cn
186      nDelete( &xx );
187      xx= nCopy(x[i]);
188      xx= nNeg(xx);               // xx= -x[i]
189
190      for ( j= (cn-i-1); j <= (cn-2); j++) { // j=(cn+1-i); j <= (cn-1)
191        nDelete( &tmp1 );
192        tmp1= nMult( xx, c[j+1] );           // c[j]= c[j] + (xx * c[j+1])
193        newnum= nAdd( c[j], tmp1 );
194        nDelete( &c[j] );
195        c[j]= newnum;
196      }
197
198      newnum= nAdd( xx, c[cn-1] );           // c[cn-1]= c[cn-1] + xx
199      nDelete( &c[cn-1] );
200      c[cn-1]= newnum;
201    }
202
203    for ( i= 0; i < cn; i++ ) {              // i=1; i <= cn
204      nDelete( &xx );
205      xx= nCopy(x[i]);                     // xx= x[i]
206
207      nDelete( &t );
208      t= nInit( 1 );                         // t= b= 1
209      nDelete( &b );
210      b= nInit( 1 );
211      nDelete( &s );                         // s= q[cn-1]
212      s= nCopy( q[cn-1] );
213
214      for ( k= cn-1; k >= 1; k-- ) {         // k=cn; k >= 2
215        nDelete( &tmp1 );
216        tmp1= nMult( xx, b );                // b= c[k] + (xx * b)
217        nDelete( &b );
218        b= nAdd( c[k], tmp1 );
219
220        nDelete( &tmp1 );
221        tmp1= nMult( q[k-1], b );            // s= s + (q[k-1] * b)
222        newnum= nAdd( s, tmp1 );
223        nDelete( &s );
224        s= newnum;
225
226        nDelete( &tmp1 );
227        tmp1= nMult( xx, t );                // t= (t * xx) + b
228        newnum= nAdd( tmp1, b );
229        nDelete( &t );
230        t= newnum;
231      }
232
233      nDelete( &w[i] );                      // w[i]= s/t
234      w[i]= nDiv( s, t );
235      nNormalize( w[i] );
236
237      mprSTICKYPROT(ST_VANDER_STEP);
238    }
239  }
240  mprSTICKYPROT("\n");
241
242  // free mem
243  for ( j= 0; j < cn; j++ ) nDelete( c+j );
244  omFreeSize( (ADDRESS)c, cn * sizeof( number ) );
245
246  nDelete( &tmp1 );
247  nDelete( &s );
248  nDelete( &t );
249  nDelete( &b );
250  nDelete( &xx );
251
252  // makes quotiens smaller
253  for ( j= 0; j < cn; j++ ) nNormalize( w[j] );
254
255  return w;
256}
257//<-
258
259//-----------------------------------------------------------------------------
260//-------------- rootContainer ------------------------------------------------
261//-----------------------------------------------------------------------------
262
263//-> definitions
264#define MR       8        // never change this value
265#define MT      20
266#define MAXIT   (MT*MR)   // max number of iterations in laguer root finder
267
268// set these values according to gmp_default_prec_bits and gmp_equalupto_bits!
269#define EPS     2.0e-34   // used by rootContainer::laguer_driver(), imag() == 0.0 ???
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  switch (polishmode)
463  {
464  case PM_NONE:
465  case PM_POLISH:
466    found_roots= laguer_driver( ad, theroots, polishmode == PM_POLISH );
467    if (!found_roots)
468    {
469      WarnS("rootContainer::solver: No roots found!");
470      goto solverend;
471    }
472    break;
473  case PM_CORRUPT:
474    found_roots= laguer_driver( ad, theroots, false );
475    // corrupt the roots
476    for ( i= 0; i < tdg; i++ )
477      *theroots[i]= *theroots[i] * (gmp_float)(1.0+0.01*(mprfloat)i);
478    // and interpolate again
479    found_roots= laguer_driver( ad, theroots, true );
480    if (!found_roots)
481    {
482      WarnS("rootContainer::solver: No roots found!");
483      goto solverend;
484    }
485    break;
486  default:
487    Warn("rootContainer::solver: Unsupported polish mode %d! Valid are [0,1,2].",polishmode);
488    found_roots= false;
489  } // switch
490
491 solverend:
492  for ( i=0; i <= tdg; i++ ) delete ad[i];
493  omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof(gmp_complex*) );
494
495  return found_roots;
496}
497//<-
498
499//-> gmp_complex* rootContainer::laguer_driver( bool polish )
500bool rootContainer::laguer_driver(gmp_complex ** a, gmp_complex ** roots, bool polish )
501{
502  int i,j,jj;
503  int its;
504  gmp_complex x,b,c;
505  bool ret= true;
506
507  gmp_complex ** ad= (gmp_complex**)omAlloc( (tdg+1)*sizeof(gmp_complex*) );
508  for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] );
509
510  for ( j= tdg; j >= 1; j-- )
511  {
512    x= gmp_complex();
513
514    // run laguer alg
515    laguer(ad, j, &x, &its);
516
517    mprSTICKYPROT(ST_ROOTS_LGSTEP);
518    if ( its > MAXIT ) {  // error
519      WarnS("rootContainer::laguer_driver: To many iterations!");
520      ret= false;
521      goto theend;
522    }
523    if ( abs(x.imag()) <= (gmp_float)(2.0*EPS)*abs(x.real()))
524    {
525      x= gmp_complex( x.real() );
526    }
527    *roots[j-1]= x;
528    b= *ad[j];
529    for ( jj= j-1; jj >= 0; jj-- )
530    {
531      c= *ad[jj];
532      *ad[jj]= b;
533      b= ( x * b ) + c;
534    }
535  }
536
537  if ( polish )
538  {
539    mprSTICKYPROT(ST_ROOTS_LGPOLISH);
540    for ( i=0; i <= tdg; i++ ) *ad[i]=*a[i];
541
542    for ( j= 1; j <= tdg; j++ )
543    {
544      // run laguer alg with corrupted roots
545      laguer( ad, tdg, roots[j-1], &its );
546
547      mprSTICKYPROT(ST_ROOTS_LGSTEP);
548      if ( its > MAXIT )
549      {  // error
550        WarnS("rootContainer::laguer_driver: To many iterations!");
551        ret= false;
552        goto theend;
553      }
554    }
555    for ( j= 2; j <= tdg; j++ )
556    {
557      // sort root by their absolute real parts by straight insertion
558      x= *roots[j-1];
559      for ( i= j-1; i >= 1; i-- )
560      {
561        if ( abs(roots[i-1]->real()) <= abs(x.real()) ) break;
562        *roots[i]= *roots[i-1];
563      }
564      *roots[i]= x;
565    }
566  }
567
568 theend:
569  mprSTICKYPROT("\n");
570  for ( i=0; i <= tdg; i++ ) delete ad[i];
571  omFreeSize( (ADDRESS) ad, (tdg+1)*sizeof( gmp_complex* ));
572
573  return ret;
574}
575//<-
576
577//-> void rootContainer::laguer(...)
578void rootContainer::laguer(gmp_complex ** a, int m, gmp_complex *x, int *its)
579{
580  int iter,j;
581  gmp_float abx_g, abp_g, abm_g, err_g;
582  gmp_complex dx, x1, b, d, f, g, h, sq, gp, gm, g2;
583  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 };
584
585  gmp_float epss(1.0/pow(10.0,(int)(gmp_output_digits+gmp_output_digits/4)));
586
587  for ( iter= 1; iter <= MAXIT; iter++ )
588  {
589    mprSTICKYPROT(ST_ROOTS_LG);
590
591    *its=iter;
592
593    b= *a[m];
594    err_g= abs(b);
595    d= gmp_complex();
596    f= gmp_complex();
597    abx_g= abs(*x);
598
599// gcc 2.95.2 on the dec alpha chokes on this
600#if defined(__GNUC__)
601#if ! (defined(__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 95)
602    for ( j= m-1; j >= 0; j-- )
603    {
604      f= ( *x * f ) + d;
605      d= ( *x * d ) + b;
606      b= ( *x * b ) + *a[j];
607      err_g= abs( b ) + ( abx_g * err_g );
608    }
609
610    err_g *= epss; // EPSS;
611
612    if ( abs(b) <= err_g ) return;
613
614    g= d / b;
615    g2 = g * g;
616    h= g2 - ( ( f / b ) * (gmp_float)2.0 );
617    sq= sqrt(( ( h * (gmp_float)m ) - g2 ) * (gmp_float)(m - 1));
618    gp= g + sq;
619    gm= g - sq;
620    abp_g= abs( gp );
621    abm_g= abs( gm );
622
623    if ( abp_g < abm_g ) gp= gm;
624
625    dx = ( (max(abp_g,abm_g) > (gmp_float)0.0)
626           ? ( gmp_complex( (mprfloat)m ) / gp )
627           : ( gmp_complex( cos((mprfloat)iter),sin((mprfloat)iter))
628               * exp(log((gmp_float)1.0+abx_g))) );
629
630    x1= *x - dx;
631
632    if ( (*x == x1) ) return;
633
634    if ( iter % MT ) *x= x1;
635    else *x -= ( dx * frac_g[ iter / MT ] );
636#endif
637#endif
638  }
639
640  *its= MAXIT+1;
641  return;
642}
643//<-
644
645//-----------------------------------------------------------------------------
646//-------------- rootArranger -------------------------------------------------
647//-----------------------------------------------------------------------------
648
649//-> rootArranger::rootArranger(...)
650rootArranger::rootArranger( rootContainer ** _roots,
651                            rootContainer ** _mu,
652                            const int _howclean )
653  : roots(_roots), mu(_mu), howclean(_howclean)
654{
655  found_roots=false;
656}
657//<-
658
659//-> void rootArranger::solve_all()
660void rootArranger::solve_all()
661{
662  int i;
663  found_roots= true;
664
665  // find roots of polys given by coeffs in roots
666  rc= roots[0]->getAnzElems();
667  for ( i= 0; i < rc; i++ )
668    if ( !roots[i]->solver( howclean ) )
669    {
670      found_roots= false;
671      return;
672    }
673  // find roots of polys given by coeffs in mu
674  mc= mu[0]->getAnzElems();
675  for ( i= 0; i < mc; i++ )
676    if ( ! mu[i]->solver( howclean ) )
677    {
678      found_roots= false;
679      return;
680    }
681}
682//<-
683
684//-> void rootArranger::arrange()
685void rootArranger::arrange()
686{
687  gmp_complex tmp,zwerg;
688  int anzm= mu[0]->getAnzElems();
689  int anzr= roots[0]->getAnzRoots();
690  int xkoord, r, rtest, xk, mtest;
691  bool found;
692  //gmp_complex mprec(1.0/pow(10,gmp_output_digits-5),1.0/pow(10,gmp_output_digits-5));
693  gmp_float mprec(1.0/pow(10.0,(int)(gmp_output_digits/3)));
694
695  for ( xkoord= 0; xkoord < anzm; xkoord++ ) {    // für x1,x2, x1,x2,x3, x1,x2,...,xn
696    for ( r= 0; r < anzr; r++ ) {                 // für jede Nullstelle
697      // (x1-koordinate) * evp[1] + (x2-koordinate) * evp[2] +
698      //                                  ... + (xkoord-koordinate) * evp[xkoord]
699      tmp= gmp_complex();
700      for ( xk =0; xk <= xkoord; xk++ )
701      {
702        tmp -= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); //xk+1
703      }
704      found= false;
705      for ( rtest= r; rtest < anzr; rtest++ ) {   // für jede Nullstelle
706        zwerg = tmp - (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xk+1, xkoord+2
707        for ( mtest= 0; mtest < anzr; mtest++ )
708        {
709          //          if ( tmp == (*mu[xkoord])[mtest] )
710          //          {
711          if ( ((zwerg.real() <= (*mu[xkoord])[mtest].real() + mprec) &&
712                (zwerg.real() >= (*mu[xkoord])[mtest].real() - mprec)) &&
713               ((zwerg.imag() <= (*mu[xkoord])[mtest].imag() + mprec) &&
714                (zwerg.imag() >= (*mu[xkoord])[mtest].imag() - mprec)) )
715           {
716             roots[xk]->swapRoots( r, rtest );
717             found= true;
718             break;
719           }
720        }
721      } // rtest
722      if ( !found )
723      {
724        Warn("rootArranger::arrange: No match? coord %d, root %d.",xkoord,r);
725#ifdef mprDEBUG_PROT
726        WarnS("One of these ...");
727        for ( rtest= r; rtest < anzr; rtest++ )
728        {
729          tmp= gmp_complex();
730          for ( xk =0; xk <= xkoord; xk++ )
731          {
732            tmp-= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1);
733          }
734          tmp-= (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xkoord+2
735          Warn("  %s",complexToStr(tmp,gmp_output_digits+1),rtest);
736        }
737        WarnS(" ... must match to one of these:");
738        for ( mtest= 0; mtest < anzr; mtest++ )
739        {
740          Warn("                  %s",complexToStr((*mu[xkoord])[mtest],gmp_output_digits+1));
741        }
742#endif
743      }
744    } // r
745  } // xkoord
746}
747//<-
748
749//-> lists rootArranger::listOfRoots( int oprec )
750lists rootArranger::listOfRoots( const unsigned int oprec )
751{
752  int i,j,tr;
753  int count= roots[0]->getAnzRoots(); // number of roots
754  int elem= roots[0]->getAnzElems();  // number of koordinates per root
755
756  lists listofroots= (lists)omAlloc( sizeof(slists) ); // must be done this way!
757
758  if ( found_roots )
759  {
760    listofroots->Init( count );
761
762    for (i=0; i < count; i++)
763    {
764      lists onepoint= (lists)omAlloc(sizeof(slists)); // must be done this way!
765      onepoint->Init(elem);
766      for ( j= 0; j < elem; j++ )
767      {
768        if ( !rField_is_long_C() )
769        {
770          onepoint->m[j].rtyp=STRING_CMD;
771          onepoint->m[j].data=(void *)complexToStr((*roots[j])[i],oprec);
772        }
773        else
774        {
775          onepoint->m[j].rtyp=NUMBER_CMD;
776          onepoint->m[j].data=(void *)nCopy((number)(roots[j]->getRoot(i)));
777        }
778        onepoint->m[j].next= NULL;
779        onepoint->m[j].name= NULL;
780      }
781      listofroots->m[i].rtyp=LIST_CMD;
782      listofroots->m[i].data=(void *)onepoint;
783      listofroots->m[j].next= NULL;
784      listofroots->m[j].name= NULL;
785    }
786
787  }
788  else
789  {
790    listofroots->Init( 0 );
791  }
792
793  return listofroots;
794}
795//<-
796
797//-----------------------------------------------------------------------------
798//-------------- simplex ----- ------------------------------------------------
799//-----------------------------------------------------------------------------
800
801//  #ifdef mprDEBUG_PROT
802//  #define error(a) a
803//  #else
804//  #define error(a)
805//  #endif
806
807#define error(a) a
808
809#define MAXPOINTS      1000
810
811//-> simplex::*
812//
813simplex::simplex( int rows, int cols )
814   : LiPM_cols(cols), LiPM_rows(rows)
815{
816  int i;
817
818  LiPM_rows=LiPM_rows+3;
819  LiPM_cols=LiPM_cols+2;
820
821  LiPM = (mprfloat **)omAlloc( LiPM_rows * sizeof(mprfloat *) );  // LP matrix
822  for( i= 0; i < LiPM_rows; i++ )
823  {
824    // Mem must be allocated aligned, also for type double!
825    LiPM[i] = (mprfloat *)omAlloc0Aligned( LiPM_cols * sizeof(mprfloat) );
826  }
827
828  iposv = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
829  izrov = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
830
831  m=n=m1=m2=m3=icase=0;
832
833#ifdef mprDEBUG_ALL
834  Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols);
835#endif
836}
837
838simplex::~simplex()
839{
840  // clean up
841  int i;
842  for( i= 0; i < LiPM_rows; i++ )
843  {
844    omFreeSize( (ADDRESS) LiPM[i], LiPM_cols * sizeof(mprfloat) );
845  }
846  omFreeSize( (ADDRESS) LiPM, LiPM_rows * sizeof(mprfloat *) );
847
848  omFreeSize( (ADDRESS) iposv, 2*LiPM_rows*sizeof(int) );
849  omFreeSize( (ADDRESS) izrov, 2*LiPM_rows*sizeof(int) );
850}
851
852BOOLEAN simplex::mapFromMatrix( matrix m )
853{
854  int i,j;
855//    if ( MATROWS( m ) > LiPM_rows ||  MATCOLS( m ) > LiPM_cols ) {
856//      WarnS("");
857//      return FALSE;
858//    }
859
860  number coef;
861  for ( i= 1; i <= MATROWS( m ); i++ )
862  {
863     for ( j= 1; j <= MATCOLS( m ); j++ )
864     {
865        if ( MATELEM(m,i,j) != NULL )
866        {
867           coef= pGetCoeff( MATELEM(m,i,j) );
868           if ( coef != NULL && !nIsZero(coef) )
869              LiPM[i][j]= (double)(*(gmp_float*)coef);
870           //#ifdef mpr_DEBUG_PROT
871           //Print("%f ",LiPM[i][j]);
872           //#endif
873        }
874     }
875     //     PrintLn();
876  }
877
878  return TRUE;
879}
880
881matrix simplex::mapToMatrix( matrix m )
882{
883  int i,j;
884//    if ( MATROWS( m ) < LiPM_rows-3 ||  MATCOLS( m ) < LiPM_cols-2 ) {
885//      WarnS("");
886//      return NULL;
887//    }
888
889//Print(" %d x %d\n",MATROWS( m ),MATCOLS( m ));
890
891  number coef;
892  gmp_float * bla;
893  for ( i= 1; i <= MATROWS( m ); i++ )
894  {
895    for ( j= 1; j <= MATCOLS( m ); j++ )
896    {
897       pDelete( &(MATELEM(m,i,j)) );
898       MATELEM(m,i,j)= NULL;
899//Print(" %3.0f ",LiPM[i][j]);
900       if ( LiPM[i][j] != 0.0 )
901       {
902          bla= new gmp_float(LiPM[i][j]);
903          coef= (number)bla;
904          MATELEM(m,i,j)= pOne();
905          pSetCoeff( MATELEM(m,i,j), coef );
906       }
907    }
908//PrintLn();
909  }
910
911  return m;
912}
913
914intvec * simplex::posvToIV()
915{
916   int i;
917   intvec * iv = new intvec( m );
918   for ( i= 1; i <= m; i++ )
919   {
920      IMATELEM(*iv,i,1)= iposv[i];
921   }
922   return iv;
923}
924
925intvec * simplex::zrovToIV()
926{
927   int i;
928   intvec * iv = new intvec( n );
929   for ( i= 1; i <= n; i++ )
930   {
931      IMATELEM(*iv,i,1)= izrov[i];
932   }
933   return iv;
934}
935
936void simplex::compute()
937{
938  int i,ip,ir,is,k,kh,kp,m12,nl1,nl2;
939  int *l1,*l2,*l3;
940  mprfloat q1, bmax;
941
942  if ( m != (m1+m2+m3) )
943  {
944    // error: bad input
945    error(WarnS("simplex::compute: Bad input constraint counts!");)
946    icase=-2;
947    return;
948  }
949
950  l1= (int *) omAlloc0( (n+1) * sizeof(int) );
951  l2= (int *) omAlloc0( (m+1) * sizeof(int) );
952  l3= (int *) omAlloc0( (m+1) * sizeof(int) );
953
954  nl1= n;
955  for ( k=1; k<=n; k++ ) l1[k]=izrov[k]=k;
956  nl2=m;
957  for ( i=1; i<=m; i++ )
958  {
959    if ( LiPM[i+1][1] < 0.0 )
960    {
961      // error: bad input
962      error(WarnS("simplex::compute: Bad input tableau!");)
963      error(Warn("simplex::compute: in input Matrix row %d, column 1, value %f",i+1,LiPM[i+1][1]);)
964      icase=-2;
965      // free mem l1,l2,l3;
966      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
967      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
968      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
969      return;
970    }
971    l2[i]= i;
972    iposv[i]= n+i;
973  }
974  for ( i=1; i<=m2; i++) l3[i]= 1;
975  ir= 0;
976  if (m2+m3)
977  {
978    ir=1;
979    for ( k=1; k <= (n+1); k++ )
980    {
981      q1=0.0;
982      for ( i=m1+1; i <= m; i++ ) q1+= LiPM[i+1][k];
983      LiPM[m+2][k]= -q1;
984    }
985
986    do
987    {
988      simp1(LiPM,m+1,l1,nl1,0,&kp,&bmax);
989      if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] < -SIMPLEX_EPS )
990      {
991        icase= -1; // no solution found
992        // free mem l1,l2,l3;
993        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
994        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
995        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
996        return;
997      }
998      else if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] <= SIMPLEX_EPS )
999      {
1000        m12= m1+m2+1;
1001        if ( m12 <= m )
1002        {
1003          for ( ip= m12; ip <= m; ip++ )
1004          {
1005            if ( iposv[ip] == (ip+n) )
1006            {
1007              simp1(LiPM,ip,l1,nl1,1,&kp,&bmax);
1008              if ( fabs(bmax) >= SIMPLEX_EPS)
1009                goto one;
1010            }
1011          }
1012        }
1013        ir= 0;
1014        --m12;
1015        if ( m1+1 <= m12 )
1016          for ( i=m1+1; i <= m12; i++ )
1017            if ( l3[i-m1] == 1 )
1018              for ( k=1; k <= n+1; k++ )
1019                LiPM[i+1][k] = -(LiPM[i+1][k]);
1020        break;
1021      }
1022      //#if DEBUG
1023      //print_bmat( a, m+2, n+3);
1024      //#endif
1025      simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1026      if ( ip == 0 )
1027      {
1028        icase = -1; // no solution found
1029        // free mem l1,l2,l3;
1030        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1031        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1032        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1033        return;
1034      }
1035    one: simp3(LiPM,m+1,n,ip,kp);
1036      // #if DEBUG
1037      // print_bmat(a,m+2,n+3);
1038      // #endif
1039      if ( iposv[ip] >= (n+m1+m2+1))
1040      {
1041        for ( k= 1; k <= nl1; k++ )
1042          if ( l1[k] == kp ) break;
1043        --nl1;
1044        for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1];
1045        ++(LiPM[m+2][kp+1]);
1046        for ( i= 1; i <= m+2; i++ ) LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1047      }
1048      else
1049      {
1050        if ( iposv[ip] >= (n+m1+1) )
1051        {
1052          kh= iposv[ip]-m1-n;
1053          if ( l3[kh] )
1054          {
1055            l3[kh]= 0;
1056            ++(LiPM[m+2][kp+1]);
1057            for ( i=1; i<= m+2; i++ )
1058              LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1059          }
1060        }
1061      }
1062      is= izrov[kp];
1063      izrov[kp]= iposv[ip];
1064      iposv[ip]= is;
1065    } while (ir);
1066  }
1067  /* end of phase 1, have feasible sol, now optimize it */
1068  loop
1069  {
1070    // #if DEBUG
1071    // print_bmat( a, m+1, n+5);
1072    // #endif
1073    simp1(LiPM,0,l1,nl1,0,&kp,&bmax);
1074    if (bmax <= /*SIMPLEX_EPS*/0.0)
1075    {
1076      icase=0; // finite solution found
1077      // free mem l1,l2,l3
1078      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1079      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1080      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1081      return;
1082    }
1083    simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1084    if (ip == 0)
1085    {
1086      //printf("Unbounded:");
1087      // #if DEBUG
1088      //       print_bmat( a, m+1, n+1);
1089      // #endif
1090      icase=1;                /* unbounded */
1091      // free mem
1092      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1093      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1094      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1095      return;
1096    }
1097    simp3(LiPM,m,n,ip,kp);
1098    is= izrov[kp];
1099    izrov[kp]= iposv[ip];
1100    iposv[ip]= is;
1101  }/*for ;;*/
1102}
1103
1104void simplex::simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax )
1105{
1106  int k;
1107  mprfloat test;
1108
1109  if( nll <= 0)
1110  {                        /* init'tion: fixed */
1111    *bmax = 0.0;
1112    return;
1113  }
1114  *kp=ll[1];
1115  *bmax=a[mm+1][*kp+1];
1116  for (k=2;k<=nll;k++)
1117  {
1118    if (iabf == 0)
1119    {
1120      test=a[mm+1][ll[k]+1]-(*bmax);
1121      if (test > 0.0)
1122      {
1123        *bmax=a[mm+1][ll[k]+1];
1124        *kp=ll[k];
1125      }
1126    }
1127    else
1128    {                        /* abs values: have fixed it */
1129      test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
1130      if (test > 0.0)
1131      {
1132        *bmax=a[mm+1][ll[k]+1];
1133        *kp=ll[k];
1134      }
1135    }
1136  }
1137}
1138
1139void simplex::simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 )
1140{
1141  int k,ii,i;
1142  mprfloat qp,q0,q;
1143
1144  *ip= 0;
1145  for ( i=1; i <= nl2; i++ )
1146  {
1147    if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS )
1148    {
1149      *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1];
1150      *ip= l2[i];
1151      for ( i= i+1; i <= nl2; i++ )
1152      {
1153        ii= l2[i];
1154        if (a[ii+1][kp+1] < -SIMPLEX_EPS)
1155        {
1156          q= -a[ii+1][1] / a[ii+1][kp+1];
1157          if (q - *q1 < -SIMPLEX_EPS)
1158          {
1159            *ip=ii;
1160            *q1=q;
1161          }
1162          else if (q - *q1 < SIMPLEX_EPS)
1163          {
1164            for ( k=1; k<= n; k++ )
1165            {
1166              qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
1167              q0= -a[ii+1][k+1]/a[ii+1][kp+1];
1168              if ( q0 != qp ) break;
1169            }
1170            if ( q0 < qp ) *ip= ii;
1171          }
1172        }
1173      }
1174    }
1175  }
1176}
1177
1178void simplex::simp3( mprfloat **a, int i1, int k1, int ip, int kp )
1179{
1180  int kk,ii;
1181  mprfloat piv;
1182
1183  piv= 1.0 / a[ip+1][kp+1];
1184  for ( ii=1; ii <= i1+1; ii++ )
1185  {
1186    if ( ii -1 != ip )
1187    {
1188      a[ii][kp+1] *= piv;
1189      for ( kk=1; kk <= k1+1; kk++ )
1190        if ( kk-1 != kp )
1191          a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
1192    }
1193  }
1194  for ( kk=1; kk<= k1+1; kk++ )
1195    if ( kk-1 != kp ) a[ip+1][kk] *= -piv;
1196  a[ip+1][kp+1]= piv;
1197}
1198//<-
1199
1200//-----------------------------------------------------------------------------
1201
1202//#endif // HAVE_MPR
1203
1204// local Variables: ***
1205// folded-file: t ***
1206// compile-command-1: "make installg" ***
1207// compile-command-2: "make install" ***
1208// End: ***
Note: See TracBrowser for help on using the repository browser.