diff --git a/check/features.frm b/check/features.frm index 04be478b..03d40d51 100644 --- a/check/features.frm +++ b/check/features.frm @@ -617,15 +617,15 @@ Print; .end #pend_if wordsize == 2 assert result("PI0") =~ expr("1.0e+00*a*b") -assert result("PI1") =~ expr("3.14159265358979323846e+00") -assert result("PI2") =~ expr("9.42477796076937971538e+00") -assert result("PI3") =~ expr("5.4413980927026535518e+00") -assert result("PI4") =~ expr("3.10062766802998201755e+01") -assert result("PI5") =~ expr("3.14159265358979323846e+00*a*b") -assert result("PI6") =~ expr("1.77245385090551602731e+00") -assert result("EE1") =~ expr("2.71828182845904523537e+00") -assert result("EE2") =~ expr("8.53973422267356706549e+00") -assert result("EM1") =~ expr("5.77215664901532860607e-01") +assert result("PI1") =~ expr("3.141592653589793238e+00") +assert result("PI2") =~ expr("9.424777960769379715e+00") +assert result("PI3") =~ expr("5.441398092702653552e+00") +assert result("PI4") =~ expr("3.100627668029982018e+01") +assert result("PI5") =~ expr("3.141592653589793238e+00*a*b") +assert result("PI6") =~ expr("1.772453850905516027e+00") +assert result("EE1") =~ expr("2.718281828459045235e+00") +assert result("EE2") =~ expr("8.539734222673567065e+00") +assert result("EM1") =~ expr("5.772156649015328606e-01") *--#] evaluate_symbol : *--#[ evaluate_symbol_pi : #- @@ -641,7 +641,7 @@ Evaluate pi_; Print; .end #pend_if wordsize == 2 -assert result("PI") =~ expr("3.141592653589793238462643383279502884198e+00") +assert result("PI") =~ expr("3.1415926535897932384626433832795028842e+00") assert result("EE") =~ expr("1.0e+00*ee_") assert result("EM") =~ expr("1.0e+00*em_") *--#] evaluate_symbol_pi : @@ -660,7 +660,7 @@ Print; .end #pend_if wordsize == 2 assert result("PI") =~ expr("1.0e+00*pi_") -assert result("EE") =~ expr("2.718281828459045235360287471352662497757247093699959574967e+00") +assert result("EE") =~ expr("2.7182818284590452353602874713526624977572470937e+00") assert result("EM") =~ expr("1.0e+00*em_") *--#] evaluate_symbol_ee : *--#[ evaluate_symbol_em : @@ -679,8 +679,115 @@ Print; #pend_if wordsize == 2 assert result("PI") =~ expr("1.0e+00*pi_") assert result("EE") =~ expr("1.0e+00*ee_") -assert result("EM") =~ expr("5.7721566490153286060651209008240243104215933593992359880577e-01") +assert result("EM") =~ expr("5.77215664901532860606512090082402431042159335939923598806e-01") *--#] evaluate_symbol_em : +*--#[ evaluate_mzv_2 : +#- +L F = mzv_(2); +.sort + +Hide; +#do bits=5,65,3 + #StartFloat `bits',2 + Local F`bits' = F; + Evaluate mzv_; + Print; + .sort + Hide; + #endfloat +#enddo +.end +#pend_if wordsize == 2 +assert result("F5") =~ expr("2e+00") +assert result("F8") =~ expr("1.6e+00") +assert result("F11") =~ expr("1.64e+00") +assert result("F14") =~ expr("1.645e+00") +assert result("F17") =~ expr("1.6449e+00") +assert result("F20") =~ expr("1.64493e+00") +assert result("F23") =~ expr("1.64493e+00") +assert result("F26") =~ expr("1.644934e+00") +assert result("F29") =~ expr("1.6449341e+00") +assert result("F32") =~ expr("1.64493407e+00") +assert result("F35") =~ expr("1.644934067e+00") +assert result("F38") =~ expr("1.6449340668e+00") +assert result("F41") =~ expr("1.64493406685e+00") +assert result("F44") =~ expr("1.644934066848e+00") +assert result("F47") =~ expr("1.6449340668482e+00") +assert result("F50") =~ expr("1.64493406684823e+00") +assert result("F53") =~ expr("1.64493406684823e+00") +assert result("F56") =~ expr("1.644934066848226e+00") +assert result("F59") =~ expr("1.6449340668482264e+00") +assert result("F62") =~ expr("1.64493406684822644e+00") +assert result("F65") =~ expr("1.644934066848226436e+00") +*--#] evaluate_mzv_2 : +*--#[ evaluate_all_mzv_2-6 : +#- +Symbol a,n,x,jj; +CFunction mzv; +#do weight=2,6 + L F`weight' = x^`weight'*mzv(); +#enddo +* Generate all possible arguments +repeat id x^n?{>0}*mzv(?a) = sum_(jj,1,n, x^(n-jj)*mzv(?a,jj)); +* Only keep convergent MZVs +id mzv(1,?a) = 0; +.sort + +Hide; +#do weight=2,6 + #StartFloat 74,`weight' + Local MZV`weight' = F`weight'; + id mzv(?a) = mzv(?a)*mzv_(?a); + Evaluate mzv_; + Print +s; + .sort + Hide; + #endfloat +#enddo +.end +#pend_if wordsize == 2 +assert result("MZV2") =~ expr(" + + 1.644934066848226436472e+00*mzv(2) +") +assert result("MZV3") =~ expr(" + + 1.2020569031595942854e+00*mzv(2,1) + + 1.2020569031595942854e+00*mzv(3) +") +assert result("MZV4") =~ expr(" + + 1.082323233711138191516e+00*mzv(2,1,1) + + 8.11742425283353643637e-01*mzv(2,2) + + 2.70580808427784547879e-01*mzv(3,1) + + 1.082323233711138191516e+00*mzv(4) +") +assert result("MZV5") =~ expr(" + + 1.036927755143369926331e+00*mzv(2,1,1,1) + + 7.11566197550572432097e-01*mzv(2,1,2) + + 2.288103976033537597687e-01*mzv(2,2,1) + + 7.11566197550572432097e-01*mzv(2,3) + + 9.655115998944373446565e-02*mzv(3,1,1) + + 2.288103976033537597687e-01*mzv(3,2) + + 9.655115998944373446565e-02*mzv(4,1) + + 1.036927755143369926331e+00*mzv(5) +") +assert result("MZV6") =~ expr(" + + 1.017343061984449139715e+00*mzv(2,1,1,1,1) + + 6.745239140339681404916e-01*mzv(2,1,1,2) + + 2.137988682245925470996e-01*mzv(2,1,2,1) + + 6.183495605712693078956e-01*mzv(2,1,3) + + 8.848338245436871429433e-02*mzv(2,2,1,1) + + 1.907518241220842136965e-01*mzv(2,2,2) + + 7.922139756520716599903e-02*mzv(2,3,1) + + 6.745239140339681404916e-01*mzv(2,4) + + 4.053689727151973782905e-02*mzv(3,1,1,1) + + 7.922139756520716599903e-02*mzv(3,1,2) + + 3.230902899166988169841e-02*mzv(3,2,1) + + 2.137988682245925470996e-01*mzv(3,3) + + 1.748985316901140442593e-02*mzv(4,1,1) + + 8.848338245436871429433e-02*mzv(4,2) + + 4.053689727151973782905e-02*mzv(5,1) + + 1.017343061984449139715e+00*mzv(6) +") +*--#] evaluate_all_mzv_2-6 : *--#[ Issue49 : * Add mul_ function for polynomial multiplications Symbols x,y,z; diff --git a/sources/declare.h b/sources/declare.h index 951def5f..60d9dd80 100644 --- a/sources/declare.h +++ b/sources/declare.h @@ -1755,6 +1755,7 @@ void PathOutput(PHEAD WORD *,WORD,WORD *,WORD); #ifdef WITHFLOAT int DoStartFloat(UBYTE *); int DoEndFloat(UBYTE *); +int SetFloatPrecision(WORD); void SetupMZVTables(void); void SetupMPFTables(void); void ClearMZVTables(void); diff --git a/sources/float.c b/sources/float.c index 9f2ad055..f6d17a6b 100644 --- a/sources/float.c +++ b/sources/float.c @@ -54,7 +54,6 @@ void ZtoForm(UWORD *a,WORD *na,mpz_t z); long FloatToInteger(UWORD *out, mpf_t floatin, long *bitsused); void IntegerToFloat(mpf_t result, UWORD *formlong, int longsize); int FloatToRat(UWORD *ratout, WORD *nratout, mpf_t floatin); -int SetFloatPrecision(WORD prec); int AddFloats(PHEAD WORD *fun3, WORD *fun1, WORD *fun2); int MulFloats(PHEAD WORD *fun3, WORD *fun1, WORD *fun2); int DivFloats(PHEAD WORD *fun3, WORD *fun1, WORD *fun2); @@ -946,8 +945,12 @@ int PrintFloat(WORD *fun,int numdigits) if ( numdigits > prec || numdigits == 0 ) { numdigits = prec; } +/* + GMP's gmp_snprintf always prints a non-zero number before the decimal point, so we + ask for one digit less. +*/ if ( UnpackFloat(aux4,fun) == 0 ) - n = gmp_snprintf((char *)(AO.floatspace),AO.floatsize,"%.*Fe",numdigits,aux4); + n = gmp_snprintf((char *)(AO.floatspace),AO.floatsize,"%.*Fe",numdigits-1,aux4); if ( numdigits == prec ) { UBYTE *s1, *s2; int n1, n2 = n; @@ -1069,85 +1072,79 @@ void SetupMZVTables(void) and each deeper sum goes one higher, we make the tablesize a bit bigger. This may not be needed if we fiddle with the sum boundaries. */ - size_t nt = sizeof(mp_limb_t); - int prec; #ifdef WITHPTHREADS int i, Nw, id, totnum; - size_t fullsize, N, sumsize, j; - mp_limb_t *d; - Nw = AC.DefaultPrecision+AC.MaxWeight+1; - SetFloatPrecision(Nw); -/* prec = (AC.DefaultPrecision + 8*nt-1)/(8*nt); */ + size_t N; + mpf_t *a; + Nw = AC.DefaultPrecision; N = (size_t)Nw; - for ( i = 0, sumsize = 0; i <= Nw+1; i++ ) sumsize += (Nw-i+8*nt)/(8*nt)+1; - fullsize = (N+2)*sizeof(mpf_t)+sumsize*nt; totnum = AM.totalnumberofthreads; for ( id = 0; id < totnum; id++ ) { - if ( AB[id]->T.mpf_tab1 ) M_free(AB[id]->T.mpf_tab1,"mpftab1"); - AB[id]->T.mpf_tab1 = (void *)Malloc1(fullsize,"mpftab1"); - d = (mp_limb_t *)(((mpf_t *)(AB[id]->T.mpf_tab1))+N+2); - for ( j = 0; j < sumsize; j++ ) d[j] = 0; - for ( i = 0; i <= Nw; i++ ) { - ((mpf_t *)(AB[id]->T.mpf_tab1))[i]->_mp_prec = (Nw-i+8*nt)/(8*nt); - ((mpf_t *)(AB[id]->T.mpf_tab1))[i]->_mp_size = 0; - ((mpf_t *)(AB[id]->T.mpf_tab1))[i]->_mp_exp = 0; - ((mpf_t *)(AB[id]->T.mpf_tab1))[i]->_mp_d = d; - d += (Nw-i+8*nt)/(8*nt)+1; - } - if ( AB[id]->T.mpf_tab2 ) M_free(AB[id]->T.mpf_tab2,"mpftab2"); - AB[id]->T.mpf_tab2 = (void *)Malloc1(fullsize,"mpftab2"); - d = (mp_limb_t *)(((mpf_t *)(AB[id]->T.mpf_tab2))+N+1); - for ( j = 0; j < sumsize; j++ ) d[j] = 0; - for ( i = 0; i <= Nw; i++ ) { - ((mpf_t *)(AB[id]->T.mpf_tab2))[i]->_mp_prec = (Nw-i+8*nt)/(8*nt); - ((mpf_t *)(AB[id]->T.mpf_tab2))[i]->_mp_size = 0; - ((mpf_t *)(AB[id]->T.mpf_tab2))[i]->_mp_exp = 0; - ((mpf_t *)(AB[id]->T.mpf_tab2))[i]->_mp_d = d; - d += (Nw-i+8*nt)/(8*nt)+1; - } + if ( AB[id]->T.mpf_tab1 ) { + a = (mpf_t *)AB[id]->T.mpf_tab1; + for ( i = 0; i <=Nw; i++ ) { + mpf_clear(a[i]); + } + M_free(AB[id]->T.mpf_tab1,"mpftab1"); + } + AB[id]->T.mpf_tab1 = (void *)Malloc1((N+2)*sizeof(mpf_t),"mpftab1"); + a = (mpf_t *)AB[id]->T.mpf_tab1; + for ( i = 0; i <=Nw; i++ ) { +/* + As explained in the comment above, we could make this variable precision + using mpf_init2. +*/ + mpf_init(a[i]); + } + if ( AB[id]->T.mpf_tab2 ) { + a = (mpf_t *)AB[id]->T.mpf_tab2; + for ( i = 0; i <=Nw; i++ ) { + mpf_clear(a[i]); + } + M_free(AB[id]->T.mpf_tab2,"mpftab2"); + } + AB[id]->T.mpf_tab2 = (void *)Malloc1((N+2)*sizeof(mpf_t),"mpftab2"); + a = (mpf_t *)AB[id]->T.mpf_tab2; + for ( i = 0; i <=Nw; i++ ) { + mpf_init(a[i]); + } } #else int i, Nw; - size_t fullsize, N, sumsize, j; - mp_limb_t *d; -/* Nw = AC.DefaultPrecision+AC.MaxWeight+sizeof(mp_limb_t); */ - Nw = AC.DefaultPrecision+AC.MaxWeight+1; - SetFloatPrecision(Nw); -/* prec = (AC.DefaultPrecision + 8*nt-1)/(8*nt); */ + size_t N; + Nw = AC.DefaultPrecision; N = (size_t)Nw; - for ( i = 0, sumsize = 0; i <= Nw+1; i++ ) sumsize += (Nw-i+8*nt)/(8*nt)+1; - fullsize = (N+2)*sizeof(mpf_t)+sumsize*nt; - if ( mpftab1 ) M_free(mpftab1,"mpftab1"); - AT.mpf_tab1 = (void *)Malloc1(fullsize,"mpftab1"); - d = (mp_limb_t *)(((mpf_t *)(AT.mpf_tab1))+N+1); - for ( j = 0; j < sumsize; j++ ) d[j] = 0; + if ( AT.mpf_tab1 ) { + for ( i = 0; i <= Nw; i++ ) { + mpf_clear(mpftab1[i]); + } + M_free(AT.mpf_tab1,"mpftab1"); + } + AT.mpf_tab1 = (void *)Malloc1((N+2)*sizeof(mpf_t),"mpftab1"); for ( i = 0; i <= Nw; i++ ) { - mpftab1[i]->_mp_prec = (Nw-i+8*nt)/(8*nt); - mpftab1[i]->_mp_size = 0; - mpftab1[i]->_mp_exp = 0; - mpftab1[i]->_mp_d = d; - d += (Nw-i+8*nt)/(8*nt)+1; - } - if ( mpftab2 ) M_free(mpftab2,"mpftab2"); - AT.mpf_tab2 = (void *)Malloc1(fullsize,"mpftab2"); - d = (mp_limb_t *)(((mpf_t *)(AT.mpf_tab2))+N+1); - for ( j = 0; j < sumsize; j++ ) d[j] = 0; +/* + As explained in the comment above, we could make this variable precision + using mpf_init2. +*/ + mpf_init(mpftab1[i]); + } + if ( AT.mpf_tab2 ) { + for ( i = 0; i <= Nw; i++ ) { + mpf_clear(mpftab2[i]); + } + M_free(AT.mpf_tab2,"mpftab2"); + } + AT.mpf_tab2 = (void *)Malloc1((N+2)*sizeof(mpf_t),"mpftab2"); for ( i = 0; i <= Nw; i++ ) { - mpftab2[i]->_mp_prec = (Nw-i+8*nt)/(8*nt); - mpftab2[i]->_mp_size = 0; - mpftab2[i]->_mp_exp = 0; - mpftab2[i]->_mp_d = d; - d += (Nw-i+8*nt)/(8*nt)+1; + mpf_init(mpftab2[i]); } #endif - if ( AS.delta_1 ) M_free(AS.delta_1,"delta1"); - prec = (AC.DefaultPrecision + 8*nt-1)/(8*nt); - AS.delta_1 = (void *)Malloc1((prec+1)*sizeof(mp_limb_t)+sizeof(mpf_t),"delta1"); - d = (mp_limb_t *)(((mpf_t *)(AS.delta_1))+1); - mpfdelta1->_mp_prec = prec; - mpfdelta1->_mp_size = 0; - mpfdelta1->_mp_exp = 0; - mpfdelta1->_mp_d = d; + if ( AS.delta_1 ) { + mpf_clear(mpfdelta1); + M_free(AS.delta_1,"delta1"); + } + AS.delta_1 = (void *)Malloc1(sizeof(mpf_t),"delta1"); + mpf_init(mpfdelta1); SimpleDelta(mpfdelta1,1); /* this can speed up things. delta1 = ln(2) */ /* Finally the character buffer for printing @@ -1163,16 +1160,9 @@ void SetupMZVTables(void) void SetupMPFTables(void) { - size_t nt = sizeof(mp_limb_t); - int prec, prec1, i; #ifdef WITHPTHREADS - int Nw, id, totnum; - size_t j; - mp_limb_t *d; - Nw = AC.DefaultPrecision+AC.MaxWeight+1; - SetFloatPrecision(Nw); - prec = (AC.DefaultPrecision + 8*nt-1)/(8*nt); - prec1 = prec+1; + int id, totnum; + mpf_t *a; /* Now the aux variables */ @@ -1180,44 +1170,32 @@ void SetupMPFTables(void) totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads); #endif for ( id = 0; id < totnum; id++ ) { - if ( AB[id]->T.aux_ ) M_free(AB[id]->T.aux_,"aux_"); - AB[id]->T.aux_ = (void *)Malloc1(8*(prec1*sizeof(mp_limb_t)+sizeof(mpf_t)),"aux-mp"); - d = (mp_limb_t *)(((mpf_t *)(AB[id]->T.aux_))+8); - for ( j = 0; j < (size_t)(8*prec); j++ ) d[j] = 0; - for ( i = 0; i < 8; i++ ) { - ((mpf_t *)(AB[id]->T.aux_))[i]->_mp_prec = prec; - ((mpf_t *)(AB[id]->T.aux_))[i]->_mp_size = 0; - ((mpf_t *)(AB[id]->T.aux_))[i]->_mp_exp = 0; - ((mpf_t *)(AB[id]->T.aux_))[i]->_mp_d = d; - d += prec1; + if ( AB[id]->T.aux_ ) { +/* + We work here with a[0] etc because the aux1 etc contain B which + in the current routine would be AB[0] only +*/ + a = (mpf_t *)AB[id]->T.aux_; + mpf_clears(a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7],(mpf_ptr)0); + M_free(AB[id]->T.aux_,"AB[id]->T.aux_"); } + AB[id]->T.aux_ = (void *)Malloc1(sizeof(mpf_t)*8,"AB[id]->T.aux_"); + a = (mpf_t *)AB[id]->T.aux_; + mpf_inits(a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7],(mpf_ptr)0); if ( AB[id]->T.indi1 ) M_free(AB[id]->T.indi1,"indi1"); AB[id]->T.indi1 = (int *)Malloc1(sizeof(int)*AC.MaxWeight*2,"indi1"); AB[id]->T.indi2 = AB[id]->T.indi1 + AC.MaxWeight; } #else - int Nw; - size_t j; - mp_limb_t *d; -/* Nw = AC.DefaultPrecision+AC.MaxWeight+sizeof(mp_limb_t); */ - Nw = AC.DefaultPrecision+AC.MaxWeight+1; - SetFloatPrecision(Nw); - prec = (AC.DefaultPrecision + 8*nt-1)/(8*nt); - prec1 = prec+1; /* Now the aux variables */ - if ( mpfaux_ ) M_free(mpfaux_,"aux_"); - AT.aux_ = (void *)Malloc1(8*(prec1*sizeof(mp_limb_t)+sizeof(mpf_t)),"aux_"); - d = (mp_limb_t *)(((mpf_t *)(AT.aux_))+8); - for ( j = 0; j < (size_t)(8*prec); j++ ) d[j] = 0; - for ( i = 0; i < 8; i++ ) { - ((mpf_t *)(AT.aux_))[i]->_mp_prec = prec; - ((mpf_t *)(AT.aux_))[i]->_mp_size = 0; - ((mpf_t *)(AT.aux_))[i]->_mp_exp = 0; - ((mpf_t *)(AT.aux_))[i]->_mp_d = d; - d += prec1; + if ( AT.aux_ ) { + mpf_clears(aux1,aux2,aux3,aux4,aux5,auxjm,auxjjm,auxsum,(mpf_ptr)0); + M_free(AT.aux_,"AT.aux"); } + AT.aux_ = (void *)Malloc1(sizeof(mpf_t)*8,"AT.aux_"); + mpf_inits(aux1,aux2,aux3,aux4,aux5,auxjm,auxjjm,auxsum,(mpf_ptr)0); if ( AT.indi1 ) M_free(AT.indi1,"indi1"); AT.indi1 = (int *)Malloc1(sizeof(int)*AC.MaxWeight*2,"indi1"); AT.indi2 = AT.indi1 + AC.MaxWeight; @@ -1237,28 +1215,73 @@ void SetupMPFTables(void) void ClearMZVTables(void) { #ifdef WITHPTHREADS - int id, totnum; + int i, id, totnum; + mpf_t *a; totnum = AM.totalnumberofthreads; for ( id = 0; id < totnum; id++ ) { - if ( AB[id]->T.mpf_tab1 ) { M_free(AB[id]->T.mpf_tab1,"mpftab1"); AB[id]->T.mpf_tab1 = 0; } - if ( AB[id]->T.mpf_tab2 ) { M_free(AB[id]->T.mpf_tab2,"mpftab2"); AB[id]->T.mpf_tab2 = 0; } + if ( AB[id]->T.mpf_tab1 ) { +/* + We work here with a[0] etc because the aux1, mpftab1 etc contain B + which in the current routine would be AB[0] only +*/ + a = (mpf_t *)AB[id]->T.mpf_tab1; + for ( i = 0; i <=AC.DefaultPrecision; i++ ) { + mpf_clear(a[i]); + } + M_free(AB[id]->T.mpf_tab1,"mpftab1"); + AB[id]->T.mpf_tab1 = 0; + } + if ( AB[id]->T.mpf_tab2 ) { + a = (mpf_t *)AB[id]->T.mpf_tab2; + for ( i = 0; i <=AC.DefaultPrecision; i++ ) { + mpf_clear(a[i]); + } + M_free(AB[id]->T.mpf_tab2,"mpftab2"); + AB[id]->T.mpf_tab2 = 0; + } } #ifdef WITHSORTBOTS totnum = MaX(2*AM.totalnumberofthreads-3,AM.totalnumberofthreads); #endif for ( id = 0; id < totnum; id++ ) { - if ( AB[id]->T.aux_ ) { M_free(AB[id]->T.aux_,"aux-mp"); AB[id]->T.aux_ = 0; } + if ( AB[id]->T.aux_ ) { + a = (mpf_t *)AB[id]->T.aux_; + mpf_clears(a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7],(mpf_ptr)0); + M_free(AB[id]->T.aux_,"AB[id]->T.aux_"); + AB[id]->T.aux_ = 0; + } if ( AB[id]->T.indi1 ) { M_free(AB[id]->T.indi1,"indi1"); AB[id]->T.indi1 = 0; } } #else - if ( AT.mpf_tab1 ) { M_free(AT.mpf_tab1,"mpftab1"); AT.mpf_tab1 = 0; } - if ( AT.mpf_tab2 ) { M_free(AT.mpf_tab2,"mpftab2"); AT.mpf_tab2 = 0; } - if ( AT.aux_ ) { M_free(AT.aux_,"aux-mp"); AT.aux_ = 0; } + int i; + if ( AT.mpf_tab1 ) { + for ( i = 0; i <= AC.DefaultPrecision; i++ ) { + mpf_clear(mpftab1[i]); + } + M_free(AT.mpf_tab1,"mpftab1"); + AT.mpf_tab1 = 0; + } + if ( AT.mpf_tab2 ) { + for ( i = 0; i <= AC.DefaultPrecision; i++ ) { + mpf_clear(mpftab2[i]); + } + M_free(AT.mpf_tab2,"mpftab2"); + AT.mpf_tab2 = 0; + } + if ( AT.aux_ ) { + mpf_clears(aux1,aux2,aux3,aux4,aux5,auxjm,auxjjm,auxsum,(mpf_ptr)0); + M_free(AT.aux_,"AT.aux_"); + AT.aux_ = 0; + } if ( AT.indi1 ) { M_free(AT.indi1,"indi1"); AT.indi1 = 0; } #endif if ( AO.floatspace ) { M_free(AO.floatspace,"floatspace"); AO.floatspace = 0; AO.floatsize = 0; } - if ( AS.delta_1 ) { mpfdelta1->_mp_d = 0; M_free(AS.delta_1,"delta1"); } + if ( AS.delta_1 ) { + mpf_clear(mpfdelta1); + M_free(AS.delta_1,"delta1"); + AS.delta_1 = 0; + } } /* @@ -1668,15 +1691,19 @@ void SimpleDelta(mpf_t sum, int m) We will loop until 1/2^j/j^m is smaller than the default precision. Just running to prec is however overkill, specially when m is large. We try to estimate a better value. - jmax = prec - ((2log(prec)-1)*m). + jmax = prec - ((log2(prec)-1)*m). Hence we need the leading bit of prec. We are still overshooting a bit. */ n = 0; x = xprec; while ( x ) { x >>= 1; n++; } - jmax = (int)((int)xprec - n*m); +/* + We have now n = floor(log2(x))+1. +*/ + n--; + jmax = (int)((int)xprec - (n-1)*m); mpf_set_ui(sum,0); - for ( j = 1; j < jmax; j++ ) { + for ( j = 1; j <= jmax; j++ ) { #ifdef WITHCUTOFF xprec--; mpf_set_prec_raw(jm,xprec); @@ -1710,16 +1737,20 @@ void SimpleDeltaC(mpf_t sum, int m) We will loop until 1/2^j/j^m is smaller than the default precision. Just running to prec is however overkill, specially when m is large. We try to estimate a better value. - jmax = prec - ((2log(prec)-1)*m). + jmax = prec - ((log2(prec)-1)*m). Hence we need the leading bit of prec. We are still overshooting a bit. */ n = 0; x = xprec; while ( x ) { x >>= 1; n++; } - jmax = (int)((int)xprec - n*m); +/* + We have now n = floor(log2(x))+1. +*/ + n--; + jmax = (int)((int)xprec - (n-1)*m); if ( s < 0 ) jmax /= 2; mpf_set_si(sum,0L); - for ( j = 1; j < jmax; j++ ) { + for ( j = 1; j <= jmax; j++ ) { #ifdef WITHCUTOFF xprec--; mpf_set_prec_raw(jm,xprec); diff --git a/sources/fsizes.h b/sources/fsizes.h index 6fa740ac..8e8f4af5 100644 --- a/sources/fsizes.h +++ b/sources/fsizes.h @@ -194,6 +194,6 @@ #define MAXLEGS 20 #ifdef WITHFLOAT -#define MAXWEIGHT 40 +#define MAXWEIGHT 0 #define DEFAULTPRECISION 1000 #endif diff --git a/sources/pre.c b/sources/pre.c index 9c6a15af..5c7c0448 100644 --- a/sources/pre.c +++ b/sources/pre.c @@ -7727,6 +7727,7 @@ int DoStartFloat(UBYTE *s) AC.MaxWeight = AC.tMaxWeight; AC.tMaxWeight = 0; } + SetFloatPrecision(AC.DefaultPrecision+AC.MaxWeight+1); SetupMPFTables(); if ( AC.MaxWeight > 0 ) SetupMZVTables(); SetfFloatPrecision(AC.DefaultPrecision);