Index: pobysoPythonSage/src/sageSLZ/sageSLZ.sage
===================================================================
--- pobysoPythonSage/src/sageSLZ/sageSLZ.sage	(revision 100)
+++ pobysoPythonSage/src/sageSLZ/sageSLZ.sage	(revision 101)
@@ -36,14 +36,14 @@
                                                     absoluteErrorTypeSo)
     maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
     while maxErrorSa > approxPrecSa:
+        #print "++Approximation error:", maxErrorSa
         sollya_lib_clear_obj(maxErrorSo)
         sollya_lib_clear_obj(polySo)
         sollya_lib_clear_obj(intervalCenterSo)
-        shrinkFactorSa = RRR('5.0')/(maxErrorSa/approxPrecSa).log2().abs()
+        shrinkFactorSa = RRR('5')/(maxErrorSa/approxPrecSa).log2().abs()
         #shrinkFactorSa = 1.5/(maxErrorSa/approxPrecSa)
         #errorRatioSa = approxPrecSa/maxErrorSa
         #print "Error ratio: ", errorRatioSa
-        
         if shrinkFactorSa > intervalShrinkConstFactorSa:
             actualShrinkFactorSa = intervalShrinkConstFactorSa
             #print "Fixed"
@@ -51,13 +51,15 @@
             actualShrinkFactorSa = shrinkFactorSa
             #print "Computed",shrinkFactorSa,maxErrorSa
             #print shrinkFactorSa, maxErrorSa
+        #print "Shrink factor", actualShrinkFactorSa
         currentUpperBoundSa = currentLowerBoundSa + \
                                   (currentUpperBoundSa - currentLowerBoundSa) * \
                                    actualShrinkFactorSa
         #print "Current upper bound:", currentUpperBoundSa
         sollya_lib_clear_obj(currentRangeSo)
         sollya_lib_clear_obj(polySo)
-        if currentUpperBoundSa <= currentLowerBoundSa:
+        if currentUpperBoundSa <= currentLowerBoundSa or \
+              currentUpperBoundSa == currentLowerBoundSa.nextabove():
             sollya_lib_clear_obj(absoluteErrorTypeSo)
             print "Can't find an interval."
             print "Use either or both a higher polynomial degree or a higher",
