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 2684 of file mpr_base.cc.

2685  : rmt( _rmt )
2686 {
2687  if ( extIdeal )
2688  {
2689  // extend given ideal by linear poly F0=u0x0 + u1x1 +...+ unxn
2690  gls= extendIdeal( _gls, linearPoly( rmt ), rmt );
2691  n= IDELEMS( gls );
2692  }
2693  else
2694  gls= idCopy( _gls );
2695 
2696  switch ( rmt )
2697  {
2698  case sparseResMat:
2699  resMat= new resMatrixSparse( gls );
2700  break;
2701  case denseResMat:
2702  resMat= new resMatrixDense( gls );
2703  break;
2704  default:
2705  WerrorS("uResultant::uResultant: Unknown chosen resultant matrix type!");
2706  }
2707 }
poly linearPoly(const resMatType rmt)
Definition: mpr_base.cc:2742
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:2714
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 2709 of file mpr_base.cc.

2710 {
2711  delete resMat;
2712 }

◆ 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 2714 of file mpr_base.cc.

2715 {
2716  ideal newGls= idCopy( igls );
2717  newGls->m= (poly *)omReallocSize( newGls->m,
2718  IDELEMS(igls) * sizeof(poly),
2719  (IDELEMS(igls) + 1) * sizeof(poly) );
2720  IDELEMS(newGls)++;
2721 
2722  switch ( rrmt )
2723  {
2724  case sparseResMat:
2725  case denseResMat:
2726  {
2727  int i;
2728  for ( i= IDELEMS(newGls)-1; i > 0; i-- )
2729  {
2730  newGls->m[i]= newGls->m[i-1];
2731  }
2732  newGls->m[0]= linPoly;
2733  }
2734  break;
2735  default:
2736  WerrorS("uResultant::extendIdeal: Unknown chosen resultant matrix type!");
2737  }
2738 
2739  return( newGls );
2740 }
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 2769 of file mpr_base.cc.

