commit 04d64614b764878649198d0b2cad50bb1bc8a9da
Author: Serge Torres <Serge.Torres@ens-lyon.fr>
Date:   Thu Dec 31 16:52:35 2015 +0100

    Added: more functions for points and coefficients bounds exponents.

diff --git a/src/functions_for_cvp.sage b/src/functions_for_cvp.sage
index ad7cfff..1db1106 100644
--- a/src/functions_for_cvp.sage
+++ b/src/functions_for_cvp.sage
@@ -3,6 +3,43 @@ Auxiliary functions used for CVP.
 """
 print "Functions for CVP loading..."
 #
+def cvp_affine_from_chebyshev(lowerBound, upperBound):
+    """
+    Compute the affine parameters to map [-1, 1] to the interval given
+    as argument.
+    @param lowerBound (of the target interval)
+    @param upperBound (of the target interval)
+    @return the (scalingCoefficient, offset) tuple.
+    """
+    ## Check bounds consistency. Bounds must differ.
+    if lowerBound >= upperBound:
+        return None
+    scalingCoefficient = (upperBound - lowerBound) / 2
+    offset             = (lowerBound + upperBound) / 2
+    return (scalingCoefficient, offset)
+# End cvp_affine_from_chebyshev
+#
+def cvp_affine_to_chebyshev(lowerBound, upperBound):
+    """
+    Compute the affine parameters to map the interval given
+    as argument to  [-1, 1].
+    @param lowerBound (of the target interval)
+    @param upperBound (of the target interval)
+    @return the (scalingCoefficient, offset) tuple.
+    """
+    ## Check bounds consistency. Bounds must differ.
+    if lowerBound >= upperBound:
+        return None
+    scalingCoefficient = 2 / (upperBound - lowerBound)
+    ## If bounds are symmetrical with relation to 0, return 0 straight before
+    #  attempting a division by 0.
+    if lowerBound == -upperBound:
+        offset         = 0
+    else:
+        offset         =  (lowerBound + upperBound) / (lowerBound - upperBound)
+    return (scalingCoefficient, offset)
+# End cvp_affine_to_chebyshev
+#
 def cvp_babai(redBasis, redBasisGso, vect):
     """
     Closest plane vector implementation as per Babaï.
@@ -34,6 +71,141 @@ def cvp_bkz_reduce(intBasis, blockSize):
     return fplllIntBasis._sage_()
 # End cvp_hkz_reduce.
 #
