fix: the estimate for jmax in SimpleDelta and SimpleDeltaC could be wrong and even be negative.

- It is now increased iteratively until it meets the required precision.
- This also fixes the remaining rounding issues in the euler_ and mzvhalf_ tests (see PR699).
- furthermore, the raw form in the Float_1 test has been slightly adjusted.

Author: cbmarini

check/features.frm

@@ -662,13 +662,13 @@ assert stdout =~ exact_pattern(<<'EOF')
    G1 =
       float_(10,10,1,225649930152063087544280124519603661924016376904153173222\
       961313295752588365049496190204609112667217257624737508104376810254484309\
-      3830597279756558899541412580219255737022507716035);
+      3830597279756558899541412586941038636793415360118);
 EOF
 assert stdout =~ exact_pattern(<<'EOF')
    G2 =
       float_(10,10,1,247337966873870703631573653423368526098272821023387457752\
       643030990710929556575503489173572446024607642191456056189912528713904095\
-      704444384800390830056125051605311704862654126);
+      704444384800390830056125051897088334731440434);
 EOF
 *--#] Float_1 :
 *--#[ evaluate_symbol :
@@ -919,7 +919,7 @@ assert result("EULER3") =~ expr("
        - 3.888958461681063290997e-01*euler(-1,-2)
        + 2.140723708667062274342e-01*euler(-1,-1,-1)
        - 5.372131936080402009406e-01*euler(-1,-1,1)
-       + 9.475300423012770572182e-02*euler(-1,1,-1)
+       + 9.475300423012770572183e-02*euler(-1,1,-1)
        - 5.550410866482157995314e-02*euler(-1,1,1)
        + 2.695764795315278073874e-01*euler(-1,2)
        - 5.082152128046848508121e-01*euler(2,-1)
@@ -993,7 +993,7 @@ assert result("MZVHALF1") =~ expr("
         +6.931471805599453094172e-01*mzvhalf(1)
 ")
 assert result("MZVHALF2") =~ expr("
-        +2.402265069591007123335e-01*mzvhalf(1,1)
+        +2.402265069591007123336e-01*mzvhalf(1,1)
        + 5.822405264650125059027e-01*mzvhalf(2)
 ")
 assert result("MZVHALF3") =~ expr("

sources/float.c

@@ -1920,6 +1920,18 @@ void SimpleDelta(mpf_t sum, int m)
 */
 	n--;
 	jmax = (int)((int)xprec - (n-1)*m);
+/*
+	For small prec and large m, the estimate can be wrong and even be negative, 
+	so we increase jmax until jmax + m*log2(jmax) > prec
+*/
+	if ( jmax < 0 ) jmax = 1;
+	do {
+		n = 0;
+		x = (unsigned long)jmax;
+		while (x) { x >>= 1; n++; }
+		n--; // floor(log2(jmax))
+		jmax++; 
+	} while ( jmax + m * n <= prec );
 	mpf_set_ui(sum,0);
 	for ( j = 1; j <= jmax; j++ ) {
 #ifdef WITHCUTOFF
@@ -1966,6 +1978,18 @@ void SimpleDeltaC(mpf_t sum, int m)
 */
 	n--;
 	jmax = (int)((int)xprec - (n-1)*m);
+/*
+	For small prec and large m, the estimate can be wrong and even be negative, 
+	so we increase jmax until jmax + m*log2(jmax) > prec
+*/
+	if ( jmax < 0 ) jmax = 1;
+	do {
+		n = 0;
+		x = (unsigned long)jmax;
+		while (x) { x >>= 1; n++; }
+		n--; // floor(log2(jmax))
+		jmax++; 
+	} while ( jmax + m * n <= prec );
 	if ( s < 0 ) jmax /= 2;
 	mpf_set_si(sum,0L);
 	for ( j = 1; j <= jmax; j++ ) {