Index: pobysoPythonSage/src/sageSLZ/sageSLZ.sage
===================================================================
--- pobysoPythonSage/src/sageSLZ/sageSLZ.sage	(revision 0)
+++ pobysoPythonSage/src/sageSLZ/sageSLZ.sage	(revision 66)
@@ -0,0 +1,208 @@
+
+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 = (approxPrecSa/maxErrorSa).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, variableNameSa, degreeSa, 
+                                      lowerBoundSa, 
+                                      upperBoundSa, floatingPointPrecSa, 
+                                      internalSollyaPrecSa, approxPrecSa):
+    """
+    Under the assumption that:
+    - functionSa is monotonic on the [lowerBoundSa, upperBoundSa] interval;
+    - lowerBound and upperBound belong to the same binade.
+    from a:
+    - function;
+    - a degree
+    - a pair of bounds;
+    - the floating-point precision we work on;
+    - the internal Sollya precision;
+    - the requested approximation error
+    compute a list of tuples made of:
+    - a polynomial approximating the function (a Sollya object);
+    - the range for which the polynomial approximates the function 
+      (a Sollya object);
+    - the center of the interval (a Sollya object);
+    - the actual approximation error (a Sage object). 
+    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.
+    """
+    x = var(variableNameSa)
+    # Scalling the domain -> [1,2[.
+    boundsIntervalRifSa = RealIntervalField(floatingPointPrecSa)
+    domainBoundsIntervalSa = boundsIntervalRifSa(lowerBoundSa, upperBoundSa)
+    (domainScalingExpressionSa, invDomainScalingExpressionSa) = \
+        slz_interval_scaling_expression(domainBoundsIntervalSa, variableNameSa)
+    print "domainScalingExpression for argument :", domainScalingExpressionSa
+    print "f: ", f
+    ff = f.subs({x : domainScalingExpressionSa})
+    #ff = f.subs_expr(x==domainScalingExpressionSa)
+    scaledLowerBoundSa = invDomainScalingExpressionSa(lowerBoundSa).n()
+    scaledUpperBoundSa = invDomainScalingExpressionSa(upperBoundSa).n()
+    print 'ff:', ff, "- Domain:", scaledLowerBoundSa, scaledUpperBoundSa
+    #
+    # Scalling the image -> [1,2[.
+    flbSa = f(lowerBoundSa).n()
+    fubSa = f(upperBoundSa).n()
+    if flbSa <= fubSa: # Increasing
+        imageBinadeBottomSa = floor(flbSa.log2())
+    else: # Decreasing
+        imageBinadeBottomSa = floor(fubSa.log2())
+    print 'ff:', ff, '- Image:', flbSa, fubSa, imageBinadeBottomSa
+    imageBoundsIntervalSa = boundsIntervalRifSa(flbSa, fubSa)
+    (imageScalingExpressionSa, invImageScalingExpressionSa) = \
+        slz_interval_scaling_expression(imageBoundsIntervalSa, variableNameSa)
+    iis = invImageScalingExpressionSa.function(x)
+    fff = iis.subs({x:ff})
+    print "fff:", fff,
+    print " - Image:", fff(scaledLowerBoundSa), fff(scaledUpperBoundSa)
+    #
+    resultArray = []
+    #
+    print "Approximation precision: ", RR(approxPrecSa)
+    # Prepare the arguments for the Taylor expansion computation with Sollya.
+    functionSo = pobyso_parse_string_sa_so(fff._assume_str())
+    degreeSo = pobyso_constant_from_int_sa_so(degreeSa)
+    scaledBoundsSo = pobyso_bounds_to_range_sa_so(scaledLowerBoundSa, 
+                                                  scaledUpperBoundSa)
+    # Compute the first Taylor expansion.
+    (polySo, boundsSo, intervalCenterSo, maxErrorSo) = \
+        slz_compute_polynomial_and_interval(functionSo, degreeSo,
+                                        scaledLowerBoundSa, scaledUpperBoundSa,
+                                        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))
+    realIntervalField = RealIntervalField(max(lowerBoundSa.parent().precision(),
+                                              upperBoundSa.parent().precision()))
+    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) = \
+            slz_compute_polynomial_and_interval(functionSo, degreeSo,
+                                            currentScaledLowerBoundSa, 
+                                            scaledUpperBoundSa, 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))
+        boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)
+    sollya_lib_clear_obj(functionSo)
+    sollya_lib_clear_obj(degreeSo)
+    sollya_lib_clear_obj(scaledBoundsSo)
+    return(resultArray)
+    # End slz_get_intervals_and_polynomials
+
+def slz_interval_scaling_expression(boundsInterval, varName):
+    """
+    Compute the scaling expression to map an interval that span only
+    a binade to [1, 2) and the inverse expression as well.
+    Not very sure that the transformation makes sense for negative numbers.
+    """
+    # The scaling offset is only used for negative numbers.
+    if abs(boundsInterval.endpoints()[0]) < 1:
+        if boundsInterval.endpoints()[0] >= 0:
+            scalingCoeff = 2^floor(boundsInterval.endpoints()[0].log2())
+            invScalingCoeff = 1/scalingCoeff
+            return((scalingCoeff * eval(varName), 
+                    invScalingCoeff * eval(varName)))
+        else:
+            scalingCoeff = \
+                2^(floor((-boundsInterval.endpoints()[0]).log2()) - 1)
+            scalingOffset = -3 * scalingCoeff
+            return((scalingCoeff * eval(varName) + scalingOffset,
+                    1/scalingCoeff * eval(varName) + 3))
+    else: 
+        if boundsInterval.endpoints()[0] >= 0:
+            scalingCoeff = 2^floor(boundsInterval.endpoints()[0].log2())
+            scalingOffset = 0
+            return((scalingCoeff * eval(varName), 
+                    1/scalingCoeff * eval(varName)))
+        else:
+            scalingCoeff = \
+                2^(floor((-boundsInterval.endpoints()[1]).log2()))
+            scalingOffset =  -3 * scalingCoeff
+            #scalingOffset = 0
+            return((scalingCoeff * eval(varName) + scalingOffset,
+                    1/scalingCoeff * eval(varName) + 3))
+
+   
+def slz_polynomial_and_interval_to_sage(polyRangeCenterErrorSo):
+    polynomialSa = pobyso_get_poly_so_sa(polyRangeCenterErrorSo[0])
+    polynomialChangedVarSa = pobyso_get_poly_so_sa(polyRangeCenterErrorSo[1])
+    intervalSa = \
+            pobyso_get_interval_from_range_so_sa(polyRangeCenterErrorSo[2])
+    centerSa = \
+            pobyso_get_constant_as_rn_with_rf_so_sa(polyRangeCenterErrorSo[3])
+    errorSa = \
+            pobyso_get_constant_as_rn_with_rf_so_sa(polyRangeCenterErrorSo[4])
+    return((polynomialSa, polynomialChangedVarSa, intervalSa, centerSa, errorSa))
+    # End slz_polynomial_and_interval_to_sage
+
+print "sageSLZ loaded..."