+def cvp_coefficients_bounds_projection(executablePath, arguments):
+    """
+    Compute the min and max of the coefficients with linear
+    programming projection.
+    @param executablePath: the path to the binary program;
+    @param arguments: a list of arguments to be givent to the binary
+    @return: the min and max coefficients value arrays (in this order).
+    """
+    from subprocess import check_output
+    commandLine = [executablePath] + arguments
+    ## Get rid of output on stderr.
+    devNull = open("/dev/null", "rw")
+    ## Run ther program.
+    otp = check_output(commandLine, stderr=devNull)
+    devNull.close()
+    ## Recover the results
+    bounds = sage_eval(otp)
+    minCoeffsExpList = []
+    maxCoeffsExpList = [] 
+    print "bounds:", bounds
+    for boundsPair in bounds:
+        minCoeffsExpList.append(boundsPair[0])
+        maxCoeffsExpList.append(boundsPair[1])
+    print "minCoeffsExpList:", minCoeffsExpList
+    print "maxCoeffsExpList:", maxCoeffsExpList
+    return (minCoeffsExpList, maxCoeffsExpList)
+# End cvp_coefficients_bounds_projection
+
+def cvp_coefficients_bounds_projection_exps(executablePath, 
+                                            arguments, 
+                                            realField = None):
+    """
+    Compute the min and max exponents of the coefficients with linear
+    programming projection.
+    @param executablePath: the path to the binary program;
+    @param arguments: a list of arguments to be givent to the binary
+    @param realField: the realField to use for floating-point conversion.
+    @return: the min and max exponents arrays (in this order).
+    """
+    ## If no realField is givne, fall back on RR.
+    if realField is None:
+        realField = RR
+    (minCoeffsExpList, maxCoeffsExpList) = \
+        cvp_coefficients_bounds_projection(executablePath,arguments)
+    for index in xrange(0, len(minCoeffsExpList)):
+        minCoeffsExpList[index] = \
+            realField(minCoeffsExpList[index]).log2().floor()
+        maxCoeffsExpList[index] = \
+            realField(maxCoeffsExpList[index]).log2().floor()
+    return (minCoeffsExpList, maxCoeffsExpList)
+# End cvp_coefficients_bounds_projection_exps
+
+def cvp_chebyshev_zeros_1k(numPoints, realField, contFracMaxErr = None):
+    """
+    Compute the Chebyshev zeros for some polynomial degree (numPoints, for the
+    zeros) scaled to the [lowerBound, upperBound] interval.
+    The list is returned as row floating-point numbers is contFracMaxErr is None.
+    Otherwise elements are transformed into rational numbers. 
+    """
+    if numPoints < 1:
+        return None
+    if numPoints == 0:
+        return [0]
+    ## Check that realField has a "prec()" member.
+    try:
+        realField.prec()
+    except:
+        return None
+    #
+    zerosList = []
+    ## formulas: cos (((2*i - 1) * pi) / (2*numPoints)) for i = 1,...,numPoints.
+    #  If i is -1-shifted, as in the following loop, formula is 
+    #  cos (((2*i + 1) * pi) / (2*numPoints)) for i = 0,...,numPoints-1.
+    for index in xrange(0, numPoints):
+        ## When number of points is odd (then, the Chebyshev degree is odd two), 
+        #  the central point is 0. */
+        if (numPoints % 2 == 1) and ((2 * index + 1) == numPoints):
+            if contFracMaxErr is None:
+                zerosList.append(realField(0))
+            else:
+                zerosList.append(0)
+        else:
+            currentCheby = \
+                ((realField(2*index+1) * realField(pi)) / 
+                 realField(2*numPoints)).cos()
+            #print  "Current Cheby:", currentCheby
+            ## Result is negated to have an increasing order list.
+            if contFracMaxErr is None:
+                zerosList.append(-currentCheby)
+            else:
+                zerosList.append(-currentCheby.nearby_rational(contFracMaxErr))
+            #print "Relative error:", (currentCheby/zerosList[index])
+    return zerosList
+# End cvp_chebyshev_zeros_1k.
+#
+def cvp_chebyshev_zeros_for_interval(lowerBound, 
+                                     upperBound, 
+                                     numPoints,
+                                     realField = None, 
+                                     contFracMaxErr = None):
+    """
+    Compute the Chebyshev zeros for some polynomial degree (numPoints, for the
+    zeros) scaled to the [lowerBound, upperBound] interval.
+    If no contFracMaxErr argument is given, we return the list as "raw"
+    floating-points. Otherwise, rational numbers are returned.
+    """
+    ## Arguments check.
+    if lowerBound >= upperBound:
+        return None
+    if numPoints < 1:
+        return None
+    ## If no realField argument is given try to retrieve it from the 
+    #  bounds. If impossible (e.g. rational bound), fall back on RR. 
+    if realField is None:
+        try:
+            ### Force failure if parent does not have prec() member.
+            lowerBound.parent().prec()
+            ### If no exception is thrown, set realField.
+            realField = lowerBound.parent()
+        except:
+            realField = RR
+    #print "Real field:", realField
+    ## We want "raw"floating-point nodes to only have one final rounding.
+    chebyshevZerosList = \
+        cvp_chebyshev_zeros_1k(numPoints, realField)
+    scalingFactor, offset = cvp_affine_from_chebyshev(lowerBound, upperBound)
+    for index in xrange(0, len(chebyshevZerosList)):
+        chebyshevZerosList[index] = \
+            chebyshevZerosList[index] * scalingFactor + offset
+        if not contFracMaxErr is None:
+            chebyshevZerosList[index] = \
+                chebyshevZerosList[index].nearby_rational(contFracMaxErr)
+    return chebyshevZerosList                                                             
+# End cvp_chebyshev_zeros_for_interval.
+#
 def cvp_common_scaling_exp(precision, 
                            precisionFraction = None):
     """
@@ -244,19 +416,30 @@ def cvp_polynomial_coeffs_from_vect(cvpVectInBasis,
     return polynomialCoeffsList
 # En polynomial_coeffs_from_vect.
 #
-
 def cvp_remez_all_poly_error_func_zeros(funct, 
                                         maxDegree, 
                                         lowerBound, 
                                         upperBound,
                                         precsList, 
-                                        realField = None):
+                                        realField = None,
+                                        contFracMaxErr = None):
     """
     For a given function f, a degree and an interval, compute the 
