My Project
Loading...
Searching...
No Matches
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.

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
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: