Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

""" 

Enumeration of Primitive Totally Real Fields 

  

This module contains functions for enumerating all primitive 

totally real number fields of given degree and small discriminant. 

Here a number field is called *primitive* if it contains no proper 

subfields except `\QQ`. 

  

See also :mod:`sage.rings.number_field.totallyreal_rel`, which handles the non-primitive 

case using relative extensions. 

  

Algorithm 

--------- 

  

We use Hunter's algorithm ([Cohen2000]_, Section 9.3) with modifications 

due to Takeuchi [Takeuchi1999]_ and the author [Voight2008]_. 

  

We enumerate polynomials `f(x) = x^n + a_{n-1} x^{n-1} + \dots + a_0`. 

Hunter's theorem gives bounds on `a_{n-1}` and `a_{n-2}`; then given 

`a_{n-1}` and `a_{n-2}`, one can recursively compute bounds on `a_{n-3}, 

\dots, a_0`, using the fact that the polynomial is totally real by 

looking at the zeros of successive derivatives and applying 

Rolle's theorem. See [Takeuchi1999]_ for more details. 

  

Examples 

-------- 

  

In this first simple example, we compute the totally real quadratic 

fields of discriminant `\le 50`. 

  

:: 

  

sage: enumerate_totallyreal_fields_prim(2,50) 

[[5, x^2 - x - 1], 

[8, x^2 - 2], 

[12, x^2 - 3], 

[13, x^2 - x - 3], 

[17, x^2 - x - 4], 

[21, x^2 - x - 5], 

[24, x^2 - 6], 

[28, x^2 - 7], 

[29, x^2 - x - 7], 

[33, x^2 - x - 8], 

[37, x^2 - x - 9], 

[40, x^2 - 10], 

[41, x^2 - x - 10], 

[44, x^2 - 11]] 

sage: [ d for d in range(5,50) if (is_squarefree(d) and d%4 == 1) or (d%4 == 0 and is_squarefree(d/4)) ] 

[5, 8, 12, 13, 17, 20, 21, 24, 28, 29, 33, 37, 40, 41, 44] 

  

Next, we compute all totally real quintic fields of discriminant `\le 10^5`:: 

  

sage: ls = enumerate_totallyreal_fields_prim(5,10^5) ; ls 

[[14641, x^5 - x^4 - 4*x^3 + 3*x^2 + 3*x - 1], 

[24217, x^5 - 5*x^3 - x^2 + 3*x + 1], 

[36497, x^5 - 2*x^4 - 3*x^3 + 5*x^2 + x - 1], 

[38569, x^5 - 5*x^3 + 4*x - 1], 

[65657, x^5 - x^4 - 5*x^3 + 2*x^2 + 5*x + 1], 

[70601, x^5 - x^4 - 5*x^3 + 2*x^2 + 3*x - 1], 

[81509, x^5 - x^4 - 5*x^3 + 3*x^2 + 5*x - 2], 

[81589, x^5 - 6*x^3 + 8*x - 1], 

[89417, x^5 - 6*x^3 - x^2 + 8*x + 3]] 

sage: len(ls) 

9 

  

We see that there are 9 such fields (up to isomorphism!). 

  

References 

---------- 

  

.. [Cohen2000] Henri Cohen, Advanced topics in computational number 

theory, Graduate Texts in Mathematics, vol. 193, 

Springer-Verlag, New York, 2000. 

  

.. [Martinet1980] Jacques Martinet, Petits discriminants des corps de nombres, Journ. Arithm. 1980, 

Cambridge Univ. Press, 1982, 151--193. 

  

.. [Takeuchi1999] Kisao Takeuchi, Totally real algebraic number fields of 

degree 9 with small discriminant, Saitama Math. J. 

17 (1999), 63--85 (2000). 

  

.. [Voight2008] John Voight, Enumeration of totally real number fields of bounded root 

discriminant, Lect. Notes in Comp. Sci. 5011 (2008). 

  

Authors 

------- 

  

- John Voight (2007-09-01): Initial version. 

- John Voight (2007-09-19): Various optimization tweaks. 

- John Voight (2007-10-09): Added DSage module. 

- John Voight (2007-10-17): Added pari functions to avoid recomputations. 

- John Voight (2007-10-27): Separated DSage component. 

- Craig Citro and John Voight (2007-11-04): Additional doctests and type checking. 

- Craig Citro and John Voight (2008-02-10): Final modifications for submission. 

  

------ 

""" 

  

