Révision 101

pobysoPythonSage/src/sageSLZ/sageSLZ.sage (revision 101)
36 36 
                                                    absoluteErrorTypeSo)
37 37 
    maxErrorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
38 38 
    while maxErrorSa > approxPrecSa:
39 
        #print "++Approximation error:", maxErrorSa
39 40 
        sollya_lib_clear_obj(maxErrorSo)
40 41 
        sollya_lib_clear_obj(polySo)
41 42 
        sollya_lib_clear_obj(intervalCenterSo)
42 
        shrinkFactorSa = RRR('5.0')/(maxErrorSa/approxPrecSa).log2().abs()
43 
        shrinkFactorSa = RRR('5')/(maxErrorSa/approxPrecSa).log2().abs()
43 44 
        #shrinkFactorSa = 1.5/(maxErrorSa/approxPrecSa)
44 45 
        #errorRatioSa = approxPrecSa/maxErrorSa
45 46 
        #print "Error ratio: ", errorRatioSa
46
47 47 
        if shrinkFactorSa > intervalShrinkConstFactorSa:
48 48 
            actualShrinkFactorSa = intervalShrinkConstFactorSa
49 49 
            #print "Fixed"
......
51 51 
            actualShrinkFactorSa = shrinkFactorSa
52 52 
            #print "Computed",shrinkFactorSa,maxErrorSa
53 53 
            #print shrinkFactorSa, maxErrorSa
54 
        #print "Shrink factor", actualShrinkFactorSa
54 55 
        currentUpperBoundSa = currentLowerBoundSa + \
55 56 
                                  (currentUpperBoundSa - currentLowerBoundSa) * \
56 57 
                                   actualShrinkFactorSa
57 58 
        #print "Current upper bound:", currentUpperBoundSa
58 59 
        sollya_lib_clear_obj(currentRangeSo)
59 60 
        sollya_lib_clear_obj(polySo)
60 
        if currentUpperBoundSa <= currentLowerBoundSa:
61 
        if currentUpperBoundSa <= currentLowerBoundSa or \
62 
              currentUpperBoundSa == currentLowerBoundSa.nextabove():
61 63 
            sollya_lib_clear_obj(absoluteErrorTypeSo)
62 64 
            print "Can't find an interval."
63 65 
            print "Use either or both a higher polynomial degree or a higher",
......
210 212 
    - the corresponding approximation error.
211 213 
    TODO: fix endless looping for some parameters sets.
