Laboratoire de l'Informatique et du Parallélisme » Pobyso

def pobyso_polynomial_coefficients_progressive_round_so_so(polySo,

                                                           funcSo,

                                                           icSo,

                                                           intervalSo,

                                                           itpSo,

                                                           ftpSo,

                                                           maxPrecSo,

                                                           maxErrSo):

    ## Higher order function see: 

    #  http://effbot.org/pyfaq/how-do-you-make-a-higher-order-function-in-python.htm

    def precision_decay_ratio_function(degreeSa):

        def outer(x):

            def inner(x):

                we = 3/8

                wq = 2/8

                a  = 2.2

                b  = 2 

                return we*(exp(x/a)-1) +  wq*((b*x)**2) + (1-we-wq)*x

            return  inner(x)/inner(degreeSa)

        return outer

    #

    degreeSa        = pobyso_polynomial_degree_so_sa(polySo)

    ratio           = precision_decay_ratio_function(degreeSa)

    itpSa           = pobyso_constant_from_int_so_sa(itpSo)

    ftpSa           = pobyso_constant_from_int_so_sa(ftpSo)

    maxPrecSa       = pobyso_constant_from_int_so_sa(maxPrecSo)

    maxErrSa        = pobyso_get_constant_as_rn_so_sa(maxErrSo)

    resPolySo       = pobyso_constant_0_sa_so()

    lastResPolySo   = None

    while True: 

        pDeltaSa    = ftpSa - itpSa

        for indexSa in reversed(xrange(0,degreeSa+1)):

            print "Index:", indexSa

            indexSo = pobyso_constant_from_int_sa_so(indexSa)

            #pobyso_autoprint(indexSo)

            #print ratio(indexSa)

            ctpSa = floor(ftpSa - (pDeltaSa * ratio(indexSa)))

            ctpSo = pobyso_constant_from_int_sa_so(ctpSa)

            print "Index:", indexSa, " - Target precision:", 

            pobyso_autoprint(ctpSo)

            cmonSo  = \

                sollya_lib_build_function_mul(sollya_lib_coeff(polySo, indexSo),

                                      sollya_lib_build_function_pow( \

                                          sollya_lib_build_function_free_variable(), \

                                          indexSo))

            #pobyso_autoprint(cmonSo)

            cmonrSo = pobyso_round_coefficients_single_so_so(cmonSo, ctpSo)

            sollya_lib_clear_obj(cmonSo)

            #pobyso_autoprint(cmonrSo)

            resPolySo = sollya_lib_build_function_add(resPolySo,

                                                      cmonrSo)

            pobyso_autoprint(resPolySo)

        # End for index

        freeVarSo     = sollya_lib_build_function_free_variable()

        changeVarSo   = sollya_lib_sub(freeVarSo, icSo)

        resPolyCvSo   = sollya_lib_evaluate(resPolySo, changeVarSo)

        errFuncSo = sollya_lib_build_function_sub(funcSo, resPolyCvSo)

        infNormSo = sollya_lib_dirtyinform(errFuncSo, intervalSo)

        print "Infnorm (Sollya):", ; pobyso_autoprint(infNormSo)

        cerrSa    = pobyso_get_constant_as_rn_so_sa(infNormSo)

        sollya_lib_clear_obj(resPolyCvSo)

        sollya_lib_clear_obj(infNormSo)

        print "Infnorm  (Sage):", cerrSa

        if (cerrSa > maxErrSa):

            print "Error is too large."

            if not lastResPolySo is None:

                print "Enlarging prec."

                ntpSa = floor(ftpSa + ftpSa/50)

                ## Can't enlarge (numerical)

                if ntpSa == ftpSa:

                    sollya_lib_clear_obj(lastResPolySo)

                    sollya_lib_clear_obj(resPolySo)

                    return None

                ## Can't enlarge (not enough precision left)

                if ntpSa > maxPrecSa:

                    sollya_lib_clear_obj(lastResPolySo)

                    sollya_lib_clear_obj(resPolySo)

                    return None

                ftpSa = ntpSa

                lastResPolySo = resPolySo

                continue

            ## One enlargement took place.

            else:

                sollya_lib_clear_obj(resPolySo)

                return lastResPolySo

        # cerrSa <= maxErrSa: scrap more bits.

        else:

            print "Shrinking prec."

            ntpSa = floor(ftpSa - ftpSa/50)

            ## Can't shrink (numerical)

            if ntpSa == ftpSa:

                if not lastResPolySo is None:

                    sollya_lib_clear_obj(lastResPolySo)

                return resPolySo

            ## Can't enlarge (not enough precision left)

            if ntpSa <= itpSa:

                if not lastResPolySo is None:

                    sollya_lib_clear_obj(lastResPolySo)

                return resPolySo

            ftpSa = ntpSa

            continue

    # End wile True

# End pobyso_polynomial_coefficients_progressive_truncate_so_so.

def pobyso_polynomial_degree_so_sa(polySo):

    """

    Return the degree of a Sollya polynomial as a Sage int.

    """

    degreeSo = sollya_lib_degree(polySo)

    return pobyso_constant_from_int_so_sa(degreeSo)

# End pobyso_polynomial_degree_so_sa

def pobyso_polynomial_degree_so_so(polySo):

    """

    Thin wrapper around lib_sollya_degree().

    """

    return sollya_lib_degree(polySo)

# End pobyso_polynomial_degree_so_so

    sollya_lib_clear_obj(endEllipListSo)