#***************************************************************************** 

# Copyright (C) 2007 William Stein and John Voight 

# 

# This program is free software: you can redistribute it and/or modify 

# it under the terms of the GNU General Public License as published by 

# the Free Software Foundation, either version 2 of the License, or 

# (at your option) any later version. 

# http://www.gnu.org/licenses/ 

#***************************************************************************** 

  

from __future__ import absolute_import, print_function 

  

from cysignals.memory cimport check_calloc, sig_free 

  

import math 

import sys 

  

from sage.libs.gmp.mpz cimport * 

from sage.libs.pari.all import pari 

from cypari2.gen cimport Gen as pari_gen 

from sage.libs.pari.misc cimport new_t_POL_from_int_star 

  

from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing 

from sage.rings.integer import Integer 

from sage.rings.integer cimport Integer 

from sage.rings.integer_ring import IntegerRing 

from sage.rings.all import ZZ, QQ 

from sage.misc.misc import cputime 

  

from sage.rings.number_field.totallyreal_data import tr_data, int_has_small_square_divisor 

from sage.rings.number_field.totallyreal_data cimport tr_data 

  

  

cpdef double odlyzko_bound_totallyreal(int n): 

r""" 

This function returns the unconditional Odlyzko bound for the root 

discriminant of a totally real number field of degree n. 

  

.. NOTE:: 

  

The bounds for n > 50 are not necessarily optimal. 

  

INPUT: 

  

- n (integer) the degree 

  

OUTPUT: 

  

a lower bound on the root discriminant (as a real number) 

  

EXAMPLES:: 

  

sage: [sage.rings.number_field.totallyreal.odlyzko_bound_totallyreal(n) for n in range(1,5)] 

[1.0, 2.223, 3.61, 5.067] 

  

AUTHORS: 

  

- John Voight (2007-09-03) 

  

NOTES: 

The values are calculated by Martinet [Martinet1980]_. 

""" 

  

if n <= 10: 

dB = [1.,2.223,3.610,5.067,6.523,7.941,9.301,10.596,11.823,12.985][n-1] 

elif n <= 20: 

dB = [14.0831,15.1217,16.1047,17.0359,17.9192,18.7580,19.5555,20.3148,21.0386,21.7294][n-11] 

elif n <= 30: 

dB = [22.3896,23.0212,23.6261,24.2061,24.7628,25.2976,25.2976,26.3071,26.3071,27.2440][n-21] 

elif n <= 40: 

dB = [27.2440,28.1165,28.1165,28.9315,28.9315,29.6948,29.6948,30.4117,30.4117,31.0865][n-31] 

elif n <= 50: 

dB = [31.0865,31.7232,31.7232,32.3252,32.3252,32.8954,32.8954,33.4365,33.4365,33.9508][n-41] 

else: 

dB = 33.9508 

return dB 

  

  

def enumerate_totallyreal_fields_prim(n, B, a = [], verbose=0, return_seqs=False, 

phc=False, keep_fields=False, t_2=False, 

just_print=False, 

return_pari_objects=True): 

