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

#! /opt/sage/sage

#from scipy.constants.codata import precision

def initialize_env():

    """

    Load all necessary modules.

    """

    if not 'mpfi' in sage.misc.cython.standard_libs:

        sage.misc.cython.standard_libs.append('mpfi')

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sollya_lib.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageMpfr.spyx")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageGMP.spyx")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/pobyso.py")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageSLZ.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageNumericalOperations.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageRationalOperations.sage")

    # Matrix operations are loaded by polynomial operations.

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sagePolynomialOperations.sage")

    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageRunSLZ.sage")

print "Running HTRN..."

initialize_env()

#

f(x)                  = exp(x)

precision             = 53

htrAccuracyShift      = 50#58

precisionRF           = RealField(precision)

precPlusOneRF         = RealField(precisionRF.prec()+1)

quasiExactRF          = RealField(1024)

targetAccuracy        = 2^-(precision + htrAccuracyShift)

binade                = 2^-31

binadeAbsLowerBound   = precisionRF(binade)

binadeUlp             = binadeAbsLowerBound.ulp() 

binadeLowerBound      = binadeAbsLowerBound

htrnc                 = precisionRF("1.9E9CBBFD6080B",16)  * 2^-31

print binade.n()

print binadeAbsLowerBound.n()

htrncOffset  = htrnc - binadeLowerBound

iterStep     = binadeUlp * 2^22

iterVar      = binadeLowerBound + floor(htrncOffset/iterStep) * iterStep

iterLimit    = iterVar + iterStep

print iterVar.n()

print htrnc.n()

print iterLimit.n()

iterCount = 0

iterationsStartTime = cputime()

kiloIterStartTime = cputime()

kiloIterCount     = 0

kiloIterCountStep = 2^20

print "Starting iterations..."

while iterVar < iterLimit:

    if iterVar == htrnc:

        print "Found at iter count:", iterCount

        print "After:", cputime(iterationsStartTime)

    if slz_is_htrn(argument=iterVar, 

                  function=f, 

                  targetAccuracy=targetAccuracy, 

                  targetRF=precisionRF, 

                  targetPlusOnePrecRF=precPlusOneRF, 

                  quasiExactRF=quasiExactRF):

        print "HTRN candidate:", iterVar.n().str(base=2)

        print "HTRN candidate:", iterVar.n().str(base=10)

        print "HTRN candidate:", iterVar.n().str(base=16)

        print

    if kiloIterCount == kiloIterCountStep:

        print "1024 iterations computed in:", 

        print cputime(kiloIterStartTime), "s"

        kiloIterCount = 0

        kiloIterStartTime = cputime()

    # End if

    iterVar         += binadeUlp

    iterCount       += 1

    kiloIterCount   += 1

iterationsFullTime = cputime(iterationsStartTime)

print "Iterations terminates in", iterationsStartTime, "s."