212 214 
    """
215 
    # Set Sollya to the necessary internal precision.
213 216 
    currentSollyaPrecSo = pobyso_get_prec_so()
214 217 
    currentSollyaPrecSa = pobyso_constant_from_int_so_sa(currentSollyaPrecSo)
215 218 
    if internalSollyaPrecSa > currentSollyaPrecSa:
216 219 
        pobyso_set_prec_sa_so(internalSollyaPrecSa)
220 
    #
217 221 
    x = functionSa.variables()[0] # Actual variable name can be anything.
222 
    # Scaled function: [1=,2] -> [1,2].
218 223 
    (fff, scaledLowerBoundSa, scaledUpperBoundSa, \
219 224 
            scaledLowerBoundImageSa, scaledUpperBoundImageSa) = \
220 225 
                slz_compute_scaled_function(functionSa, \
......
236 241 
                                        scaledLowerBoundSa, scaledUpperBoundSa,
237 242 
                                        approxPrecSa, internalSollyaPrecSa)
238 243 
    if polySo is None:
239 
        print "Aborting"
244 
        print "slz_get_intervals_and_polynomials: Aborting and returning None!"
245 
        if internalSollyaPrecSa != currentSollyaPrecSa:
246 
            pobyso_set_prec_sa_so(currentSollyaPrecSa)
240 247 
        return None
241 
    # Change variable stuff
242 
    changeVarExpressionSo = sollya_lib_build_function_sub(
243 
                              sollya_lib_build_function_free_variable(),
244 
                              sollya_lib_copy_obj(intervalCenterSo))
245 
    polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo)
246 
    resultArray.append((polySo, polyVarChangedSo, boundsSo, intervalCenterSo,\
247 
                         maxErrorSo))
248 248 
    realIntervalField = RealIntervalField(max(lowerBoundSa.parent().precision(),
249 249 
                                              upperBoundSa.parent().precision()))
250 250 
    boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)
251 
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
252 
    #print "First approximation error:", errorSa
253 
    # If the error and interval are OK a the first try, just return.
254 
    if boundsSa.endpoints()[1] >= scaledUpperBoundSa:
255 
        # Change variable stuff in Sollya x -> x0-x.
256 
        changeVarExpressionSo = sollya_lib_build_function_sub( \
257 
                              sollya_lib_build_function_free_variable(), \
258 
                              sollya_lib_copy_obj(intervalCenterSo))
259 
        polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo)
260 
        resultArray.append((polySo, polyVarChangedSo, boundsSo, \
261 
                            intervalCenterSo, maxErrorSo))
262 
        if internalSollyaPrecSa != currentSollyaPrecSa:
263 
            pobyso_set_prec_sa_so(currentSollyaPrecSa)
264 
        sollya_lib_clear_obj(functionSo)
265 
        sollya_lib_clear_obj(degreeSo)
266 
        sollya_lib_clear_obj(scaledBoundsSo)
267 
        #print "Approximation error:", errorSa
268 
        return resultArray
251 269 
    # Compute the next upper bound.
252 
    # If the error of approximation is more than half of the target,
253 
    # use the same interval.
254 
    # If it is less, increase it a bit.
255 
    errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
270 
    # The following ratio is always >= 1
256 271 
    currentErrorRatio = approxPrecSa / errorSa
272 
    # Starting point for the next upper bound.
257 273 
    currentScaledUpperBoundSa = boundsSa.endpoints()[1]
258 
    if currentErrorRatio  < 2 :
274 
    boundsWidthSa = boundsSa.endpoints()[1] - boundsSa.endpoints()[0]
275 
    # Compute the increment.
276 
    if currentErrorRatio > RR('1000'): # ]1.5, infinity[
259 277 
        currentScaledUpperBoundSa += \
260 
            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0])
261 
    else:
278 
                        currentErrorRatio * boundsWidthSa * 2
279 
    else:  # [1, 1.5]
262 280 
        currentScaledUpperBoundSa += \
263 
            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0]) \
264 
                         * currentErrorRatio.log2() * 2
281 
                        (RR('1.0') + currentErrorRatio.log() / 500) * boundsWidthSa
282 
    # Take into account the original interval upper bound.
265 283 
    if currentScaledUpperBoundSa > scaledUpperBoundSa:
266 284 
        currentScaledUpperBoundSa = scaledUpperBoundSa
267 285 
    # Compute the other expansions.
......
273 291 
                                            currentScaledUpperBoundSa,
274 292 
                                            approxPrecSa,
275 293 
                                            internalSollyaPrecSa)
276 
        # Change variable stuff
277 
        changeVarExpressionSo = sollya_lib_build_function_sub(
278 
                                  sollya_lib_build_function_free_variable(),
279 
                                  sollya_lib_copy_obj(intervalCenterSo))
280 
        polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo)
281 
        resultArray.append((polySo, polyVarChangedSo, boundsSo, \
282 
                            intervalCenterSo, maxErrorSo))
294 
        errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
295 
        if  errorSa < approxPrecSa:
296 
            # Change variable stuff
297 
            #print "Approximation error:", errorSa
298 
            changeVarExpressionSo = sollya_lib_build_function_sub(
299 
                                    sollya_lib_build_function_free_variable(),
300 
                                    sollya_lib_copy_obj(intervalCenterSo))
301 
            polyVarChangedSo = sollya_lib_evaluate(polySo, changeVarExpressionSo)
302 
            resultArray.append((polySo, polyVarChangedSo, boundsSo, \
303 
                                intervalCenterSo, maxErrorSo))
283 304 
        boundsSa = pobyso_range_to_interval_so_sa(boundsSo, realIntervalField)
284 305 
        # Compute the next upper bound.
285 
        # If the error of approximation is more than half of the target,
286 
        # use the same interval.
287 
        # If it is less, increase it a bit.
288 
        errorSa = pobyso_get_constant_as_rn_with_rf_so_sa(maxErrorSo)
306 
        # The following ratio is always >= 1
289 307 
        currentErrorRatio = approxPrecSa / errorSa
290 
        if currentErrorRatio  < RR('1.5') :
291 
            currentScaledUpperBoundSa = \
292 
                            boundsSa.endpoints()[1] + \
293 
                            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0])
294 
        elif currentErrorRatio < 2:
295 
            currentScaledUpperBoundSa = \
296 
                            boundsSa.endpoints()[1] + \
297 
                            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0]) \
298 
                             * currentErrorRatio.log2()
299 
        else:
300 
            currentScaledUpperBoundSa = \
301 
                            boundsSa.endpoints()[1] + \
302 
                            (boundsSa.endpoints()[1] - boundsSa.endpoints()[0]) \
303 
                             * currentErrorRatio.log2() * 2
308 
        # Starting point for the next upper bound.
309 
        currentScaledUpperBoundSa = boundsSa.endpoints()[1]
310 
        boundsWidthSa = boundsSa.endpoints()[1] - boundsSa.endpoints()[0]
311 
        # Compute the increment.
312 
        if currentErrorRatio > RR('1000'): # ]1.5, infinity[
313 
            currentScaledUpperBoundSa += \
314 
                            currentErrorRatio * boundsWidthSa * 2
315 
        else:  # [1, 1.5]
316 
            currentScaledUpperBoundSa += \
317 
                            (RR('1.0') + currentErrorRatio.log()/500) * boundsWidthSa
318 
        #print "currentErrorRatio:", currentErrorRatio
319 
        #print "currentScaledUpperBoundSa", currentScaledUpperBoundSa
304 320 
        # Test for insufficient precision.
305 321 
        if currentScaledUpperBoundSa == scaledLowerBoundSa:
306 322 
            print "Can't shrink the interval anymore!"
307 323 
            print "You should consider increasing the Sollya internal precision"
308 324 
            print "or the polynomial degree."
309 325 
            print "Giving up!"
326 
            if internalSollyaPrecSa != currentSollyaPrecSa:
327 
                pobyso_set_prec_sa_so(currentSollyaPrecSa)
310 328 
            sollya_lib_clear_obj(functionSo)
311 329 
            sollya_lib_clear_obj(degreeSo)
312 330 
            sollya_lib_clear_obj(scaledBoundsSo)

Formats disponibles : Unified diff