Index: sageSLZ.sage
===================================================================
--- sageSLZ.sage	(revision 60)
+++ sageSLZ.sage	(revision 61)
@@ -1,3 +1,61 @@
+
+def slz_compute_polynomial_and_interval(functionSo, degreeSo, lowerBoundSa, 
+                                        upperBoundSa, approxPrecSa, 
+                                        sollyaPrecSa=None):
+    """
+    Under the assumptions listed for slz_get_intervals_and_polynomials, compute
+    a polynomial that approximates the function on a an interval starting
+    at lowerBoundSa and finishing at a value that guarantees that the polynomial
+    approximates with the expected precision.
+    The interval upper bound is lowered until the expected approximation 
+    precision is reached.
+    The polynomial, the bounds, the center of the interval and the error
+    are returned.
+    """
+    RRR = lowerBoundSa.parent()
+    #goldenRatioSa = RRR(5.sqrt() / 2 - 1/2)
+    #intervalShrinkConstFactorSa = goldenRatioSa
+    intervalShrinkConstFactorSa = RRR('0.5')
+    absoluteErrorTypeSo = pobyso_absolute_so_so()
+    currentRangeSo = pobyso_bounds_to_range_sa_so(lowerBoundSa, upperBoundSa)
+    currentUpperBoundSa = upperBoundSa
+    currentLowerBoundSa = lowerBoundSa
+    # What we want here is the polynomial without the variable change, 
+    # since our actual variable will be x-intervalCenter defined over the 
+    # domain [lowerBound-intervalCenter , upperBound-intervalCenter]. 
+    (polySo, intervalCenterSo, maxErrorSo) = \
+        pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
+                                                    currentRangeSo, 
+                                                    absoluteErrorTypeSo)
+    maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
+    while maxErrorSa > approxPrecSa:
+        sollya_lib_clear_obj(maxErrorSo)
+        errorRatioSa = 1/(maxErrorSa/approxPrecSa).log2()
+        #print "Error ratio: ", errorRatioSa
+        if errorRatioSa < intervalShrinkConstFactorSa:
+            #currentUpperBoundSa = currentLowerBoundSa + (currentUpperBoundSa - currentLowerBoundSa) * errorRatioSa
+            currentUpperBoundSa = currentLowerBoundSa + \
+                                  (currentUpperBoundSa - currentLowerBoundSa) * \
+                                   intervalShrinkConstFactorSa
+        else:
+            currentUpperBoundSa = currentLowerBoundSa + \
+                                  (currentUpperBoundSa - currentLowerBoundSa) * \
+                                  intervalShrinkConstFactorSa
+        #print lowerBoundSa, currentUpperBoundSa
+        sollya_lib_clear_obj(currentRangeSo)
+        sollya_lib_clear_obj(polySo)
+        currentRangeSo = pobyso_bounds_to_range_sa_so(currentLowerBoundSa, 
+                                                      currentUpperBoundSa)
+        (polySo, intervalCenterSo, maxErrorSo) = \
+            pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
+                                                        currentRangeSo, 
+                                                        absoluteErrorTypeSo)
+        #maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo, RRR)
+        maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
+    sollya_lib_clear_obj(absoluteErrorTypeSo)
+    return((polySo, currentRangeSo, intervalCenterSo, maxErrorSo))
+    # End slz_compute_polynomial_and_interval
+
 def slz_get_intervals_and_polynomials(functionSa, degreeSa, lowerBoundSa, 
                                       upperBoundSa, floatingPointPrecSa, 
                                       internalSollyaPrecSa):
@@ -13,14 +71,17 @@
     - the internal Sollya precision;
     compute a list of tuples made of:
     - a polynomial approximating the function (a Sollya object);
-    - the bounds for which the polynomial approximates the function 
-      (a Sage object);
-    - the center of the interval;
+    - the range for which the polynomial approximates the function 
+      (a Sollya object);
+    - the center of the interval (a Sollya object);
     - the approximation error (a Sage object). 
       with the error given as the last element (a Sage object);