r""" 

This function enumerates primitive totally real fields of degree 

`n>1` with discriminant `d \leq B`; optionally one can specify the 

first few coefficients, where the sequence `a` corresponds to 

  

:: 

  

a[d]*x^n + ... + a[0]*x^(n-d) 

  

where ``length(a) = d+1``, so in particular always ``a[d] = 1``. 

  

.. NOTE:: 

  

This is guaranteed to give all primitive such fields, and 

seems in practice to give many imprimitive ones. 

  

INPUT: 

  

- ``n`` -- (integer) the degree 

- ``B`` -- (integer) the discriminant bound 

- ``a`` -- (list, default: []) the coefficient list to begin with 

- ``verbose`` -- (integer or string, default: 0) if ``verbose == 1`` 

(or ``2``), then print to the screen (really) verbosely; if verbose is 

a string, then print verbosely to the file specified by verbose. 

- ``return_seqs`` -- (boolean, default False) If ``True``, then return 

the polynomials as sequences (for easier exporting to a file). 

- ``phc`` -- boolean or integer (default: False) 

- ``keep_fields`` -- (boolean or integer, default: False) If 

``keep_fields`` is True, then keep fields up to ``B*log(B)``; if 

``keep_fields`` is an integer, then keep fields up to that integer. 

- ``t_2`` -- (boolean or integer, default: False) If ``t_2 = T``, then 

keep only polynomials with t_2 norm >= T. 

- ``just_print`` -- (boolean, default: False): if ``just_print`` is not 

False, instead of creating a sorted list of totally real number 

fields, we simply write each totally real field we find to the file 

whose filename is given by ``just_print``. In this case, we don't 

return anything. 

- ``return_pari_objects`` -- (boolean, default: True) if 

both ``return_seqs`` and ``return_pari_objects`` are ``False`` then 

it returns the elements as Sage objects; otherwise it returns pari 

objects. 

  

OUTPUT: 

  

the list of fields with entries ``[d,f]``, where ``d`` is the 

discriminant and ``f`` is a defining polynomial, sorted by 

discriminant. 

  

  

AUTHORS: 

  

- John Voight (2007-09-03) 

- Craig Citro (2008-09-19): moved to Cython for speed improvement 

  

TESTS:: 

  

sage: len(enumerate_totallyreal_fields_prim(2,10**4)) 

3043 

sage: len(enumerate_totallyreal_fields_prim(3,3**8)) 

237 

sage: len(enumerate_totallyreal_fields_prim(5,5**7)) 

6 

sage: len(enumerate_totallyreal_fields_prim(2,2**15)) # long time 

9957 

sage: len(enumerate_totallyreal_fields_prim(3,3**10)) # long time 

2720 

sage: len(enumerate_totallyreal_fields_prim(5,5**8)) # long time 

103 

  

Each of the outputs must be elements of Sage if ``return_pari_objects`` 

is set to ``False``:: 

  

sage: enumerate_totallyreal_fields_prim(2, 10) 

[[5, x^2 - x - 1], [8, x^2 - 2]] 

sage: type(enumerate_totallyreal_fields_prim(2, 10)[0][1]) 

<type 'cypari2.gen.Gen'> 

sage: enumerate_totallyreal_fields_prim(2, 10, return_pari_objects=False)[0][0].parent() 

Integer Ring 

sage: enumerate_totallyreal_fields_prim(2, 10, return_pari_objects=False)[0][1].parent() 

Univariate Polynomial Ring in x over Rational Field 

sage: enumerate_totallyreal_fields_prim(2, 10, return_seqs=True)[1][0][1][0].parent() 

Rational Field 

  

""" 

  

cdef pari_gen B_pari, d, d_poly, keepB, nf, t2val, ngt2, ng 

cdef int *f_out 

cdef int counts[4] 

cdef int i, n_int, j, skip_jp 

cdef bint found, use_t2, phc_flag, verb_int, temp_bint 

cdef Py_ssize_t k0, lenS 

cdef tr_data T 

cdef Integer dB 

cdef double db_odlyzko 

  

if not isinstance(n, Integer): 

try: 

n = Integer(n) 

except TypeError: 

raise TypeError("cannot coerce n (= %s) to an integer" % n) 

if (n < 1): 

raise ValueError("n must be at least 1.") 

  

# Initialize 

n_int = int(n) 

S = set() # set of pairs (d, f) 

lenS = 0 

  

# This is just to quiet valgrind down 

B_pari = pari(0) 

d = B_pari 

d_poly = B_pari 

keepB = B_pari 

nf = B_pari 

t2val = B_pari 

ngt2 = B_pari 

ng = B_pari 

pari_tmp1 = B_pari 

  

dB = Integer.__new__(Integer) 

dB_odlyzko = odlyzko_bound_totallyreal(n_int) 

mpz_set_d(dB.value, dB_odlyzko) 

dB = 40000*((dB+1)**n_int) 

for i from 0 <= i < 4: 

counts[i] = 0 

  

B_pari = pari(B) 

