My Project
Public Types | Public Member Functions | Private Member Functions | Private Attributes
uResultant Class Reference

Base class for solving 0-dim poly systems using u-resultant. More...

#include <mpr_base.h>

Public Types

enum  resMatType { none , sparseResMat , denseResMat }
 

Public Member Functions

 uResultant (const ideal _gls, const resMatType _rmt=sparseResMat, BOOLEAN extIdeal=true)
 
 ~uResultant ()
 
poly interpolateDense (const number subDetVal=NULL)
 
rootContainer ** interpolateDenseSP (BOOLEAN matchUp=false, const number subDetVal=NULL)
 
rootContainer ** specializeInU (BOOLEAN matchUp=false, const number subDetVal=NULL)
 
resMatrixBaseaccessResMat ()
 

Private Member Functions

 uResultant (const uResultant &)
 
ideal extendIdeal (const ideal gls, poly linPoly, const resMatType rmt)
 
poly linearPoly (const resMatType rmt)
 
int nextPrime (const int p)
 

Private Attributes

ideal gls
 
int n
 
resMatType rmt
 
resMatrixBaseresMat
 

Detailed Description

Base class for solving 0-dim poly systems using u-resultant.

Definition at line 62 of file mpr_base.h.

Member Enumeration Documentation

◆ resMatType

Enumerator
none 
sparseResMat 
denseResMat 

Definition at line 65 of file mpr_base.h.

@ denseResMat
Definition: mpr_base.h:65
@ sparseResMat
Definition: mpr_base.h:65

Constructor & Destructor Documentation

◆ uResultant() [1/2]

uResultant::uResultant ( const ideal  _gls,
const resMatType  _rmt = sparseResMat,
BOOLEAN  extIdeal = true 
)

Definition at line 2685 of file mpr_base.cc.

2686  : rmt( _rmt )
2687 {
2688  if ( extIdeal )
2689  {
2690  // extend given ideal by linear poly F0=u0x0 + u1x1 +...+ unxn
2691  gls= extendIdeal( _gls, linearPoly( rmt ), rmt );
2692  n= IDELEMS( gls );
2693  }
2694  else
2695  gls= idCopy( _gls );
2696 
2697  switch ( rmt )
2698  {
2699  case sparseResMat:
2700  resMat= new resMatrixSparse( gls );
2701  break;
2702  case denseResMat:
2703  resMat= new resMatrixDense( gls );
2704  break;
2705  default:
2706  WerrorS("uResultant::uResultant: Unknown chosen resultant matrix type!");
2707  }
2708 }
poly linearPoly(const resMatType rmt)
Definition: mpr_base.cc:2743
resMatType rmt
Definition: mpr_base.h:91
resMatrixBase * resMat
Definition: mpr_base.h:92
ideal extendIdeal(const ideal gls, poly linPoly, const resMatType rmt)
Definition: mpr_base.cc:2715
ideal gls
Definition: mpr_base.h:88
void WerrorS(const char *s)
Definition: feFopen.cc:24
ideal idCopy(ideal A)
Definition: ideals.h:60
#define IDELEMS(i)
Definition: simpleideals.h:23

◆ ~uResultant()

uResultant::~uResultant ( )

Definition at line 2710 of file mpr_base.cc.

2711 {
2712  delete resMat;
2713 }

◆ uResultant() [2/2]

uResultant::uResultant ( const uResultant )
private

Member Function Documentation

◆ accessResMat()

resMatrixBase* uResultant::accessResMat ( )
inline

Definition at line 78 of file mpr_base.h.

78 { return resMat; }

◆ extendIdeal()

ideal uResultant::extendIdeal ( const ideal  gls,
poly  linPoly,
const resMatType  rmt 
)
private

Definition at line 2715 of file mpr_base.cc.