-    The initial interval is, possibly, splitted into a list of smaller interval
-    each associated to an approximation polynomial and a the corresponding 
-    approximation error.
+    The initial interval is, possibly, splitted into smaller intervals.
+    It return a list of tuples, each made of:
+    - a polynomial;
+    - the approximation interval;
+    - the center, x0,  of the interval (the polynomial is defined as p(x-x0));
+    - the corresponding approximation error.
     """
     # Scalling the domain -> [1,2[.
     # Notice the clumsy notation for log2.
@@ -44,12 +105,13 @@
     #
     approxPrecSa = 1/(2^(floatingPointPrecSa + 1))
     print "Approximation precision: ", RR(approxPrecSa)
-    # Prepare the arguments for the Taylor expansion computation.
+    # Prepare the arguments for the Taylor expansion computation with Sollya.
     functionSo = pobyso_parse_string_sa_so(functionSa._assume_str())
     degreeSo = pobyso_constant_from_int_sa_so(degreeSa)
-    scaledBoundsSo = pobyso_range_sa_so(scaledLowerBoundSa, scaledUpperBoundSa)
+    scaledBoundsSo = pobyso_bounds_to_range_sa_so(scaledLowerBoundSa, 
+                                                  scaledUpperBoundSa)
     absoluteErrorTypeSo = pobyso_absolute_so_so()
-    #Compute the first Taylor expansion.
+    # Compute the first Taylor expansion.
     (polySo, boundsSo, intervalCenterSo, maxErrorSo) = \
         slz_compute_polynomial_and_interval(functionSo, degreeSo,
                                         scaledLowerBoundSa, scaledUpperBoundSa,
@@ -57,7 +119,8 @@
     resultArray.append((polySo, boundsSo, intervalCenterSo, maxErrorSo))
     realIntervalField = RealIntervalField(max(lowerBoundSa.parent().precision(),
                                               upperBoundSa.parent().precision()))
-    boundsSa = pobyso_range_so_sa(boundsSo, realIntervalField)
+    boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)
+    # Compute the other expansions.
     while boundsSa.endpoints()[1] < scaledUpperBoundSa:
         currentScaledLowerBoundSa = boundsSa.endpoints()[1]
         (polySo, boundsSo, intervalCenterSo, maxErrorSo) = \
@@ -66,7 +129,7 @@
                                             scaledUpperBoundSa, approxPrecSa, 
                                             internalSollyaPrecSa)
         resultArray.append((polySo, boundsSo, intervalCenterSo, maxErrorSo))
-        boundsSa = pobyso_range_so_sa(boundsSo, realIntervalField)
+        boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)
     sollya_lib_clear_obj(functionSo)
     sollya_lib_clear_obj(degreeSo)
     sollya_lib_clear_obj(scaledBoundsSo)
@@ -74,53 +137,30 @@
     return(resultArray)
     # End slz_get_intervals_and_polynomials
 
-def slz_compute_polynomial_and_interval(functionSo, degreeSo, lowerBoundSa, 
-                                        upperBoundSa, approxPrecSa, 
-                                        sollyaPrecSa=None):
-    RRR = lowerBoundSa.parent()
-    goldenRatioSa = RRR(5.sqrt() / 2 - 1/2)
-    #intervalShrinkConstFactorSa = goldenRatioSa
-    intervalShrinkConstFactorSa = RRR('0.5')
-    absoluteErrorTypeSo = pobyso_absolute_so_so()
-    currentRangeSo = pobyso_range_sa_so(lowerBoundSa, upperBoundSa)
-    currentUpperBoundSa = upperBoundSa
-    currentLowerBoundSa = lowerBoundSa
-    # What we want here is the polynomial without the variable change, 
-    # since our actual variable will be x-intervalCenter defined over the 
-    # domain [lowerBound-intervalCenter , upperBound-intervalCenter]. 
-    (polySo, intervalCenterSo, maxErrorSo) = \
-        pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
-                                                    currentRangeSo, 
-                                                    absoluteErrorTypeSo)
-    maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
-    while maxErrorSa > approxPrecSa:
-        sollya_lib_clear_obj(maxErrorSo)
-        errorRatioSa = 1/(maxErrorSa/approxPrecSa).log2()
-        #print "Error ratio: ", errorRatioSa
-        if errorRatioSa < intervalShrinkConstFactorSa:
-            #currentUpperBoundSa = currentLowerBoundSa + (currentUpperBoundSa - currentLowerBoundSa) * errorRatioSa
-            currentUpperBoundSa = currentLowerBoundSa + \
-                                  (currentUpperBoundSa - currentLowerBoundSa) * \
-                                   intervalShrinkConstFactorSa
+def slz_interval_scaling_expression(boundsInterval, varName):
+    """
+    Compute the scaling expression to map an interval that span only
+    a binade to [1, 2)
+    """
+    if abs(boundsInterval.endpoints()[0]) < 1:
+        if boundsInterval.endpoints()[0] >= 0:
+            scalingCoeff = 2^(-floor(boundsInterval.endpoints()[0].log2()))
+            return(scalingCoeff * eval(varName))
         else:
-            currentUpperBoundSa = currentLowerBoundSa + \
-                                  (currentUpperBoundSa - currentLowerBoundSa) * \
-                                  intervalShrinkConstFactorSa
-        #print lowerBoundSa, currentUpperBoundSa
-        sollya_lib_clear_obj(currentRangeSo)
-        sollya_lib_clear_obj(polySo)
-        currentRangeSo = pobyso_range_sa_so(currentLowerBoundSa, 
-                                            currentUpperBoundSa)
-        (polySo, intervalCenterSo, maxErrorSo) = \
-            pobyso_taylor_expansion_no_change_var_so_so(functionSo, degreeSo, 
-                                                        currentRangeSo, 
-                                                        absoluteErrorTypeSo)
-        #maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo, RRR)
-        maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
-    sollya_lib_clear_obj(absoluteErrorTypeSo)
-    return((polySo, currentRangeSo, intervalCenterSo, maxErrorSo))
-    # End slz_compute_polynomial_and_interval
-    
+            scalingCoeff = 2^(-floor((-boundsInterval.endpoints()[1]).log2()))
+            scalingOffset = -ceil(scalingCoeff * boundsInterval.endpoints()[0])
+            return(scalingCoeff * eval(varName) + scalingOffset)
+    else: 
+        if boundsInterval.endpoints()[0] >= 0:
+            scalingCoeff = 2^(-floor(boundsInterval.endpoints()[0].log2()))
+            scalingOffset = 0
+            return(scalingCoeff * eval(varName))
+        else:
+            scalingCoeff = 2^(-floor((-boundsInterval.endpoints()[1]).log2()))
+            scalingOffset = floor(-(scalingCoeff * boundsInterval.endpoints()[1]) + 2)
+            return(scalingCoeff * eval(varName) + scalingOffset)
+
+   
 def slz_polynomial_and_interval_to_sage(polyRangeCenterErrorSo):
     polynomialSa = pobyso_get_poly_so_sa(polyRangeCenterErrorSo[0])
     intervalSa = \