f_out = <int *>check_calloc(n_int + 1, sizeof(int)) 

f_out[n_int] = 1 

  

if keep_fields: 

if type(keep_fields) == bool: 

keepB = pari(int(math.floor(B*math.log(B)))) 

else: 

keepB = pari(keep_fields) 

else: 

keepB = pari(0) 

  

if B > keepB: 

keepB = B_pari 

  

if t_2: 

k0 = len(a) 

if isinstance(t_2, Integer): 

t2val = pari(t_2) 

else: 

t2val = pari(a[k0-2]**2-2*a[k0-3]) 

use_t2 = 1 

else: 

use_t2 = 0 

  

if phc: 

phc_flag = 1 

else: 

phc_flag = 0 

  

if just_print: 

skip_jp = 0 

jp_file = open(just_print, "w") 

else: 

skip_jp = 1 

  

# Trivial case 

if n == 1: 

sig_free(f_out) 

if return_seqs: 

return [[0,0,0,0],[[1,[-1,1]]]] 

elif return_pari_objects: 

return [[1,pari('x-1')]] 

else: 

Px = PolynomialRing(QQ, 'x') 

return [[ZZ(1), Px.gen()-1]] 

  

if verbose: 

verb_int = 1 

saveout = sys.stdout 

if type(verbose) == str: 

fsock = open(verbose, 'w') 

sys.stdout = fsock 

# Else, print to screen 

else: 

verb_int = 0 

  

T = tr_data(n_int,B,a) 

if verb_int == 2: 

T.incr(f_out,verb_int,0,phc_flag) 

else: 

T.incr(f_out,0,0,phc_flag) 

  

while f_out[n]: 

nf = new_t_POL_from_int_star(f_out, n_int+1, 0) 

if verb_int: 

print("==>", nf, "[") 

for j from 0 <= j < n-1: 

print("%s " % f_out[j]) 

print("%s]" % f_out[n - 1]) 

d_poly = nf.poldisc() 

counts[0] += 1 

if d_poly > 0 and nf.polsturm() == n: 

da = int_has_small_square_divisor(Integer(d_poly)) 

if d_poly > dB or d_poly <= B*da: 

counts[1] += 1 

if nf.polisirreducible(): 

counts[2] += 1 

zk, d = nf.nfbasis_d() 

  

if d <= keepB: 

if verb_int: 

print("has discriminant", d, end='') 

  

# Find a minimal lattice element 

counts[3] += 1 

ng = <pari_gen>((<pari_gen>(pari([nf,zk]))).polredabs()) 

  

dng = (d, ng) 

  

if skip_jp: 

# Check if K is contained in the list. 

found = dng in S 

if found and verb_int: 

print("but is not new") 

  

ngt2 = <pari_gen>(ng[n_int-1]**2-2*ng[n_int-2]) 

if not found: 

temp_bint = ngt2 >= t2val 

if ((not use_t2) or temp_bint): 

if verb_int: 

print("and is new!") 

S.add(dng) 

lenS += 1 

else: 

if ((not use_t2) or ngt2 >= t2val): 

jp_file.write(str([d, ng.Vecrev()]) + "\n") 

  

else: 

if verb_int: 

print("has discriminant", abs(d), "> B") 

else: 

if verb_int: 

print("is not irreducible") 

else: 

if verb_int: 

print("has discriminant", abs(d), "with no large enough square divisor") 

else: 

if verb_int: 

if d == 0: 

print("is not squarefree") 

else: 

print("is not totally real") 

  

if verb_int == 2: 

T.incr(f_out,verb_int,0,phc_flag) 

else: 

T.incr(f_out,0,0,phc_flag) 

  

if not skip_jp: 

if n_int == 2 and B >= 5 and ((not use_t2) or t2val <= 5): 

jp_file.write(str([2,[-1,-1,1]]) + "\n") 

if B >= 8 and B < 32: 

jp_file.write(str([2,[-2,0,1]]) + "\n") 

elif n_int == 3 and B >= 49 and ((not use_t2) or 5 >= t2val): 

jp_file.write(str([3,[1,-2,-1,1]]) + "\n") 

jp_file.close() 

sig_free(f_out) 

return 

  

# Convert S to a sorted list of pairs [d, f] 

