source: git/Singular/mpr_inout.cc @ 9a6c4a

spielwiese
Last change on this file since 9a6c4a was 9a6c4a, checked in by Moritz Wenk <wenk@…>, 24 years ago
* wenk: fixed bug in ngfRead git-svn-id: file:///usr/local/Singular/svn/trunk@3182 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 14.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5/* $Id: mpr_inout.cc,v 1.3 1999-06-29 09:03:45 wenk Exp $ */
6
7/*
8* ABSTRACT - multipolynomial resultant
9*/
10
11#include "mod2.h"
12//#ifdef HAVE_MPR
13
14//-> includes
15#include "tok.h"
16#include "structs.h"
17#include "subexpr.h"
18#include "polys.h"
19#include "ideals.h"
20#include "ring.h"
21#include "ipid.h"
22#include "ipshell.h"
23#include "febase.h"
24#include "mmemory.h"
25#include "numbers.h"
26#include "lists.h"
27
28#include <math.h>
29
30#include "mpr_global.h"
31#include "mpr_inout.h"
32#include "mpr_base.h"
33#include "mpr_numeric.h"
34
35// to get detailed timigs, define MPR_TIMING
36#ifdef MPR_TIMING
37#define TIMING
38#endif
39#include "../factory/timing.h"
40TIMING_DEFINE_PRINT(mpr_overall)
41TIMING_DEFINE_PRINT(mpr_check)
42TIMING_DEFINE_PRINT(mpr_constr)
43TIMING_DEFINE_PRINT(mpr_ures)
44TIMING_DEFINE_PRINT(mpr_mures)
45TIMING_DEFINE_PRINT(mpr_arrange)
46TIMING_DEFINE_PRINT(mpr_solver)
47
48#define TIMING_EPR(t,msg) TIMING_END_AND_PRINT(t,msg);TIMING_RESET(t);
49
50enum mprState
51{
52    mprOk,
53    mprWrongRType,
54    mprHasOne,
55    mprInfNumOfVars,
56    mprNotReduced,
57    mprNotZeroDim,
58    mprNotHomog,
59    mprUnSupField
60};
61
62extern size_t gmp_output_digits;
63//<-
64
65//-> nPrint(number n)
66void nPrint(number n)
67{
68  poly o=pOne();
69  pSetCoeff(o, nCopy(n) );
70  pWrite0( o );
71  pDelete( &o );
72}
73//<-
74
75//------------------------------------------------------------------------------
76
77//-> void mprPrintError( mprState state )
78void mprPrintError( mprState state, const char * name )
79{
80  switch (state)
81  {
82  case mprWrongRType:
83    WerrorS("Unknown resultant matrix type choosen!");
84    break;
85  case mprHasOne:
86    Werror("One element of the ideal %s is constant!",name);
87    break;
88  case mprInfNumOfVars:
89    Werror("Numer of elements in given ideal %s must be equal to %d!",
90           name,pVariables+1);
91    break;
92  case mprNotZeroDim:
93    Werror("The given ideal %s must 0-dimensional!",name);
94    break;
95  case mprNotHomog:
96    Werror("The given ideal %s has to be homogeneous in"
97           " the first ring variable!",name);
98    break;
99  case mprNotReduced:
100    Werror("The given ideal %s has to reduced!",name);
101    break;
102  case mprUnSupField:
103    WerrorS("Ground field not implemented!");
104    break;
105  default:
106    break;
107  }
108}
109//<-
110
111//-> mprState mprIdealCheck()
112mprState mprIdealCheck( const ideal theIdeal,
113                        const char * name,
114                        uResultant::resMatType mtype,
115                        BOOLEAN rmatrix= false )
116{
117  mprState state = mprOk;
118  int power;
119  int k;
120
121  int numOfVars= mtype == uResultant::denseResMat?pVariables-1:pVariables;
122  if ( rmatrix ) numOfVars++;
123
124  if ( mtype == uResultant::none )
125    state= mprWrongRType;
126
127  if ( IDELEMS(theIdeal) != numOfVars )
128    state= mprInfNumOfVars;
129
130  for ( k= IDELEMS(theIdeal) - 1; (state == mprOk) && (k >= 0); k-- )
131  {
132    poly p = (theIdeal->m)[k];
133    if ( pIsConstant(p) ) state= mprHasOne;
134    else
135    if ( (mtype == uResultant::denseResMat) && !pIsHomogeneous(p) )
136      state=mprNotHomog;
137  }
138
139  if ( !(rField_is_R()||
140         rField_is_Q()||
141         rField_is_long_R()||
142         rField_is_long_C()||
143         (rmatrix && rField_is_Q_a())) )
144    state= mprUnSupField;
145
146  if ( state != mprOk ) mprPrintError( state, "" /* name */ );
147
148  return state;
149}
150//<-
151
152//-> uResultant::resMatType determineMType( int imtype )
153uResultant::resMatType determineMType( int imtype )
154{
155  switch ( imtype )
156  {
157  case MPR_DENSE:
158    return uResultant::denseResMat;
159  case 0:
160  case MPR_SPARSE:
161    return uResultant::sparseResMat;
162  default:
163    return uResultant::none;
164  }
165}
166//<-
167
168//-> BOOLEAN nuUResSolve( leftv res, leftv arg1, leftv arg2, leftv arg3 )
169BOOLEAN nuUResSolve( leftv res, leftv arg1, leftv arg2, leftv arg3 )
170{
171  ideal gls;
172  gls= (ideal)(arg1->Data());
173  int imtype= (int)arg2->Data();
174  int howclean= (int)arg3->Data();
175
176  uResultant::resMatType mtype= determineMType( imtype );
177  int i,c,count;
178  lists listofroots= NULL;
179  lists emptylist;
180  number smv= NULL;
181  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
182
183  emptylist= (lists)Alloc( sizeof(slists) );
184  emptylist->Init( 0 );
185
186  res->rtyp = LIST_CMD;
187  res->data= (void *)emptylist;
188
189  TIMING_START(mpr_overall);
190
191  // check input ideal ( = polynomial system )
192  if ( mprIdealCheck( gls, arg1->Name(), mtype ) != mprOk )
193  {
194    return TRUE;
195  }
196
197  uResultant * ures;
198  rootContainer ** iproots;
199  rootContainer ** muiproots;
200  rootArranger * arranger;
201
202  // main task 1: setup of resultant matrix
203  TIMING_START(mpr_constr);
204  ures= new uResultant( gls, mtype );
205  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
206  {
207    WerrorS("Error occurred during matrix setup!");
208    return TRUE;
209  }
210  TIMING_EPR(mpr_constr, "construction\t\t")
211
212  // if dense resultant, check if minor nonsingular
213  TIMING_START(mpr_check);
214  if ( mtype == uResultant::denseResMat )
215  {
216    smv= ures->accessResMat()->getSubDet();
217#ifdef mprDEBUG_PROT
218    Print("// Determinant of submatrix: ");nPrint(smv);PrintLn();
219#endif
220    if ( nIsZero(smv) )
221    {
222      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
223      return TRUE;
224    }
225  }
226  TIMING_EPR(mpr_check,  "input check\t\t")
227
228  // main task 2: Interpolate specialized resultant polynomials
229  TIMING_START(mpr_ures);
230  if ( interpolate_det )
231    iproots= ures->interpolateDenseSP( false, smv );
232  else
233    iproots= ures->specializeInU( false, smv );
234  TIMING_EPR(mpr_ures,   "interpolation ures\t")
235
236  // main task 3: Interpolate specialized resultant polynomials
237  TIMING_START(mpr_mures);
238  if ( interpolate_det )
239    muiproots= ures->interpolateDenseSP( true, smv );
240  else
241    muiproots= ures->specializeInU( true, smv );
242  TIMING_EPR(mpr_mures,  "interpolation mures\t")
243
244#ifdef mprDEBUG_PROT
245  c= iproots[0]->getAnzElems();
246  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
247  c= muiproots[0]->getAnzElems();
248  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
249#endif
250
251  // main task 4: Compute roots of specialized polys and match them up
252  arranger= new rootArranger( iproots, muiproots, howclean );
253  TIMING_START(mpr_solver);
254  arranger->solve_all();
255  TIMING_EPR(mpr_solver, "solver time\t\t");
256
257  // get list of roots
258  if ( arranger->success() )
259  {
260    TIMING_START(mpr_arrange);
261    arranger->arrange();
262    TIMING_EPR(mpr_arrange, "arrange time\t\t");
263    listofroots= arranger->listOfRoots( gmp_output_digits );
264  }
265  else
266  {
267    WerrorS("Solver was unable to find any root!");
268    return TRUE;
269  }
270
271  // free everything
272  count= iproots[0]->getAnzElems();
273  for (i=0; i < count; i++) delete iproots[i];
274  Free( (ADDRESS) iproots, count * sizeof(rootContainer*) );
275  count= muiproots[0]->getAnzElems();
276  for (i=0; i < count; i++) delete muiproots[i];
277  Free( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
278
279  delete ures;
280  delete arranger;
281  nDelete( &smv );
282
283  res->data= (void *)listofroots;
284
285  emptylist->Clean();
286  //Free( (ADDRESS) emptylist, sizeof(slists) );
287
288  TIMING_EPR(mpr_overall,"overall time\t\t")
289
290  return FALSE;
291}
292//<-
293
294//-> BOOLEAN nuMPResMat( leftv res, leftv arg1, leftv arg2 )
295BOOLEAN nuMPResMat( leftv res, leftv arg1, leftv arg2 )
296{
297  ideal gls;
298  gls= (ideal)(arg1->Data());
299  int imtype= (int)arg2->Data();
300
301  uResultant::resMatType mtype= determineMType( imtype );
302
303  // check input ideal ( = polynomial system )
304  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
305  {
306    return TRUE;
307  }
308
309  uResultant *resMat= new uResultant( gls, mtype, false );
310
311  res->rtyp = MODUL_CMD;
312  res->data= (void*)resMat->accessResMat()->getMatrix();
313
314  delete resMat;
315
316  return FALSE;
317}
318//<-
319
320//-> BOOLEAN nuLagSolve( leftv res, leftv arg1, leftv arg2 )
321BOOLEAN nuLagSolve( leftv res, leftv arg1, leftv arg2 )
322{
323
324  poly gls;
325  gls= (poly)(arg1->Data());
326  int howclean= (int)arg2->Data();
327
328  int deg= pTotaldegree( gls );
329  //  int deg= pDeg( gls );
330  int len= pLength( gls );
331  int i,vpos;
332  poly piter;
333  lists elist;
334  lists rlist;
335
336  elist= (lists)Alloc( sizeof(slists) );
337  elist->Init( 0 );
338
339  if ( !(rField_is_R() ||
340         rField_is_Q() ||
341         rField_is_long_R() ||
342         rField_is_long_C()) )
343         {
344    WerrorS("Ground field not implemented!");
345    return TRUE;
346  }
347
348  if ( pVariables > 1 )
349  {
350    piter= gls;
351    for ( i= 1; i <= pVariables; i++ )
352      if ( pGetExp( piter, i ) )
353      {
354        vpos= i;
355        break;
356      }
357    while ( piter )
358    {
359      for ( i= 1; i <= pVariables; i++ )
360        if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
361        {
362          Werror("The polynomial %s must be univariate!",arg1->Name());
363          return TRUE;
364        }
365      pIter( piter );
366    }
367  }
368
369  rootContainer * roots= new rootContainer();
370  number * pcoeffs= (number *)Alloc( (deg+1) * sizeof( number ) );
371  piter= gls;
372  for ( i= deg; i >= 0; i-- )
373  {
374    //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
375    if ( piter && pTotaldegree(piter) == i )
376    {
377      pcoeffs[i]= nCopy( pGetCoeff( piter ) );
378      //nPrint( pcoeffs[i] );Print("  ");
379      pIter( piter );
380    }
381    else
382    {
383      pcoeffs[i]= nInit(0);
384    }
385  }
386
387#ifdef mprDEBUG_PROT
388  for (i=deg; i >= 0; i--)
389  {
390    nPrint( pcoeffs[i] );Print("  ");
391  }
392  PrintLn();
393#endif
394
395  TIMING_START(mpr_solver);
396  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
397  roots->solver( howclean );
398  TIMING_EPR(mpr_solver, "solver time\t\t");
399
400  int elem= roots->getAnzRoots();
401  char *out;
402  char *dummy;
403  int j;
404
405  rlist= (lists)Alloc( sizeof(slists) );
406  rlist->Init( elem );
407
408  if (rField_is_long_C())
409  {
410    for ( j= 0; j < elem; j++ )
411    {
412      rlist->m[j].rtyp=NUMBER_CMD;
413      rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
414      //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
415    }
416  }
417  else
418  {
419    for ( j= 0; j < elem; j++ )
420    {
421      dummy = complexToStr( (*roots)[j], gmp_output_digits );
422      rlist->m[j].rtyp=STRING_CMD;
423      rlist->m[j].data=(void *)dummy;
424    }
425  }
426
427  elist->Clean();
428  //Free( (ADDRESS) elist, sizeof(slists) );
429
430  for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
431  Free( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
432
433  res->rtyp= LIST_CMD;
434  res->data= (void*)rlist;
435
436  return FALSE;
437}
438//<-
439
440//-> BOOLEAN nuVanderSys( leftv res, leftv arg1, leftv arg2 )
441BOOLEAN nuVanderSys( leftv res, leftv arg1, leftv arg2, leftv arg3)
442{
443  int i;
444  ideal p,w;
445  p= (ideal)arg1->Data();
446  w= (ideal)arg2->Data();
447
448  // w[0] = f(p^0)
449  // w[1] = f(p^1)
450  // ...
451  // p can be a vector of numbers (multivariate polynom)
452  //   or one number (univariate polynom)
453  // tdg = deg(f)
454
455  int n= IDELEMS( p );
456  int m= IDELEMS( w );
457  int tdg= (int)arg3->Data();
458
459  res->data= (void*)NULL;
460
461  // check the input
462  if ( tdg < 1 )
463  {
464    WerrorS("Last input parameter must be > 0!");
465    return TRUE;
466  }
467  if ( n != pVariables )
468  {
469    Werror("Size of first input ideal must be equal to %d!",pVariables);
470    return TRUE;
471  }
472  if ( m != (int)pow((double)tdg+1,(int)n) )
473  {
474    Werror("Size of second input ideal must be equal to %d!",
475      (int)pow((double)tdg+1,(int)n));
476    return TRUE;
477  }
478  if ( !(rField_is_Q() /* ||
479         rField_is_R() || rField_is_long_R() ||
480         rField_is_long_C()*/ ) )
481         {
482    WerrorS("Ground field not implemented!");
483    return TRUE;
484  }
485
486  number tmp;
487  number *pevpoint= (number *)Alloc( n * sizeof( number ) );
488  for ( i= 0; i < n; i++ )
489  {
490    pevpoint[i]=nInit(0);
491    if (  (p->m)[i] )
492    {
493      tmp = pGetCoeff( (p->m)[i] );
494      if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
495      {
496        Free( (ADDRESS)pevpoint, n * sizeof( number ) );
497        WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
498        return TRUE;
499      }
500    } else tmp= NULL;
501    if ( !nIsZero(tmp) )
502    {
503      if ( !pIsConstant((p->m)[i]))
504      {
505        Free( (ADDRESS)pevpoint, n * sizeof( number ) );
506        WerrorS("Elements of first input ideal must be numbers!");
507        return TRUE;
508      }
509      pevpoint[i]= nCopy( tmp );
510    }
511  }
512
513  number *wresults= (number *)Alloc( m * sizeof( number ) );
514  for ( i= 0; i < m; i++ )
515  {
516    wresults[i]= nInit(0);
517    if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
518    {
519      if ( !pIsConstant((w->m)[i]))
520      {
521        Free( (ADDRESS)pevpoint, n * sizeof( number ) );
522        Free( (ADDRESS)wresults, m * sizeof( number ) );
523        WerrorS("Elements of second input ideal must be numbers!");
524        return TRUE;
525      }
526      wresults[i]= nCopy(pGetCoeff((w->m)[i]));
527    }
528  }
529
530  vandermonde vm( m, n, tdg, pevpoint, FALSE );
531  number *ncpoly= vm.interpolateDense( wresults );
532  // do not free ncpoly[]!!
533  poly rpoly= vm.numvec2poly( ncpoly );
534
535  Free( (ADDRESS)pevpoint, n * sizeof( number ) );
536  Free( (ADDRESS)wresults, m * sizeof( number ) );
537
538  res->data= (void*)rpoly;
539  return FALSE;
540}
541//<-
542
543//-> function u_resultant_det
544poly u_resultant_det( ideal gls, int imtype )
545{
546  uResultant::resMatType mtype= determineMType( imtype );
547  poly resdet;
548  poly emptypoly= pInit();
549  number smv= NULL;
550
551  TIMING_START(mpr_overall);
552
553  // check input ideal ( = polynomial system )
554  if ( mprIdealCheck( gls, "", mtype ) != mprOk )
555  {
556    return emptypoly;
557  }
558
559  uResultant *ures;
560
561  // main task 1: setup of resultant matrix
562  TIMING_START(mpr_constr);
563  ures= new uResultant( gls, mtype );
564  TIMING_EPR(mpr_constr,"construction");
565
566  // if dense resultant, check if minor nonsingular
567  if ( mtype == uResultant::denseResMat )
568  {
569    smv= ures->accessResMat()->getSubDet();
570#ifdef mprDEBUG_PROT
571    Print("// Determinant of submatrix: ");nPrint(smv); PrintLn();
572#endif
573    if ( nIsZero(smv) )
574    {
575      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
576      return emptypoly;
577    }
578  }
579
580  // main task 2: Interpolate resultant polynomial
581  TIMING_START(mpr_ures);
582  resdet= ures->interpolateDense( smv );
583  TIMING_EPR(mpr_ures,"ures");
584
585  // free mem
586  delete ures;
587  nDelete( &smv );
588  pDelete( &emptypoly );
589
590  TIMING_EPR(mpr_overall,"overall");
591
592  return ( resdet );
593}
594//<-
595
596//-----------------------------------------------------------------------------
597
598//#endif // HAVE_MPR
599
600// local Variables: ***
601// folded-file: t ***
602// compile-command-1: "make installg" ***
603// compile-command-2: "make install" ***
604// End: ***
605
606// in folding: C-c x
607// leave fold: C-c y
608//   foldmode: F10
609
610
611
612
613
614
615
616
617
Note: See TracBrowser for help on using the repository browser.