Ticket #10532: mu2.sage

File mu2.sage, 1.6 KB (added by pernici, 9 years ago)

benchmark

Line 
1from time import time
2import os
3import sys
4
5def lhsums(a):
6  if len(a) %2 == 1:
7    a = a[1:]
8  n = len(a)/2
9  #print a[:n]
10  #print a[n:]
11  return float(sum(a[n:]))/sum(a[:n])
12
13verbose=1
14
15try:
16  F = sys.argv[1]
17  n = int(sys.argv[2])
18  h = int(sys.argv[3])
19  B = int(sys.argv[4])
20  N = int(sys.argv[5])
21  assert n > 1
22except:
23  print '''prog F n h N
24  F = field: QQ, RR, GF7
25  n = number of variables in PolynomialRing > 1
26  h = precision
27  B = bound
28  N = number of times square is taken
29  '''
30  sys.exit()
31
32if F == 'QQ':
33    field = QQ
34elif F == 'RR':
35    field = RR
36elif F == 'GF7':
37    field = GF(7)
38else:
39   print 'field not allowed'
40   sys.exit()
41
42print 'using', field
43
44if field == QQ:
45  field_ext = ''
46else:
47  field_ext = '_' + F
48
49R = PowerSeriesRing(field,n,'t')
50d = R.gens_dict()
51d1 = {}
52for i in range(n):
53  d1[d.keys()[i]] = d.values()[i] + 1
54#print d
55print d1
56fn = 'dat/dat_%s%s_%d_%d_%d' %(sys.argv[0][:-3],field_ext,n,h,B)
57print fn
58
59if os.path.exists(fn):
60  print 'use polynomial in', fn
61  s = open(fn).read()
62  ns = s.find('+ O')
63  s = s[:ns]
64  p = sage_eval(s,d) + R(0).O(h)
65else:
66  print 'choose a random polynomial'
67  p = R.random_element(prec=h,bound=B)
68  #print 'p=',p
69  pp = p.polynomial()
70  pp = pp.subs(**d1)
71  p = pp + R(0).O(h)
72  #print 'new p=',p
73  f = open(fn, 'w')
74  f.write('%s' % p)
75
76pw =1
77for k in range(N):
78  if verbose:
79    a = [len(x.coefficients()) for x in p._bg_value.coefficients()]
80    r = lhsums(a)
81    if verbose == 2:
82      print a
83    print 'pw=%d ratio=%.2f' % (pw, r)
84  t0 = time()
85  p = p*(p+1)
86  #print pw, p
87  pw*=2
88  t1 = time()
89  print '%d %.4f' %(k, t1-t0)