Miroirs/newrpl
1
0
Fork 0
mirror of https://git.code.sf.net/p/newrpl/sources synced 2024-11-16 19:51:25 +01:00
Code Issues Projects Releases Packages Wiki Activity Actions

Added PROOT multiplicity, better error management.

Browse source
This commit is contained in:
claudiol 2017-10-20 17:48:14 -04:00
parent 2eb26c0a67
commit 704742c66a
2 changed files with 199 additions and 54 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

158
newrpl/lib-64-arithmetic.c
View file

@ -1692,16 +1692,17 @@ case IPPOST:
WORDPTR vect_val=rplPeekData(2);
WORDPTR real_val=rplPeekData(1);
/*
if(ISLIST(*real_val)) {
rplListBinaryDoCmd();
return;
}
else if((ISIDENT(*vect_val) || ISSYMBOLIC(*vect_val)) || (ISIDENT(*real_val) || ISSYMBOLIC(*real_val))){
else if((ISIDENT(*vect_val) || ISSYMBOLIC(*vect_val)))
{
rplSymbApplyOperator(CurOpcode,2);
return;
}
else */ if(ISMATRIX(*vect_val) && ISNUMBER(*real_val) ){
else if(ISMATRIX(*vect_val) ){
BINT rows=MATROWS(vect_val[1]),cols=MATCOLS(vect_val[1]);
@ -1709,48 +1710,24 @@ case IPPOST:
rplError(ERR_VECTOREXPECTED);
return;
}
BINT f;
for(f=0;f<cols;++f) {
WORDPTR entry=rplMatrixFastGet(vect_val,1,f+1);
if(!ISNUMBER(*entry)) {
rplError(ERR_VECTOROFNUMBERSEXPECTED);
return;
}
WORDPTR *savestk=DSTop;
rplPushData(rplPeekData(2));
if(Exceptions) { DSTop=savestk; return; }
WORDPTR *first=rplMatrixExplode();
if(!first || Exceptions) { DSTop=savestk; return; }
WORDPTR result=rplPolyEvalEx(first,cols-1,savestk-1);
if(!result || Exceptions) { DSTop=savestk; return; }
if(ISSYMBOLIC(*result)) {
result=rplSymbNumericReduce(result);
if(!result || Exceptions) { DSTop=savestk; return; }
}
// DO IT ALL WITH REALS
BINT saveprec=Context.precdigits;
BINT argdigits=(cols+7)&~7;
if(argdigits>MAX_USERPRECISION) {
argdigits=MAX_USERPRECISION;
}
argdigits=(argdigits>Context.precdigits)? argdigits:Context.precdigits;
// AUTOMATICALLY INCREASE PRECISION TEMPORARILY
Context.precdigits=argdigits;
// use Horner scheme
rplNumberToRReg(0,real_val);
rplNumberToRReg(1,(WORDPTR)zero_bint);
for(f=0;f<cols;++f) {
WORDPTR entry=rplMatrixFastGet(vect_val,1,f+1);
rplNumberToRReg(3,entry);
mulReal(&RReg[2], &RReg[1], &RReg[0]);
addReal(&RReg[1], &RReg[2], &RReg[3]);
}
Context.precdigits=saveprec;
WORDPTR newnumber=rplNewReal(&RReg[1]);
if(!newnumber) return;
// drop one value and replace level 1 value
rplDropData(1);
rplOverwriteData(1,newnumber);
DSTop=savestk-1;
rplOverwriteData(1,result);
return;
@ -3078,20 +3055,20 @@ case PROOT:
WORDPTR vect_val=rplPeekData(1);
/*
if(ISLIST(*real_val)) {
rplListBinaryDoCmd();
if(ISLIST(*vect_val)) {
rplListUnaryDoCmd();
return;
}
else if((ISIDENT(*vect_val) || ISSYMBOLIC(*vect_val)) || (ISIDENT(*real_val) || ISSYMBOLIC(*real_val))){
rplSymbApplyOperator(CurOpcode,2);
return;
}
else */ if(ISMATRIX(*vect_val)){
else if(ISMATRIX(*vect_val)){
BINT cplxmode=rplTestSystemFlag(FL_COMPLEXMODE);
rplSetSystemFlag(FL_COMPLEXMODE);
BINT rows=MATROWS(vect_val[1]),cols=MATCOLS(vect_val[1]);
if(rows) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
rplError(ERR_VECTOREXPECTED);
return;
}
@ -3101,6 +3078,7 @@ case PROOT:
for(f=0;f<cols;++f) {
WORDPTR entry=rplMatrixFastGet(vect_val,1,f+1);
if(!ISNUMBERCPLX(*entry)) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
rplError(ERR_VECTOROFNUMBERSEXPECTED);
return;
@ -3109,10 +3087,53 @@ case PROOT:
}
WORDPTR solution;
REAL re;
for(f=1;f<cols-1;++f) {
if(f>1) {
// TEST IF PREVIOUS SOLUTION IS ALSO A SOLUTION OF THE DEFLATED POLYNOMIAL
solution=rplPolyEvalEx(savestk,cols-f,DSTop-1);
if(!solution || Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
rplPushData(solution);
rplCallOvrOperator(CMD_OVR_ABS);
if(Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
rplReadNumberAsReal(rplPopData(),&re);
if(Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
if(iszeroReal(&re) || (intdigitsReal(&re)<-Context.precdigits)) {
// THE CURRENT ROOT HAS MULTIPLICITY
rplPushData(rplPeekData(1));
// DEFLATE THE POLYNOMIAL
rplPolyDeflateEx(savestk,cols-f,DSTop-1);
if(Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
continue;
}
}
solution=rplPolyRootEx(savestk,cols-f);
if(Exceptions) {
if(!solution || Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
@ -3120,12 +3141,14 @@ case PROOT:
// WE HAVE ONE SOLUTION!
rplPushData(solution);
if(Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
// DEFLATE THE POLYNOMIAL
rplPolyDeflateEx(savestk,cols-f,DSTop-1);
if(Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
@ -3138,17 +3161,46 @@ case PROOT:
// THEREFORE THE LAST ROOT IS X=-a1/a0
rplPushData(savestk[1]);
rplPushData(savestk[0]);
rplCallOvrOperator(CMD_OVR_DIV);
rplCallOvrOperator(CMD_OVR_NEG);
if(Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
rplCallOvrOperator(CMD_OVR_DIV);
if(Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
rplCallOvrOperator(CMD_OVR_NEG);
if(Exceptions) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
BINT doerror=0;
if(!cplxmode) {
// ISSUE AN ERROR IF ANY OF THE ROOTS ARE COMPLEX
for(f=1;f<cols;++f) {
if(ISCOMPLEX(*rplPeekData(f))) { doerror=1; break; }
}
}
solution=rplMatrixCompose(0,cols-1);
if(!solution) {
if(!cplxmode) rplClrSystemFlag(FL_COMPLEXMODE);
DSTop=savestk;
return;
}
DSTop=savestk;
rplOverwriteData(1,solution);
if(!cplxmode) {
rplClrSystemFlag(FL_COMPLEXMODE);
if(doerror) rplError(ERR_COMPLEXRESULT);
}
return;

95
newrpl/solvers.c
View file

@ -254,8 +254,100 @@ WORDPTR rplPolyRootEx(WORDPTR *first,BINT degree)
} while(1); // ADD SOME LOOP LIMIT HERE IN CASE IT DOESN'T CONVERGE!
Context.precdigits=oldprec;
// FOUND A ROOT, CLEAN IT UP BEFORE SENDING
pk=rplPeekData(1);
if(ISCOMPLEX(*pk)) {
BINT cclass=rplComplexClass(pk);
if(cclass==CPLX_ZERO) rplOverwriteData(1,(WORDPTR)zero_bint);
else if(cclass==CPLX_NORMAL) {
REAL re,im;
rplReadCNumberAsReal(pk,&re);
rplReadCNumberAsImag(pk,&im);
BINT digre,digim;
digre=intdigitsReal(&re);
digim=intdigitsReal(&im);
if( (digim<-10) && (digre>digim+10)) {
// IF IM(r)<1E-10 AND RE(r)>10^10*IM(r)
// TRY A REAL ROOT AND SEE IF IT'S STILL GOOD
rplNewRealPush(&re);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
pk=rplPolyEvalEx(first,degree,DSTop-1);
if(!pk) { Context.precdigits=oldprec; return 0; }
rplPushData(pk);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
rplCallOvrOperator(CMD_OVR_ABS);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
rplReadNumberAsReal(rplPopData(),&err);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
if(iszeroReal(&err) || (intdigitsReal(&err)<-(2*oldprec)))
{ pk=rplPopData(); rplOverwriteData(1,pk); } // REAL ROOT ALONE IS STILL GOOD, USE IT
}
else
if( (digre<-10) && (digim>digre+10)) {
// IF RE(r)<1E-10 AND IM()>10^10*RE()
// TRY A PURE IMAGINARY ROOT AND SEE IF IT'S STILL GOOD
rplZeroToRReg(0);
rplNewComplexPush(&RReg[0],&im,ANGLENONE);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
pk=rplPolyEvalEx(first,degree,DSTop-1);
if(!pk) { Context.precdigits=oldprec; return 0; }
rplPushData(pk);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
rplCallOvrOperator(CMD_OVR_ABS);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
rplReadNumberAsReal(rplPopData(),&err);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
if(iszeroReal(&err) || (intdigitsReal(&err)<-(2*oldprec)))
{ pk=rplPopData(); rplOverwriteData(1,pk); } // IMAG. ROOT ALONE IS STILL GOOD, USE IT
}
}
}
if(ISNUMBER(*pk)) {
// ADJUST ROOTS TO BE INTEGER IF THEY ARE CLOSE TO IT
REAL re;
rplReadNumberAsReal(pk,&re);
fracReal(&RReg[0],&re);
if(!iszeroReal(&RReg[0]) && (intdigitsReal(&RReg[0])<-10)) {
ipReal(&RReg[1],&re,0);
rplNewRealPush(&RReg[1]);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
pk=rplPolyEvalEx(first,degree,DSTop-1);
if(!pk) { Context.precdigits=oldprec; return 0; }
rplPushData(pk);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
rplCallOvrOperator(CMD_OVR_ABS);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
rplReadNumberAsReal(rplPopData(),&err);
if(Exceptions) { Context.precdigits=oldprec; return 0; }
if(iszeroReal(&err) || (intdigitsReal(&err)<-(2*oldprec)))
{ pk=rplPopData(); rplOverwriteData(1,pk); } // INTEGER ROOT ALONE IS STILL GOOD, USE IT
}
}
Context.precdigits=oldprec;
rplDropData(2);
return pk;
}
@ -292,3 +384,4 @@ WORDPTR rplPolyDeflateEx(WORDPTR *first,BINT degree,WORDPTR *value)
return rplPopData();
}