Statistiques
| Révision :

root / pobysoPythonSage / src / sageSLZ / runSLZ-53proj.sage @ 268

Historique | Voir | Annoter | Télécharger (5,59 ko)

1
#! /opt/sage/sage
2
# @file runSLZ-113proj.sage
3
#
4
#@par Changes from runSLZ-113.sage
5
# LLL reduction is not performed on the matrix itself but rather on the
6
# product of the matrix with a uniform random matrix.
7
# The reduced matrix obtained is discarded but the transformation matrix 
8
# obtained is used to multiply the original matrix in order to reduced it.
9
# If a sufficient level of reduction is obtained, we stop here. If not
10
# the product matrix obtained above is LLL reduced. But as it has been
11
# pre-reduced at the above step, reduction is supposed to be much faster.
12
#
13
# Both reductions combined should hopefully be faster than a straight single
14
# reduction.
15
#
16
# Run SLZ for p=113
17
#from scipy.constants.codata import precision
18
def initialize_env():
19
    """
20
    Load all necessary modules.
21
    """
22
    compiledSpyxDir = "/home/storres/recherche/arithmetique/pobysoPythonSage/compiledSpyx"
23
    if compiledSpyxDir not in sys.path:
24
        sys.path.append(compiledSpyxDir)
25
    if not 'mpfi' in sage.misc.cython.standard_libs:
26
        sage.misc.cython.standard_libs.append('mpfi')
27
    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sollya_lib.sage")
28
#    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageMpfr.spyx")
29
#    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageGMP.spyx")
30
    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/pobyso.py")
31
    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageSLZ.sage")
32
    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageNumericalOperations.sage")
33
    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageRationalOperations.sage")
34
    # Matrix operations are loaded by polynomial operations.
35
    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sagePolynomialOperations.sage")
36
    load("/home/storres/recherche/arithmetique/pobysoPythonSage/src/sageSLZ/sageRunSLZ.sage")
37

    
38

    
39
print "Running SLZ..."
40
initialize_env()
41
from sageMpfr import *
42
from sageGMP  import *
43
import sys
44
from subprocess import call
45
#
46
## Main variables and parameters.
47
x         = var('x')
48
func(x)   = exp(x)
49
precision = 53
50
emin      = -1022
51
emax      = 1023 
52
RRR = RealField(precision)
53
degree              = 0
54
alpha               = 0
55
htrn                = 0
56
intervalCenter      = 0
57
intervalRadius      = 0
58
debugMode           = False          
59
## Local functions
60
#
61
def usage():
62
    write = sys.stderr.write
63
    write("\nUsage:\n")
64
    write("  " + scriptName + " <degree> <alpha> <htrn> <intervalCenter>\n")
65
    write("               <numberOfNumbers> [debug]\n")
66
    write("\nArguments:\n")
67
    write("  degree          the degree of the polynomial (integer)\n")
68
    write("  alpha           alpha (integer)\n")
69
    write("  htrn            hardness-to-round - a number of bits (integer)\n")
70
    write("  intervalCenter  the interval center (a floating-point number)\n")
71
    write("  numberOfNumbers the number of floating-point numbers in the interval\n")
72
    write("                  as a positive integral expression\n")
73
    write("  debug           debug mode (\"debug\", in any case)\n\n")
74
    sys.exit(2)
75
# End usage.
76
#
77
argsCount = len(sys.argv)
78
scriptName = os.path.basename(__file__)
79
if argsCount < 5:
80
    usage()
81
for index in xrange(1,argsCount):
82
    if index == 1:
83
        degree = int(sys.argv[index])
84
    elif index == 2:
85
        alpha = int(sys.argv[index])
86
    elif index == 3:
87
        htrn = int(eval(sys.argv[index]))
88
    elif index == 4:
89
        try:
90
            intervalCenter = QQ(sage_eval(sys.argv[index]))
91
        except:
92
            intervalCenter = RRR(sys.argv[index])
93
        intervalCenter = RRR(intervalCenter)
94
    elif index == 5:
95
        ## Can be read as rational number but must end up as an integer.
96
        numberOfNumbers = QQ(sage_eval(sys.argv[index]))
97
        if numberOfNumbers != numberOfNumbers.round():
98
            raise Exception("Invalid number of numbers: " + sys.argv[index] + ".")
99
        numberOfNumbers = numberOfNumbers.round()
100
        ## The number must be strictly positive.
101
        if numberOfNumbers <= 0:
102
            raise Exception("Invalid number of numbers: " + sys.argv[index] + ".")
103
    elif index == 6:
104
        debugMode = sys.argv[index].upper()
105
        debugMode = (debugMode == "DEBUG")
106
# Done with command line arguments collection.
107
#
108
## Debug printing
109
print "degree         :", degree
110
print "alpha          :", alpha
111
print "htrn           :", htrn
112
print "interval center:", intervalCenter.n(prec=10).str(truncate=False)
113
print "num of nums    :", RR(numberOfNumbers).log2().n(prec=10).str(truncate=False)
114
print "debug mode     :", debugMode
115
print
116
#
117
## Set the terminal window title.
118
terminalWindowTitle = ['stt', str(degree), str(alpha), str(htrn), 
119
                       intervalCenter.n(prec=10).str(truncate=False), 
120
                       RRR(numberOfNumbers).log2().n(prec=10).str(truncate=False)]
121
call(terminalWindowTitle)
122
#
123
intervalCenterBinade = slz_compute_binade(intervalCenter)
124
intervalRadius       = \
125
    2^intervalCenterBinade * 2^(-precision + 1) * numberOfNumbers / 2
126
srs_run_SLZ_v05_proj(inputFunction=func, 
127
                     inputLowerBound         = intervalCenter - intervalRadius, 
128
                     inputUpperBound         = intervalCenter + intervalRadius, 
129
                     alpha                   = alpha, 
130
                     degree                  = degree, 
131
                     precision               = precision, 
132
                     emin                    = emin, 
133
                     emax                    = emax, 
134
                     targetHardnessToRound   = htrn, 
135
                     debug                   = debugMode)
136