source: git/Singular/mpr_numeric.cc @ 512a2b

spielwiese
Last change on this file since 512a2b was 512a2b, checked in by Olaf Bachmann <obachman@…>, 24 years ago
p_polys.h git-svn-id: file:///usr/local/Singular/svn/trunk@4606 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 28.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5/* $Id: mpr_numeric.cc,v 1.10 2000-09-18 09:19:20 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    for ( j= m-1; j >= 0; j-- )
600    {
601      f= ( *x * f ) + d;
602      d= ( *x * d ) + b;
603      b= ( *x * b ) + *a[j];
604      err_g= abs( b ) + ( abx_g * err_g );
605    }
606
607    err_g *= epss; // EPSS;
608
609    if ( abs(b) <= err_g ) return;
610
611    g= d / b;
612    g2 = g * g;
613    h= g2 - ( ( f / b ) * (gmp_float)2.0 );
614    sq= sqrt(( ( h * (gmp_float)m ) - g2 ) * (gmp_float)(m - 1));
615    gp= g + sq;
616    gm= g - sq;
617    abp_g= abs( gp );
618    abm_g= abs( gm );
619
620    if ( abp_g < abm_g ) gp= gm;
621
622    dx = ( (max(abp_g,abm_g) > (gmp_float)0.0)
623           ? ( gmp_complex( (mprfloat)m ) / gp )
624           : ( gmp_complex( cos((mprfloat)iter),sin((mprfloat)iter))
625               * exp(log((gmp_float)1.0+abx_g))) );
626
627    x1= *x - dx;
628
629    if ( (*x == x1) ) return;
630
631    if ( iter % MT ) *x= x1;
632    else *x -= ( dx * frac_g[ iter / MT ] );
633  }
634
635  *its= MAXIT+1;
636
637  return;
638}
639//<-
640
641//-----------------------------------------------------------------------------
642//-------------- rootArranger -------------------------------------------------
643//-----------------------------------------------------------------------------
644
645//-> rootArranger::rootArranger(...)
646rootArranger::rootArranger( rootContainer ** _roots,
647                            rootContainer ** _mu,
648                            const int _howclean )
649  : roots(_roots), mu(_mu), howclean(_howclean)
650{
651  found_roots=false;
652}
653//<-
654
655//-> void rootArranger::solve_all()
656void rootArranger::solve_all()
657{
658  int i;
659  found_roots= true;
660
661  // find roots of polys given by coeffs in roots
662  rc= roots[0]->getAnzElems();
663  for ( i= 0; i < rc; i++ )
664    if ( !roots[i]->solver( howclean ) )
665    {
666      found_roots= false;
667      return;
668    }
669  // find roots of polys given by coeffs in mu
670  mc= mu[0]->getAnzElems();
671  for ( i= 0; i < mc; i++ )
672    if ( ! mu[i]->solver( howclean ) )
673    {
674      found_roots= false;
675      return;
676    }
677}
678//<-
679
680//-> void rootArranger::arrange()
681void rootArranger::arrange()
682{
683  gmp_complex tmp,zwerg;
684  int anzm= mu[0]->getAnzElems();
685  int anzr= roots[0]->getAnzRoots();
686  int xkoord, r, rtest, xk, mtest;
687  bool found;
688  //gmp_complex mprec(1.0/pow(10,gmp_output_digits-5),1.0/pow(10,gmp_output_digits-5));
689  gmp_float mprec(1.0/pow(10.0,(int)(gmp_output_digits/3)));
690
691  for ( xkoord= 0; xkoord < anzm; xkoord++ ) {    // für x1,x2, x1,x2,x3, x1,x2,...,xn
692    for ( r= 0; r < anzr; r++ ) {                 // für jede Nullstelle
693      // (x1-koordinate) * evp[1] + (x2-koordinate) * evp[2] +
694      //                                  ... + (xkoord-koordinate) * evp[xkoord]
695      tmp= gmp_complex();
696      for ( xk =0; xk <= xkoord; xk++ )
697      {
698        tmp -= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); //xk+1
699      }
700      found= false;
701      for ( rtest= r; rtest < anzr; rtest++ ) {   // für jede Nullstelle
702        zwerg = tmp - (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xk+1, xkoord+2
703        for ( mtest= 0; mtest < anzr; mtest++ )
704        {
705          //          if ( tmp == (*mu[xkoord])[mtest] )
706          //          {
707          if ( ((zwerg.real() <= (*mu[xkoord])[mtest].real() + mprec) &&
708                (zwerg.real() >= (*mu[xkoord])[mtest].real() - mprec)) &&
709               ((zwerg.imag() <= (*mu[xkoord])[mtest].imag() + mprec) &&
710                (zwerg.imag() >= (*mu[xkoord])[mtest].imag() - mprec)) )
711           {
712             roots[xk]->swapRoots( r, rtest );
713             found= true;
714             break;
715           }
716        }
717      } // rtest
718      if ( !found )
719      {
720        Warn("rootArranger::arrange: No match? coord %d, root %d.",xkoord,r);
721#ifdef mprDEBUG_PROT
722        WarnS("One of these ...");
723        for ( rtest= r; rtest < anzr; rtest++ )
724        {
725          tmp= gmp_complex();
726          for ( xk =0; xk <= xkoord; xk++ )
727          {
728            tmp-= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1);
729          }
730          tmp-= (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xkoord+2
731          Warn("  %s",complexToStr(tmp,gmp_output_digits+1),rtest);
732        }
733        WarnS(" ... must match to one of these:");
734        for ( mtest= 0; mtest < anzr; mtest++ )
735        {
736          Warn("                  %s",complexToStr((*mu[xkoord])[mtest],gmp_output_digits+1));
737        }
738#endif
739      }
740    } // r
741  } // xkoord
742}
743//<-
744
745//-> lists rootArranger::listOfRoots( int oprec )
746lists rootArranger::listOfRoots( const unsigned int oprec )
747{
748  int i,j,tr;
749  int count= roots[0]->getAnzRoots(); // number of roots
750  int elem= roots[0]->getAnzElems();  // number of koordinates per root
751
752  lists listofroots= (lists)omAlloc( sizeof(slists) ); // must be done this way!
753
754  if ( found_roots )
755  {
756    listofroots->Init( count );
757
758    for (i=0; i < count; i++)
759    {
760      lists onepoint= (lists)omAlloc(sizeof(slists)); // must be done this way!
761      onepoint->Init(elem);
762      for ( j= 0; j < elem; j++ )
763      {
764        if ( !rField_is_long_C() )
765        {
766          onepoint->m[j].rtyp=STRING_CMD;
767          onepoint->m[j].data=(void *)complexToStr((*roots[j])[i],oprec);
768        }
769        else
770        {
771          onepoint->m[j].rtyp=NUMBER_CMD;
772          onepoint->m[j].data=(void *)nCopy((number)(roots[j]->getRoot(i)));
773        }
774        onepoint->m[j].next= NULL;
775        onepoint->m[j].name= NULL;
776      }
777      listofroots->m[i].rtyp=LIST_CMD;
778      listofroots->m[i].data=(void *)onepoint;
779      listofroots->m[j].next= NULL;
780      listofroots->m[j].name= NULL;
781    }
782
783  }
784  else
785  {
786    listofroots->Init( 0 );
787  }
788
789  return listofroots;
790}
791//<-
792
793//-----------------------------------------------------------------------------
794//-------------- simplex ----- ------------------------------------------------
795//-----------------------------------------------------------------------------
796
797//  #ifdef mprDEBUG_PROT
798//  #define error(a) a
799//  #else
800//  #define error(a)
801//  #endif
802
803#define error(a) a
804
805#define MAXPOINTS      1000
806
807//-> simplex::*
808//
809simplex::simplex( int rows, int cols )
810   : LiPM_cols(cols), LiPM_rows(rows)
811{
812  int i;
813
814  LiPM_rows=LiPM_rows+3;
815  LiPM_cols=LiPM_cols+2;
816
817  LiPM = (mprfloat **)omAlloc( LiPM_rows * sizeof(mprfloat *) );  // LP matrix
818  for( i= 0; i < LiPM_rows; i++ )
819  {
820    // Mem must be allocated aligned, also for type double!
821    LiPM[i] = (mprfloat *)omAlloc0Aligned( LiPM_cols * sizeof(mprfloat) );
822  }
823
824  iposv = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
825  izrov = (int *)omAlloc0( 2*LiPM_rows*sizeof(int) );
826
827  m=n=m1=m2=m3=icase=0;
828
829#ifdef mprDEBUG_ALL
830  Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols);
831#endif
832}
833
834simplex::~simplex()
835{
836  // clean up
837  int i;
838  for( i= 0; i < LiPM_rows; i++ )
839  {
840    omFreeSize( (ADDRESS) LiPM[i], LiPM_cols * sizeof(mprfloat) );
841  }
842  omFreeSize( (ADDRESS) LiPM, LiPM_rows * sizeof(mprfloat *) );
843
844  omFreeSize( (ADDRESS) iposv, 2*LiPM_rows*sizeof(int) );
845  omFreeSize( (ADDRESS) izrov, 2*LiPM_rows*sizeof(int) );
846}
847
848BOOLEAN simplex::mapFromMatrix( matrix m )
849{
850  int i,j;
851//    if ( MATROWS( m ) > LiPM_rows ||  MATCOLS( m ) > LiPM_cols ) {
852//      WarnS("");
853//      return FALSE;
854//    }
855
856  number coef;
857  for ( i= 1; i <= MATROWS( m ); i++ )
858  {
859     for ( j= 1; j <= MATCOLS( m ); j++ ) 
860     {
861        if ( MATELEM(m,i,j) != NULL ) 
862        {
863           coef= pGetCoeff( MATELEM(m,i,j) );
864           if ( coef != NULL && !nIsZero(coef) )
865              LiPM[i][j]= (double)(*(gmp_float*)coef);
866           //#ifdef mpr_DEBUG_PROT
867           //Print("%f ",LiPM[i][j]);
868           //#endif
869        }
870     }
871     //     PrintLn();
872  }
873
874  return TRUE;
875}
876
877matrix simplex::mapToMatrix( matrix m )
878{
879  int i,j;
880//    if ( MATROWS( m ) < LiPM_rows-3 ||  MATCOLS( m ) < LiPM_cols-2 ) {
881//      WarnS("");
882//      return NULL;
883//    }
884
885//Print(" %d x %d\n",MATROWS( m ),MATCOLS( m ));
886
887  number coef;
888  gmp_float * bla;
889  for ( i= 1; i <= MATROWS( m ); i++ )
890  {
891    for ( j= 1; j <= MATCOLS( m ); j++ )
892    {
893       pDelete( &(MATELEM(m,i,j)) );
894       MATELEM(m,i,j)= NULL;
895//Print(" %3.0f ",LiPM[i][j]);
896       if ( LiPM[i][j] != 0.0 ) 
897       {
898          bla= new gmp_float(LiPM[i][j]);
899          coef= (number)bla;
900          MATELEM(m,i,j)= pOne();
901          pSetCoeff( MATELEM(m,i,j), coef );
902       }
903    }
904//PrintLn();
905  }
906 
907  return m;
908}
909
910intvec * simplex::posvToIV()
911{
912   int i;
913   intvec * iv = new intvec( m );
914   for ( i= 1; i <= m; i++ )
915   {
916      IMATELEM(*iv,i,1)= iposv[i];
917   }
918   return iv;
919}
920
921intvec * simplex::zrovToIV()
922{
923   int i;
924   intvec * iv = new intvec( n );
925   for ( i= 1; i <= n; i++ )
926   {
927      IMATELEM(*iv,i,1)= izrov[i];
928   }
929   return iv; 
930}
931
932void simplex::compute()
933{
934  int i,ip,ir,is,k,kh,kp,m12,nl1,nl2;
935  int *l1,*l2,*l3;
936  mprfloat q1, bmax;
937
938  if ( m != (m1+m2+m3) )
939  {
940    // error: bad input
941    error(WarnS("simplex::compute: Bad input constraint counts!");)
942    icase=-2;
943    return;
944  }
945
946  l1= (int *) omAlloc0( (n+1) * sizeof(int) );
947  l2= (int *) omAlloc0( (m+1) * sizeof(int) );
948  l3= (int *) omAlloc0( (m+1) * sizeof(int) );
949
950  nl1= n;
951  for ( k=1; k<=n; k++ ) l1[k]=izrov[k]=k;
952  nl2=m;
953  for ( i=1; i<=m; i++ )
954  {
955    if ( LiPM[i+1][1] < 0.0 )
956    {
957      // error: bad input
958      error(WarnS("simplex::compute: Bad input tableau!");)
959      error(Warn("simplex::compute: in input Matrix row %d, column 1, value %f",i+1,LiPM[i+1][1]);)
960      icase=-2;
961      // free mem l1,l2,l3;
962      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
963      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
964      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
965      return;
966    }
967    l2[i]= i;
968    iposv[i]= n+i;
969  }
970  for ( i=1; i<=m2; i++) l3[i]= 1;
971  ir= 0;
972  if (m2+m3)
973  {
974    ir=1;
975    for ( k=1; k <= (n+1); k++ )
976    {
977      q1=0.0;
978      for ( i=m1+1; i <= m; i++ ) q1+= LiPM[i+1][k];
979      LiPM[m+2][k]= -q1;
980    }
981
982    do
983    {
984      simp1(LiPM,m+1,l1,nl1,0,&kp,&bmax);
985      if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] < -SIMPLEX_EPS )
986      {
987        icase= -1; // no solution found
988        // free mem l1,l2,l3;
989        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
990        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
991        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
992        return;
993      }
994      else if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] <= SIMPLEX_EPS )
995      {
996        m12= m1+m2+1;
997        if ( m12 <= m )
998        {
999          for ( ip= m12; ip <= m; ip++ )
1000          {
1001            if ( iposv[ip] == (ip+n) )
1002            {
1003              simp1(LiPM,ip,l1,nl1,1,&kp,&bmax);
1004              if ( fabs(bmax) >= SIMPLEX_EPS)
1005                goto one;
1006            }
1007          }
1008        }
1009        ir= 0;
1010        --m12;
1011        if ( m1+1 <= m12 )
1012          for ( i=m1+1; i <= m12; i++ )
1013            if ( l3[i-m1] == 1 )
1014              for ( k=1; k <= n+1; k++ )
1015                LiPM[i+1][k] = -(LiPM[i+1][k]);
1016        break;
1017      }
1018      //#if DEBUG
1019      //print_bmat( a, m+2, n+3);
1020      //#endif
1021      simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1022      if ( ip == 0 )
1023      {
1024        icase = -1; // no solution found
1025        // free mem l1,l2,l3;
1026        omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1027        omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1028        omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1029        return;
1030      }
1031    one: simp3(LiPM,m+1,n,ip,kp);
1032      // #if DEBUG
1033      // print_bmat(a,m+2,n+3);
1034      // #endif
1035      if ( iposv[ip] >= (n+m1+m2+1))
1036      {
1037        for ( k= 1; k <= nl1; k++ )
1038          if ( l1[k] == kp ) break;
1039        --nl1;
1040        for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1];
1041        ++(LiPM[m+2][kp+1]);
1042        for ( i= 1; i <= m+2; i++ ) LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1043      }
1044      else
1045      {
1046        if ( iposv[ip] >= (n+m1+1) )
1047        {
1048          kh= iposv[ip]-m1-n;
1049          if ( l3[kh] )
1050          {
1051            l3[kh]= 0;
1052            ++(LiPM[m+2][kp+1]);
1053            for ( i=1; i<= m+2; i++ )
1054              LiPM[i][kp+1] = -(LiPM[i][kp+1]);
1055          }
1056        }
1057      }
1058      is= izrov[kp];
1059      izrov[kp]= iposv[ip];
1060      iposv[ip]= is;
1061    } while (ir);
1062  }
1063  /* end of phase 1, have feasible sol, now optimize it */
1064  loop
1065  {
1066    // #if DEBUG
1067    // print_bmat( a, m+1, n+5);
1068    // #endif
1069    simp1(LiPM,0,l1,nl1,0,&kp,&bmax);
1070    if (bmax <= /*SIMPLEX_EPS*/0.0)
1071    {
1072      icase=0; // finite solution found
1073      // free mem l1,l2,l3
1074      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1075      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1076      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1077      return;
1078    }
1079    simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
1080    if (ip == 0)
1081    {
1082      //printf("Unbounded:");
1083      // #if DEBUG
1084      //       print_bmat( a, m+1, n+1);
1085      // #endif
1086      icase=1;                /* unbounded */
1087      // free mem
1088      omFreeSize( (ADDRESS) l3, (m+1) * sizeof(int) );
1089      omFreeSize( (ADDRESS) l2, (m+1) * sizeof(int) );
1090      omFreeSize( (ADDRESS) l1, (n+1) * sizeof(int) );
1091      return;
1092    }
1093    simp3(LiPM,m,n,ip,kp);
1094    is= izrov[kp];
1095    izrov[kp]= iposv[ip];
1096    iposv[ip]= is;
1097  }/*for ;;*/
1098}
1099
1100void simplex::simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax )
1101{
1102  int k;
1103  mprfloat test;
1104
1105  if( nll <= 0)
1106  {                        /* init'tion: fixed */
1107    *bmax = 0.0;
1108    return;
1109  }
1110  *kp=ll[1];
1111  *bmax=a[mm+1][*kp+1];
1112  for (k=2;k<=nll;k++)
1113  {
1114    if (iabf == 0)
1115    {
1116      test=a[mm+1][ll[k]+1]-(*bmax);
1117      if (test > 0.0)
1118      {
1119        *bmax=a[mm+1][ll[k]+1];
1120        *kp=ll[k];
1121      }
1122    }
1123    else
1124    {                        /* abs values: have fixed it */
1125      test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
1126      if (test > 0.0)
1127      {
1128        *bmax=a[mm+1][ll[k]+1];
1129        *kp=ll[k];
1130      }
1131    }
1132  }
1133}
1134
1135void simplex::simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 )
1136{
1137  int k,ii,i;
1138  mprfloat qp,q0,q;
1139
1140  *ip= 0;
1141  for ( i=1; i <= nl2; i++ )
1142  {
1143    if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS )
1144    {
1145      *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1];
1146      *ip= l2[i];
1147      for ( i= i+1; i <= nl2; i++ )
1148      {
1149        ii= l2[i];
1150        if (a[ii+1][kp+1] < -SIMPLEX_EPS)
1151        {
1152          q= -a[ii+1][1] / a[ii+1][kp+1];
1153          if (q - *q1 < -SIMPLEX_EPS)
1154          {
1155            *ip=ii;
1156            *q1=q;
1157          }
1158          else if (q - *q1 < SIMPLEX_EPS)
1159          {
1160            for ( k=1; k<= n; k++ )
1161            {
1162              qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
1163              q0= -a[ii+1][k+1]/a[ii+1][kp+1];
1164              if ( q0 != qp ) break;
1165            }
1166            if ( q0 < qp ) *ip= ii;
1167          }
1168        }
1169      }
1170    }
1171  }
1172}
1173
1174void simplex::simp3( mprfloat **a, int i1, int k1, int ip, int kp )
1175{
1176  int kk,ii;
1177  mprfloat piv;
1178
1179  piv= 1.0 / a[ip+1][kp+1];
1180  for ( ii=1; ii <= i1+1; ii++ )
1181  {
1182    if ( ii -1 != ip )
1183    {
1184      a[ii][kp+1] *= piv;
1185      for ( kk=1; kk <= k1+1; kk++ )
1186        if ( kk-1 != kp )
1187          a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
1188    }
1189  }
1190  for ( kk=1; kk<= k1+1; kk++ )
1191    if ( kk-1 != kp ) a[ip+1][kk] *= -piv;
1192  a[ip+1][kp+1]= piv;
1193}
1194//<-
1195
1196//-----------------------------------------------------------------------------
1197
1198//#endif // HAVE_MPR
1199
1200// local Variables: ***
1201// folded-file: t ***
1202// compile-command-1: "make installg" ***
1203// compile-command-2: "make install" ***
1204// End: ***
Note: See TracBrowser for help on using the repository browser.