2716 {
2717  ideal newGls= idCopy( igls );
2718  newGls->m= (poly *)omReallocSize( newGls->m,
2719  IDELEMS(igls) * sizeof(poly),
2720  (IDELEMS(igls) + 1) * sizeof(poly) );
2721  IDELEMS(newGls)++;
2722 
2723  switch ( rrmt )
2724  {
2725  case sparseResMat:
2726  case denseResMat:
2727  {
2728  int i;
2729  for ( i= IDELEMS(newGls)-1; i > 0; i-- )
2730  {
2731  newGls->m[i]= newGls->m[i-1];
2732  }
2733  newGls->m[0]= linPoly;
2734  }
2735  break;
2736  default:
2737  WerrorS("uResultant::extendIdeal: Unknown chosen resultant matrix type!");
2738  }
2739 
2740  return( newGls );
2741 }
int i
Definition: cfEzgcd.cc:132
#define omReallocSize(addr, o_size, size)
Definition: omAllocDecl.h:220

◆ interpolateDense()

poly uResultant::interpolateDense ( const number  subDetVal = NULL)

Definition at line 2770 of file mpr_base.cc.

2771 {
2772  int i,j,p;
2773  long tdg;
2774 
2775  // D is a Polynom homogeneous in the coeffs of F0 and of degree tdg = d0*d1*...*dn
2776  tdg= resMat->getDetDeg();
2777 
2778  // maximum number of terms in polynom D (homogeneous, of degree tdg)
2779  // long mdg= (facul(tdg+n-1) / facul( tdg )) / facul( n - 1 );
2780  long mdg= over( n-1, tdg );
2781 
2782  // maximal number of terms in a polynom of degree tdg
2783  long l=(long)pow( (double)(tdg+1), n );
2784 
2785 #ifdef mprDEBUG_PROT
2786  Print("// total deg of D: tdg %ld\n",tdg);
2787  Print("// maximum number of terms in D: mdg: %ld\n",mdg);
2788  Print("// maximum number of terms in polynom of deg tdg: l %ld\n",l);
2789 #endif
2790 
2791  // we need mdg results of D(p0,p1,...,pn)
2792  number *presults;
2793  presults= (number *)omAlloc( mdg * sizeof( number ) );
2794  for (i=0; i < mdg; i++) presults[i]= nInit(0);
2795 
2796  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
2797  number *pev= (number *)omAlloc( n * sizeof( number ) );
2798  for (i=0; i < n; i++) pev[i]= nInit(0);
2799 
2800  mprPROTnl("// initial evaluation point: ");
2801  // initial evaluatoin point
2802  p=1;
2803  for (i=0; i < n; i++)
2804  {
2805  // init pevpoint with primes 3,5,7,11, ...
2806  p= nextPrime( p );
2807  pevpoint[i]=nInit( p );
2808  nTest(pevpoint[i]);
2809  mprPROTNnl(" ",pevpoint[i]);
2810  }
2811 
2812  // evaluate the determinant in the points pev^0, pev^1, ..., pev^mdg
2813  mprPROTnl("// evaluating:");
2814  for ( i=0; i < mdg; i++ )
2815  {
2816  for (j=0; j < n; j++)
2817  {
2818  nDelete( &pev[j] );
2819  nPower(pevpoint[j],i,&pev[j]);
2820  mprPROTN(" ",pev[j]);
2821  }
2822  mprPROTnl("");
2823 
2824  nDelete( &presults[i] );
2825  presults[i]=resMat->getDetAt( pev );
2826 
2828  }
2829  mprSTICKYPROT("\n");
2830 
2831  // now interpolate using vandermode interpolation
2832  mprPROTnl("// interpolating:");
2833  number *ncpoly;
2834  {
2835  vandermonde vm( mdg, n, tdg, pevpoint );
2836  ncpoly= vm.interpolateDense( presults );
2837  }
2838 
2839  if ( subDetVal != NULL )
2840  { // divide by common factor
2841  number detdiv;
2842  for ( i= 0; i <= mdg; i++ )
2843  {
2844  detdiv= nDiv( ncpoly[i], subDetVal );
2845  nNormalize( detdiv );
2846  nDelete( &ncpoly[i] );
2847  ncpoly[i]= detdiv;
2848  }
2849  }
2850 
2851 #ifdef mprDEBUG_ALL
2852  PrintLn();
2853  for ( i=0; i < mdg; i++ )
2854  {
2855  nPrint(ncpoly[i]); PrintS(" --- ");
2856  }
2857  PrintLn();
2858 #endif
2859 
2860  // prepare ncpoly for later use
2861  number nn=nInit(0);
2862  for ( i=0; i < mdg; i++ )
2863  {
2864  if ( nEqual(ncpoly[i],nn) )
2865  {
2866  nDelete( &ncpoly[i] );
2867  ncpoly[i]=NULL;
2868  }
2869  }
2870  nDelete( &nn );
2871 
2872  // create poly presenting the determinat of the uResultant
2873  intvec exp( n );
2874  for ( i= 0; i < n; i++ ) exp[i]=0;
2875 
2876  poly result= NULL;
2877 
2878  long sum=0;
2879  long c=0;
2880 
2881  for ( i=0; i < l; i++ )
2882  {
2883  if ( sum == tdg )
2884  {
2885  if ( !nIsZero(ncpoly[c]) )
2886  {
2887  poly p= pOne();
2888  if ( rmt == denseResMat )
2889  {
2890  for ( j= 0; j < n; j++ ) pSetExp( p, j+1, exp[j] );
2891  }
2892  else if ( rmt == sparseResMat )
2893  {
2894  for ( j= 1; j < n; j++ ) pSetExp( p, j, exp[j] );
2895  }
2896  pSetCoeff( p, ncpoly[c] );
2897  pSetm( p );
2898  if (result!=NULL) result= pAdd( result, p );
2899  else result= p;
2900  }
2901  c++;
2902  }
2903  sum=0;
2904  exp[0]++;
2905  for ( j= 0; j < n - 1; j++ )
2906  {
2907  if ( exp[j] > tdg )
2908  {
2909  exp[j]= 0;
2910  exp[j + 1]++;
2911  }
2912  sum+=exp[j];
2913  }
2914  sum+=exp[n-1];
2915  }
2916 
2917  pTest( result );
2918 
2919  return result;
2920 }
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
int l
Definition: cfEzgcd.cc:100
int p
Definition: cfModGcd.cc:4078
Definition: intvec.h:23
virtual number getDetAt(const number *)
Definition: mpr_base.h:36
virtual long getDetDeg()
Definition: mpr_base.h:39
int nextPrime(const int p)
Definition: mpr_base.cc:3173
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:29
#define Print
Definition: emacs.cc:80
return result
Definition: facAbsBiFact.cc:75
int j
Definition: facHensel.cc:110
unsigned long over(const unsigned long n, const unsigned long d)
Definition: mpr_base.cc:2659
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:357
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
#define mprPROTN(msg, nval)
Definition: mpr_global.h:48
#define ST_BASE_EV
Definition: mpr_global.h:62
#define mprPROTnl(msg)
Definition: mpr_global.h:42
#define mprPROTNnl(msg, nval)
Definition: mpr_global.h:49
#define nDiv(a, b)
Definition: numbers.h:32
#define nDelete(n)
Definition: numbers.h:16
#define nIsZero(n)
Definition: numbers.h:19
#define nEqual(n1, n2)
Definition: numbers.h:20
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:46
#define nNormalize(n)
Definition: numbers.h:30
#define nInit(i)
Definition: numbers.h:24
#define nTest(a)
Definition: numbers.h:35
#define nPower(a, b, res)
Definition: numbers.h:38
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define NULL
Definition: omList.c:12
#define pAdd(p, q)
Definition: polys.h:203
#define pTest(p)
Definition: polys.h:414
#define pSetm(p)
Definition: polys.h:271
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pOne()
Definition: polys.h:315
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310