2770 {
2771  int i,j,p;
2772  long tdg;
2773 
2774  // D is a Polynom homogeneous in the coeffs of F0 and of degree tdg = d0*d1*...*dn
2775  tdg= resMat->getDetDeg();
2776 
2777  // maximum number of terms in polynom D (homogeneous, of degree tdg)
2778  // long mdg= (facul(tdg+n-1) / facul( tdg )) / facul( n - 1 );
2779  long mdg= over( n-1, tdg );
2780 
2781  // maximal number of terms in a polynom of degree tdg
2782  long l=(long)pow( (double)(tdg+1), n );
2783 
2784 #ifdef mprDEBUG_PROT
2785  Print("// total deg of D: tdg %ld\n",tdg);
2786  Print("// maximum number of terms in D: mdg: %ld\n",mdg);
2787  Print("// maximum number of terms in polynom of deg tdg: l %ld\n",l);
2788 #endif
2789 
2790  // we need mdg results of D(p0,p1,...,pn)
2791  number *presults;
2792  presults= (number *)omAlloc( mdg * sizeof( number ) );
2793  for (i=0; i < mdg; i++) presults[i]= nInit(0);
2794 
2795  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
2796  number *pev= (number *)omAlloc( n * sizeof( number ) );
2797  for (i=0; i < n; i++) pev[i]= nInit(0);
2798 
2799  mprPROTnl("// initial evaluation point: ");
2800  // initial evaluatoin point
2801  p=1;
2802  for (i=0; i < n; i++)
2803  {
2804  // init pevpoint with primes 3,5,7,11, ...
2805  p= nextPrime( p );
2806  pevpoint[i]=nInit( p );
2807  nTest(pevpoint[i]);
2808  mprPROTNnl(" ",pevpoint[i]);
2809  }
2810 
2811  // evaluate the determinant in the points pev^0, pev^1, ..., pev^mdg
2812  mprPROTnl("// evaluating:");
2813  for ( i=0; i < mdg; i++ )
2814  {
2815  for (j=0; j < n; j++)
2816  {
2817  nDelete( &pev[j] );
2818  nPower(pevpoint[j],i,&pev[j]);
2819  mprPROTN(" ",pev[j]);
2820  }
2821  mprPROTnl("");
2822 
2823  nDelete( &presults[i] );
2824  presults[i]=resMat->getDetAt( pev );
2825 
2827  }
2828  mprSTICKYPROT("\n");
2829 
2830  // now interpolate using vandermode interpolation
2831  mprPROTnl("// interpolating:");
2832  number *ncpoly;
2833  {
2834  vandermonde vm( mdg, n, tdg, pevpoint );
2835  ncpoly= vm.interpolateDense( presults );
2836  }
2837 
2838  if ( subDetVal != NULL )
2839  { // divide by common factor
2840  number detdiv;
2841  for ( i= 0; i <= mdg; i++ )
2842  {
2843  detdiv= nDiv( ncpoly[i], subDetVal );
2844  nNormalize( detdiv );
2845  nDelete( &ncpoly[i] );
2846  ncpoly[i]= detdiv;
2847  }
2848  }
2849 
2850 #ifdef mprDEBUG_ALL
2851  PrintLn();
2852  for ( i=0; i < mdg; i++ )
2853  {
2854  nPrint(ncpoly[i]); PrintS(" --- ");
2855  }
2856  PrintLn();
2857 #endif
2858 
2859  // prepare ncpoly for later use
2860  number nn=nInit(0);
2861  for ( i=0; i < mdg; i++ )
2862  {
2863  if ( nEqual(ncpoly[i],nn) )
2864  {
2865  nDelete( &ncpoly[i] );
2866  ncpoly[i]=NULL;
2867  }
2868  }
2869  nDelete( &nn );
2870 
2871  // create poly presenting the determinat of the uResultant
2872  intvec exp( n );
2873  for ( i= 0; i < n; i++ ) exp[i]=0;
2874 
2875  poly result= NULL;
2876 
2877  long sum=0;
2878  long c=0;
2879 
2880  for ( i=0; i < l; i++ )
2881  {
2882  if ( sum == tdg )
2883  {
2884  if ( !nIsZero(ncpoly[c]) )
2885  {
2886  poly p= pOne();
2887  if ( rmt == denseResMat )
2888  {
2889  for ( j= 0; j < n; j++ ) pSetExp( p, j+1, exp[j] );
2890  }
2891  else if ( rmt == sparseResMat )
2892  {
2893  for ( j= 1; j < n; j++ ) pSetExp( p, j, exp[j] );
2894  }
2895  pSetCoeff( p, ncpoly[c] );
2896  pSetm( p );
2897  if (result!=NULL) result= pAdd( result, p );
2898  else result= p;
2899  }
2900  c++;
2901  }
2902  sum=0;
2903  exp[0]++;
2904  for ( j= 0; j < n - 1; j++ )
2905  {
2906  if ( exp[j] > tdg )
2907  {
2908  exp[j]= 0;
2909  exp[j + 1]++;
2910  }
2911  sum+=exp[j];
2912  }
2913  sum+=exp[n-1];
2914  }
2915 
2916  pTest( result );
2917 
2918  return result;
2919 }
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:3172
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:2658
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:415
#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 2921 of file mpr_base.cc.

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

2743 {
2744  int i;
2745 
2746  poly newlp= pOne();
2747  poly actlp, rootlp= newlp;
2748 
2749  for ( i= 1; i <= (currRing->N); i++ )
2750  {
2751  actlp= newlp;
2752  pSetExp( actlp, i, 1 );
2753  pSetm( actlp );
2754  newlp= pOne();
2755  actlp->next= newlp;
2756  }
2757  actlp->next= NULL;
2758  pDelete( &newlp );
2759 
2760  if ( rrmt == sparseResMat )
2761  {
2762  newlp= pOne();
2763  actlp->next= newlp;
2764  newlp->next= NULL;
2765  }
2766  return ( rootlp );
2767 }
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 3172 of file mpr_base.cc.

3173 {
3174  int init=i;
3175  int ii=i+2;
3176  extern int IsPrime(int p); // from Singular/ipshell.{h,cc}
3177  int j= IsPrime( ii );
3178  while ( j <= init )
3179  {
3180  ii+=2;
3181  j= IsPrime( ii );
3182  }
3183  return j;
3184 }
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 3059 of file mpr_base.cc.

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