-    zeros of the ||f-remez_d(f)|| function and those of the truncated Remez 
-    ovet the interval of the interval
+    zeros of the (f-remez_d(f)) function and those of the (f-truncRemez_d(f)) 
+    function over the interval. If the (f-truncRemez_d(f)) function has very 
+    few zeros, compute the zeros of the derivative instead.
     """
-    
+    ## If no realField argument is given try to retrieve it from the 
+    #  bounds. If impossible (e.g. rational bound), fall back on RR. 
+    if realField is None:
+        try:
+            ### Force failure if parent does not have prec() member.
+            lowerBound.parent().prec()
+            ### If no exception is thrown, set realField.
+            realField = lowerBound.parent()
+        except:
+            realField = RR
+    #print "Real field:", realField
     funcSa         = func
     funcAsStringSa = func._assume_str().replace("_SAGE_VAR_","",100)
     print "Function as a string:", funcAsStringSa
@@ -267,20 +450,11 @@ def cvp_remez_all_poly_error_func_zeros(funct,
                                            maxDegree, 
                                            lowerBound,
                                            upperBound)
-    ## Compute the infnorm of the error functions.
-    absoluteSo = pobyso_absolute_so_so()
-    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)
-    infNormSo  = pobyso_range_max_abs_so_so(pobyso_supnorm_so_so(pStarSo, 
-                                                                 funcSo, 
-                                                                 intervalSo, 
-                                                                 absoluteSo))
-    print "\nError function infnorm: ", ;pobyso_autoprint(infNormSo)
-    pobyso_clear_obj(infNormSo)
-    #
     ## Compute the zeroes of the error function.
     ### First create the error function with copies since they are needed later.
     errorFuncSo = sollya_lib_build_function_sub(sollya_lib_copy_obj(pStarSo), 
                                                 sollya_lib_copy_obj(funcSo))
+    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)
     errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncSo, intervalSo)
     pobyso_clear_obj(errorFuncSo) # Clears funcSo and pStarSo as well.
     #print "Zeroes of the error function from Sollya:"
@@ -320,13 +494,66 @@ def cvp_remez_all_poly_error_func_zeros(funct,
             pobyso_clear_obj(errorFuncDiffSo) # Clears funcSo and pTruncSo as well.
     pobyso_clear_obj(errorFuncSo) # Clears funcSo and pTruncSo as well.
     ## Sollya objects clean up.
-    pobyso_clear_obj(absoluteSo)
     pobyso_clear_obj(intervalSo)
     errorFuncZerosSa += errorFuncTruncZerosSa
     errorFuncZerosSa.sort()
-    ## Convert the numbers to rationals.
-    for index in xrange(0, len(errorFuncZerosSa)):
-        errorFuncZerosSa[index] = errorFuncZerosSa[index].nearby_rational(1/100000)
+    ## If required, convert the numbers to rational numbers.
+    if not contFracMaxErr is None:
+        for index in xrange(0, len(errorFuncZerosSa)):
+            errorFuncZerosSa[index] = \
+                errorFuncZerosSa[index].nearby_rational(contFracMaxErr)
+    return errorFuncZerosSa
+    #
+# End remez_all_poly_error_func_zeros
+#
+def cvp_remez_poly_error_func_zeros(funct, 
+                                    maxDegree, 
+                                    lowerBound, 
+                                    upperBound,
+                                    realField = None,
+                                    contFracMaxErr = None):
+    """
+    For a given function f, a degree and an interval, compute the 
+    zeros of the (f-remez_d(f)) function.
+    """
+    ## If no realField argument is given try to retrieve it from the 
+    #  bounds. If impossible (e.g. rational bound), fall back on RR. 
+    if realField is None:
+        try:
+            ### Force failure if parent does not have prec() member.
+            lowerBound.parent().prec()
+            ### If no exception is thrown, set realField.
+            realField = lowerBound.parent()
+        except:
+            realField = RR
+    #print "Real field:", realField
+    funcSa         = func
+    funcAsStringSa = func._assume_str().replace("_SAGE_VAR_","",100)
+    #print "Function as a string:", funcAsStringSa
+    ## Create the Sollya version of the function and compute the
+    #  Remez approximant
+    funcSo  = pobyso_parse_string(funcAsStringSa)
+    pStarSo = pobyso_remez_canonical_sa_so(funcAsStringSa, 
+                                           maxDegree, 
+                                           lowerBound,
+                                           upperBound)
+    ## Compute the zeroes of the error function.
+    ### Error function creation consumes both pStarSo and funcSo.
+    errorFuncSo = sollya_lib_build_function_sub(pStarSo, funcSo)
+    intervalSo = pobyso_range_from_bounds_sa_so(lowerBound, upperBound)
+    errorFuncZerosSo = pobyso_dirty_find_zeros_so_so(errorFuncSo, intervalSo)
+    pobyso_clear_obj(errorFuncSo) # Clears funcSo and pStarSo as well.
+    #print "Zeroes of the error function from Sollya:"
+    #pobyso_autoprint(errorFuncZerosSo)
+    errorFuncZerosSa = pobyso_float_list_so_sa(errorFuncZerosSo)
+    pobyso_clear_obj(errorFuncZerosSo)
+    ## Sollya objects clean up.
+    pobyso_clear_obj(intervalSo)
+    ## Converting to rational numbers, if required.
+    if not contFracMaxErr is None:
+        for index in xrange(0, len(errorFuncZerosSa)):
+            errorFuncZerosSa[index] = \
+                errorFuncZerosSa[index].nearby_rational(contFracMaxErr)
     return errorFuncZerosSa
     #
 # End remez_poly_error_func_zeros