# we sort only according to d 

# because we cannot compare t_POL objects (PARI polynomials) 

S = sorted([list(s) for s in S], key=lambda x: x[0]) 

  

# In the application of Smyth's theorem above (and easy 

# irreducibility test), we exclude finitely many possibilities 

# which we must now throw back in. 

if n_int == 2 and B >= 5 and ((not use_t2) or t2val <= 5): 

S = [[5, pari('x^2-x-1')]] + S 

lenS += 1 

if B >= 8 and B < 32: 

S.insert(1, [8, pari('x^2-2')]) 

lenS += 1 

elif n_int == 3 and B >= 49 and ((not use_t2) or 5 >= t2val): 

S = [[49, pari('x^3-x^2-2*x+1')]] + S 

lenS += 1 

# The polynomials with n = 4 define imprimitive number fields. 

  

# Now check for isomorphic fields 

lenS = weed_fields(S, lenS) 

  

# Output. 

if verb_int: 

print("=" * 80) 

print("Polynomials tested:", counts[0]) 

print("Polynomials with sssd poldisc:", counts[1]) 

print("Irreducible polynomials:", counts[2]) 

print("Polynomials with nfdisc <= B:", counts[3]) 

for i from 0 <= i < lenS: 

print(S[i]) 

if type(verbose) == str: 

fsock.close() 

sys.stdout = saveout 

  

sig_free(f_out) 

# Make sure to return elements that belong to Sage 

if return_seqs: 

return [[ZZ(counts[i]) for i in range(4)], 

[[ZZ(s[0]), map(QQ, s[1].polrecip().Vec())] for s in S]] 

elif return_pari_objects: 

return S 

else: 

Px = PolynomialRing(QQ, 'x') 

return [[ZZ(s[0]), Px(map(QQ, s[1].list()))] 

for s in S] 

  

def weed_fields(S, Py_ssize_t lenS=0): 

r""" 

Function used internally by the :func:`~enumerate_totallyreal_fields_prim` 

routine. (Weeds the fields listed by [discriminant, polynomial] 

for isomorphism classes.) Returns the size of the resulting list. 

  

EXAMPLES:: 

  

sage: ls = [[5,pari('x^2-3*x+1')],[5,pari('x^2-5')]] 

sage: sage.rings.number_field.totallyreal.weed_fields(ls) 

1 

sage: ls 

[[5, x^2 - 3*x + 1]] 

""" 

cdef Py_ssize_t i, j, n 

if lenS == 0: 

lenS = len(S) 

i = 0 

if not lenS: 

return lenS 

n = len(S[0][1])-1 

while i < lenS-1: 

j = i+1 

while j < lenS and S[i][0] == S[j][0]: 

if S[i][1].nfisisom(S[j][1]): 

# Keep the one with a smallest T_2 

T_2i = S[i][1][n-1]**2 - 2*S[i][1][n-2] 

T_2j = S[j][1][n-1]**2 - 2*S[j][1][n-2] 

if T_2i <= T_2j: 

S.pop(j) 

lenS -= 1 

else: 

t = S.pop(j) 

S.pop(i) 

S.insert(i, t) 

lenS -= 1 

else: 

j += 1 

i += 1 

  

return lenS 

  

def timestr(m): 

r""" 

Converts seconds to a human-readable time string. 

  

INPUT: 

  

- m -- integer, number of seconds 

  

OUTPUT: 

  

The time in days, hours, etc. 

  

EXAMPLES:: 

  

sage: sage.rings.number_field.totallyreal.timestr(3765) 

'1h 2m 45.0s' 

""" 

  

n = math.floor(m) 

p = m-n 

outstr = '' 

if m >= 60*60*24: 

t = n//(60*60*24) 

outstr += str(t)[:len(str(t))-2] + 'd ' 

n -= t*(60*60*24) 

if m >= 60*60: 

t = n//(60*60) 

outstr += str(t)[:len(str(t))-2] + 'h ' 

n -= t*(60*60) 

if m >= 60: 

t = n//60 

outstr += str(t)[:len(str(t))-2] + 'm ' 

n -= t*60 

n += p 

outstr += '%.1f' % n + 's' 

  

return outstr