@@ -210,11 +212,14 @@
     - the corresponding approximation error.
     TODO: fix endless looping for some parameters sets.
     """
+    # Set Sollya to the necessary internal precision.
     currentSollyaPrecSo = pobyso_get_prec_so()
     currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)
     if internalSollyaPrecSa > currentSollyaPrecSa:
         pobyso_set_prec_sa_so(internalSollyaPrecSa)
+    #
     x = functionSa.variables()[0] # Actual variable name can be anything.
+    # Scaled function: [1=,2] -> [1,2].
     (fff, scaledLowerBoundSa, scaledUpperBoundSa, \
             scaledLowerBoundImageSa, scaledUpperBoundImageSa) = \
                 slz_compute_scaled_function(functionSa, \
@@ -236,32 +241,45 @@
                                         scaledLowerBoundSa, scaledUpperBoundSa,
                                         approxPrecSa, internalSollyaPrecSa)
     if polySo is None:
-        print "Aborting"
+        print "slz_get_intervals_and_polynomials: Aborting and returning None!"
+        if internalSollyaPrecSa != currentSollyaPrecSa:
+            pobyso_set_prec_sa_so(currentSollyaPrecSa)
         return None
-    # Change variable stuff
-    changeVarExpressionSo = sollya_lib_build_function_sub(
-                              sollya_lib_build_function_free_variable(), 
-                              sollya_lib_copy_obj(intervalCenterSo))
-    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo) 
-    resultArray.append((polySo, polyVarChangedSo, boundsSo, intervalCenterSo,\
-                         maxErrorSo))
     realIntervalField = RealIntervalField(max(lowerBoundSa.parent().precision(),
                                               upperBoundSa.parent().precision()))
     boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)
+    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
+    #print "First approximation error:", errorSa
+    # If the error and interval are OK a the first try, just return.
+    if boundsSa.endpoints()[1] >= scaledUpperBoundSa:
+        # Change variable stuff in Sollya x -> x0-x.
+        changeVarExpressionSo = sollya_lib_build_function_sub( \
+                              sollya_lib_build_function_free_variable(), \
+                              sollya_lib_copy_obj(intervalCenterSo))
+        polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo) 
+        resultArray.append((polySo, polyVarChangedSo, boundsSo, \
+                            intervalCenterSo, maxErrorSo))
+        if internalSollyaPrecSa != currentSollyaPrecSa:
+            pobyso_set_prec_sa_so(currentSollyaPrecSa)
+        sollya_lib_clear_obj(functionSo)
+        sollya_lib_clear_obj(degreeSo)
+        sollya_lib_clear_obj(scaledBoundsSo)
+        #print "Approximation error:", errorSa
+        return resultArray
     # Compute the next upper bound.
-    # If the error of approximation is more than half of the target,
-    # use the same interval.
-    # If it is less, increase it a bit.
-    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
+    # The following ratio is always >= 1
     currentErrorRatio = approxPrecSa / errorSa
+    # Starting point for the next upper bound.
     currentScaledUpperBoundSa = boundsSa.endpoints()[1]
-    if currentErrorRatio  < 2 :
+    boundsWidthSa = boundsSa.endpoints()[1] - boundsSa.endpoints()[0]
+    # Compute the increment. 
+    if currentErrorRatio > RR('1000'): # ]1.5, infinity[
         currentScaledUpperBoundSa += \
-            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0])
-    else:
+                        currentErrorRatio * boundsWidthSa * 2
+    else:  # [1, 1.5]
         currentScaledUpperBoundSa += \
-            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0]) \
-                         * currentErrorRatio.log2() * 2
+                        (RR('1.0') + currentErrorRatio.log() / 500) * boundsWidthSa
+    # Take into account the original interval upper bound.                     
     if currentScaledUpperBoundSa > scaledUpperBoundSa:
         currentScaledUpperBoundSa = scaledUpperBoundSa
     # Compute the other expansions.
@@ -273,40 +291,40 @@
                                             currentScaledUpperBoundSa, 
                                             approxPrecSa, 
                                             internalSollyaPrecSa)
-        # Change variable stuff
-        changeVarExpressionSo = sollya_lib_build_function_sub(
-                                  sollya_lib_build_function_free_variable(), 
-                                  sollya_lib_copy_obj(intervalCenterSo))
-        polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo) 
-        resultArray.append((polySo, polyVarChangedSo, boundsSo, \
-                            intervalCenterSo, maxErrorSo))
+        errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
+        if  errorSa < approxPrecSa:
+            # Change variable stuff
+            #print "Approximation error:", errorSa
+            changeVarExpressionSo = sollya_lib_build_function_sub(
+                                    sollya_lib_build_function_free_variable(), 
+                                    sollya_lib_copy_obj(intervalCenterSo))
+            polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo) 
+            resultArray.append((polySo, polyVarChangedSo, boundsSo, \
+                                intervalCenterSo, maxErrorSo))
         boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)
         # Compute the next upper bound.
-        # If the error of approximation is more than half of the target,
-        # use the same interval.
-        # If it is less, increase it a bit.
-        errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
+        # The following ratio is always >= 1
         currentErrorRatio = approxPrecSa / errorSa
-        if currentErrorRatio  < RR('1.5') :
-            currentScaledUpperBoundSa = \
-                            boundsSa.endpoints()[1] + \
-                            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0])
-        elif currentErrorRatio < 2:
-            currentScaledUpperBoundSa = \
-                            boundsSa.endpoints()[1] + \
-                            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0]) \
-                             * currentErrorRatio.log2()
-        else:
-            currentScaledUpperBoundSa = \
-                            boundsSa.endpoints()[1] + \
-                            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0]) \
-                             * currentErrorRatio.log2() * 2
+        # Starting point for the next upper bound.
+        currentScaledUpperBoundSa = boundsSa.endpoints()[1]
+        boundsWidthSa = boundsSa.endpoints()[1] - boundsSa.endpoints()[0]
+        # Compute the increment. 
+        if currentErrorRatio > RR('1000'): # ]1.5, infinity[
+            currentScaledUpperBoundSa += \
+                            currentErrorRatio * boundsWidthSa * 2
+        else:  # [1, 1.5]
+            currentScaledUpperBoundSa += \
+                            (RR('1.0') + currentErrorRatio.log()/500) * boundsWidthSa
+        #print "currentErrorRatio:", currentErrorRatio
+        #print "currentScaledUpperBoundSa", currentScaledUpperBoundSa
         # Test for insufficient precision.
         if currentScaledUpperBoundSa == scaledLowerBoundSa:
             print "Can't shrink the interval anymore!"
             print "You should consider increasing the Sollya internal precision"
             print "or the polynomial degree."
             print "Giving up!"
+            if internalSollyaPrecSa != currentSollyaPrecSa:
+                pobyso_set_prec_sa_so(currentSollyaPrecSa)
             sollya_lib_clear_obj(functionSo)
             sollya_lib_clear_obj(degreeSo)
             sollya_lib_clear_obj(scaledBoundsSo)