◆ interpolateDenseSP()

rootContainer ** uResultant::interpolateDenseSP ( BOOLEAN  matchUp = false,
const number  subDetVal = NULL 
)

Definition at line 2922 of file mpr_base.cc.

2923 {
2924  int i,p,uvar;
2925  long tdg;
2926  int loops= (matchUp?n-2:n-1);
2927 
2928  mprPROTnl("uResultant::interpolateDenseSP");
2929 
2930  tdg= resMat->getDetDeg();
2931 
2932  // evaluate D in tdg+1 distinct points, so
2933  // we need tdg+1 results of D(p0,1,0,...,0) =
2934  // c(0)*u0^tdg + c(1)*u0^tdg-1 + ... + c(tdg-1)*u0 + c(tdg)
2935  number *presults;
2936  presults= (number *)omAlloc( (tdg + 1) * sizeof( number ) );
2937  for ( i=0; i <= tdg; i++ ) presults[i]= nInit(0);
2938 
2939  rootContainer ** roots;
2940  roots= (rootContainer **) omAlloc( loops * sizeof(rootContainer*) );
2941  for ( i=0; i < loops; i++ ) roots[i]= new rootContainer(); // 0..n-2
2942 
2943  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
2944  for (i=0; i < n; i++) pevpoint[i]= nInit(0);
2945 
2946  number *pev= (number *)omAlloc( n * sizeof( number ) );
2947  for (i=0; i < n; i++) pev[i]= nInit(0);
2948 
2949  // now we evaluate D(u0,-1,0,...0), D(u0,0,-1,0,...,0), ..., D(u0,0,..,0,-1)
2950  // or D(u0,k1,k2,0,...,0), D(u0,k1,k2,k3,0,...,0), ..., D(u0,k1,k2,k3,...,kn)
2951  // this gives us n-1 evaluations
2952  p=3;
2953  for ( uvar= 0; uvar < loops; uvar++ )
2954  {
2955  // generate initial evaluation point
2956  if ( matchUp )
2957  {
2958  for (i=0; i < n; i++)
2959  {
2960  // prime(random number) between 1 and MAXEVPOINT
2961  nDelete( &pevpoint[i] );
2962  if ( i == 0 )
2963  {
2964  //p= nextPrime( p );
2965  pevpoint[i]= nInit( p );
2966  }
2967  else if ( i <= uvar + 2 )
2968  {
2969  pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
2970  //pevpoint[i]=nInit(383);
2971  }
2972  else
2973  pevpoint[i]=nInit(0);
2974  mprPROTNnl(" ",pevpoint[i]);
2975  }
2976  }
2977  else
2978  {
2979  for (i=0; i < n; i++)
2980  {
2981  // init pevpoint with prime,0,...0,1,0,...,0
2982  nDelete( &pevpoint[i] );
2983  if ( i == 0 )
2984  {
2985  //p=nextPrime( p );
2986  pevpoint[i]=nInit( p );
2987  }
2988  else
2989  {
2990  if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
2991  else pevpoint[i]= nInit(0);
2992  }
2993  mprPROTNnl(" ",pevpoint[i]);
2994  }
2995  }
2996 
2997  // prepare aktual evaluation point
2998  for (i=0; i < n; i++)
2999  {
3000  nDelete( &pev[i] );
3001  pev[i]= nCopy( pevpoint[i] );
3002  }
3003  // evaluate the determinant in the points pev^0, pev^1, ..., pev^tdg
3004  for ( i=0; i <= tdg; i++ )
3005  {
3006  nDelete( &pev[0] );
3007  nPower(pevpoint[0],i,&pev[0]); // new evpoint
3008 
3009  nDelete( &presults[i] );
3010  presults[i]=resMat->getDetAt( pev ); // evaluate det at point evpoint
3011 
3012  mprPROTNnl("",presults[i]);
3013 
3015  mprPROTL("",tdg-i);
3016  }
3017  mprSTICKYPROT("\n");
3018 
3019  // now interpolate
3020  vandermonde vm( tdg + 1, 1, tdg, pevpoint, FALSE );
3021  number *ncpoly= vm.interpolateDense( presults );
3022 
3023  if ( subDetVal != NULL )
3024  { // divide by common factor
3025  number detdiv;
3026  for ( i= 0; i <= tdg; i++ )
3027  {
3028  detdiv= nDiv( ncpoly[i], subDetVal );
3029  nNormalize( detdiv );
3030  nDelete( &ncpoly[i] );
3031  ncpoly[i]= detdiv;
3032  }
3033  }
3034 
3035 #ifdef mprDEBUG_ALL
3036  PrintLn();
3037  for ( i=0; i <= tdg; i++ )
3038  {
3039  nPrint(ncpoly[i]); PrintS(" --- ");
3040  }
3041  PrintLn();
3042 #endif
3043 
3044  // save results
3045  roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
3047  loops );
3048  }
3049 
3050  // free some stuff: pev, presult
3051  for ( i=0; i < n; i++ ) nDelete( pev + i );
3052  omFreeSize( (void *)pev, n * sizeof( number ) );
3053 
3054  for ( i=0; i <= tdg; i++ ) nDelete( presults+i );
3055  omFreeSize( (void *)presults, (tdg + 1) * sizeof( number ) );
3056 
3057  return roots;
3058 }
#define FALSE
Definition: auxiliary.h:96
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:66
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:300
#define MAXEVPOINT
Definition: mpr_base.cc:2653
#define mprPROTL(msg, intval)
Definition: mpr_global.h:46
#define nCopy(n)
Definition: numbers.h:15
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
int siRand()
Definition: sirandom.c:42

◆ linearPoly()

poly uResultant::linearPoly ( const resMatType  rmt)
private

Definition at line 2743 of file mpr_base.cc.

2744 {
2745  int i;
2746 
2747  poly newlp= pOne();
2748  poly actlp, rootlp= newlp;
2749 
2750  for ( i= 1; i <= (currRing->N); i++ )
2751  {
2752  actlp= newlp;
2753  pSetExp( actlp, i, 1 );
2754  pSetm( actlp );
2755  newlp= pOne();
2756  actlp->next= newlp;
2757  }
2758  actlp->next= NULL;
2759  pDelete( &newlp );
2760 
2761  if ( rrmt == sparseResMat )
2762  {
2763  newlp= pOne();
2764  actlp->next= newlp;
2765  newlp->next= NULL;
2766  }
2767  return ( rootlp );
2768 }
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define pDelete(p_ptr)
Definition: polys.h:186

◆ nextPrime()

int uResultant::nextPrime ( const int  p)
private

Definition at line 3173 of file mpr_base.cc.

3174 {
3175  int init=i;
3176  int ii=i+2;
3177  extern int IsPrime(int p); // from Singular/ipshell.{h,cc}
3178  int j= IsPrime( ii );
3179  while ( j <= init )
3180  {
3181  ii+=2;
3182  j= IsPrime( ii );
3183  }
3184  return j;
3185 }
void init()
Definition: lintree.cc:864
int IsPrime(int p)
Definition: prime.cc:61

◆ specializeInU()

rootContainer ** uResultant::specializeInU ( BOOLEAN  matchUp = false,
const number  subDetVal = NULL 
)

Definition at line 3060 of file mpr_base.cc.

3061 {
3062  int i,/*p,*/uvar;
3063  long tdg;
3064  poly pures,piter;
3065  int loops=(matchUp?n-2:n-1);
3066  int nn=n;
3067  if (loops==0) { loops=1;nn++;}
3068 
3069  mprPROTnl("uResultant::specializeInU");
3070 
3071  tdg= resMat->getDetDeg();
3072 
3073  rootContainer ** roots;
3074  roots= (rootContainer **) omAlloc( loops * sizeof(rootContainer*) );
3075  for ( i=0; i < loops; i++ ) roots[i]= new rootContainer(); // 0..n-2
3076 
3077  number *pevpoint= (number *)omAlloc( nn * sizeof( number ) );
3078  for (i=0; i < nn; i++) pevpoint[i]= nInit(0);
3079 
3080  // now we evaluate D(u0,-1,0,...0), D(u0,0,-1,0,...,0), ..., D(u0,0,..,0,-1)
3081  // or D(u0,k1,k2,0,...,0), D(u0,k1,k2,k3,0,...,0), ..., D(u0,k1,k2,k3,...,kn)
3082  // p=3;
3083  for ( uvar= 0; uvar < loops; uvar++ )
3084  {
3085  // generate initial evaluation point
3086  if ( matchUp )
3087  {
3088  for (i=0; i < n; i++)
3089  {
3090  // prime(random number) between 1 and MAXEVPOINT
3091  nDelete( &pevpoint[i] );
3092  if ( i <= uvar + 2 )
3093  {
3094  pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
3095  //pevpoint[i]=nInit(383);
3096  }
3097  else pevpoint[i]=nInit(0);
3098  mprPROTNnl(" ",pevpoint[i]);
3099  }
3100  }
3101  else
3102  {
3103  for (i=0; i < n; i++)
3104  {
3105  // init pevpoint with prime,0,...0,-1,0,...,0
3106  nDelete( &(pevpoint[i]) );
3107  if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
3108  else pevpoint[i]= nInit(0);
3109  mprPROTNnl(" ",pevpoint[i]);
3110  }
3111  }
3112 
3113  pures= resMat->getUDet( pevpoint );
3114 
3115  number *ncpoly= (number *)omAlloc( (tdg+1) * sizeof(number) );
3116 
3117 #ifdef MPR_MASI
3118  BOOLEAN masi=true;
3119 #endif
3120 
3121  piter= pures;
3122  for ( i= tdg; i >= 0; i-- )
3123  {
3124  if ( piter && pTotaldegree(piter) == i )
3125  {
3126  ncpoly[i]= nCopy( pGetCoeff( piter ) );
3127  pIter( piter );
3128 #ifdef MPR_MASI
3129  masi=false;
3130 #endif
3131  }
3132  else
3133  {
3134  ncpoly[i]= nInit(0);
3135  }
3136  mprPROTNnl("", ncpoly[i] );
3137  }
3138 #ifdef MPR_MASI
3139  if ( masi ) mprSTICKYPROT("MASI");
3140 #endif
3141 
3142  mprSTICKYPROT(ST_BASE_EV); // .
3143 
3144  if ( subDetVal != NULL ) // divide by common factor
3145  {
3146  number detdiv;
3147  for ( i= 0; i <= tdg; i++ )
3148  {
3149  detdiv= nDiv( ncpoly[i], subDetVal );
3150  nNormalize( detdiv );
3151  nDelete( &ncpoly[i] );
3152  ncpoly[i]= detdiv;
3153  }
3154  }
3155 
3156  pDelete( &pures );
3157 
3158  // save results
3159  roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
3161  loops );
3162  }
3163 
3164  mprSTICKYPROT("\n");
3165 
3166  // free some stuff: pev, presult
3167  for ( i=0; i < n; i++ ) nDelete( pevpoint + i );
3168  omFreeSize( (void *)pevpoint, n * sizeof( number ) );
3169 
3170  return roots;
3171 }
int BOOLEAN
Definition: auxiliary.h:87
virtual poly getUDet(const number *)
Definition: mpr_base.h:34
#define pIter(p)
Definition: monomials.h:37
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
static long pTotaldegree(poly p)
Definition: polys.h:282

Field Documentation

◆ gls

ideal uResultant::gls
private

Definition at line 88 of file mpr_base.h.

◆ n

int uResultant::n
private

Definition at line 89 of file mpr_base.h.

◆ resMat

resMatrixBase* uResultant::resMat
private

Definition at line 92 of file mpr_base.h.

◆ rmt

resMatType uResultant::rmt
private

Definition at line 91 of file mpr_base.h.


The documentation for this class was generated from the following files: