CoCalc -- absalgtext.sage

📚 The CoCalc Library - books, templates and other resources
cocalc-examples / paulsen-abstract-algebra / absalgtext.sage
¹³²⁹²² views
License: OTHER
1
#######################################################
2
##                                                   ##
3
##  absalgtext.sage                                  ##
4
##                                                   ##
5
##  These routines are used with the sage workbooks  ##
6
##  that are provided with the book "Abstract        ##
7
##  Algebra, an Interactive Approach" by William     ##
8
##  Paulsen.                                         ##
9
##                                                   ##
10
##  Users of this file are encouraged to buy a copy  ##
11
##  of the textbook.                                 ##
12
##                                                   ##
13
#######################################################
14

15
## First of all, we need to see that the Small Groups package is loaded with Sage.  It is not automatically loaded when Sage is installed, so
16
## we probably have to load it the first time this package is run.  
17

18
try:
19
   gap.eval("SmallGroup(1,1)")
20
except:   
21
   print "Additional databases must be downloaded from the internet.  Please wait."
22
   try:
23
       install_package('database_gap')
24
   except:
25
       print "Error in installing database_gap."
26
       print "For Linux users, this will have to be installed manually.  Type:"
27
       print "    sudo sage -i database_gap"
28
       print "in the terminal window."
29
   else:    
30
       print "These new databases will be available the next time Sage is started."
31

32

33
######################################################
34
##                                                  ##
35
##  GAP routines                                    ##
36
##                                                  ##
37
## These routines are written in GAP, not Sage.     ##
38
## They are evaluated using gap.eval( ), in order   ##
39
## that we may use these routines later.            ##
40
##                                                  ##
41
######################################################
42

43
gap.eval("TempFreeGroup := FreeGroup(0);")
44
gap.eval("CurrentGroup := TempFreeGroup/[];")
45
gap.eval("CurrentPCGroup := TempFreeGroup/[];")
46

47
gap.eval("ITER_POLY_WARN := false;");  # So we can define polynomials over a ring without a warning message.
48

49
# These functions will be introduced later, but to prevent 
50
# "Variable must have a value" errors, we set them to 0 here
51
gap.eval("Mult := 0")
52
gap.eval("AddC := 0")
53
gap.eval("PermToInt := 0")
54

55
gap.eval("QuickTable := function(L) " +
56
            "local i, j, k, LL, G, n, R, M, p; " +
57
            "LL := List(L); " +
58
            "G := Flat(LL); " + # find the base group
59
            "n := Size(LL); " +
60
            "M := NullMat(n,n); " + # matrix of zeros, mutable
61
            "for i in [1..n] do " +
62
               "for j in [1..n] do " +
63
                  "if IsList(LL[i]) then " +
64
                     "p := Mult(G,LL[i],LL[j]); " + # Multiply two cosets
65
                  "else " +
66
                     "p := LL[i]*LL[j]; " +
67
                  "fi; " +
68
                  "for k in [1..n] do " +
69
                     "if IsList(p) then " +
70
                        "if( Set(p) = Set(LL[k]) ) then " +
71
                           "M[i][j] := k; " +
72
                           "break; " +
73
                        "fi; " +
74
                     "elif (p = LL[k]) then " +
75
                        "M[i][j] := k; " +
76
                        "break; " +
77
                     "fi; " +
78
                  "od; " +
79
               "od; " +
80
             "od; " +
81
             "return M; " +
82
          "end;")
83

84
gap.eval("QuickAdd := function(L) " +
85
            "local i, j, k, LL, G, n, R, M, p; " +
86
            "LL := List(L); " +
87
            "G := Flat(LL); " + # find the base group
88
            "n := Size(LL); " +
89
            "M := NullMat(n,n); " + # matrix of zeros, mutable
90
            "for i in [1..n] do " +
91
               "for j in [1..n] do " +
92
                  "if IsList(LL[i]) then " +
93
                     "p := AddC(G,LL[i],LL[j]); " + # Add two cosets
94
                  "else " +
95
                     "p := LL[i] + LL[j]; " +
96
                  "fi; " +
97
                  "for k in [1..n] do " +
98
                     "if IsList(p) then " +
99
                        "if( Set(p) = Set(LL[k]) ) then " +
100
                           "M[i][j] := k; " +
101
                           "break; " +
102
                        "fi; " +
103
                     "elif (p = LL[k]) then " +
104
                        "M[i][j] := k; " +
105
                        "break; " +
106
                     "fi; " +
107
                  "od; " +
108
               "od; " +
109
            "od; " +
110
            "return M; " +
111
         "end;")
112

113
gap.eval("ListGroup := function(g) " +
114
            "local L, Ord, n, iList, pointer, Cont, GList, ele; " +
115
            "if (IsGroup(g) = false) then " +
116
               'Print("Not a group\n"); ' +
117
               "return fail; " +
118
            "fi; " +
119
            "if (Size(g) = infinity) then " +
120
               'Print("Infinite group\n"); ' +
121
               "return fail; " +
122
            "fi; " +
123
            "L := GeneratorsOfGroup(g); " +
124
            "Ord := List(L, x->Order(x)); " +
125
             # do a arbitrarily nested for loop
126
            "n := Size(L); " +
127
            "GList := [ ]; " +
128
            "iList := List([1..n],x->0); " + # Vector of zeros, of length n 
129
            "Cont := true; " +
130
            "while Cont do " +
131
               "ele := Product([1..n],x->L[x]^iList[x]); " +
132
               "if not (ele in GList) then " +
133
                  "Add(GList,ele); " +
134
               "fi; " +
135
                # increment the iList
136
               "pointer := 1; " +
137
               "iList[pointer] := iList[pointer] + 1; " +
138
               "while (iList[pointer] >= Ord[pointer]) do " +
139
                  "iList[pointer] := 0; " +
140
                  "pointer := pointer + 1; " +
141
                  "if (pointer > n) then " +
142
                     "Cont := false; " +
143
                     "pointer := 1; " +
144
                  "fi; " +
145
                  "iList[pointer] := iList[pointer] + 1; " +
146
               "od; " +
147
            "od; " +
148
            "return GList; " +
149
         "end;")
150

151
gap.eval("ListRing := function(g) " +
152
            "local L, Ord, n, iList, pointer, Cont, GList, ele; " +
153
            "if (IsRing(g) = false) then " +
154
               'Print("Not a ring\n"); ' +
155
               "return fail; " +
156
            "fi; " +
157
            "if (Size(g) = infinity) then " +
158
               'Print("Infinite ring\n"); ' +
159
               "return fail; " +
160
            "fi; " +
161
            "L := GeneratorsOfRing(g); " +
162
            "Ord := List(L, x->Size(Ring(x))); " +
163
             # do a arbitrarily nested for loop 
164
            "n := Size(L); " +
165
            "GList := [ ]; " +
166
            "iList := List([1..n],x->0); " + # Vector of zeros, of length n
167
            "Cont := true; " +
168
            "while Cont do " +
169
               "ele := Sum([1..n],x->iList[x]*L[x]); " +
170
               "if not (ele in GList) then " +
171
                  "Add(GList,ele); " +
172
               "fi; " +
173
                # increment the iList
174
               "pointer := 1; " +
175
               "iList[pointer] := iList[pointer] + 1; " +
176
               "while (iList[pointer] >= Ord[pointer]) do " +
177
                  "iList[pointer] := 0; " +
178
                  "pointer := pointer + 1; " +
179
                  "if (pointer > n) then " +
180
                     "Cont := false; " +
181
                     "pointer := 1; " +
182
                  "fi; " +
183
                  "iList[pointer] := iList[pointer] + 1; " +
184
               "od; " +
185
            "od; " +
186
            "return GList; " +
187
         "end;")
188

189
gap.eval("CheckRing := function(R) " +
190
            "local i,j,k,T,L; " +
191
            "if Size(R) = infinity then " +
192
               "L := GeneratorsOfRing(R); " +
193
            "else " +
194
               "L := List(R); " +
195
            "fi; " +
196
            "T := true; " +
197
            "for i in L do " +
198
               "for j in L do " +
199
                  "for k in L do " +
200
                     "if ( (i+j)*k <> (i*k) + (j*k)) then " +
201
                        "T := false; " +
202
                     "fi; " +
203
                  "od; " +
204
               "od; " +
205
            "od; " +
206
            "if (not T) then " +
207
               'return "Ring is not left distributive."; ' +
208
            "else " +
209
               "for i in L do " +
210
                  "for j in L do " +
211
                     "for k in L do " +
212
                        "if ( i*(j+k) <> (i*j) + (i*k)) then " +
213
                           "T := false; " +
214
                        "fi; " +
215
                     "od; " +
216
                  "od; " +
217
               "od; " +
218
               "if (not T) then " +
219
                  'return "Ring is not right distributive."; ' +
220
               "else " +
221
                  "for i in L do " +
222
                     "for j in L do " +
223
                        "for k in L do " +
224
                           "if ( (i*j)*k <> i*(j*k) ) then " +
225
                              "T := false; " +
226
                           "fi; " +
227
                        "od; " +
228
                     "od; " +
229
                  "od; " +
230
                  "if T then " +
231
                      'return "This is a ring."; ' +
232
                  "else " +
233
                      'return "Associative law does not hold."; ' +
234
                  "fi; " +
235
               "fi; " +
236
            "fi; " +
237
         "end;")    
238

239
gap.eval("Mult := function( G, H, Y ) " +
240
            "local GG, XX, YY, i, j, FlatG, SizeG, p, base, k, L; " +
241
             # Find the base group 
242
            "GG := List(G); " +
243
            "FlatG := Flat(GG); " +
244
             # Determine the operation
245
            "SizeG := Size(FlatG); " +
246
             # base = 0 means G is a group, base > 0 means multiply mod base,
247
             # base < 0 means add mod base.
248
            "base := 0; " +
249
            "if (Number(FlatG, IsInt) = SizeG) then " +
250
               "base := Maximum(FlatG) + 1; " +
251
               "if (0 in FlatG) then " +
252
                  "base := -base; " +
253
               "fi; " +
254
            "fi; " +
255
            "L := []; " +
256
            "if (IsList(H) or IsGroup(H) or IsRing(H)) then " +
257
               "XX := List(H); " +
258
            "else " +
259
               "XX := [H]; " +
260
            "fi; " +
261
            "if (IsList(Y) or IsGroup(Y) or IsRing(Y)) then " +
262
               "YY := List(Y); " +
263
            "else " +
264
               "YY := [Y]; " +
265
            "fi; " +
266
             # assume first that XX and YY are simple lists. " +
267
            "for i in [1..Size(XX)] do " +
268
               "for j in [1..Size(YY)] do " +
269
                  "if(IsList(XX[i]) or IsList(YY[j]) ) then " +
270
                     "Add(L, Mult(FlatG, XX[i], YY[j])); " + # recursively defined for sets of sets
271
                  "else " +
272
                      # find the product p " +
273
                     "if (base > 0) then " +
274
                        "p := (XX[i] * YY[j]) mod base; " +
275
                     "elif (base < 0) then " +
276
                        "p := (XX[i] + YY[j]) mod (- base); " +
277
                     "else " +
278
                        "p := XX[i]*YY[j]; " +
279
                     "fi; " +
280
                     "for k in [1..Size(GG)] do " +
281
                        "if (p = GG[k]) then " +
282
                           "Add(L, GG[k]); " +
283
                           "break; " +
284
                        "fi; " +
285
                     "od; " +
286
                  "fi; " +
287
               "od; " +
288
            "od; " +
289
            "Sort(L); " +
290
             # remove duplicates from the list
291
            "for i in Reversed([2..Size(L)]) do " +
292
               "for j in Reversed([1..i-1]) do " +
293
                  "if (L[i] = L[j]) then " +
294
                     "Remove(L,i); " +
295
                     "break; " +
296
                  "fi; " +
297
               "od; " +
298
            "od; " +
299
            "return L; " +
300
         "end;")
301

302
gap.eval("LSide := 0;")
303
gap.eval("RSide := 0;")
304

305
gap.eval("LeftSide := function(x) " +
306
            "LSide := x; " +
307
         "end;")   
308
 
309
gap.eval("RightSide := function(x) " +
310
            "RSide := x; " +
311
         "end;")  
312

313
gap.eval("Prod := function(G) " +
314
            "local i, j, k, p, L; " +
315
            "L := []; " +
316
            "for i in [1..Size(LSide)] do " +
317
               "for j in [1..Size(RSide)] do " +
318
                   # find the product p
319
                  "p := LSide[i]*RSide[j]; " +
320
                  "if (G = []) then " +
321
                     "Add(L, p); " +
322
                  "else " +
323
                     "for k in [1..Size(G)] do " +
324
                        "if (p = G[k]) then " +
325
                           "Add(L, G[k]); " +
326
                           "break; " +
327
                        "fi; " +
328
                     "od; " +
329
                  "fi; " +
330
               "od; " +
331
            "od; " +
332
            "Sort(L); " +
333
             # remove duplicates from the list (which will be together if sorted)
334
            "for i in Reversed([2..Size(L)]) do " +
335
               "if (L[i] = L[i-1]) then " +
336
                  "Remove(L,i); " +
337
               "fi; " +
338
            "od; " +
339
            "return L; " +
340
         "end;")
341

342
gap.eval("AddCoset := function(R) " +
343
            "local i, j, k, p, L; " +
344
            "L := []; " +
345
            "for i in [1..Size(LSide)] do " +
346
               "for j in [1..Size(R)] do " +
347
                   # find the sum p " +
348
                  "p := LSide[i] + R[j]; " +
349
                  "Add(L, p); " +
350
               "od; " +
351
            "od; " +
352
            "Sort(L); " +
353
             # remove duplicates from the list (which will be together if sorted)
354
            "for i in Reversed([2..Size(L)]) do " +
355
               "if (L[i] = L[i-1]) then " +
356
                  "Remove(L,i); " +
357
               "fi; " +
358
            "od; " +
359
            "return L; " +
360
         "end;")
361

362
gap.eval("OrigGroup := 0;")
363

364
gap.eval("OriginalGroup := function(x) " +
365
            "OrigGroup := x; " +
366
         "end;")  
367

368
gap.eval("Conform := function(sub) " +
369
            "local i, k, L; " +
370
            "L := []; " +
371
            "for i in [1..Size(sub)] do " +
372
               "for k in [1..Size(OrigGroup)] do " +
373
                   "if (sub[i] = OrigGroup[k]) then " +
374
                      "Add(L, OrigGroup[k]); " +
375
                      "break; " +
376
                   "fi; " +
377
               "od; " +
378
            "od; " +
379
            "return L; " +
380
         "end;")   
381

382
gap.eval("LftCoset := function(H) " +
383
            "local GG, HH, i, j, k, q, base, p, FlatQ, SizeG, t; " +
384
            "GG := OrigGroup; " +
385
            "SizeG := Size(GG); " +
386
            "HH := Conform(H); " +
387
            "Sort(HH); " +
388
            "q := [HH]; " +
389
            "Print(q); " +
390
            "FlatQ := Flat(HH); " +
391
            "for i in [1..SizeG] do " +
392
               "if not(GG[i] in FlatQ) then " +
393
                   "t := []; " +
394
                   "for j in [1..Size(HH)] do " +
395
                      "p := GG[i]*HH[j]; " +
396
                      "for k in [1..SizeG] do " +
397
                         "if p = GG[k] then " +
398
                            "Add(t, GG[k]); " +
399
                            "Add(FlatQ, GG[k]); " +
400
                         "fi; " +
401
                      "od; " +
402
                   "od; " +
403
                   "Sort(t); " +
404
                   "Add(q, t); " +
405
                   "Print(q); " +
406
               "fi; " +
407
            "od; " +
408
            "return q; " +
409
         "end;")   
410

411
gap.eval("RtCoset := function(H) " +
412
            "local GG, HH, i, j, k, q, p, FlatQ, SizeG, t; " +
413
            "GG := OrigGroup; " +
414
            "SizeG := Size(GG); " +
415
            "HH := Conform(H); " +
416
            "Sort(HH); " +
417
            "q := [HH]; " +
418
            "FlatQ := Flat(HH); " +
419
            "for i in [1..SizeG] do " +
420
               "if not(GG[i] in FlatQ) then " +
421
                   "t := []; " +
422
                   "for j in [1..Size(HH)] do " +
423
                      "p := HH[j]*GG[i]; " +
424
                      "for k in [1..SizeG] do " +
425
                         "if p = GG[k] then " +
426
                            "Add(t, GG[k]); " +
427
                            "Add(FlatQ, GG[k]); " +
428
                         "fi; " +
429
                      "od; " +
430
                   "od; " +
431
                   "Sort(t); " +
432
                   "Add(q, t); " +
433
               "fi; " +
434
            "od; " +
435
            "return q; " +
436
         "end;") 
437

438
gap.eval("RingCoset := function(H) " +
439
            "local GG, HH, i, j, k, q, p, FlatQ, SizeG, t; " +
440
            "GG := OrigGroup; " +
441
            "SizeG := Size(GG); " +
442
            "HH := H; " +
443
            "Sort(HH); " +
444
            "q := [HH]; " +
445
            "FlatQ := Flat(HH); " +
446
            "for i in [1..SizeG] do " +
447
               "if not(GG[i] in FlatQ) then " +
448
                   "t := []; " +
449
                   "for j in [1..Size(HH)] do " +
450
                      "p := HH[j] + GG[i]; " +
451
                      "Add(t, p); " +
452
                      "Add(FlatQ, p); " +
453
                   "od; " +
454
                   "Sort(t); " +
455
                   "Add(q, t); " +
456
               "fi; " +
457
            "od; " +
458
            "return q; " +
459
         "end;")     
460

461
gap.eval("NormClosure := function(S) " +
462
            "local GG, SS, NN, LL; " +
463
            "SS := Group(S); " + # make both S and G into groups.
464
            "GG := Group(OrigGroup); " +
465
            "NN := NormalClosure(GG, SS); " +
466
            "LL := List(NN); " +
467
            "return Conform(LL); " +
468
         "end;")   
469

470
gap.eval("IdealClosure := function(S) " +
471
            "local RR, SS, II; " +
472
            "SS := List(S); " +         # make S into a list.
473
            "RR := Ring(OrigGroup); " + # make GG into a ring
474
            "II := Ideal(RR, SS); " +
475
            "return List(II); " +
476
         "end;")   
477

478
gap.eval("ConjClasses := function(G) " +
479
             "local GG, SS, LL; " +
480
             "OriginalGroup(G); " +
481
             "GG := Group(G); " +
482
             "SS := ConjugacyClasses(GG); " +
483
             "LL := List(SS, x -> Conform(List(x))); " +
484
             "return LL; " +
485
         "end;")
486

487
gap.eval("CommSubgroup := function(K) " +
488
             "local HH, KK, MM; " +
489
             "HH := Group(OrigGroup); " +
490
             "KK := Group(K); " +
491
             "MM := CommutatorSubgroup(HH, KK); " +
492
             "if IsSubset(HH, MM) then " +        # we can use OrigGroup to Conform the list
493
                "return Conform(List(MM)); " +
494
             "else " +
495
                "return List(MM); " +
496
             "fi; " +
497
         "end;")   
498
 
499
gap.eval("AddC := function( G, H, Y ) " +
500
            "local GG, XX, YY, i, j, FlatG, SizeG, p, base, k, L; " +
501
             # Find the base ring.
502
            "GG := List(G); " +
503
            "FlatG := Flat(GG); " +
504
             # Find the operation 
505
            "SizeG := Size(FlatG); " +
506
             # base = 0 means G is a ring, base > 0 means add mod base.
507
            "base := 0; " +
508
            "if (Number(FlatG, IsInt) = SizeG) then " +
509
               "base := Maximum(FlatG) + 1; " +
510
            "fi; " +
511
            "L := []; " +
512
            "if (IsList(H) or IsRing(H)) then " +
513
               "XX := List(H); " +
514
            "else " +
515
               "XX := [H]; " + # list with 1 element
516
            "fi; " +
517
            "if (IsList(Y) or IsRing(Y)) then " +
518
               "YY := List(Y); " +
519
            "else " +
520
               "YY := [Y]; " + # list with 1 element
521
            "fi; " +
522
             # assume first that XX and YY are simple lists.
523
            "for i in [1..Size(XX)] do " +
524
               "for j in [1..Size(YY)] do " +
525
                  "if(IsList(XX[i]) or IsList(YY[j]) ) then " +
526
                     "Add(L, AddC(FlatG, XX[i], YY[j])); " + # recursively defined for sets of sets
527
                  "else " +
528
                      # find the sum p
529
                     "if (base > 0) then " +
530
                        "p := (XX[i] + YY[j]) mod base; " +
531
                     "else " +
532
                        "p := XX[i]+YY[j]; " +
533
                     "fi; " +
534
                     "for k in [1..Size(GG)] do " +
535
                        "if (p = GG[k]) then " +
536
                           "Add(L, GG[k]); " +
537
                           "break; " +
538
                        "fi; " +
539
                     "od; " +
540
                  "fi; " +
541
               "od; " +
542
            "od; " +
543
            "Sort(L); " +
544
             # remove duplicates from the list (which will be together if sorted)
545
            "for i in Reversed([2..Size(L)]) do " +
546
               "if (L[i] = L[i-1]) then " +
547
                  "Remove(L,i); " +
548
               "fi; " +
549
            "od; " +
550
            "return L; " +
551
         "end;")
552
  
553
gap.eval("ListCurrentGroup := function(x) " +
554
            "return List(CurrentGroup); " +
555
         "end;")
556

557
gap.eval("ListCurrentPCGroup := function(x) " +
558
            "return List(CurrentPCGroup); " +
559
         "end;")
560

561

562

563

564

565

566

567

568

569
## Set up type of elements to test whether an element is of a certain type 
570

571
GapInstance = type(gap("(1,2)")) # used for isinstance to detect Gap elements
572
IntModType = type(Integers(2)(1)) # used for isinstance to detect multiplication mod n
573
IntType = type(1)  
574
GFType = type(GF(2))
575
QuatType = type(QuaternionAlgebra(SR, -1, -1).gens()[0])
576

577
## Set flags to their default position
578

579
reset("e")
580
E = e # since e is often used for the identity element of a group, we need a way to access 2.718281828...
581

582
LastInit = "None"
583
UnivFieldVar = ""
584
LastFieldVar = ""
585
#_QuickMult_ = False
586

587
## Allow alias of several Sage functions to be consistent with the CamelCase convention of Mathematica
588
    
589
def EulerPhi(x):
590
    return euler_phi(x)
591
    
592
def PowerMod(x, r, n): 
593
    return pow(x, r, n)    
594
    
595
def PartitionsP(x):
596
    return Partitions(x).cardinality()    
597

598
def ConwayPolynomial(p,d):
599
    return conway_polynomial(p, d)    
600

601
def List(x):
602
    if str(parent(x)) == 'Gap':
603
        return list(x.List())
604
    else:
605
        return list(x)
606
        
607
def PolynomialQuotient(dividend, divisor):
608
    return (dividend._maxima_().divide(divisor).sage())[0]   
609

610
def PolynomialRemainder(dividend, divisor):
611
    return (dividend._maxima_().divide(divisor).sage())[1]
612
   
613
maxima_calculus('algebraic: true;') # allows Together to also rationalize the denominator.
614
        
615
def Together(x):
616
    if isinstance(x, QuatType):
617
        return x[0].simplify_rational() + x[1].simplify_rational() * i + x[2].simplify_rational() * j + x[3].simplify_rational() * k
618
    return x.simplify_rational()
619
    
620
def AddMod(x):
621
    return ZGroup(x)
622
    
623
def MultMod(x):
624
    return ZRing(x)  
625

626
def TrigReduce(x):
627
    y = expand(x)
628
    try:
629
        y = y.trig_reduce()
630
    except:
631
        return y    
632
    y = y.subs(cos(4*pi/9) == cos(pi/9) - cos(2*pi/9))
633
    y = y.subs(cos(5*pi/9) == cos(2*pi/9) - cos(pi/9))
634
    y = y.subs(cos(3*pi/7) == cos(2*pi/7) - cos(pi/7) + 1/2) 
635
    y = y.subs(cos(4*pi/7) == cos(pi/7) - cos(2*pi/7) - 1/2)   
636
    y = expand(y)
637
    return y #y.trig_reduce()  
638

639
########################################################
640
##                                                    ##
641
##  Graphical Routines                                ##
642
##                                                    ##
643
##  These routines take advantage of Sage's graphics  ##
644
##  capabilities.  They include Terry the triangle,   ##
645
##  swapping books, the octahedron, and the Pyriminx  ##
646
##  puzzle.                                           ##
647
##                                                    ##
648
########################################################
649

650
ColorTable = {0:'#000000', 1:'#FFFF00', 2:'#FF00FF', 3:'#00FFFF', 4:'#FF0000', 5:'#00FF00', 6:'#9966CC', 7:'#FF9900', 8:'#B3CCB3',
651
9:'#CC9980', 10:'#8099CC', 11:'#CC3380', 12:'#CC9933', 13:'#999933', 14:'#FF6666', 15:'#33CCB3', 16:'#FF99FF', 17:'#B3FF00',
652
18:'#0000FF', 19:'#996666', 20:'#0099FF', 21:'#FFCCB3', 22:'#6633FF', 23:'#FFCC66', 24:'#339966', 25:'#80CCCC', 26:'#66CC80',
653
27:'#00FF99', 28:'#CC8033', 29:'#3333E5'}
654

655
       
656
TerryColor = ['#FF0000', '#00FF00', '#00FFFF']
657

658
def ShowTerry():
659
    global TerryColor
660
    global Stay
661
    global RotLft
662
    global RotRt
663
    global Spin
664
    global FlipLft
665
    global FlipRt
666
    global IdentityElement
667
    if(IdentityElement != "Stay"):
668
        # If Terry's group is not already defined, set up the variable names.
669
        Stay = var('Stay')
670
        RotLft = var('RotLft')
671
        RotRt = var('RotRt')
672
        Spin = var('Spin')
673
        FlipLft = var('FlipLft')
674
        FlipRt = var('FlipRt')
675
    sqrt3 = sqrt(3).numerical_approx()
676
    Graff = Graphics()
677
    Graff += polygon([(0,0),(0,2),(sqrt3,-1)], color = TerryColor[0])
678
    Graff += polygon([(0,0),(0,2),(-sqrt3,-1)], color = TerryColor[1]) 
679
    Graff += polygon([(0,0),(sqrt3,-1),(-sqrt3,-1)], color = TerryColor[2])
680
    return Graff.show(axes = false, aspect_ratio = 1, xmin = -2, xmax = 2, ymin = -2, ymax = 2)
681
       
682
def Terry(*args):
683
    global TerryColor
684
    sqrt3 = sqrt(3).numerical_approx()
685
    v = [] 
686
    if(IdentityElement == 'Stay'):
687
        ## If InitTerry has been defined, then like Mathematica, we want Terry to take the shortcut.
688
        Prod = Stay
689
        for arg in args:
690
            Prod = Prod * arg
691
        newList = []
692
        newList.append(Prod)
693
    else:
694
        newList = []
695
        for arg in args:
696
            newList.append(arg)    
697
    for arg in newList:
698
        if(arg == Stay): 
699
            Graff = Graphics()
700
            Graff += polygon([(0,0),(0,2),(sqrt3,-1)], color = TerryColor[0])
701
            Graff += polygon([(0,0),(0,2),(-sqrt3,-1)], color = TerryColor[1]) 
702
            Graff += polygon([(0,0),(sqrt3,-1),(-sqrt3,-1)], color = TerryColor[2])
703
            for i in range(5):    
704
                v.append(Graff)
705
        elif(arg == RotRt):
706
            for i in range(6): 
707
                th = 2*pi*i/3/5.0  # theta ranges from 0 to 2pi/3, inclusive   
708
                Graff = Graphics()
709
                Graff += polygon([(0,0),(2*sin(th),2*cos(th)), (2*sin(th + 2*pi/3),2*cos(th + 2*pi/3))], color = TerryColor[0])
710
                Graff += polygon([(0,0), (2*sin(th + 2*pi/3),2*cos(th + 2*pi/3)), (2*sin(th - 2*pi/3),2*cos(th - 2*pi/3))], color = TerryColor[2])
711
                Graff += polygon([(0,0), (2*sin(th - 2*pi/3),2*cos(th - 2*pi/3)), (2*sin(th),2*cos(th))], color = TerryColor[1])
712
                v.append(Graff)
713
            Temp = TerryColor[0]
714
            TerryColor[0] = TerryColor[1]
715
            TerryColor[1] = TerryColor[2]
716
            TerryColor[2] = Temp        
717
        elif(arg == RotLft):
718
            for i in range(6): 
719
                th = 2*pi*i/3/5.0  # theta ranges from 0 to 2pi/3, inclusive   
720
                Graff = Graphics()
721
                Graff += polygon([(0,0),(-2*sin(th),2*cos(th)), (2*sin(2*pi/3 - th),2*cos(2*pi/3 - th))], color = TerryColor[0])
722
                Graff += polygon([(0,0), (2*sin(2*pi/3-th),2*cos(2*pi/3-th)), (-2*sin(th + 2*pi/3),2*cos(th + 2*pi/3))], color = TerryColor[2])
723
                Graff += polygon([(0,0), (-2*sin(th + 2*pi/3),2*cos(th + 2*pi/3)), (-2*sin(th),2*cos(th))], color = TerryColor[1])
724
                v.append(Graff)
725
            Temp = TerryColor[0]
726
            TerryColor[0] = TerryColor[2]
727
            TerryColor[2] = TerryColor[1]
728
            TerryColor[1] = Temp  
729
        elif(arg == Spin):
730
            for i in range(6): 
731
                th = pi*i/5.0  # theta ranges from 0 to pi, inclusive   
732
                Graff = Graphics()
733
                Graff += polygon([(0,0),(0,2),(sqrt3*cos(th),sin(th)/2-1)], color = TerryColor[0])
734
                Graff += polygon([(0,0), (sqrt3*cos(th),sin(th)/2-1), (-sqrt3*cos(th),-sin(th)/2-1)], color = TerryColor[2])
735
                Graff += polygon([(0,0), (-sqrt3*cos(th),-sin(th)/2-1), (0,2)], color = TerryColor[1])
736
                v.append(Graff)
737
            Temp = TerryColor[0]
738
            TerryColor[0] = TerryColor[1]
739
            TerryColor[1] = Temp  
740
        elif(arg == FlipRt):
741
            for i in range(6): 
742
                th = pi*i/5.0  # theta ranges from 0 to pi, inclusive   
743
                Graff = Graphics()
744
                Graff += polygon([(0,0),(-sqrt3,-1),(sqrt3/2 - sqrt3*cos(th)/2 - sqrt3*sin(th)/4, 1/2 + 3*cos(th)/2 - sin(th)/4)], color = TerryColor[1])
745
                Graff += polygon([(0,0), (sqrt3/2 - sqrt3*cos(th)/2 - sqrt3*sin(th)/4, 1/2 + 3*cos(th)/2 - sin(th)/4),
746
                                         (sqrt3/2 + sqrt3*cos(th)/2 + sqrt3*sin(th)/4, 1/2 - 3*cos(th)/2 + sin(th)/4)], color = TerryColor[0])
747
                Graff += polygon([(0,0), (sqrt3/2 + sqrt3*cos(th)/2 + sqrt3*sin(th)/4, 1/2 - 3*cos(th)/2 + sin(th)/4), 
748
                                         (-sqrt3,-1)], color = TerryColor[2])
749
                v.append(Graff)
750
            Temp = TerryColor[1]
751
            TerryColor[1] = TerryColor[2]
752
            TerryColor[2] = Temp 
753
        elif(arg == FlipLft):
754
            for i in range(6): 
755
                th = pi*i/5.0  # theta ranges from 0 to pi, inclusive   
756
                Graff = Graphics()
757
                Graff += polygon([(0,0),(sqrt3,-1),(-sqrt3/2 - sqrt3*cos(th)/2 + sqrt3*sin(th)/4, 1/2 - 3*cos(th)/2 - sin(th)/4)], color = TerryColor[2])
758
                Graff += polygon([(0,0), (-sqrt3/2 - sqrt3*cos(th)/2 + sqrt3*sin(th)/4, 1/2 - 3*cos(th)/2 - sin(th)/4),
759
                                         (-sqrt3/2 + sqrt3*cos(th)/2 - sqrt3*sin(th)/4, 1/2 + 3*cos(th)/2 + sin(th)/4)], color = TerryColor[1])
760
                Graff += polygon([(0,0), (-sqrt3/2 + sqrt3*cos(th)/2 - sqrt3*sin(th)/4, 1/2 + 3*cos(th)/2 + sin(th)/4), 
761
                                         (sqrt3,-1)], color = TerryColor[0])
762
                v.append(Graff)
763
            Temp = TerryColor[0]
764
            TerryColor[0] = TerryColor[2]
765
            TerryColor[2] = Temp                     
766
        else:
767
            print("Invalid rotation")
768
    aniTerry = animate(v, xmin = -2, xmax = 2, ymin = -2, ymax = 2, aspect_ratio = 1, axes = false)
769
    return aniTerry.show(delay = 50, iterations = 1)       
770
       
771
BookColor = []  # Initialize so we can use it later
772
 
773
def ShowBooks():
774
    global BookColor
775
    Graff = Graphics()
776
    for i in range(len(BookColor)):
777
        Graff += polygon([(2*i,-i),(2*i+.16,-.34-i),(2*i+.382,-.618-i),(2*i+.66,-.84-i),(2*i+1,-1-i),(2*i+1.421,-1.079-i),(2*i+2,-1-i), (2*i+2,-.875-i),(2*i+2.25,-.75-i),(2*i+5.333,2.333-i),(2*i+5.333,9.333-i),(2*i+5.167,9.417-i),(2*i+3.5,7.75-i),(2*i+3.5,10.25-i), (2*i+3.333,10.333-i),(2*i+.25,7.25-i),(2*i+.125,7-i),(2*i,7-i)], color = BookColor[i])
778
        Graff += polygon([(2*i+4.5,8.75-i),(2*i+2,6.25-i),(2*i+2,6-i),(2*i+1.421,5.921-i),(2*i+1,6-i), (2*i+.66,6.16-i),(2*i+.382,6.382-i),(2*i+.271,6.521-i),(2*i+3.167,9.417-i)], rgbcolor=(1,1,1))
779
        Graff += line([(2*i,-i),(2*i+.16,-.34-i),(2*i+.382,-.618-i),(2*i+.66,-.84-i),(2*i+1,-1-i),(2*i+1.421,-1.079-i),(2*i+2,-1-i),(2*i+2,6-i), (2*i+1.421,5.921-i),(2*i+1,6-i),(2*i+.66,6.16-i),(2*i+.382,6.382-i),(2*i+.16,6.66-i),(2*i,7-i),(2*i,-i)], rgbcolor=(0,0,0))
780
        Graff += line([(2*i+3.167,9.417-i),(2*i+4.5,8.75-i)], rgbcolor=(0,0,0))
781
        Graff += line([(2*i+5.333,9.333-i),(2*i+2.25,6.25-i),(2*i+2,6.125-i)], rgbcolor=(0,0,0))
782
        Graff += line([(2*i,7-i),(2*i+.25,7-i),(2*i+3.5,10.25-i)], rgbcolor=(0,0,0))
783
        Graff += line([(2*i+2,6-i),(2*i+2,6.25-i),(2*i+5.167,9.417-i),(2*i+5.333,9.333-i)], rgbcolor=(0,0,0))
784
        Graff += line([(2*i+.25,6.547-i),(2*i+.25,7-i)], rgbcolor=(0,0,0))
785
        Graff += line([(2*i+.125,7-i),(2*i+.25,7.25-i),(2*i+3.333,10.333-i),(2*i+3.5,10.25-i),(2*i+3.5,9.25-i)], rgbcolor=(0,0,0))
786
    i = len(BookColor)-1    
787
    Graff += line([(2*i+2.25,6.25-i),(2*i+2.25,-.75-i),(2*i+5.333,2.333-i),(2*i+5.333,9.333-i)], rgbcolor=(0,0,0))   
788
    return Graff.show(axes = false)
789
   
790
def InitBooks(n):
791
    global BookColor
792
    global Stay
793
    global Left
794
    global Right
795
    global First
796
    global Last
797
    global Rev
798
    Stay = var('Stay')
799
    Left = var('Left')
800
    Right = var('Right')
801
    First = var('First')
802
    Last = var('Last')
803
    Rev = var('Rev')
804
    BookColor = []
805
    for i in range(n):
806
        BookColor.append(ColorTable[i+4])
807
    return ShowBooks()
808
            
809
def MoveBooks(*args):
810
    global BookColor
811
    n = len(BookColor)
812
    if n > 1:
813
        for arg in args:
814
            if (arg == Stay):
815
                ShowBooks()
816
            elif (arg == Left):
817
                BookColor.append(BookColor.pop(0))
818
                ShowBooks()
819
            elif (arg == Right):
820
                BookColor.insert(0,BookColor.pop())
821
                ShowBooks()
822
            elif (arg == First):
823
                Temp = BookColor[0]
824
                BookColor[0] = BookColor[1]
825
                BookColor[1] = Temp
826
                ShowBooks()
827
            elif (arg == Last):
828
                Temp = BookColor[n-1]
829
                BookColor[n-1] = BookColor[n-2]
830
                BookColor[n-2] = Temp
831
                ShowBooks()
832
            elif (arg == Rev):
833
                BookColor.reverse()
834
                ShowBooks()
835
            else:
836
                print("Invalid step")
837
    return            
838
                
839
OctColor = [(1,0,0), (0,1,0), (0,0,1), (1,1,0), (1,0,1), (0,1,1), (1,.6,0), (0,0,0)]
840

841
def OctProj(x,y,z):  # projects a 3D point into 2D, using the projection angle that was used in TeX
842
    return (-0.3*x + y, -0.375*x -0.15*y + z) 
843
    
844
def ShowOctahedron():
845
    global OctColor
846
    Graff = Graphics()
847
    Graff += polygon([OctProj(0,0,1), OctProj(0,1,0), OctProj(1,0,0)], rgbcolor = OctColor[0])
848
    Graff += polygon([OctProj(0,0,1), OctProj(0,-1,0), OctProj(1,0,0)], rgbcolor = OctColor[2])
849
    Graff += polygon([OctProj(0,0,-1), OctProj(0,1,0), OctProj(1,0,0)], rgbcolor = OctColor[4])
850
    Graff += polygon([OctProj(0,0,-1), OctProj(0,-1,0), OctProj(1,0,0)], rgbcolor = OctColor[6])
851
    return Graff.show(axes = false, aspect_ratio = 1)
852

853
def ShowOctahedronWithQuaternions():
854
    global OctColor
855
    Graff = Graphics()
856
    Graff += polygon([OctProj(0,0,1), OctProj(0,1,0), OctProj(1,0,0)], rgbcolor = OctColor[0])
857
    Graff += polygon([OctProj(0,0,1), OctProj(0,-1,0), OctProj(1,0,0)], rgbcolor = OctColor[2])
858
    Graff += polygon([OctProj(0,0,-1), OctProj(0,1,0), OctProj(1,0,0)], rgbcolor = OctColor[4])
859
    Graff += polygon([OctProj(0,0,-1), OctProj(0,-1,0), OctProj(1,0,0)], rgbcolor = OctColor[6])
860
    Graff += text("i", OctProj(0,1.1,0), fontsize=20, rgbcolor=(0,0,0))
861
    Graff += text("-i", OctProj(0,-1.1,0), fontsize=20, rgbcolor=(0,0,0))
862
    Graff += text("j", OctProj(0.7,0,0), fontsize=20, rgbcolor=(0,0,0))
863
    Graff += text("k", OctProj(0,0,1.1), fontsize=20, rgbcolor=(0,0,0))
864
    Graff += text("-k", OctProj(0,0,-1.1), fontsize=20, rgbcolor=(0,0,0))
865
    Graff += text("<--------  -j", OctProj(-1,0,0), horizontal_alignment='left', fontsize=20, rgbcolor=(0,0,0))
866
    return Graff.show(axes = false, aspect_ratio = 1)
867
    
868
def InitOctahedron():
869
    global OctColor
870
    global a
871
    global b
872
    global c
873
    a = var('a')
874
    b = var('b')
875
    c = var('c')    
876
    OctColor = [(1,0,0), (0,1,0), (0,0,1), (1,1,0), (1,0,1), (0,1,1), (1,.6,0), (0,0,0)]
877
    return ShowOctahedron()
878
    
879
def RotateOctahedron(*args):
880
    global OctColor
881
    v = [] 
882
    for arg in args:
883
        if(arg == a):
884
            for i in range(11): 
885
                th = pi*i/10.0  # theta ranges from 0 to pi, inclusive  
886
                X = OctProj(numerical_approx((1 + cos(th))/2), numerical_approx((-1 + cos(th))/2), numerical_approx(-sqrt(2)*sin(th)/2))
887
                nX = OctProj(numerical_approx((-1 - cos(th))/2), numerical_approx((1 - cos(th))/2), numerical_approx(sqrt(2)*sin(th)/2))
888
                Y = OctProj(numerical_approx((-1 + cos(th))/2), numerical_approx((1 + cos(th))/2), numerical_approx(-sqrt(2)*sin(th)/2))
889
                nY = OctProj(numerical_approx((1 - cos(th))/2), numerical_approx((-1 - cos(th))/2), numerical_approx(sqrt(2)*sin(th)/2))
890
                Z = OctProj(numerical_approx(sqrt(2)*sin(th)/2), numerical_approx(sqrt(2)*sin(th)/2), numerical_approx(cos(th)))
891
                nZ = OctProj(numerical_approx(-sqrt(2)*sin(th)/2), numerical_approx(-sqrt(2)*sin(th)/2), numerical_approx(-cos(th)))
892
                Graff = Graphics()
893
                if th < 0.3455: # cutoff for this particular projection function
894
                    Graff += polygon([X,Y,Z], rgbcolor = OctColor[0])
895
                    Graff += polygon([X,nY,Z], rgbcolor = OctColor[2])
896
                    Graff += polygon([X,Y,nZ], rgbcolor = OctColor[4])
897
                    Graff += polygon([X,nY,nZ], rgbcolor = OctColor[6])
898
                elif th < 0.5278:
899
                    Graff += polygon([X,Y,Z], rgbcolor = OctColor[0])
900
                    Graff += polygon([X,nY,Z], rgbcolor = OctColor[2])
901
                    Graff += polygon([X,Y,nZ], rgbcolor = OctColor[4])
902
                    Graff += polygon([nX,Y,Z], rgbcolor = OctColor[1])
903
                elif th < 1.7593:
904
                    Graff += polygon([X,Y,Z], rgbcolor = OctColor[0])
905
                    Graff += polygon([X,nY,Z], rgbcolor = OctColor[2])
906
                    Graff += polygon([nX,nY,Z], rgbcolor = OctColor[3])
907
                    Graff += polygon([nX,Y,Z], rgbcolor = OctColor[1])
908
                elif th < 1.9478:    
909
                    Graff += polygon([nX,nY,nZ], rgbcolor = OctColor[7])
910
                    Graff += polygon([X,nY,Z], rgbcolor = OctColor[2])
911
                    Graff += polygon([nX,nY,Z], rgbcolor = OctColor[3])
912
                    Graff += polygon([nX,Y,Z], rgbcolor = OctColor[1])
913
                else: 
914
                    Graff += polygon([nX,nY,nZ], rgbcolor = OctColor[7])
915
                    Graff += polygon([X,nY,Z], rgbcolor = OctColor[2])
916
                    Graff += polygon([nX,nY,Z], rgbcolor = OctColor[3])
917
                    Graff += polygon([X,nY,nZ], rgbcolor = OctColor[6])   
918
                v.append(Graff)    
919
            Temp = OctColor[0]
920
            OctColor[0] = OctColor[7]
921
            OctColor[7] = Temp
922
            Temp = OctColor[3]
923
            OctColor[3] = OctColor[4]
924
            OctColor[4] = Temp
925
            Temp = OctColor[1]
926
            OctColor[1] = OctColor[5]
927
            OctColor[5] = Temp
928
            Temp = OctColor[2]
929
            OctColor[2] = OctColor[6]
930
            OctColor[6] = Temp      
931
        elif(arg == b):
932
            for i in range(11):
933
                th = 2*pi*i/3/10.0  # theta ranges from 0 to 2pi/3, inclusive 
934
                X = OctProj(numerical_approx(1/3 + 2*cos(th)/3), numerical_approx(1/3 - cos(th)/3 + sin(th)/sqrt(3)), numerical_approx(1/3 - cos(th)/3 - sin(th)/sqrt(3)))
935
                nX = OctProj(numerical_approx(-1/3 - 2*cos(th)/3), numerical_approx(-1/3 + cos(th)/3 - sin(th)/sqrt(3)), numerical_approx(-1/3 + cos(th)/3 + sin(th)/sqrt(3)))
936
                Y = OctProj(numerical_approx(1/3 - cos(th)/3 - sin(th)/sqrt(3)), numerical_approx(1/3 + 2*cos(th)/3), numerical_approx(1/3 - cos(th)/3 + sin(th)/sqrt(3))) 
937
                nY = OctProj(numerical_approx(-1/3 + cos(th)/3 + sin(th)/sqrt(3)), numerical_approx(-1/3 - 2*cos(th)/3), numerical_approx(-1/3 + cos(th)/3 - sin(th)/sqrt(3)))    
938
                Z = OctProj(numerical_approx(1/3 - cos(th)/3 + sin(th)/sqrt(3)), numerical_approx(1/3 - cos(th)/3 - sin(th)/sqrt(3)), numerical_approx(1/3 + 2*cos(th)/3)) 
939
                nZ = OctProj(numerical_approx(-1/3 + cos(th)/3 - sin(th)/sqrt(3)), numerical_approx(-1/3 + cos(th)/3 + sin(th)/sqrt(3)), numerical_approx(-1/3 - 2*cos(th)/3)) 
940
                Graff = Graphics()
941
                if th < .68:
942
                    Graff += polygon([X,Y,Z], rgbcolor = OctColor[0])
943
                    Graff += polygon([X,nY,Z], rgbcolor = OctColor[2])
944
                    Graff += polygon([X,Y,nZ], rgbcolor = OctColor[4])
945
                    Graff += polygon([X,nY,nZ], rgbcolor = OctColor[6])
946
                elif th < 1.09:
947
                    Graff += polygon([X,Y,Z], rgbcolor = OctColor[0])
948
                    Graff += polygon([X,nY,Z], rgbcolor = OctColor[2])
949
                    Graff += polygon([X,Y,nZ], rgbcolor = OctColor[4])
950
                    Graff += polygon([nX,Y,Z], rgbcolor = OctColor[1])  
951
                else:
952
                    Graff += polygon([X,Y,Z], rgbcolor = OctColor[0])
953
                    Graff += polygon([X,nY,Z], rgbcolor = OctColor[2])
954
                    Graff += polygon([nX,nY,Z], rgbcolor = OctColor[3])
955
                    Graff += polygon([nX,Y,Z], rgbcolor = OctColor[1])    
956
                v.append(Graff)
957
            Temp = OctColor[1]
958
            OctColor[1] = OctColor[4]
959
            OctColor[4] = OctColor[2]
960
            OctColor[2] = Temp
961
            Temp = OctColor[3]
962
            OctColor[3] = OctColor[5]
963
            OctColor[5] = OctColor[6]
964
            OctColor[6] = Temp
965
        elif(arg == c):
966
            for i in range(11):
967
                th = pi*i/2/10.0  # theta ranges from 0 to pi/2, inclusive
968
                X = OctProj(1.0,0.0,0.0)
969
                nX = OctProj(-1.0,0.0,0.0)
970
                Y = OctProj(0.0, numerical_approx(cos(th)), numerical_approx(-sin(th)))
971
                nY = OctProj(0.0, numerical_approx(-cos(th)), numerical_approx(sin(th)))
972
                Z = OctProj(0.0, numerical_approx(sin(th)), numerical_approx(cos(th)))
973
                nZ = OctProj(0.0, numerical_approx(-sin(th)), numerical_approx(-cos(th)))
974
                Graff = Graphics()
975
                Graff += polygon([X,Y,Z], rgbcolor = OctColor[0])
976
                Graff += polygon([X,nY,Z], rgbcolor = OctColor[2])
977
                Graff += polygon([X,Y,nZ], rgbcolor = OctColor[4])
978
                Graff += polygon([X,nY,nZ], rgbcolor = OctColor[6])
979
                v.append(Graff)
980
            Temp = OctColor[0]
981
            OctColor[0] = OctColor[2]
982
            OctColor[2] = OctColor[6]
983
            OctColor[6] = OctColor[4]
984
            OctColor[4] = Temp
985
            Temp = OctColor[1]
986
            OctColor[1] = OctColor[3]
987
            OctColor[3] = OctColor[7]
988
            OctColor[7] = OctColor[5]
989
            OctColor[5] = Temp
990
        else:
991
            print("Invalid rotation")
992
    aniOct = animate(v, xmin = -1, xmax = 1, ymin = -1.03, ymax = 1.03, aspect_ratio = 1, axes = false)
993
    return aniOct.show(delay = 50, iterations = 1)
994
    
995
PuzColor = [(1,0,0), (1,0,0), (1,0,0), (1,0,0), (1,0,0), (1,0,0), (1,1,0), (1,1,0), (1,1,0), (1,1,0), (1,1,0), (1,1,0), (0,1,0), (0,1,0), (0,1,0), (0,1,0), (0,1,0), (0,1,0), (0,0,1), (0,0,1), (0,0,1), (0,0,1), (0,0,1), (0,0,1)]    
996
    
997
def PuzProj(x,y,z):  # projects a 3D point into 2D, using the projection angle that was used in TeX
998
    return (-0.45*x + 0.55*y+0.05*z, -0.3125*x -0.3125*y + .5625*z)    
999
    
1000
def ShowPuzzle():
1001
    global PuzColor
1002
    Graff = Graphics()
1003
    if PuzColor[0] != 0:
1004
        Graff += polygon([PuzProj(-1,-1,-1),PuzProj(-3,1,-1),PuzProj(-3,-1,1)], rgbcolor = PuzColor[18])
1005
        Graff += polygon([PuzProj(-1,-1,-1),PuzProj(-3,-1,1),PuzProj(-1,-3,1)], rgbcolor = PuzColor[19])
1006
        Graff += polygon([PuzProj(-1,-1,-1),PuzProj(-1,-3,1),PuzProj(1,-3,-1)], rgbcolor = PuzColor[20])
1007
        Graff += polygon([PuzProj(-1,-1,-1),PuzProj(1,-3,-1),PuzProj(1,-1,-3)], rgbcolor = PuzColor[21])
1008
        Graff += polygon([PuzProj(-1,-1,-1),PuzProj(1,-1,-3),PuzProj(-1,1,-3)], rgbcolor = PuzColor[22])
1009
        Graff += polygon([PuzProj(-1,-1,-1),PuzProj(-1,1,-3),PuzProj(-3,1,-1)], rgbcolor = PuzColor[23])
1010
        Graff += polygon([PuzProj(-1,1,1),PuzProj(1,3,1),PuzProj(1,1,3)], rgbcolor = PuzColor[0])
1011
        Graff += polygon([PuzProj(-1,1,1),PuzProj(1,1,3),PuzProj(-1,-1,3)], rgbcolor = PuzColor[1])
1012
        Graff += polygon([PuzProj(-1,1,1),PuzProj(-1,-1,3),PuzProj(-3,-1,1)], rgbcolor = PuzColor[2])
1013
        Graff += polygon([PuzProj(-1,1,1),PuzProj(-3,-1,1),PuzProj(-3,1,-1)], rgbcolor = PuzColor[3])
1014
        Graff += polygon([PuzProj(-1,1,1),PuzProj(-3,1,-1),PuzProj(-1,3,-1)], rgbcolor = PuzColor[4])
1015
        Graff += polygon([PuzProj(-1,1,1),PuzProj(-1,3,-1),PuzProj(1,3,1)], rgbcolor = PuzColor[5])
1016
        Graff += polygon([PuzProj(1,-1,1),PuzProj(1,1,3),PuzProj(3,1,1)], rgbcolor = PuzColor[6])
1017
        Graff += polygon([PuzProj(1,-1,1),PuzProj(3,1,1),PuzProj(3,-1,-1)], rgbcolor = PuzColor[7])
1018
        Graff += polygon([PuzProj(1,-1,1),PuzProj(3,-1,-1),PuzProj(1,-3,-1)], rgbcolor = PuzColor[8])
1019
        Graff += polygon([PuzProj(1,-1,1),PuzProj(1,-3,-1),PuzProj(-1,-3,1)], rgbcolor = PuzColor[9])
1020
        Graff += polygon([PuzProj(1,-1,1),PuzProj(-1,-3,1),PuzProj(-1,-1,3)], rgbcolor = PuzColor[10])
1021
        Graff += polygon([PuzProj(1,-1,1),PuzProj(-1,-1,3),PuzProj(1,1,3)], rgbcolor = PuzColor[11])
1022
        Graff += polygon([PuzProj(1,1,-1),PuzProj(3,1,1),PuzProj(1,3,1)], rgbcolor = PuzColor[12])
1023
        Graff += polygon([PuzProj(1,1,-1),PuzProj(1,3,1),PuzProj(-1,3,-1)], rgbcolor = PuzColor[13])
1024
        Graff += polygon([PuzProj(1,1,-1),PuzProj(-1,3,-1),PuzProj(-1,1,-3)], rgbcolor = PuzColor[14])
1025
        Graff += polygon([PuzProj(1,1,-1),PuzProj(-1,1,-3),PuzProj(1,-1,-3)], rgbcolor = PuzColor[15])
1026
        Graff += polygon([PuzProj(1,1,-1),PuzProj(1,-1,-3),PuzProj(3,-1,-1)], rgbcolor = PuzColor[16])
1027
        Graff += polygon([PuzProj(1,1,-1),PuzProj(3,-1,-1),PuzProj(3,1,1)], rgbcolor = PuzColor[17])
1028
        Graff += line([(-0.3214, 0.0625),(0.8214, 0.0625)], rgbcolor = (0,0,0))
1029
        Graff += line([(-0.57857, -0.6875),(0.33571, 0.9125)], rgbcolor = (0,0,0))
1030
        Graff += line([(0.27857, -0.6875),(-0.23571, 0.2125)], rgbcolor = (0,0,0))
1031
        Graff += line([PuzProj(3,1,1), PuzProj(1,3,1), PuzProj(1,1,3), PuzProj(3,1,1), PuzProj(3,-1,-1), PuzProj(-1,-1,3), PuzProj(-1,3,-1), PuzProj(3,-1,-1), PuzProj(1,-3,-1), PuzProj(1,1,3), PuzProj(-3,1,-1), PuzProj(-3,-1,1), PuzProj(1,3,1), PuzProj(1,-1,-3), PuzProj(-1,1,-3), PuzProj(3,1,1), PuzProj(-1,-3,1), PuzProj(-1,-1,3), PuzProj(1,1,3)], rgbcolor = (0,0,0))
1032
        Graff += line([PuzProj(1,-3,-1),PuzProj(-1,-3,1)], rgbcolor = (0,0,0))
1033
        Graff += line([PuzProj(3,-1,-1),PuzProj(1,-1,-3)], rgbcolor = (0,0,0))
1034
        Graff += line([PuzProj(-1,-1,3),PuzProj(-3,-1,1)], rgbcolor = (0,0,0))
1035
        Graff += line([PuzProj(-1,3,-1),PuzProj(1,3,1)], rgbcolor = (0,0,0))
1036
        Graff += line([PuzProj(-1,3,-1),PuzProj(-3,1,-1)], rgbcolor = (0,0,0))
1037
        Graff += line([PuzProj(-1,3,-1),PuzProj(-1,1,-3)], rgbcolor = (0,0,0))        
1038
    else:    
1039
        Graff += polygon([PuzProj(-1,-1,-1),PuzProj(-3,1,-1),PuzProj(-3,-1,1)], rgbcolor = PuzColor[18])
1040
        Graff += polygon([PuzProj(-1,-1,-1),PuzProj(-1,-3,1),PuzProj(1,-3,-1)], rgbcolor = PuzColor[20])
1041
        Graff += polygon([PuzProj(-1,-1,-1),PuzProj(1,-1,-3),PuzProj(-1,1,-3)], rgbcolor = PuzColor[22])
1042
        Graff += polygon([PuzProj(-1,1,1),PuzProj(1,1,3),PuzProj(-1,-1,3)], rgbcolor = PuzColor[1])
1043
        Graff += polygon([PuzProj(-1,1,1),PuzProj(-3,-1,1),PuzProj(-3,1,-1)], rgbcolor = PuzColor[3])
1044
        Graff += polygon([PuzProj(-1,1,1),PuzProj(-1,3,-1),PuzProj(1,3,1)], rgbcolor = PuzColor[5])
1045
        Graff += polygon([PuzProj(1,-1,1),PuzProj(3,1,1),PuzProj(3,-1,-1)], rgbcolor = PuzColor[7])
1046
        Graff += polygon([PuzProj(1,-1,1),PuzProj(1,-3,-1),PuzProj(-1,-3,1)], rgbcolor = PuzColor[9])
1047
        Graff += polygon([PuzProj(1,-1,1),PuzProj(-1,-1,3),PuzProj(1,1,3)], rgbcolor = PuzColor[11])
1048
        Graff += polygon([PuzProj(1,1,-1),PuzProj(1,3,1),PuzProj(-1,3,-1)], rgbcolor = PuzColor[13])
1049
        Graff += polygon([PuzProj(1,1,-1),PuzProj(-1,1,-3),PuzProj(1,-1,-3)], rgbcolor = PuzColor[15])
1050
        Graff += polygon([PuzProj(1,1,-1),PuzProj(3,-1,-1),PuzProj(3,1,1)], rgbcolor = PuzColor[17])
1051
        Graff += line([PuzProj(3,1,1), PuzProj(-1,-3,1), PuzProj(1,-3,-1), PuzProj(1,1,3), PuzProj(-1,-1,3), PuzProj(3,-1,-1), PuzProj(3,1,1), PuzProj(-1,1,-3), PuzProj(1,-1,-3), PuzProj(1,3,1), PuzProj(-1,3,-1), PuzProj(3,-1,-1)], rgbcolor = (0,0,0)) 
1052
        Graff += line([PuzProj(1,3,1), PuzProj(-3,-1,1), PuzProj(-3,1,-1), PuzProj(1,1,3)], rgbcolor = (0,0,0)) 
1053
        Graff += line([PuzProj(-1,-1,3), PuzProj(-1,3,-1)], rgbcolor = (0,0,0)) 
1054
        Graff += line([(-2.15, 0.0625), (-1.2457,  0.0625)], rgbcolor = (0,0,0)) 
1055
        Graff += line([(1.85, 0.0625), (1.3357,  0.0625)], rgbcolor = (0,0,0))   
1056
        Graff += line([(0.85, -1.6875), (0.5643, -1.1875)], rgbcolor = (0,0,0))  
1057
        Graff += line([(-1.15, -1.6875), (-0.8543, -1.1875)], rgbcolor = (0,0,0))   
1058
        Graff += line([(-1.15, 1.8125), (-0.6929, 1.0125)], rgbcolor = (0,0,0)) 
1059
        Graff += line([(0.85, 1.8125), (0.5929, 1.3625)], rgbcolor = (0,0,0))   
1060
        Graff += line([(-0.6237, -0.7664), (0.3763, 0.9836)], rgbcolor = (0,0,0))  
1061
        Graff += line([(-0.87, 0.0625), (1.13,  0.0625)], rgbcolor = (0,0,0))   
1062
        Graff += line([(0.4654, -1.0144), (-0.5346,  0.7356)], rgbcolor = (0,0,0))   
1063
    return Graff.show(axes = false, aspect_ratio = 1)
1064
    
1065
def InitPuzzle():
1066
    global PuzColor
1067
    global f
1068
    global b
1069
    global r
1070
    global l   
1071
    f = var('f')
1072
    b = var('b')
1073
    r = var('r')
1074
    l = var('l')
1075
    PuzColor = [(1,0,0), (1,0,0), (1,0,0), (1,0,0), (1,0,0), (1,0,0), (1,1,0), (1,1,0), (1,1,0), (1,1,0), (1,1,0), (1,1,0), (0,1,0), (0,1,0), (0,1,0), (0,1,0), (0,1,0), (0,1,0), (0,0,1), (0,0,1), (0,0,1), (0,0,1), (0,0,1), (0,0,1)]    
1076
    return ShowPuzzle()   
1077
    
1078
def RotatePuzzle(*args):
1079
    global PuzColor   
1080
    for arg in args:
1081
        if arg == f:
1082
            Temp = PuzColor[0]
1083
            PuzColor[0] = PuzColor[6]
1084
            PuzColor[6] = PuzColor[12]
1085
            PuzColor[12] = Temp
1086
            Temp = PuzColor[1]
1087
            PuzColor[1] = PuzColor[7]
1088
            PuzColor[7] = PuzColor[13]
1089
            PuzColor[13] = Temp
1090
            Temp = PuzColor[5]
1091
            PuzColor[5] = PuzColor[11]
1092
            PuzColor[11] = PuzColor[17]
1093
            PuzColor[17] = Temp
1094
        elif arg == b:    
1095
            Temp = PuzColor[1]
1096
            PuzColor[1] = PuzColor[18]
1097
            PuzColor[18] = PuzColor[9]
1098
            PuzColor[9] = Temp
1099
            Temp = PuzColor[2]
1100
            PuzColor[2] = PuzColor[19]
1101
            PuzColor[19] = PuzColor[10]
1102
            PuzColor[10] = Temp
1103
            Temp = PuzColor[3]
1104
            PuzColor[3] = PuzColor[20]
1105
            PuzColor[20] = PuzColor[11]
1106
            PuzColor[11] = Temp            
1107
        elif arg == r:    
1108
            Temp = PuzColor[3]
1109
            PuzColor[3] = PuzColor[13]
1110
            PuzColor[13] = PuzColor[22]
1111
            PuzColor[22] = Temp
1112
            Temp = PuzColor[4]
1113
            PuzColor[4] = PuzColor[14]
1114
            PuzColor[14] = PuzColor[23]
1115
            PuzColor[23] = Temp
1116
            Temp = PuzColor[5]
1117
            PuzColor[5] = PuzColor[15]
1118
            PuzColor[15] = PuzColor[18]
1119
            PuzColor[18] = Temp              
1120
        elif arg == l:    
1121
            Temp = PuzColor[7]
1122
            PuzColor[7] = PuzColor[20]
1123
            PuzColor[20] = PuzColor[15]
1124
            PuzColor[15] = Temp
1125
            Temp = PuzColor[8]
1126
            PuzColor[8] = PuzColor[21]
1127
            PuzColor[21] = PuzColor[16]
1128
            PuzColor[16] = Temp
1129
            Temp = PuzColor[9]
1130
            PuzColor[9] = PuzColor[22]
1131
            PuzColor[22] = PuzColor[17]
1132
            PuzColor[17] = Temp 
1133
        else:
1134
            print("Invalid rotation")
1135
        ShowPuzzle()
1136
    return None    
1137

1138
def PuzzlePosition1():
1139
    global PuzColor
1140
    PuzColor = [(1,1,0), (1,0,0), (1,0,0), (1,0,0), (1,0,0), (1,0,0), (0,1,0), (1,1,0), (1,1,0), (1,1,0), (1,1,0), (1,1,0), (1,0,0), (0,1,0), (0,1,0), (0,1,0), (0,1,0), (0,1,0), (0,0,1), (0,0,1), (0,0,1), (0,0,1), (0,0,1), (0,0,1)]
1141
    return ShowPuzzle()        
1142
    
1143
def PuzzlePosition2():
1144
    global PuzColor
1145
    PuzColor = [(1,0,0), (1,1,0), (1,0,0), (0,0,1), (1,0,0), (0,1,0), (1,1,0), (0,1,0), (1,1,0), (0,0,1), (1,1,0), (1,0,0), (0,1,0), (1,0,0), (0,1,0), (0,0,1), (0,1,0), (1,1,0), (1,0,0), (0,0,1), (1,1,0), (0,0,1), (0,1,0), (0,0,1)]
1146
    return ShowPuzzle()     
1147
    
1148
def HideCorners():
1149
    global PuzColor
1150
    for i in [0,2,4,6,8,10,12,14,16,19,21,23]:
1151
        PuzColor[i] = 0
1152
    return ShowPuzzle()   
1153

1154
    
1155
def ShowRationals(a,b):
1156
    import fractions
1157
    A = numerical_approx(a)
1158
    B = numerical_approx(b)
1159
    if A == B:
1160
        print "Endpoints must be different"
1161
        return None   
1162
    if A > B:
1163
        q = A
1164
        A = B
1165
        B = q
1166
    QuoList = []
1167
    r = 1
1168
    q = 1
1169
    while r < 18:
1170
        PrintPoint = False
1171
        for p in range(ceil(A*q), B*q +1):
1172
            if fractions.gcd(p, q) == 1:
1173
                QuoList.append((numerical_approx(p)/numerical_approx(q), 0.75^(r-1)))
1174
                PrintPoint = True
1175
        if PrintPoint:
1176
            r = r+1  
1177
        q = q+1
1178
    P = list_plot(QuoList, aspect_ratio = (B-A)/2, frame = True, axes = False)        
1179
#    G = Graphics()
1180
#    G += point(QuoList)
1181
#    G.show(aspect_ratio = (B-A)/2, frame = True, axes = False) 
1182
#    return G
1183
    return P
1184

1185
def CountableQ(n):
1186
    import fractions
1187
    global QuoList
1188
    global G
1189
    G = Graphics()
1190
    NN = abs(floor(n))
1191
    if NN == 0:
1192
        NN = 1
1193
    QuoList = []    
1194
    LineList = []
1195
    for q in range(1,17):
1196
        for p in range(-NN*q, NN*q +1):
1197
            if fractions.gcd(p, q) == 1:
1198
                QuoList.append((numerical_approx(p)/numerical_approx(q), 0.75^(q-1)))
1199
    G += point(QuoList)            
1200
    for q in range(-floor(NN/2),(NN+1)/2):
1201
        G += line([(2*q,1),(2*q+1,1)])
1202
    for q in range(1, NN+1):
1203
        G += line([(q-.5,.75),(q,1)])    
1204
        G += line([(-q+.5,.75),(-q,1)])    
1205
        G += line([(-numerical_approx(1)/(q+1),(.75)^q),(numerical_approx(1)/(q+1),(.75)^q)])
1206
    for p in range(1,NN):
1207
        for q in range(2, NN-p+2):
1208
            G += line([(p+numerical_approx(1)/(q+1)-1,(.75)^q),  (p-numerical_approx(1)/(q+1),(.75)^q),  (p+numerical_approx(1)/q,(.75)^(q-1))])    
1209
            G += line([(1-p-numerical_approx(1)/(q+1),(.75)^q),  (numerical_approx(1)/(q+1)-p,(.75)^q),  (-p-numerical_approx(1)/q,(.75)^(q-1))])
1210
    return G.show(aspect_ratio = NN, frame = True, axes = False)           
1211
        
1212
def PolarPlot():
1213
    Graff = Graphics()
1214
    Graff.set_aspect_ratio(1.0)
1215
    Graff += line([(0,0), (3,4)], rgbcolor=(0,0,0))
1216
    Graff += arc((0,0), 1, sector = (0, 0.92729), rgbcolor=(0,0,0))
1217
    Graff += text("0", (0.7, 0.4), rgbcolor=(0,0,0), fontsize = 18)
1218
    Graff += text("-", (0.7, 0.4), rgbcolor=(0,0,0), fontsize = 18) # How to make a theta
1219
    Graff += text("r", (1.5, 2.5), rgbcolor=(0,0,0), fontsize = 18)
1220
    Graff += text("(x,y) = x + y i", (3.2, 4), rgbcolor=(0,0,0), horizontal_alignment = "left", fontsize = 18) 
1221
    return Graff.show() 
1222

1223
def DrawGalois(n):
1224
    Graff = Graphics()
1225
    Graff += line([(.47, .95), (0.05, .55)], rgbcolor=(0,0,0))
1226
    Graff += line([(.49, .95), (0.35, .55)], rgbcolor=(0,0,0))
1227
    Graff += line([(.51, .95), (0.65, .55)], rgbcolor=(0,0,0))
1228
    Graff += line([(.53, .95), (0.95, .55)], rgbcolor=(0,1,0)) 
1229
    Graff += line([(.47, .05), (0.05, .45)], rgbcolor=(0,1,0))
1230
    Graff += line([(.49, .05), (0.35, .45)], rgbcolor=(0,1,0))
1231
    Graff += line([(.51, .05), (0.65, .45)], rgbcolor=(0,1,0))
1232
    Graff += line([(.53, .05), (0.95, .45)], rgbcolor=(0,1,0))
1233
    Graff += line([(.50, .05), (0.50, .95)], rgbcolor=(0,1,0))
1234
    Graff += text("Z", (0.76, 0.77), rgbcolor=(0,0,1), fontsize = 12)
1235
    Graff += text("2", (0.775, 0.76), rgbcolor=(0,0,1), fontsize = 9)
1236
    Graff += text("S", (0.52, 0.77), rgbcolor=(0,0,1), fontsize = 12)
1237
    Graff += text("3", (0.535, 0.76), rgbcolor=(0,0,1), fontsize = 9)
1238
    Graff += text("Z", (0.22, 0.23), rgbcolor=(0,0,1), fontsize = 12)
1239
    Graff += text("2", (0.235, 0.22), rgbcolor=(0,0,1), fontsize = 9)
1240
    Graff += text("Z", (0.39, 0.23), rgbcolor=(0,0,1), fontsize = 12)
1241
    Graff += text("2", (0.405, 0.22), rgbcolor=(0,0,1), fontsize = 9)
1242
    Graff += text("Z", (0.59, 0.23), rgbcolor=(0,0,1), fontsize = 12)
1243
    Graff += text("2", (0.605, 0.22), rgbcolor=(0,0,1), fontsize = 9)
1244
    Graff += text("Z", (0.76, 0.23), rgbcolor=(0,0,1), fontsize = 12)
1245
    Graff += text("3", (0.775, 0.22), rgbcolor=(0,0,1), fontsize = 9)
1246
    if n == 1:    
1247
        Graff += text("Q", (0.494, 1.0), rgbcolor=(0,0,0), fontsize = 12, horizontal_alignment = "left")
1248
        Graff += text("I", (0.501, 1.0), rgbcolor=(0,0,0), fontsize = 10, horizontal_alignment = "left")
1249
        Graff += text("Q(2    )", (-0.04, 0.5), rgbcolor=(0,0,0), fontsize = 12, horizontal_alignment = "left")
1250
        Graff += text("I", (-0.033, 0.5), rgbcolor=(0,0,0), fontsize = 10, horizontal_alignment = "left")
1251
        Graff += text("1/3", (0.005, 0.51), rgbcolor=(0,0,0), fontsize = 9, horizontal_alignment = "left")
1252
        Graff += text("Q(w  2   )", (0.30, 0.5), rgbcolor=(0,0,0), fontsize = 12, horizontal_alignment = "left")
1253
        Graff += text("I", (0.307, 0.5), rgbcolor=(0,0,0), fontsize = 10, horizontal_alignment = "left")
1254
        Graff += text("1/3", (0.38, 0.51), rgbcolor=(0,0,0), fontsize = 9, horizontal_alignment = "left")
1255
        Graff += text("3", (0.35, 0.49), rgbcolor=(0,0,0), fontsize = 9, horizontal_alignment = "left")
1256
        Graff += text("Q(w  2   )", (0.60, 0.5), rgbcolor=(0,0,0), fontsize = 12, horizontal_alignment = "left")
1257
        Graff += text("I", (0.607, 0.5), rgbcolor=(0,0,0), fontsize = 10, horizontal_alignment = "left")
1258
        Graff += text("2  1/3", (0.653, 0.51), rgbcolor=(0,0,0), fontsize = 9, horizontal_alignment = "left")        
1259
        Graff += text("3", (0.65, 0.49), rgbcolor=(0,0,0), fontsize = 9, horizontal_alignment = "left")          
1260
        Graff += text("Q((-3)    )", (1.0, 0.5), rgbcolor=(0,0,0), fontsize = 12, horizontal_alignment = "left")
1261
        Graff += text("1/2", (1.08, 0.51), rgbcolor=(0,0,0), fontsize = 9, horizontal_alignment = "left") 
1262
        Graff += text("I", (1.007, 0.5), rgbcolor=(0,0,0), fontsize = 10, horizontal_alignment = "left")
1263
        Graff += text("Q(2    , w  2   )", (0.4, 0.0), rgbcolor=(0,0,0), fontsize = 12, horizontal_alignment = "left")
1264
        Graff += text("I", (0.407, 0.0), rgbcolor=(0,0,0), fontsize = 10, horizontal_alignment = "left")
1265
        Graff += text("1/3           1/3", (0.443, 0.01), rgbcolor=(0,0,0), fontsize = 9, horizontal_alignment = "left")
1266
        Graff += text("3", (0.513, -0.01), rgbcolor=(0,0,0), fontsize = 9, horizontal_alignment = "left")
1267
    if n == 2:    
1268
        Graff += text("S", (0.49, 1.0), rgbcolor=(0,0,0), fontsize = 12)
1269
        Graff += text("3", (0.505, 0.99), rgbcolor=(0,0,0), fontsize = 9)
1270
        Graff += text("{(), (1 2)}", (0.0, 0.5), rgbcolor=(0,0,0), fontsize = 12, horizontal_alignment = "right")
1271
        Graff += text("{(), (1 3)}", (0.35, 0.5), rgbcolor=(0,0,0), fontsize = 12)
1272
        Graff += text("{(), (2 3)}", (0.65, 0.5), rgbcolor=(0,0,0), fontsize = 12)
1273
        Graff += text("A", (1.0, 0.5), rgbcolor=(0,0,0), fontsize = 12)
1274
        Graff += text("3", (1.015, 0.49), rgbcolor=(0,0,0), fontsize = 9)
1275
        Graff += text("{()}", (0.5, 0.0), rgbcolor=(0,0,0), fontsize = 12)    
1276
    return Graff.show(axes = False)
1277
        
1278
######################################################
1279
##                                                  ##
1280
## Utility routines - each of these does something  ##
1281
## simple in a Python way                           ##
1282
##                                                  ##
1283
######################################################
1284

1285
def Flatten(Set):  # Fast way to flatten a list by 1 level
1286
    if isinstance(Set, list):
1287
        if len(Set) > 0:
1288
            if isinstance(Set[0], (list, GroupSet)):
1289
                return [item for sublist in Set for item in sublist]  
1290
    return Set
1291

1292
def Uniquify(Set):  # Eliminate duplicates in a list while preserving order
1293
    if isinstance(Set, list):
1294
        out = []
1295
        [out.append(item) for item in Set if item not in out]
1296
    return out
1297
    
1298
def IsSubset(shortlist, longlist): # determines if shortlist is a subset of longlist
1299
    for item in shortlist:
1300
        if not(item in longlist):
1301
            return false
1302
    return true        
1303
   
1304
def Intersection(arg1, *args):
1305
    if len(args) == 0:        #If there is only one arguement, see if it is a list of lists.
1306
        T = list(arg1)
1307
        if isinstance(T[0], list):
1308
            S = Uniquify(T[0])
1309
            for i in range(1,len(T)):
1310
                S = [item for item in S if item in T[i] ]
1311
            return sorted(S)
1312
    if isinstance(arg1, GroupSet):
1313
        S = Uniquify(arg1._List)
1314
    else:        
1315
        S = Uniquify(arg1)
1316
    for arg in args:
1317
        if isinstance(arg, GroupSet):
1318
            S = [item for item in S if item in arg._List ]
1319
        else:   
1320
            S = [item for item in S if item in arg ]
1321
    if isinstance(arg1, GroupSet):    
1322
        return GroupSet(sorted(S))
1323
    return sorted(S) 
1324
    
1325
def ConvertIdentity(str):  # Converts <Identity ...> in a string to the name of the identity
1326
    ss = str
1327
    ii = ss.find("<")         # if we find a <, we can assume that everything up to > will be converted to
1328
    if ii > -1:               # the identity element.
1329
        jj = ss.find(">",ii)
1330
        if jj > -1:                        # here is a twist.  When the identity element is displayed in a pc
1331
            if(len(ss)>jj+7):              # group, it is "<identity> of ..."  (of... is not in the brackets!)
1332
                if(ss[jj+1:jj+8] == " of ..."):
1333
                    jj = jj + 7
1334
            ss = ss[:ii] + IdentityElement + ss[jj+1:]  
1335
    return ss  
1336
    
1337
#################################################################
1338
##                                                             ##
1339
##  Many of the group elements are defined as instances of a   ##
1340
##  class.  Here we see the classes that define the different  ##
1341
##  types of groups.                                           ##
1342
##                                                             ##
1343
#################################################################
1344
    
1345
class SumModN:
1346
    def __init__(self, x, mod):
1347
        self._x = x
1348
        self._mod = mod    
1349
    def __eq__(self, other):
1350
        if isinstance(other, SumModN):
1351
            return (self._x == other._x) and (self._mod == other._mod)
1352
        else:
1353
            return NotImplemented            
1354
    def __ne__(self, other):
1355
        result = self.__eq__(other)
1356
        if result is NotImplemented:
1357
            return result
1358
        return not result          
1359
    def __mul__(self, other):
1360
        if isinstance(other, SumModN):
1361
            if self._mod == other._mod:
1362
                return SumModN((self._x + other._x) % self._mod, self._mod ) 
1363
        return NotImplemented
1364
    # We also allow addition by an integer, to allow for the CircleGraph with Add(n)    
1365
    def __add__(self, other):
1366
        if isinstance(other, (int, long, IntType)):
1367
            return SumModN((self._x + other) % self._mod, self._mod )
1368
        return NotImplemented   
1369
    def __radd__(self, other):
1370
        if isinstance(other, (int, long, IntType)):
1371
            return SumModN((self._x + other) % self._mod, self._mod )
1372
        return NotImplemented    
1373
    def __pow__(self, exp):
1374
        if isinstance(exp, (int, long, IntType)):
1375
            return SumModN((self._x * exp) % self._mod, self._mod )
1376
        return NotImplemented   
1377
    # We want subgroups to appear in order, such as [0,2,4,6,8] 
1378
    def __cmp__(self, other):
1379
        if isinstance(other, SumModN):
1380
            if self._mod == other._mod:
1381
                return cmp(self._x, other._x)                      
1382
        else:
1383
            return NotImplemented        
1384
    def __str__(self):
1385
        return str(self._x)    
1386
    def __repr__(self):
1387
        return str(self._x)   
1388

1389
####################################################################
1390
##
1391
## The GroupSet class is in essence a wrapper for a list of elements.
1392
## GroupSets differ from lists and sets in the following ways:
1393
## 1) we can add and multiply two GroupSets together, like cosets
1394
## 2) Unlike sets, GroupSets can contain mutable items, such as other GroupSets
1395
## 3) Unlike sets, GroupSets do not have to be sorted, although the sum and 
1396
##     product of GroupSets yields a _sorted_ GroupSet, for comparison
1397
## 4) Unlike lists, GroupSets are displayed with curly braces  
1398
##
1399
########################################################################
1400

1401
class GroupSet:                    
1402
    def __init__(self, L):
1403
        if isinstance(L, list):
1404
            NewList = [] 
1405
            for item in L: # Strip the wrappers off of the Gap elements
1406
                if isinstance(item, GapElement):
1407
                    NewList.append(item._x)
1408
                else:
1409
                    NewList.append(item)
1410
            self._List = NewList        
1411
        else:
1412
            return NotImplemented    
1413
    def __eq__(self, other):
1414
        if isinstance(other, GroupSet):
1415
            # Since both cosets should already be sorted, we just have to compare the lists to each other.
1416
            return(self._List == other._List)
1417
        return NotImplemented
1418
    def __ne__(self, other):
1419
        result = self.__eq__(other)
1420
        if result is NotImplemented:
1421
            return result
1422
        return not result 
1423
    def __mul__(self, other):
1424
        if isinstance(other, GroupSet):
1425
            # we are multiplying two Cosets together
1426
            if (len(self._List) == 0) or (len(other._List) == 0): # if either is the empty coset, return the empty coset.
1427
                return GroupSet([])
1428
            if isinstance(self._List[0], GapInstance) and isinstance(other._List[0], GapInstance):
1429
                # deligate the job to Gap for speed purposes
1430
                # if _QuickMult_ is set to true, we can assume that these are both cosets of the same subgroup.
1431
                # in which case, we can save time by multiplying the FIRST element of the first coset by the second coset.
1432
                # Actually, this didn't save time
1433
                #if _QuickMult_:
1434
                #    LL1 = gap([self._List[0]])
1435
                #else:    
1436
                LL1 = gap(self._List)
1437
                LL2 = gap(other._List)
1438
                G = gap([])
1439
                LL1.LeftSide()
1440
                LL2.RightSide()
1441
                Prod1 = list(G.Prod())
1442
                return GroupSet(Prod1)
1443
            else:    
1444
                # Find all possible products
1445
                Prod = []
1446
                for x1 in self._List:
1447
                    for x2 in other._List:
1448
                        ##  We could be dealing with cosets of cosets.
1449
                        ##  In which case, there is the posibility that x1 or x2 is a GroupSet
1450
                        ##  This will not cause a problem if x1 is a GroupSet, since we will just have recursion.
1451
                        ##  But if x2 is a GroupSet and x1 is NOT, we need to convert x1 to a GroupSet with a single element.
1452
                        if isinstance(x2, GroupSet) and (not isinstance(x1, GroupSet)):
1453
                            p = GroupSet([x1]) * x2
1454
                        else:    
1455
                            p = x1 * x2    
1456
                        Prod.append(p)
1457
                Prod.sort()
1458
                for i in range(len(Prod)-1,0,-1): # Count backwards
1459
                    if Prod[i] == Prod[i-1]:      # Delete duplicates - which will be together after the sort
1460
                        del Prod[i]    
1461
                return GroupSet(Prod)  
1462
        # if other is not a coset, it might be a list or single element
1463
        if isinstance(other, list):
1464
            return self * GroupSet(other)
1465
        # last chance - single element    
1466
        return self * GroupSet([other])         
1467
    def __rmul__(self, other):
1468
        # obviously other is not a GroupSet, or it would have been covered by __mul__
1469
        if isinstance(other, list):
1470
            return GroupSet(other) * self
1471
        # last chance - single element    
1472
        return GroupSet([other]) * self  
1473
    def __add__(self, other):
1474
        if isinstance(other, GroupSet):
1475
            # we are adding two Cosets together
1476
            if (len(self._List) == 0) or (len(other._List) == 0): # if either is the empty coset, return the empty coset.
1477
                return GroupSet([])
1478
            if isinstance(self._List[0], GapInstance) and isinstance(other._List[0], GapInstance):
1479
                # deligate the job to Gap for speed purposes
1480
                LL1 = gap(self._List)
1481
                LL2 = gap(other._List)
1482
                LL1.LeftSide()
1483
                Sum1 = list(LL2.AddCoset())
1484
                return GroupSet(Sum1)
1485
            else:    
1486
                # Find all possible sums
1487
                PSum = []
1488
                for x1 in self._List:
1489
                    for x2 in other._List:
1490
                        ##  We could be dealing with cosets of cosets.
1491
                        ##  In which case, there is the posibility that x1 or x2 is a GroupSet
1492
                        ##  This will not cause a problem if x1 is a GroupSet, since we will just have recursion.
1493
                        ##  But if x2 is a GroupSet and x1 is NOT, we need to convert x1 to a GroupSet with a single element.
1494
                        if isinstance(x2, GroupSet) and (not isinstance(x1, GroupSet)):
1495
                            p = GroupSet([x1]) * x2
1496
                        else:    
1497
                            p = x1 * x2    
1498
                        PSum.append(p)
1499
                PSum.sort()
1500
                for i in range(len(PSum)-1,0,-1): # Count backwards
1501
                    if PSum[i] == PSum[i-1]:      # Delete duplicates - which will be together after the sort
1502
                        del PSum[i]    
1503
                return GroupSet(PSum)  
1504
        # if other is not a coset, it might be a list or single element
1505
        if isinstance(other, list):
1506
            return self + GroupSet(other)
1507
        # last chance - single element    
1508
        return self + GroupSet([other])
1509
    def __radd__(self, other):
1510
        # obviously other is not a GroupSet, or it would have been covered by __add__
1511
        if isinstance(other, list):
1512
            return GroupSet(other) + self
1513
        # last chance - single element    
1514
        return GroupSet([other]) + self      
1515
    def __neg__(self):
1516
        # this is only used for rings.  We take the negative of all of the elements in the list.
1517
        N = []
1518
        for item in self._List:
1519
            N.append(-item)
1520
        N.sort()
1521
        return GroupSet(N)    
1522
    def __sub__(self, other):
1523
        return self + (-other) # this should work is all cases.    
1524
    def __cmp__(self, other):
1525
        if isinstance(other, GroupSet):
1526
            return cmp(sorted(self._List), sorted(other._List))
1527
        return NotImplemented     
1528
    def __len__(self):
1529
        return len(self._List)    
1530
    def __iter__(self):              # Allows iteration constructions such as
1531
        return iter(self._List)      # [order(x) for x in G] 
1532
    def __getitem__(self, index):
1533
        return self._List[index]    
1534
    def __pow__(self, exp):          # Raising a set to a power, think of it as a coset (element of a quotient group)
1535
        if isinstance(exp, (int, long, IntType)):
1536
            x = self._List[0]^(exp - 1)
1537
            return self * x 
1538
        return NotImplemented      
1539
    def Sort(self):
1540
        self._List.sort()
1541
        return None    
1542
    def Conformity(self, other):        # puts the elements in self into the format given in the larger group other
1543
        if isinstance(self._List[0], GroupSet):
1544
            ## self is a coset of a coset.  Recursively do the Conformity to each of self's members
1545
            for i in range(len(self._List)):
1546
                self._List[i].Conformity(other)
1547
            return self    
1548
        if isinstance(other, GroupSet):
1549
            if isinstance(other._List[0], GapInstance): # if not a Gap object, do nothing
1550
                Ggap = gap(other._List)
1551
                Ggap.OriginalGroup()
1552
                Lgap = gap(self._List)
1553
                L = Lgap.Conform()
1554
                self._List = list(L)
1555
            return self                
1556
    def __str__(self):
1557
        s = str(self._List)
1558
        s = "{" + s[1:-1] + "}"  # replace the square brackets (guaranteed to be on the outside) with curly braces.
1559
        s = ConvertIdentity(s)
1560
        return(s) 
1561
    def __repr__(self):
1562
        return str(self)
1563

1564
# this is a class wrapper that wraps around Gap elements, so that we can apply non-Gap operations on it.
1565

1566
class GapElement:
1567
    def __init__(self, x):
1568
        if isinstance(x, GapInstance):
1569
            self._x = x    
1570
        else:
1571
            return NotImplemented    
1572
    def __eq__(self, other):
1573
        if isinstance(other, GapElement):
1574
            return(self._x == other._x)
1575
#        if isinstance(other, GapInstance):  #Don't rely on this, since other == self would always give False
1576
#            return(self._x == other)
1577
        return NotImplemented
1578
    def __ne__(self, other):
1579
        result = self.__eq__(other)
1580
        if result is NotImplemented:
1581
            return result
1582
        return not result
1583
    def __cmp__(self, other):
1584
        if isinstance(other, GapElement):
1585
            return cmp(self._x, other._x) 
1586
        return NotImplemented   
1587
    def __mul__(self, other):
1588
        if isinstance(other, GapElement):
1589
           return GapElement(self._x * other._x)
1590
        if isinstance(other, GroupSet):
1591
           return GroupSet([self]) * other 
1592
        if isinstance(other, GapInstance):
1593
           return GapElement(self._x * other)   
1594
        if isinstance(other, (int, long, IntType)):
1595
           return GapElement(self._x * other)
1596
        return NotImplemented
1597
    def __rmul__(self, other):
1598
        if isinstance(other, (int, long, IntType)):
1599
            return GapElement(other * self._x)
1600
        # all other cases should have been covered by __mul__
1601
        return NotImplemented
1602
    def __add__(self, other):
1603
        if isinstance(other, GapElement):
1604
           return GapElement(self._x + other._x)
1605
        if isinstance(other, GroupSet):
1606
           return GroupSet([self]) + other 
1607
        return NotImplemented     
1608
    def __sub__(self, other):
1609
        if isinstance(other, GapElement):
1610
           return GapElement(self._x - other._x)
1611
        if isinstance(other, GroupSet):
1612
           return GroupSet([self]) + (-other) 
1613
        return NotImplemented         
1614
    def __neg__(self):
1615
        return GapElement(- self._x)       
1616
    def __pow__(self, exp):
1617
        if isinstance(exp, (int, long, IntType)):
1618
            return GapElement(self._x^exp) 
1619
        return NotImplemented    
1620
    def __str__(self):  
1621
        s = str(self._x) 
1622
        s = ConvertIdentity(s)  
1623
        return s
1624
    def __repr__(self):    # Make the wrapper invisible
1625
        return str(self)      
1626
        
1627

1628
         
1629
class Homomorph:
1630
    def __init__(self, domain, target):
1631
        # convert domain and target to lists, and strip off wrappers if present.
1632
        if isinstance(domain, GroupSet):
1633
            self.Domain = domain._List
1634
        elif isinstance(domain, list):
1635
            NewList = [] 
1636
            for item in domain: # Strip the wrappers off of the Gap elements
1637
                if isinstance(item, GapElement):
1638
                    NewList.append(item._x)
1639
                else:
1640
                    NewList.append(item)
1641
            self.Domain = NewList  
1642
        else:
1643
            print("Domain must be a group or list") 
1644
        if isinstance(target, GroupSet):
1645
            self.Target = target._List
1646
        elif isinstance(target, list):
1647
            NewList = [] 
1648
            for item in target: # Strip the wrappers off of the Gap elements
1649
                if isinstance(item, GapElement):
1650
                    NewList.append(item._x)
1651
                else:
1652
                    NewList.append(item)
1653
            self.Target = NewList  
1654
        else:
1655
            print("Target must be a group or list")
1656
            # TODO: we may eventually allow target to be a field
1657
        self.In = []  
1658
        self.Out = [] 
1659
        self.IsRing = false  # For RingHomo, this would be set to true after initialization            
1660
    def Set(self, input, output): ## Assumes the wrappers are stripped from input and output 
1661
        # Always conform both the input and output to the Domain and Target.
1662
        # For groups, this will put the elements in standard form (ListGroup format)
1663
        # For rings, it will allow us to enter 1 for (1 mod 6).
1664
        if input in self.In:
1665
            self.Out[self.In.index(input)] = self.Target[self.Target.index(output)] 
1666
        else:
1667
            self.In.append(self.Domain[self.Domain.index(input)])
1668
            self.Out.append(self.Target[self.Target.index(output)])            
1669
    def __call__(self, value):
1670
        valuex = value
1671
        if isinstance(value, GapElement):
1672
            valuex = value._x
1673
        try:
1674
            output = self.Out[self.In.index(valuex)] 
1675
            if isinstance(output, GapInstance):
1676
                return GapElement(output)
1677
            return output    
1678
        except ValueError:
1679
            print("Homomorphism is not defined at that value.")  
1680
    def Finish(self):
1681
        gen = copy(self.In)
1682
        grp = copy(self.In)
1683
        if isinstance(gen[0], GapInstance):  # Domain consists of Gap elements, so we use Gap to assist.
1684
            # First find the multiplication table of the domain
1685
            fDomain = gap(self.Domain)    #
1686
            M = fDomain.QuickTable() # will produce a 1 index matrix
1687
            if self.IsRing:
1688
                A = fDomain.QuickAdd() # will produce a 1 index matrix for addition
1689
            # convert the homomorphism to a translation dictionary, which uses an integer for the element
1690
            tran = [self.Domain.index(key) + 1 for key in gen]
1691
            gen = copy(tran)
1692
            grp = copy(tran)
1693
            for g1 in grp:
1694
                for g2 in gen:
1695
                    z = int(M[g1][g2])
1696
                    prod = self.Out[grp.index(g1)] * self.Out[grp.index(g2)] # multiplication is in the target group - so no tables
1697
                    if self.IsRing and isinstance(prod, GroupSet):    
1698
                        prod = prod + self.Out[0] + (-self.Out[0])   
1699
                    if z in grp:
1700
                        # z is already in the list, check that fun[g1]*fun[g2] = fun[g1*g2]
1701
                        if prod != self.Out[grp.index(z)]:
1702
                            print str(self.Out[grp.index(g1)]) + " * " + str(self.Out[grp.index(g2)] ) + " is not " + str(self.Out[grp.index(z)])
1703
                            return "Homomorphism failed"
1704
                    else:
1705
                        if not(prod in self.Target):
1706
                            print str(prod) + " is not in target."
1707
                            return "Homomorphism failed"
1708
                        prod = self.Target[self.Target.index(prod)]
1709
                        grp.append(z)   
1710
                        self.In.append(self.Domain[z-1])
1711
                        self.Out.append(prod)  
1712
                    if self.IsRing:
1713
                        z = int(A[g1][g2])
1714
                        if z in grp:
1715
                            # z is already in the list, check that fun[g1] + fun[g2] = fun[g1 + g2]  
1716
                            if self.Out[grp.index(g1)] + self.Out[grp.index(g2)] != self.Out[grp.index(z)]:
1717
                                print str(self.Out[grp.index(g1)]) + " + " + str(self.Out[grp.index(g2)] ) + " is not " + str(self.Out[grp.index(z)])
1718
                                return "Homomorphism failed" 
1719
                        else: 
1720
                            summ = self.Out[grp.index(g1)] + self.Out[grp.index(g2)]
1721
                            if not(summ in self.Target):
1722
                                print str(prod) + " is not in target."
1723
                                return "Homomorphism failed"        
1724
                            grp.append(z)   
1725
                            self.In.append(self.Domain[z-1])
1726
                            self.Out.append(summ)              
1727
        else:
1728
            for g1 in grp:
1729
                for g2 in gen:
1730
                    z = g1 * g2
1731
                    prod = self.Out[self.In.index(g1)] * self.Out[self.In.index(g2)]
1732
                    if self.IsRing and isinstance(z, GroupSet):
1733
                        z = z + g1 + (-g1) #completes the coset
1734
                    if self.IsRing and isinstance(prod, GroupSet):    
1735
                        prod = prod + self.Out[0] + (-self.Out[0])    
1736
                    if z in grp:
1737
                        # z is equivalent to an element in grp, but may have a different form
1738
                        # replace z with the form that is in grp.
1739
                        # z = grp[grp.index(z)]
1740
                        if prod != self.Out[self.In.index(z)]:
1741
                            print str(self.Out[self.In.index(g1)])+" * "+str(self.Out[self.In.index(g2)])+" is not "+str(self.Out[self.In.index(z)])
1742
                            return "Homomorphism failed"
1743
                    else:
1744
                        if not(prod in self.Target):
1745
                            print str(prod) + " is not in target."
1746
                            return "Homomorphism failed"
1747
                        prod = self.Target[self.Target.index(prod)]    
1748
                        grp.append(z) 
1749
                        self.In.append(z)
1750
                        self.Out.append(prod) 
1751
                    if self.IsRing:
1752
                        z = g1 + g2
1753
                        if z in grp:
1754
                            if self.Out[self.In.index(g1)] + self.Out[self.In.index(g2)] != self.Out[self.In.index(z)]:
1755
                                print str(self.Out[self.In.index(g1)])+" + "+str(self.Out[self.In.index(g2)])+" is not "+str(self.Out[self.In.index(z)])
1756
                                return "Homomorphism failed"
1757
                        else:
1758
                            summ = self.Out[self.In.index(g1)] + self.Out[self.In.index(g2)]
1759
                            if not(summ in self.Target):
1760
                                print str(summ) + " is not in target."
1761
                                return "Homomorphism failed"    
1762
                            grp.append(z) 
1763
                            self.In.append(z)
1764
                            self.Out.append(summ)                           
1765
        if len(self.In) == len(self.Domain):                
1766
            return "Homomorphism defined"
1767
        else:
1768
            return "Homomorphism consistent, but not defined for the whole domain."    
1769
    def GetDomain(self):
1770
        return self.Domain
1771
    def GetRange(self):
1772
        # This is what gets tricky.  We want to remove duplicates in self.Out, and put them in the order of self.Target
1773
        tempRange = []
1774
        for g1 in self.Target:
1775
            if (g1 in self.Out):
1776
                tempRange.append(g1)
1777
        return tempRange
1778
    def Inv(self, value):
1779
        # On input, value will be a list of elements (stripped of GapElement).
1780
        # returns a sorted list of elements so that Fun(x) = value
1781
        # value may be a list, and in fact, if the elements of the target are lists, add an extra list to value
1782
        n = len(self.Out)
1783
        R = []
1784
        for i in range(n):
1785
           if self.Out[i] in value:
1786
               if not self.In[i] in R:
1787
                   R.append(self.In[i])
1788
        return sorted(R)    
1789
          
1790
class FieldHomo:
1791
    def __init__(self):        
1792
        self.In = []  
1793
        self.Out = []
1794
        self.Key = [] 
1795
        self.RationalMap = false # assume it is for the current field extension.
1796
        if GenList == []:
1797
            self.RationalMap = true          
1798
    def Set(self, input, output): 
1799
        if input in self.In:
1800
            self.Out[self.In.index(input)] = output 
1801
        else:
1802
            self.In.append(input)
1803
            self.Out.append(output)     
1804
    def __call__(self, value):
1805
        x = value
1806
        if self.RationalMap:
1807
            x = TrigReduce(x)
1808
            vec = Compon(x, self.In)
1809
            total = vec[-1] # last component is the rational part
1810
            for i in range(len(self.In)):
1811
                total = total + vec[i]*self.Out[i]
1812
            return expand(total)
1813
        else:
1814
            if self.Key == []:    # CheckHomo was never done.
1815
                flag = self.Finish()
1816
                if not(flag):
1817
                    return None
1818
            try:
1819
                array = Vectorize(value)
1820
            except:
1821
                # if it did not vectorize, we must have entered a rational number
1822
                return value
1823
            n = len(array)
1824
            total = 0
1825
            for i in range(n):
1826
                total = total + array[i]*self.Key[i]    
1827
            return total      
1828
    def Finish(self): #Really CheckHomo
1829
        # This checks whether the field homomorphism works.
1830
        # In the case of a homomorphism on a field extension, it constructs 
1831
        # the vector that is used for the call function
1832
        if self.RationalMap:
1833
            n = len(self.In)
1834
            for i in range(n):
1835
                for j in range(n):
1836
                    z = TrigReduce(self.Out[i] * self.Out[j])
1837
                    w = TrigReduce(self.In[i] * self.In[j])
1838
                    if self(w) != z:
1839
                        return "f(" + str(self.In[i]) + ")*f(" + str(self.In[j]) + ") is not " + str(z) + "."
1840
            return true     
1841
        else: 
1842
            KeyList = []
1843
            n = len(GenList)
1844
            OutList = []
1845
            try:
1846
                for i in range(n):
1847
                    OutList.append(self.Out[self.In.index(GenList[i])])
1848
            except:
1849
                print "The mappings for all of the generators has not been defined."    
1850
                return false 
1851
            iList = [0 for i in DegreeList] # Generalized for loop
1852
            Cont = true
1853
            while Cont:
1854
                ele = 1
1855
                for i in range(n):
1856
                    ele = ele * OutList[i]^iList[i]
1857
                KeyList.append(ele)
1858
                # increment the iList
1859
                pointer = 0
1860
                iList[pointer] = iList[pointer] + 1
1861
                while iList[pointer] >= DegreeList[pointer]:
1862
                    iList[pointer] = 0
1863
                    pointer = pointer + 1
1864
                    if pointer == n:
1865
                        Cont = false
1866
                        pointer = 0       
1867
                    iList[pointer] = iList[pointer] + 1
1868
            self.Key = KeyList  
1869
            # Check to see if homomorphism is consistent
1870
            for i in range(n):
1871
                ele1 = GenList[i]^(DegreeList[i]-1)
1872
                ele2 = GenList[i]
1873
                ele3 = GenList[i]^DegreeList[i]
1874
                if self(ele1) * self(ele2) != self(ele3):
1875
                    print "Inconsistent definition."
1876
                    return false
1877
            return true      
1878
                          
1879
DisplayPermInt = false
1880
  
1881
class Perm:
1882
    def __init__(self, *args):
1883
        L = list(args)
1884
        if len(L) == 1: # allow for a single argument of a list
1885
            if isinstance(L, list):
1886
                L = L[0]
1887
        # test to see if this is the numbers from 1 to n
1888
        if sorted(L) != range(1,len(L)+1):
1889
            raise ValueError("Must be a rearrangement of the numbers 1 through n.")
1890
        else:    
1891
            # delete trailing number if in position
1892
            while (len(L)>0) and (L[-1] == len(L)):
1893
               L.pop() 
1894
            self._P = L
1895
    def __eq__(self, other):
1896
        if isinstance(other, Perm):
1897
            return self._P == other._P
1898
        else:
1899
            return NotImplemented
1900
    def __ne__(self, other):
1901
        result = self.__eq__(other)
1902
        if result is NotImplemented:
1903
            return result
1904
        return not result      
1905
    def __cmp__(self, other):
1906
        if isinstance(other, Perm):
1907
            if len(self._P) < len(other._P):
1908
                return -1
1909
            if len(self._P) > len(other._P):
1910
                return 1
1911
            n = len(self._P)
1912
            for i in range(n-1,0,-1):
1913
                if self._P[i] < other._P[i]:
1914
                    return 1
1915
                if self._P[i] > other._P[i]:
1916
                    return -1   
1917
            return 0                     
1918
        else:
1919
            return NotImplemented   
1920
    def __lt__(self, other):
1921
        if isinstance(other, Perm):
1922
            return (self.__cmp__(other) == -1)
1923
        else:
1924
            return NotImplemented    
1925
    def __gt__(self, other):
1926
        if isinstance(other, Perm):
1927
            return (self.__cmp__(other) == 1)
1928
        else:
1929
            return NotImplemented    
1930
    def __le__(self, other):
1931
        if isinstance(other, Perm):
1932
            return (self.__cmp__(other) < 1)
1933
        else:
1934
            return NotImplemented    
1935
    def __ge__(self, other):
1936
        if isinstance(other, Perm):
1937
            return (self.__cmp__(other) > -1)
1938
        else:
1939
            return NotImplemented     
1940
    def __call__(self, value):   
1941
        if value < 1:
1942
            return value
1943
        try:
1944
            return self._P[value-1]     
1945
        except: # probably out of range
1946
            return value    
1947
    def __mul__(self, other):
1948
        if isinstance(other, Perm):
1949
            L = copy(other._P)
1950
            M = copy(self._P)  
1951
            u = max(len(L), len(M))
1952
            # get both lists to be the same length
1953
            for i in range(len(L)+1, u+1):
1954
                L.append(i)
1955
            for i in range(len(M)+1, u+1):
1956
                M.append(i)        
1957
            ## The following line uses left to right multiplication, assuming (fog)(x) = g(f(x)). 
1958
            # T = [L[M[i]-1] for i in range(u)]
1959
            ## The following line uses right to left multiplication, assuming (fog)(x) = f(g(x)).
1960
            T = [M[L[i]-1] for i in range(u)]
1961
            return Perm(T)
1962
        elif isinstance(other, Cycle):
1963
            return self.ToCycle() * other
1964
        else:
1965
            return NotImplemented   
1966
    def __pow__(self, exp):
1967
        if isinstance(exp, (int, long, IntType)):
1968
            if exp == 0:
1969
                return Perm()
1970
            elif exp == 1:
1971
                return self
1972
            elif exp > 1:
1973
                return (self^(exp-1)) * self #Recursive definition
1974
            elif exp == -1:
1975
                L = self._P
1976
                u = len(L)
1977
                T = copy(L)
1978
                for i in range(u):
1979
                    T[L[i]-1] = i + 1
1980
                return Perm(T)
1981
            else: # exp < -1
1982
                return (self^(-1))^(-exp)  
1983
        else:
1984
            return NotImplemented      
1985
    def ToCycle(self):
1986
        T = copy(self._P)
1987
        n = len(T)
1988
        R = []
1989
        for i in range(n):
1990
            if T[i] > 0:
1991
                L = [i+1]
1992
                j = T[i]
1993
                T[i] = 0
1994
                while T[j-1] > 0:
1995
                    L.append(j)
1996
                    k = j
1997
                    j = T[k-1]
1998
                    T[k-1] = 0
1999
                if len(L) > 1:
2000
                    R.append(L)
2001
        ret = Cycle()
2002
        ret._C = R
2003
        return ret   
2004
    def ToInt(self):
2005
        total = 1
2006
        u = len(self._P)
2007
        for i in range(1,u):
2008
            for j in range(i):
2009
                if self._P[j] > self._P[i]:
2010
                    total = total + factorial(i)
2011
        return total         
2012
    def __str__(self):
2013
        if DisplayPermInt:
2014
            return str(self.ToInt())
2015
        return "P" + str(tuple(self._P))  
2016
    def __repr__(self):
2017
        return str(self)    
2018
           
2019
# This is a user defined error to handle the case where a non-integer cycle is converted to a permutation.             
2020
class CycleToPermError(Exception):
2021
    def __init__(self, value):
2022
        self.value = value
2023
    def __str__(self):
2024
        return repr(self.value)    
2025
    
2026
class Cycle:
2027
    def __init__(self, *args):
2028
        L = list(args)
2029
        if len(L) == 1: # allow for a single argument of a list
2030
            if isinstance(L, list):
2031
                L = L[0]    
2032
        # test to see if there are any duplicates in the list
2033
        if len(Uniquify(L)) != len(L):   
2034
            raise ValueError("Cannot have duplicates in cycle.")
2035
        elif len(L) == 0:
2036
            self._C = []   # Note that the identity element is stored as an empty list. 
2037
        else:
2038
            # Remove GapElement wrappers if necessary
2039
            NewL = []
2040
            for item in L:
2041
                if isinstance(item, GapElement):
2042
                    NewL.append(item._x)
2043
                else:
2044
                    NewL.append(item)    
2045
            # rotate list so that it begins with the smallest value
2046
            n = NewL.index(min(NewL))
2047
            for i in range(n):
2048
                NewL.append(NewL.pop(0))
2049
            self._C = [NewL]
2050
    def __eq__(self, other):
2051
        if isinstance(other, Cycle):
2052
            return self._C == other._C
2053
        else:
2054
            return NotImplemented            
2055
    def __ne__(self, other):
2056
        result = self.__eq__(other)
2057
        if result is NotImplemented:
2058
            return result
2059
        return not result  
2060
    def __lt__(self, other):
2061
        if isinstance(other, Cycle):
2062
            try:
2063
                return self.ToPerm() < other.ToPerm()
2064
            except CycleToPermError:
2065
                return self._C < other._C    
2066
        elif isinstance(other, Perm):
2067
            return self.ToPerm() < other
2068
        else:
2069
            return NotImplemented    
2070
    def __le__(self, other):
2071
        if isinstance(other, Cycle):
2072
            try:
2073
                return self.ToPerm() <= other.ToPerm()
2074
            except CycleToPermError:
2075
                return self._C <= other._C        
2076
        elif isinstance(other, Perm):
2077
            return self.ToPerm() <= other
2078
        else:
2079
            return NotImplemented   
2080
    def __gt__(self, other):
2081
        result = self.__le__(other)  
2082
        if result is NotImplemented:
2083
            return result
2084
        return not result  
2085
    def __ge__(self, other):
2086
        result = self.__lt__(other)  
2087
        if result is NotImplemented:
2088
            return result
2089
        return not result     
2090
    def __call__(self, value):
2091
        x = value
2092
        n = len(self._C)
2093
        for i in range(n):
2094
            if x in self._C[i]:
2095
                j = self._C[i].index(x) + 1
2096
                if j == len(self._C[i]):
2097
                    j = 0
2098
                x = self._C[i][j]
2099
        return x   
2100
    def __mul__(self, other):
2101
        if isinstance(other, Cycle):  
2102
            ## The following 2 lines uses left to right multiplication, assuming (fog)(x) = g(f(x)).  
2103
            # S = copy(self._C)
2104
            # S.extend(other._C) # combine the two lists into one list of lists
2105
            ## The following 2 lines uses right to left multiplication, assuming (fog)(x) = f(g(x)).  
2106
            S = copy(other._C)
2107
            S.extend(self._C) # combine the two lists into one list of lists
2108
            m = len(S)
2109
            T = Uniquify(Flatten(S)) # find all of the symbols in the list
2110
            T = sorted(T)
2111
            R = []
2112
            n = len(T)
2113
            for i in range(n):
2114
                if T[i] != None: # None is the indication that this symbol has been used
2115
                    L = []
2116
                    j = T[i]
2117
                    while j in T:
2118
                        L.append(j)
2119
                        T[T.index(j)] = None
2120
                        for k in range(m):  
2121
                            if j in S[k]:
2122
                                q = S[k].index(j)+1
2123
                                if q == len(S[k]):
2124
                                    q = 0
2125
                                j = S[k][q]
2126
                    if len(L) > 1:
2127
                        R.append(L)            
2128
            ret = Cycle()                
2129
            ret._C = R
2130
            return ret 
2131
        elif isinstance(other, Perm):
2132
            return self * other.ToCycle()
2133
        else:
2134
            return NotImplemented  
2135
    def __pow__(self, exp):
2136
        if isinstance(exp, (int, long, IntType)) :
2137
            if exp == 0:
2138
                return Cycle()
2139
            elif exp == 1:
2140
                return self
2141
            elif exp > 1:
2142
                return (self^(exp-1)) * self #Recursive definition
2143
            elif exp == -1:
2144
                R = []
2145
                for L in copy(self._C):
2146
                    LL = copy(L)
2147
                    LL.reverse()
2148
                    LL.insert(0, LL.pop())  # rotates last position to first
2149
                    R.append(LL)
2150
                ret = Cycle()
2151
                ret._C = R
2152
                return ret    
2153
            else: # exp < -1
2154
                return (self^(-1))^(-exp)  
2155
        else:
2156
            return NotImplemented              
2157
    def __str__(self):
2158
        if len(self._C) == 0:
2159
            return "()"
2160
        s = ""
2161
        for i in range(len(self._C)):
2162
            s = s + str(tuple(self._C[i]))
2163
        s = ConvertIdentity(s)    
2164
        return s       
2165
    def __repr__(self):
2166
        return str(self)
2167
    def ToPerm(self):
2168
        if len(self._C) == 0:
2169
            return Perm()
2170
        L = Flatten(self._C)
2171
        n = max(L)
2172
        m = min(L)
2173
        ## m should be a positive integer, and n should also be an integer, or else return an error.
2174
        if isinstance(n, (int, long, IntType)) and isinstance(m, (int, long, IntType)) and m > 0:
2175
            PP = range(1,n+1) 
2176
            for R in self._C:
2177
                m = len(R)
2178
                for i in range(m-1):
2179
                    PP[R[i]-1] = R[i+1]
2180
                PP[R[-1]-1] = R[0]    
2181
            return Perm(PP)     
2182
        else:
2183
            raise CycleToPermError(m)
2184

2185
   
2186
class Twople:       
2187
    def __init__(self, xvalue, yvalue):
2188
        # Strip away wrappers if present
2189
        if isinstance(xvalue, GapElement):    
2190
            self._x = xvalue._x
2191
        else:
2192
            self._x = xvalue
2193
        if isinstance(yvalue, GapElement): 
2194
            self._y = yvalue._x
2195
        else:    
2196
            self._y = yvalue
2197
    def __eq__(self, other):
2198
        if isinstance(other, Twople):      
2199
            return (self._x == other._x) and (self._y == other._y)
2200
        else:
2201
            return NotImplemented            
2202
    def __ne__(self, other):
2203
        result = self.__eq__(other)
2204
        if result is NotImplemented:
2205
            return result
2206
        return not result  
2207
    def __mul__(self, other):
2208
        if isinstance(other, Twople):              
2209
            return Twople( self._x * other._x, self._y * other._y )
2210
        else:
2211
            return NotImplemented  
2212
    def __pow__(self, exp):
2213
        if isinstance(exp, (int, long, IntType)) :  
2214
            return Twople(self._x^exp, self._y^exp)
2215
        else:
2216
            return NotImplemented   
2217
    def __str__(self):
2218
        return "(" + str(self._x) + ", " + str(self._y) + ")"      
2219
    def __repr__(self):
2220
        return "(" + str(self._x) + ", " + str(self._y) + ")"  
2221

2222
class Basis:
2223
    def __init__(self, *args):
2224
        if len(args) == 1:   
2225
            BasList = [1]
2226
            VecList = args[0]
2227
        else:
2228
            BasList = args[0]
2229
            VecList = args[1]   
2230
        # BasList and VecList will be lists, we need to vectorize each combination of BasList * VecList.
2231
        self._n = len(BasList)
2232
        self._m = len(VecList)
2233
        self._bas = BasList
2234
        self._worked = false
2235
        dimen = self._n * self._m
2236
        M = []
2237
        for j in range(self._m):
2238
            for i in range(self._n):
2239
                v = BasList[i]*VecList[j]
2240
                w = Vectorize(v)
2241
                if not isinstance(w, list):
2242
                    w = [w]
2243
                # pad w with zeros to make it of length dimen
2244
                while len(w) < dimen:
2245
                   w.append(0)
2246
                M.append(w)
2247
        Mat = matrix(M) # convert to a matrix
2248
        try:
2249
           Inv = ~Mat   # take the inverse matrix
2250
        except:
2251
           return None
2252
        self.inv = Inv  
2253
        self._worked = true
2254
        return None  
2255
    def Coeff(self, value):
2256
        if not(self._worked):
2257
            print "basis is not linearly independent"
2258
            return false
2259
        dimen = self._n * self._m    
2260
        w = Vectorize(value)
2261
        if not isinstance(w, list):
2262
            w = [w]
2263
        # pad w with zeros to make it of length dimen
2264
        while len(w) < dimen:
2265
            w.append(0)
2266
        T = matrix(w) * self.inv # Use inverse matrix to find new vector
2267
        # T will be a row matrix, not a list, so we have to convert it to a list.
2268
        L = []
2269
        for j in range(self._m):
2270
            tot = 0
2271
            for i in range(self._n):
2272
                k = T[0][i + j*self._n]
2273
                tot = tot + k*self._bas[i]
2274
            L.append(tot)    
2275
        return L  
2276
        
2277
##############################################################################
2278
## 
2279
##  Main group defining routines
2280
##
2281
##############################################################################
2282
    
2283
def ZGroup(n):
2284
    v = []
2285
    for i in range(n):
2286
        v.append(SumModN(i,n))
2287
    return GroupSet(v)    
2288
        
2289
def ZStar(n):
2290
    v = []
2291
    R = Integers(n) # Built in Zn ring
2292
    for i in range(1,n):
2293
        if gcd(i,n) == 1:
2294
            v.append(R(i))
2295
    return GroupSet(v)
2296
    
2297
def ZRing(n):
2298
    v = []
2299
    R = Integers(n) # Built in Zn ring
2300
    for i in range(n):
2301
        v.append(R(i))
2302
    return GroupSet(v)    
2303
    
2304
def InitQuaternions():
2305
    global i
2306
    global j
2307
    global k
2308
#  By defining Quaternions with SR, (Symbolic Ring), we allow for constructions like a + b*i + c*j + d*k
2309
    Q.<i,j,k> = QuaternionAlgebra(SR,-1,-1)
2310
    return GroupSet([1, i, j, k, -1, -i, -j, -k])    
2311

2312
def InitSkew9():
2313
    global a
2314
    global b
2315
#   This defines an unusual skew field of dimension 9 over the rational numbers.  It is defined by
2316
#   a^3 = 3a+1, b^3 = 2, and b*a =(2-a^2)*b    
2317
    A = FreeAlgebra(SR, 2, "x")
2318
    F = A.monoid()
2319
    a = F.gens()[0]
2320
    b = F.gens()[1]
2321
    mons = [F(1), a, a^2, b, a*b, a^2*b, b^2, a*b^2, a^2*b^2]
2322
    M = MatrixSpace(QQ, 9)
2323
    mats = [M([0,1,0, 0,0,0,0,0,0,  
2324
               0,0,1, 0, 0, 0, 0, 0, 0,  
2325
               1,3,0, 0, 0, 0, 0, 0, 0,   
2326
               0,0,0, 2, 0,-1, 0, 0, 0,  
2327
               0,0,0,-1,-1, 0, 0, 0, 0,  
2328
               0,0,0, 0,-1,-1, 0, 0, 0,
2329
               0,0,0, 0, 0, 0,-2,-1, 1,  
2330
               0,0,0, 0, 0, 0, 1, 1,-1,   
2331
               0,0,0, 0 ,0 ,0,-1,-2, 1]), 
2332
            M([0,0,0,1,0,0,0,0,0,   
2333
               0,0,0,0,1,0,0,0,0,  
2334
               0,0,0,0,0,1,0,0,0,
2335
               0,0,0,0,0,0,1,0,0,   
2336
               0,0,0,0,0,0,0,1,0,  
2337
               0,0,0,0,0,0,0,0,1,
2338
               2,0,0,0,0,0,0,0,0,   
2339
               0,2,0,0,0,0,0,0,0,  
2340
               0,0,2,0,0,0,0,0,0])]
2341
    Skew9.<a,b> = A.quotient(mons, mats)  
2342
    return GroupSet([a, b])
2343

2344
## The main groups are interfaced with Gap.  These routines set up a quotient group of a FreeGroup in Gap
2345

2346
IdentityElement = "e"
2347
GeneratorList = []
2348
RelationsList = []
2349
gap.eval("TempFreeGroup := FreeGroup(0)")
2350
gap.eval("CurrentGroup := TempFreeGroup/[]")
2351
CurrentGroup = gap("CurrentGroup")
2352

2353
def InitGroup(name_of_identity_element):
2354
    global IdentityElement
2355
    global GeneratorList
2356
    global RelationsList
2357
    global CurrentGroup
2358
    global LastInit
2359
    IdentityElement = name_of_identity_element
2360
    GeneratorList = []
2361
    RelationsList = []
2362
    LastInit = "InitGroup"
2363
    gap.eval("TempFreeGroup := FreeGroup(0)")
2364
    gap.eval("CurrentGroup := TempFreeGroup/[]")
2365
    CurrentGroup = gap("CurrentGroup")
2366
    gap.eval(IdentityElement + ":= Identity(CurrentGroup);")
2367
    tmpstr = "global " + IdentityElement + "; " + IdentityElement + " = GapElement(CurrentGroup.Identity())"
2368
    exec(tmpstr)
2369
    return "New group with identity " + name_of_identity_element
2370
    
2371
def AddGroupVar(*args):
2372
    global GeneratorList
2373
    global CurrentGroup
2374
    if len(args) == 0:
2375
        print("Must have at least one argument")
2376
    else:
2377
        for newvar in args:
2378
            GeneratorList.append(newvar)
2379
        tmpG = len(GeneratorList)
2380
        tmpstr = "TempFreeGroup := FreeGroup(["
2381
        for tmpI in range(tmpG-1):
2382
            tmpstr = tmpstr + '"' + GeneratorList[tmpI] + '", '
2383
        tmpstr = tmpstr + '"' + GeneratorList[tmpG-1] + '"])'   
2384
        gap.eval(tmpstr)
2385
        gap.eval("AssignGeneratorVariables(TempFreeGroup)")
2386
        tmpR = len(RelationsList)
2387
        if tmpR == 0:    
2388
            gap.eval("CurrentGroup := TempFreeGroup/[]")
2389
        else:
2390
            tmpstr = "CurrentGroup := TempFreeGroup/["  
2391
            for tmpI in range(tmpR-1):
2392
                tmpstr = tmpstr + RelationsList[tmpI] + ", " 
2393
            tmpstr = tmpstr + RelationsList[tmpR-1] + "];"
2394
            gap.eval(tmpstr)
2395
        gap.eval("AssignGeneratorVariables(CurrentGroup)")
2396
        CurrentGroup = gap("CurrentGroup")
2397
        gap.eval(IdentityElement + ":= Identity(CurrentGroup);")
2398
        tmpstr = "global " + IdentityElement + "; " + IdentityElement + " = GapElement(CurrentGroup.Identity())"
2399
        exec(tmpstr)
2400
        for tmpI in range(tmpG):
2401
            tmpstr = "global " + GeneratorList[tmpI] + "; " + GeneratorList[tmpI] + ' = GapElement(gap("' + GeneratorList[tmpI] + '"))'       
2402
            exec(tmpstr)
2403
    return None    
2404

2405
def Define(element1, element2):
2406
    global LastFieldVar
2407
    global GenList
2408
    global DegreeList
2409
    if LastInit == "InitGroup":
2410
        # when defining a Gap group, element1 and element2 will probably have the GapElement wrapper.  
2411
        if isinstance(element1, GapElement):
2412
            elem1 = element1._x
2413
        else:
2414
            elem1 = element1    
2415
        if isinstance(element2, GapElement):
2416
            elem2 = element2._x
2417
        else:
2418
            elem2 = element2     
2419
        global RelationsList
2420
        global CurrentGroup
2421
        if (elem1 in CurrentGroup) and (elem2 in CurrentGroup):
2422
            tmpstr = str(elem1 * (elem2)^(-1))
2423
            RelationsList.append(tmpstr)
2424
            tmpG = len(GeneratorList)
2425
            tmpstr = "TempFreeGroup := FreeGroup(["
2426
            for tmpI in range(tmpG-1):
2427
                tmpstr = tmpstr + '"' + GeneratorList[tmpI] + '", '
2428
            tmpstr = tmpstr + '"' + GeneratorList[tmpG-1] + '"])'   
2429
            gap.eval(tmpstr)
2430
            gap.eval("AssignGeneratorVariables(TempFreeGroup)")
2431
            tmpR = len(RelationsList)
2432
            tmpstr = "CurrentGroup := TempFreeGroup/["  
2433
            for tmpI in range(tmpR-1):
2434
                tmpstr = tmpstr + RelationsList[tmpI] + ", " 
2435
            tmpstr = tmpstr + RelationsList[tmpR-1] + "];"
2436
            gap.eval(tmpstr)
2437
            gap.eval("AssignGeneratorVariables(CurrentGroup)")
2438
            CurrentGroup = gap("CurrentGroup")
2439
            gap.eval(IdentityElement + ":= Identity(CurrentGroup);")
2440
            tmpstr = "global " + IdentityElement + "; " + IdentityElement + " = GapElement(CurrentGroup.Identity())"
2441
            exec(tmpstr)
2442
            for tmpI in range(tmpG):
2443
                tmpstr = "global " + GeneratorList[tmpI] + "; " + GeneratorList[tmpI] + ' = GapElement(gap("' + GeneratorList[tmpI] + '"))'      
2444
                exec(tmpstr)
2445
        return None   
2446
    elif LastInit == "InitDomain":
2447
        if LastFieldVar == "":
2448
            return "New variable has not been defined yet."
2449
        # we create a polynomial by subtracting the two elements
2450
        DefPoly = element1 - element2
2451
        polystr = preparse(str(DefPoly))
2452
        n = DefPoly.degree()
2453
        if BaseCharacteristic > 0 and isinstance(CurrentField, GFType):
2454
            # This is the first extension on a finite field.  
2455
            # In order to enable the factoring algorithm, we will construct 
2456
            # a Galois field using this polynomial. 
2457
            # Note that this polynomial MUST be irreducible over GF(n)
2458
            # for this to work.
2459
            m = BaseCharacteristic^n # size of new field
2460
            tmpstr = "global CurrentField; global " + LastFieldVar + "; CurrentField = GF(" + str(m) + ', name = "' 
2461
            tmpstr = tmpstr + LastFieldVar + '", modulus = ' + polystr + "); " + LastFieldVar + " = CurrentField._first_ngens(1)[0]"
2462
            exec(tmpstr)   
2463
        elif BaseCharacteristic == 0:
2464
            # extension of an infinite field
2465
            tmpstr = "global CurrentField; global " + LastFieldVar + "; CurrentField = CurrentField.extension(" + polystr + ", names=('" 
2466
            tmpstr = tmpstr + LastFieldVar + "',)); " + LastFieldVar + " = CurrentField._first_ngens(1)[0]"
2467
            exec(tmpstr)   
2468
        elif len(GenList) > 1:
2469
            return "Can't do extension of an extension of an extension of a finite field!"   
2470
        elif DegreeList[0] == 0:
2471
            # We are doing an extension of a rational function field of finite characteristic.  Proceed as an infinite field.
2472
            tmpstr = "global CurrentField; global " + LastFieldVar + "; CurrentField = CurrentField.extension(" + polystr + ", names=('" 
2473
            tmpstr = tmpstr + LastFieldVar + "',)); " + LastFieldVar + " = CurrentField._first_ngens(1)[0]"
2474
            exec(tmpstr) 
2475
        else:
2476
            # We are trying to extend an extension to a finite field.  We must first convert the GF definition of CurrentGroup to
2477
            # an extension of Z_n.
2478
            # Note: factorization will no longer be possible for this extension of an extension of a finite characteristic.
2479
            FirstVar = CurrentField.variable_names()[0]
2480
            OldPolyStr = str(CurrentField.modulus())
2481
            OldPolyStr = OldPolyStr.replace("x", FirstVar) # x may not be defined, so convert to the variable that will be defined
2482
            OldPolyStr = preparse(OldPolyStr)
2483
            # reset FirstVar to be the variable of GF(n), rather than the field.
2484
            tmpstr = "global " + FirstVar + '; TempDomain = GF(' + str(BaseCharacteristic) + ')["' + FirstVar + '"]; ' 
2485
            tmpstr = tmpstr + FirstVar + " = TempDomain._first_ngens(1)[0]"
2486
            exec(tmpstr)   
2487
            # now we redefine the field in terms of an extension of GF(p), rather than GF(p^n).
2488
            tmpstr = "global CurrentField; global " + FirstVar + "; CurrentField = GF(" + str(BaseCharacteristic) + ").extension(" 
2489
            tmpstr = tmpstr + OldPolyStr + ", names=('" 
2490
            tmpstr = tmpstr + FirstVar + "',)); " + FirstVar + " = CurrentField._first_ngens(1)[0]"
2491
            exec(tmpstr)               
2492
            # We still have to redefine the new variable to be a variable of this new field.
2493
            tmpstr = "global " + LastFieldVar + '; TempDomain = CurrentField["' + LastFieldVar + '"]; ' 
2494
            tmpstr = tmpstr + LastFieldVar + " = TempDomain._first_ngens(1)[0]"
2495
            exec(tmpstr)      
2496
            # We are now ready to define the new field as an extension of this extension
2497
            tmpstr = "global CurrentField; global " + LastFieldVar + "; CurrentField = CurrentField.extension(" 
2498
            tmpstr = tmpstr + polystr + ", names=('" 
2499
            tmpstr = tmpstr + LastFieldVar + "',)); " + LastFieldVar + " = CurrentField._first_ngens(1)[0]"
2500
            exec(tmpstr) 
2501
            # We are almost done.  Unfortunately, the FirstVar element is not set to an element of the new field, but of
2502
            # the intermediate field.  This causes a problem if we try to multiply FirstVar with LastFieldVar, giving
2503
            # a NonImplemented error.  So we need to set FirstVar to the element in the new field corresponding the the first
2504
            # generator.  This is tricky, since CurrentField.gens() is missing this information.  The only approach is a Monte Carlo method!
2505
            tmpstr = "global TempFirstVar; TempFirstVar = " + FirstVar
2506
            exec(tmpstr)
2507
            tmpstr = "global TempLastVar; TempLastVar = " + LastFieldVar
2508
            exec(tmpstr)
2509
            NewFirstVar = CurrentField.random_element()*TempLastVar
2510
            while str(NewFirstVar) != str(TempFirstVar):
2511
                NewFirstVar = CurrentField.random_element()*TempLastVar
2512
            tmpstr = "global " + FirstVar + "; " + FirstVar + " = NewFirstVar"
2513
            exec(tmpstr)  
2514
            # reset the first var of GenList so that it works.
2515
            GenList[0] = NewFirstVar         
2516
        GenList.append(CurrentField.gen())
2517
        DegreeList.append(n)
2518
        LastFieldVar = "" # So we cannot immediately do another extension, but first define the next variable.
2519
        if UnivFieldVar != "":
2520
            tmpstr = "global " + UnivFieldVar + '; TempDomain = CurrentField["' + UnivFieldVar + '"]; ' 
2521
            tmpstr = tmpstr + UnivFieldVar + " = TempDomain._first_ngens(1)[0]"
2522
            exec(tmpstr)          
2523
    else:
2524
        return "Nothing is initialized yet."  
2525
  
2526
def ToPolycyclic():
2527
    global GeneratorList
2528
    global CurrentGroup
2529
    try:
2530
        gap.eval("CurrentPCGroup := PcGroupFpGroup(CurrentGroup)")
2531
    except:
2532
        print "Group definition doesn't meet all the polycyclic requirements."
2533
    else:
2534
        gap.eval("AssignGeneratorVariables(CurrentPCGroup)")
2535
        print "Group converted to the Polycyclic format."
2536
        tmpG = len(GeneratorList)
2537
        CurrentGroup = gap("CurrentPCGroup")
2538
        gap.eval(IdentityElement + ":= Identity(CurrentPCGroup);")
2539
        tmpstr = "global " + IdentityElement + "; " + IdentityElement + " = GapElement(CurrentGroup.Identity())"
2540
        exec(tmpstr)
2541
        for tmpI in range(tmpG):
2542
            tmpstr = "global " + GeneratorList[tmpI] + "; " + GeneratorList[tmpI] + ' = GapElement(gap("' + GeneratorList[tmpI] + '"))'      
2543
            exec(tmpstr)
2544
        x = gap([]) 
2545
        L = x.ListCurrentPCGroup()   
2546
        return GroupSet(list(L))   
2547
    
2548
def GroupStringConvert(gapstr):
2549
    # StructureDescription uses GAP to identify a group, and returns a string with the name of the group.
2550
    # Unfortunately, GAP's convention is slightly different than the textbook.  GAP used C5 for Z_5, and D8 for D_4 
2551
    # This routine converts the string of the group name to the name that the user will be more familiar with.
2552
    str2 = ''
2553
    L = len(gapstr)
2554
    i = 1 # this will be the pointer
2555
    while i < L:
2556
        char = gapstr[i]
2557
        if char != '"':
2558
            # remove the quotation marks
2559
            if char == 'C':
2560
                # a 'C' followed by a number should be replaced with a 'Z'
2561
                if (gapstr[i+1] >= '0' and gapstr[i+1] <= '9'):
2562
                    str2 += 'Z'
2563
                else:
2564
                    str2 += char
2565
            elif (char == 'D' and gapstr[i-1] != 'Q'):
2566
                # a 'D' followed by a number should have the number halved.  Not applied to QD16
2567
                str2 += char
2568
                strnum = ""
2569
                j = i+1
2570
                while (gapstr[j] >= '0' and gapstr[j] <= '9'): # no fear of out of range, since the last character is '"'.
2571
                    strnum += gapstr[j]
2572
                    j = j+1
2573
                    i = i+1
2574
                str2 = str2 + str(int(strnum)/2)
2575
            else:
2576
                str2 += char 
2577
        i = i+1 
2578
    return(str2)   
2579

2580
def StructureDescription(*args):
2581
    # This uses the small groups package from GAP to identify the group CurrentGroup.  This should have been installed the first time this
2582
    # package is run.
2583
    try:
2584
        if len(args) == 0:
2585
            gapstr = gap.eval("StructureDescription(CurrentGroup)")
2586
        else:
2587
            gapstr = "Perms := Group(" + str(PermToCycle(NthPerm(args[0])))
2588
            for i in range(1, len(args)):
2589
                gapstr = gapstr + "," + str(PermToCycle(NthPerm(args[i])))
2590
            gapstr = gapstr + ")"   
2591
            gap.eval(gapstr)
2592
            gapstr = gap.eval("StructureDescription(Perms)")      
2593
    except:
2594
        print 'Error in finding the structure of the group.'
2595
    else:
2596
        print GroupStringConvert(gapstr)  
2597
        
2598
def GaloisType(poly):
2599
    # This uses GAP to find the Galois group of the polynomial.  It assumes that the polynomial has been defined in QQ[x], for 
2600
    # some variable x
2601
    varname = poly.variable_name()
2602
    deg = poly.degree()
2603
    nstr = gap.eval(varname + ' := Indeterminate(Rationals, "' + varname + '");; GaloisType(' + str(poly) + ')')
2604
    outstr = gap.eval("StructureDescription(TransitiveGroup(" + str(deg) + "," + nstr + "))")
2605
    # Although TransitiveGroup gives more information than StructureDescription, students would not understand the notation produced by
2606
    # the output of TransitiveGroup.  Thus, we find the StructureDescription to make it easier for the students.
2607
    print GroupStringConvert(outstr)     
2608
        
2609
        
2610
def InitTerry():
2611
    global IdentityElement
2612
    global GeneratorList
2613
    global RelationsList
2614
    global CurrentGroup
2615
    global Stay
2616
    global RotLft
2617
    global RotRt
2618
    global Spin
2619
    global FlipRt
2620
    global FlipLft
2621
    gap.eval('ff := FreeGroup("RotLft","RotRt","Spin","FlipLft","FlipRt")')
2622
    gap.eval('CurrentGroup := ff/[ff.1^3,ff.2^3,ff.1*ff.2,ff.3^2,ff.4^2,ff.5^2,ff.5*ff.2*ff.4,ff.3*ff.2*ff.5]')
2623
    gap.eval('AssignGeneratorVariables(CurrentGroup)')
2624
    gap.eval('SetReducedMultiplication(CurrentGroup)')
2625
    CurrentGroup = gap('CurrentGroup')
2626
    Stay = GapElement(gap('Identity(CurrentGroup)'))
2627
    RotLft = GapElement(gap('RotLft'))
2628
    RotRt = GapElement(gap('RotRt'))
2629
    Spin = GapElement(gap('Spin'))
2630
    FlipLft = GapElement(gap('FlipLft'))
2631
    FlipRt = GapElement(gap('FlipRt'))
2632
    IdentityElement = 'Stay'
2633
    return GroupSet([Stay, FlipRt, RotRt, FlipLft, RotLft, Spin])
2634
 
2635
###################################################################################
2636
##
2637
##  Once the Gap group is defined, there are several ways for it to be displayed
2638
##
2639
###################################################################################
2640
 
2641
def ListGroup():  # this will only be used for groups defined in GAP
2642
    str = gap.eval("Size(CurrentGroup)")
2643
    if str == 'infinity':
2644
        return 'Group is infinite.'
2645
    gap.eval("GList := ListGroup(CurrentGroup)")
2646
    GList = gap('GList')
2647
    return GroupSet(list(GList))
2648
    
2649
def ListRing():   # this will only be used for rings defined by InitRing
2650
    str = gap.eval("Size(CurrentRing)")
2651
    if str == 'infinity':
2652
        return 'Ring is infinite.'
2653
    gap.eval("RList := ListRing(CurrentRing)")
2654
    RList = gap('RList')
2655
    return GroupSet(list(RList))
2656
    
2657
def SetReducedMult():  # this will only be used for groups defined in GAP
2658
    gap.eval("SetReducedMultiplication(CurrentGroup)")
2659
    return None            
2660
                                                         
2661
def Conform(S_GroupSet, G_GroupSet):
2662
    if (isinstance(S_GroupSet, GroupSet) and isinstance(G_GroupSet, GroupSet)):
2663
        return S_GroupSet.Conformity(G_GroupSet) 
2664
   
2665
def Group(*args):
2666
    global gen
2667
    global grp
2668
    if len(args) == 0:  # If no arguments, do the same thing as ListGroup.
2669
        return ListGroup()    
2670
    else:
2671
        # We can assume that all arguments will be of the same type, so we can strip the GapElement wrappers
2672
        v = []
2673
        for arg in args:
2674
            if isinstance(arg, GapElement):
2675
                v.append(arg._x)
2676
            else:    
2677
                v.append(arg)
2678
        if isinstance(v[0], GapInstance):  # Delicate the job to GAP
2679
            L = gap(v)    #
2680
            Subgroup = L.Group()
2681
            GList = Subgroup.List()
2682
            return GroupSet(list(GList))
2683
        else:  
2684
            gen = Uniquify(v) # removes duplicates
2685
            grp = sorted(gen)   # sorts the list
2686
            for g1 in grp:
2687
                for g2 in gen:
2688
                    z = g1 * g2                  
2689
                    if not (z in grp):
2690
                        grp.append(z)
2691
            return GroupSet(sorted(grp)) # TODO: determine sorting algorithm                
2692
                               
2693
##################################################################################################
2694
##
2695
##  Graphical display routines
2696
##
2697
################################################################################################ 
2698

2699
FilterCommas = false
2700

2701
#
2702
#  This converts a GAP object (or similar object) to a string.  Spaces are removed, and "<identity>" is replaced with
2703
#  the current value of IdentityElement.  Mainly this is used for multiplication tables, but it can also be used for circle graphs.
2704
#
2705

2706
def ObjectToStr(obj):
2707
    global FilterCommas
2708
    ss = str(parent(obj))
2709
    if isinstance(obj, (GapInstance, list)):  # might have <identity...> that needs to be converted
2710
        ss = str(obj)
2711
        ss = ConvertIdentity(ss)
2712
        ss = ss.replace(" ","")   # Remove the spaces
2713
        if (FilterCommas == true):
2714
            ss = ss.replace(",","")   
2715
    else:                         # Easy case: do a simple convertion
2716
        ss = str(obj) 
2717
        ss = ss.replace(" ","")   # Remove the spaces    
2718
    return ss                     # TODO:  If obj is a list, it might be a list of Gap objects, or a list of list of Gap objects, etc.
2719
    
2720
def Add(x):
2721
    L = ["Add"]
2722
    L.append(x)
2723
    return(L)
2724

2725
def Mult(x):
2726
    L = ["LeftMult"]    
2727
    L.append(x)
2728
    return(L)
2729

2730
def LeftMult(x):
2731
    L = ["LeftMult"]    
2732
    L.append(x)
2733
    return(L)
2734
        
2735
def RightMult(x):
2736
    L = ["RightMult"]    
2737
    L.append(x)
2738
    return(L)
2739
    
2740
def Pow(x):
2741
    L = ["Pow"]    
2742
    L.append(x)
2743
    return(L)
2744
    
2745
Inv = ["Inv"]    
2746
    
2747
def CircleGraph(V, *args): # each args will either be a permutation, an automorphism, or an ordered pair of the form [string, element]
2748
    L = []                 # to complicate matters, there may have several arguments 
2749
    if isinstance(V, GroupSet):
2750
        ListV = V._List
2751
    else:    
2752
        TempList = list(V)   # If a list was entered, the elements may be wrapped in GapElement.
2753
        ListV = []           # In this case, we must strip off the wrapper.
2754
        for item in TempList:
2755
            if isinstance(item, GapElement):
2756
                ListV.append(item._x)
2757
            elif isinstance(item, str):
2758
                ListV.append(eval(item))
2759
            else:    
2760
                ListV.append(item)    
2761
    n = len(ListV)           
2762
    for i in range(n):
2763
       L.append(ObjectToStr(V[i]))   
2764
    pos_dict = {} # set up the dots in a circle, going clockwise from the top
2765
    for i in range(n):
2766
        x = float(-cos(pi/2 + ((2*pi)/n)*i))
2767
        y = float(sin(pi/2 + ((2*pi)/n)*i))
2768
        pos_dict[L[i]] = [x,y]              
2769
    g = DiGraph(loops = true, multiedges = true)
2770
    # find the maximum length of the strings in L (at least 2)
2771
    maxString = 2
2772
    for i in range(n):
2773
        if len(L[i])> maxString:
2774
           maxString = len(L[i])
2775
    question = false
2776
    currentcolor = 4
2777
    for op in args:
2778
        edges_covered = []
2779
        while len(edges_covered) < n:
2780
            # find first element not already covered
2781
            k = 0
2782
            while k in edges_covered:
2783
                k = k + 1
2784
            currentpoint = ListV[k]
2785
            while ((k > -1) and (not (k in edges_covered))):  # This must run at least once, since we searched for a point not in edges_covered
2786
                currentpoint = ListV[k]
2787
                edges_covered.append(k)
2788
                # Find the next point, based upon the operation in op
2789
                if isinstance(op, list):
2790
                    if op[0] == "Add":
2791
                        if isinstance(op[1], GapElement):
2792
                            nextpoint = currentpoint + op[1]._x
2793
                        else:  
2794
                            nextpoint = currentpoint + op[1]
2795
                    elif op[0] == "LeftMult":
2796
                        if isinstance(op[1], GapElement):
2797
                            nextpoint = currentpoint * op[1]._x
2798
                        else:  
2799
                            nextpoint = currentpoint * op[1]
2800
                    elif op[0] == "RightMult":
2801
                        # This is the one case that can be a headache
2802
                        # if V was a list of cosets, then currentpoint will be
2803
                        #    a GroupSet, otherwise it would be a GapElement
2804
                        if isinstance(op[1], GapElement):
2805
                            if isinstance(currentpoint, GroupSet):
2806
                                nextpoint = op[1] * currentpoint
2807
                            else:                                
2808
                                nextpoint = op[1]._x * currentpoint
2809
                        else:    
2810
                            nextpoint = op[1] * currentpoint
2811
                    elif op[0] == "Pow":
2812
                        nextpoint = currentpoint^(op[1])
2813
                    elif op[0] == "Inv":
2814
                        try:
2815
                            nextpoint = currentpoint^(-1)
2816
                        except ZeroDivisionError:
2817
                            nextpoint = "?"    
2818
                    else:    
2819
                        nextpoint = currentpoint # Catches exceptions
2820
                elif isinstance(op, (Cycle, Perm, FieldHomo)):
2821
                    nextpoint = op(currentpoint)
2822
                elif isinstance(op, Homomorph):
2823
                    nextpoint = op(currentpoint) 
2824
                    # the homomorphism might put the answer in a wrapper, which we must take off
2825
                    if isinstance(nextpoint, GapElement):
2826
                        nextpoint = nextpoint._x   
2827
                else:        
2828
                    # TODO: determine other types of operations that could happen
2829
                    nextpoint = "?"
2830
                # Find the position of the nextpoint in the list
2831
                if (nextpoint in ListV):  #  Add this edge to the graph 
2832
                    ll = ListV.index(nextpoint)
2833
                    g.add_edge((L[k], L[ll], currentcolor))
2834
                    k = ll
2835
                else:
2836
                    if(question == false):  # add a "?" to the position list
2837
                        pos_dict["?"] = [0,0]
2838
                        question = true
2839
                    g.add_edge((L[k], "?", currentcolor))
2840
                    k = -1  # get out of this while loop
2841
                currentpoint = nextpoint
2842
            if len(args) == 1:    # change the color after each cycle, only if there is only one function    
2843
                currentcolor = currentcolor + 1
2844
        currentcolor = currentcolor + 1  # update for each function       
2845
        vertsize = 70*(maxString + 1)      
2846
    return g.plot(pos = pos_dict, edge_colors = g._color_by_label(ColorTable), vertex_size = vertsize,  
2847
                                                                                     xmin = -1.025, xmax = 1.025, ymin = -1.05, ymax = 1.0)
2848

2849
QuotientRing = false
2850

2851
def MultTable(v):
2852
    Graff = Graphics()
2853
    if isinstance(v, GroupSet):
2854
        GG = v._List
2855
    else:    
2856
        TempList = list(v)   # If a list was entered, the elements may be wrapped in GapElement.
2857
        GG = []              # In this case, we must strip off the wrapper.
2858
        for item in TempList:
2859
            if isinstance(item, GapElement):
2860
                GG.append(item._x)
2861
            else:
2862
                GG.append(item)    
2863
    n = len(GG)
2864
    if(n > 28):
2865
        return "Group too big to display table."
2866
    L = []
2867
    for i in range(n):
2868
        L.append(ObjectToStr(GG[i]))
2869
    if isinstance(GG[0], GapInstance):
2870
        GG = gap(GG) # Convert list to a GAP object
2871
        M = GG.QuickTable() # will produce a 1 index matrix
2872
#    elif isinstance(GG[0], (int, long, IntType)) and isinstance(GG[-1], (int, long, IntType)) and (not InitPermMultiplication) and (min(GG) > -1):
2873
#        # for a list of integers, try to determine the operation
2874
#        base = max(GG) + 1
2875
#        M = matrix(n+1)
2876
#        if 0 in GG:
2877
#            # assume group is addition modulo base    
2878
#            for i in range(n):
2879
#                for j in range(n):
2880
#                    w = (GG[i] + GG[j]) % base
2881
#                    if w in GG:
2882
#                        M[i+1, j+1] = GG.index(w)+1
2883
#        else:    
2884
#            # assume group is multiplication modulo base    
2885
#            for i in range(n):
2886
#                for j in range(n):
2887
#                    w = (GG[i] * GG[j]) % base
2888
#                   if w in GG:
2889
#                        M[i+1, j+1] = GG.index(w)+1            
2890
    else: #defalt case - use multiplication
2891
#        GFlat = []
2892
#        if isinstance(GG[0], list):  # GG is a list of lists
2893
#            GFlat = Flatten(GG)
2894
        M = matrix(n+1)
2895
        for i in range(n):
2896
            for j in range(n):
2897
                w = GG[j] * GG[i]
2898
                if w in GG:
2899
                    M[j+1, i+1] = GG.index(w)+1
2900
                else:
2901
                    # second chance - w might be a _subset_ of an element in GG
2902
                    if isinstance(w, GroupSet):
2903
                        for k in range(n):
2904
                            if IsSubset(w._List, GG[k]):
2905
                                if QuotientRing:
2906
                                    M[j+1, i+1] = k + 1
2907
                                else:
2908
                                    M[j+1, i+1] = -(k + 1)
2909
                                    Graff += text(w,(i+1.5, n-j-0.5), rgbcolor=(0,0,0))
2910
                                break     
2911
    Graff += line([(1,0),(1,n),(n+1,n)], rgbcolor=(0,0,0))
2912
    # find the maximum length of the strings in L
2913
    maxString = 0
2914
    for i in range(n):
2915
        if len(L[i])> maxString:
2916
           maxString = len(L[i])
2917
    for i in range(n):
2918
        Graff += polygon([(i+1, n+1),(i+2, n+1),(i+2, n),(i+1, n)], color = ColorTable[i+1])
2919
        Graff += polygon([(0, n-i),(1, n-i),(1, n-i-1),(0, n-i-1)], color = ColorTable[i+1])
2920
        if (maxString*(n+1) < 64):
2921
            Graff += text(L[i],(i+1.5, n+0.5), rgbcolor=(0,0,0))
2922
            Graff += text(L[i],(0.5, n-i-0.5), rgbcolor=(0,0,0))
2923
        else:
2924
            Graff += text(L[i],(-0.01*(n+1),n-i-0.5), rgbcolor=(0,0,0), horizontal_alignment='right')   
2925
    for i in range(n):    
2926
        for j in range(n):
2927
            p = int(M[j+1][i+1]) # note: have to convert the GAP integer to a sage integer
2928
            Graff += polygon([(i+1, n-j),(i+2, n-j),(i+2, n-j-1),(i+1, n-j-1)], color = ColorTable[abs(p)])
2929
            if ((p > 0) and (maxString*(n+1) < 64)):
2930
               Graff += text(L[p-1],(i+1.5, n-j-0.5), rgbcolor=(0,0,0))
2931
    return Graff.show(axes = false)
2932
    
2933
def AddTable(v):
2934
    Graff = Graphics()
2935
    if isinstance(v, GroupSet):
2936
        GG = v._List
2937
    else:   
2938
        TempList = list(v)   # If a list was entered, the elements may be wrapped in GapElement.
2939
        GG = []              # In this case, we must strip off the wrapper.
2940
        for item in TempList:
2941
            if isinstance(item, GapElement):
2942
                GG.append(item._x)
2943
            else:
2944
                GG.append(item)   
2945
    n = len(GG)
2946
    if(n > 28):
2947
        return "Ring too big to display table."
2948
    L = []
2949
    for i in range(n):
2950
        L.append(ObjectToStr(GG[i]))
2951
    s = str(parent(GG[0]))
2952
    if s == 'Gap':  #Delegate the task to GAP
2953
        GG = gap(GG)  # Convert list to a GAP object
2954
        M = GG.QuickAdd() # will produce a 1 index matrix    
2955
    else: # probably Ring of integers modulo n   
2956
        M = matrix(n+1)
2957
        for i in range(n):
2958
            for j in range(n):
2959
                w = GG[i] + GG[j]
2960
                if w in v:
2961
                    M[i+1, j+1] = GG.index(w)+1
2962
    Graff += line([(1,0),(1,n),(n+1,n)], rgbcolor=(0,0,0))
2963
    # find the maximum length of the strings in L
2964
    maxString = 0
2965
    for i in range(n):
2966
        if len(L[i])> maxString:
2967
           maxString = len(L[i])
2968
    for i in range(n):
2969
        Graff += polygon([(i+1, n+1),(i+2, n+1),(i+2, n),(i+1, n)], color = ColorTable[i+1])
2970
        Graff += polygon([(0, n-i),(1, n-i),(1, n-i-1),(0, n-i-1)], color = ColorTable[i+1])
2971
        if (maxString*(n+1) < 64):
2972
            Graff += text(L[i],(i+1.5, n+0.5), rgbcolor=(0,0,0))
2973
            Graff += text(L[i],(0.5, n-i-0.5), rgbcolor=(0,0,0))
2974
        else:
2975
            Graff += text(L[i],(-0.01*(n+1),n-i-0.5), rgbcolor=(0,0,0), horizontal_alignment='right')   
2976
    for i in range(n):    
2977
        for j in range(n):
2978
            p = int(M[j+1][i+1]) # note: have to convert the GAP integer to a sage integer
2979
            Graff += polygon([(i+1, n-j),(i+2, n-j),(i+2, n-j-1),(i+1, n-j-1)], color = ColorTable[p])
2980
            if ((p > 0) and (maxString*(n+1) < 64)):
2981
               Graff += text(L[p-1],(i+1.5, n-j-0.5), rgbcolor=(0,0,0))
2982
    return Graff.show(axes = false)
2983
        
2984

2985
def GraphHomo(fun):
2986
    global maxString
2987
    global Lrange
2988
    global Ldomain
2989
    global pos_dict
2990
    # This graphs the function, using the notation given in the list domain.    
2991
    domain = fun.GetDomain() 
2992
    fRange = fun.GetRange()  
2993
    # Convert domain and range to strings
2994
    m = len(domain)
2995
    n = len(fRange)
2996
    Ldomain = []
2997
    Lrange = []
2998
    for i in range(m):
2999
        Ldomain.append(ObjectToStr(domain[i]))  
3000
    for i in range(n):
3001
        Lrange.append(ObjectToStr(fRange[i])) 
3002
    graff = Graphics()     
3003
    for domnum in range(m):      
3004
        key = domain[domnum]
3005
        rannum = fRange.index(fun.Out[fun.In.index(key)]) 
3006
        if m == 1:
3007
            x = 0.5
3008
        else:
3009
            x =  1 - (numerical_approx(domnum)/(m-1))    
3010
        if n == 1:
3011
            y = 0.5
3012
        else:
3013
            y =  1 - (numerical_approx(rannum)/(n-1))    
3014
        graff += arrow((0,x), (1,y), color = ColorTable[4+rannum])
3015
    if m == 1:
3016
        graff += point((0, 0.5), rgbcolor = (0,0,0), size = 50)
3017
        graff += text(Ldomain[0],(-0.02, 0.5), rgbcolor=(0,0,0), horizontal_alignment='right') 
3018
    else:   
3019
        for i in range(m):
3020
            y = 1 - (numerical_approx(i)/(m-1))
3021
            graff += point((0, y), rgbcolor = (0,0,0), size = 50)
3022
            graff += text(Ldomain[i],(-0.02, y), rgbcolor=(0,0,0), horizontal_alignment='right') 
3023
    if n == 1:
3024
        graff += point((1, 0.5), rgbcolor = (0,0,0), size = 50)
3025
        graff += text(Lrange[0],(1.02, 0.5), rgbcolor=(0,0,0), horizontal_alignment='left') 
3026
    else:   
3027
        for i in range(n):
3028
            y = 1 - (numerical_approx(i)/(n-1))
3029
            graff+= point((1, y), rgbcolor = (0,0,0), size = 50)   
3030
            graff += text(Lrange[i],(1.02, y), rgbcolor=(0,0,0), horizontal_alignment='left')   
3031
    return graff.show(axes = false)       
3032
 
3033
#######################################################################################
3034
##
3035
##   RSA routines
3036
##
3037
#######################################################################################    
3038
 
3039
def NextPrime(n):
3040
    if n%2 == 0:
3041
       n = n+1
3042
    while not(is_prime(n)):
3043
       n = n+2
3044
    return n
3045
   
3046
def CentToAscii(x):
3047
    if x < 1:
3048
        y = 32
3049
    elif x < 30:
3050
        y = x + 64
3051
    elif x < 47:
3052
        y = x + 18
3053
    elif x < 48:
3054
        y = 33
3055
    elif x < 81:
3056
        y = x + 46
3057
    elif x < 95:
3058
        y = x - 47
3059
    elif x == 95:
3060
        y = 196
3061
    elif x == 96:
3062
        y = 203
3063
    elif x == 97:
3064
        y = 214
3065
    elif x == 98:
3066
        y = 220
3067
    elif x == 99:
3068
        y = 223
3069
    else:
3070
        y = 0
3071
    return y
3072

3073
def AsciiToCent(x):
3074
    if x < 33:
3075
        y = 0
3076
    elif x < 34:
3077
        y = 47
3078
    elif x < 48:
3079
        y = x + 47
3080
    elif x < 65:
3081
        y = x - 18
3082
    elif x < 94:
3083
        y = x - 64
3084
    elif x < 127:
3085
        y = x - 46
3086
    else:
3087
        y = 0
3088
    return y
3089

3090
def MessageToNumber(s):
3091
    total = 0
3092
    for i in range(len(s)):
3093
        total = total * 100 + AsciiToCent(ord(s[i]))
3094
    return total
3095

3096
def NumberToMessage(n):
3097
    strlist = ""
3098
    Temp = int(n)
3099
    while Temp > 0:
3100
        m = Temp % 100
3101
        Temp = Temp // 100
3102
        strlist = chr(CentToAscii(m)) + strlist
3103
    return strlist    
3104
 
3105
######################################################################################       
3106
   
3107
def C(*args):      
3108
    L = [];
3109
    for arg in args:  # Convert args to a list, which will support indexing
3110
        L.append(arg)
3111
    return Cycle(L)   
3112

3113
def P(*args):      
3114
    L = [];
3115
    for arg in args:  # Convert args to a list, which will support indexing
3116
        L.append(arg)
3117
    return Perm(L)
3118

3119
def Order(x):
3120
    if isinstance(x, GapElement):
3121
        xx = x._x
3122
    else:
3123
        xx = x    
3124
    if isinstance(xx, GapInstance):    # Delicate the job to GAP
3125
        n = xx.Group().Size()
3126
        if n == gap('infinity'):
3127
            return Infinity
3128
        else:
3129
            return int(n)
3130
    elif isinstance(xx, SumModN):
3131
        return xx._mod/GCD(x._x, x._mod)    
3132
    else:
3133
        # since we do not know what type of element it is, we do not know the identity element.
3134
        # Our strategy is to find an n > 1 for which x^n = x, and return n-1.
3135
        y = x*x
3136
        n = 1
3137
        while x != y:
3138
           y = y * x
3139
           n = n + 1
3140
        return n              
3141
            
3142
def RootCount(G, k):
3143
    # Counts the number of elements for which g^k = e.
3144
    count = 0;
3145
    for x in G:
3146
       t = x
3147
       for i in range(k):
3148
          t = t*x
3149
       if t == x:
3150
          count = count + 1   
3151
    return count      
3152
            
3153
def Generators(G):
3154
    # finds all of the generators of the group, assuming it is cyclic.
3155
    L = [];
3156
    n = len(G)   
3157
    for x in G:
3158
        if (Order(x) == n):
3159
            L.append(x)
3160
    return L                     
3161
                
3162
def InitRing(*args):
3163
    # actually, all this does is var's the varables in *args, and forms a list of these variables in GeneratorsOfRing.  
3164
    # But this sets up the ability to form a ring from using the orders, and the mini-multiplication table.
3165
    global GeneratorsOfRing
3166
    GeneratorsOfRing = []
3167
    for arg in args:
3168
        var(arg)
3169
        tmpstr = 'global ' + arg
3170
        exec(tmpstr)
3171
        tmpstr = 'GeneratorsOfRing.append(' + arg + ')'
3172
        exec(tmpstr)
3173
    return None   
3174
       
3175
def DefineRing(CharR, GenTable):
3176
    # CharR is an integer vector showing the orders on the basis elements.  
3177
    # GenTable is a mini-multiplication table showing the products of two basis elements, written as a linear combination of the basis elements. 
3178
    global CurrentRing  
3179
    TempT = []
3180
    n = len(GeneratorsOfRing)  
3181
    for i in range(n):
3182
        TempT1 = []
3183
        for j in range(n):
3184
            TempVars = []
3185
            TempCk = []
3186
            for k in range(n):
3187
                m = SR(GenTable[i][j]).coefficient(GeneratorsOfRing[k],1)
3188
                m = int(m)
3189
                if m != 0:
3190
                    TempVars.append(k+1) #GAP is 1 indexed
3191
                    TempCk.append(m)
3192
            TempT1.append([TempVars, TempCk])
3193
        TempT.append(TempT1)
3194
    TempT.append(0) # Neither symmetric or anti-symmetric is assumed.
3195
    TempT.append(0) # Add zero of the defining ring   
3196
    tmpstr = 'CurrentRing := RingByStructureConstants(' + str(CharR) + ',' + str(TempT) + ',['
3197
    for i in range(n):
3198
       tmpstr = tmpstr + '"' + str(GeneratorsOfRing[i]) + '"'
3199
       if i < n-1:
3200
           tmpstr = tmpstr + ','
3201
    tmpstr = tmpstr + '])'       
3202
    gap.eval(tmpstr)
3203
    for i in range(n):
3204
        tmpstr = 'global ' + str(GeneratorsOfRing[i]) + '; ' + str(GeneratorsOfRing[i]) + ' = GapElement(gap("CurrentRing.' + str(i+1) + '"))' 
3205
        exec(tmpstr)   
3206
    CurrentRing = gap("CurrentRing") 
3207
    return None   
3208
     
3209
def NumberSmallRings(size):
3210
    return int(gap.eval("NumberSmallRings(" + str(size) + ")"))
3211

3212
def SmallRing(size, idnum):
3213
    gap.eval("CurrentRing := SmallRing(" + str(size) + "," + str(idnum) + ")")
3214
    CurrentRing = gap('CurrentRing')
3215
    L = CurrentRing.GeneratorsOfRing()  # This will be 1 indexed from GAP
3216
    n = len(L)
3217
    for i in range(n):
3218
        tmpstr = 'global ' + str(L[i+1]) + '; ' + str(L[i+1]) + ' = GapElement(gap("CurrentRing.' + str(i+1) + '"))' 
3219
        exec(tmpstr)
3220
    gap.eval("RList := ListRing(CurrentRing)")
3221
    RList = gap('RList')
3222
    return GroupSet(list(RList))
3223
         
3224
def CheckRing():
3225
    s = gap.eval("CheckRing(CurrentRing)")
3226
    ss = str(s)
3227
    ss = ss.replace('"','')
3228
    print ss
3229
    return None
3230
  
3231
def FindUnity(*args):
3232
    if len(args) == 0: 
3233
        Rring = gap(CurrentRing) 
3234
        ident = Rring.Identity()
3235
        if str(ident) == 'fail':
3236
            print "No identity element found."
3237
            return None
3238
        return GapElement(ident)
3239
    R_list = args[0]
3240
    # R must be a list or GroupSet for this to work.
3241
    # Returns the identity element, if there is one
3242
    if isinstance(R_list, list):
3243
        R_ring = GroupSet(R_list)
3244
    else:    
3245
        R_ring = R_list
3246
    if isinstance(R_ring._List[0], GapInstance):   #delegate the task to Gap
3247
        Rset = gap(R_ring._List)    
3248
        Rring = Rset.Ring()
3249
        ident = Rring.Identity()
3250
        if str(ident) == 'fail':
3251
            print "No identity element found."
3252
            return None
3253
        return GapElement(ident)
3254
    for i in R_ring._List:
3255
        TestFlag = true
3256
        for j in R_ring._List:
3257
            if ((i * j != j) or (j * i != j)):
3258
                TestFlag = false
3259
        if TestFlag:
3260
            return i
3261
    print "No identity element found."
3262
    return None         
3263

3264
def Ring(*args):
3265
    global gen
3266
    global grp
3267
    if len(args) == 0:  # If no arguments, do the same thing as ListRing.
3268
        return ListRing()    
3269
    else:
3270
        # We can assume that all arguments will be of the same type, so we can strip the GapElement wrappers
3271
        v = []
3272
        for arg in args:
3273
            if isinstance(arg, GapElement):
3274
                v.append(arg._x)
3275
            else:    
3276
                v.append(arg)
3277
        if isinstance(v[0], GapInstance):  # Delicate the job to GAP
3278
            L = gap(v)    #
3279
            Subring = L.Ring()
3280
            GList = Subring.List()
3281
            return GroupSet(list(GList))
3282
        else:  
3283
            gen = Uniquify(v) # removes duplicates
3284
            grp = sorted(gen)   # sorts the list
3285
            for g1 in grp:
3286
                for g2 in gen:
3287
                    z = g1 * g2                  
3288
                    if not (z in grp):
3289
                        grp.append(z)
3290
                    z = g1 + g2
3291
                    if not (z in grp):
3292
                        grp.append(z)    
3293
            return GroupSet(sorted(grp))     
3294
            
3295
def LftCoset(G_list, H_list):   
3296
    # both G and H must be lists or Cosets for this to work
3297
    # if they are lists, convert to Cosets (which will remove wrappers if present)
3298
    # We can assume that H is a subgroup of G
3299
    # Note:  The cosets in the output MUST be sorted in order to allow for MultTable to work
3300
    if isinstance(G_list, list):
3301
        G_coset = GroupSet(G_list)
3302
    else:
3303
        G_coset = G_list
3304
    if isinstance(H_list, list):
3305
        H_coset = GroupSet(H_list)
3306
    else:
3307
        H_coset = H_list        
3308
    if not (isinstance(G_coset, GroupSet) and isinstance(H_coset, GroupSet)):
3309
        return "Error: both parameters must be lists or subgroups" 
3310
    H_coset.Sort() 
3311
    if isinstance(G_coset._List[0], GapInstance):   # May delegate the task to Gap
3312
        # If G is a list of GapElements, then H will be too.  We can delegate the task to GAP
3313
        Ggap = gap(G_coset._List)
3314
        Ggap.OriginalGroup()    
3315
        Hgap = gap(H_coset._List)
3316
        OutGap = Hgap.LftCoset()
3317
        ## OutGap will be a gap list of lists.  We need to convert to a GroupSet of GroupSets
3318
        q = []
3319
        OutList = list(OutGap)
3320
        for item in OutList:
3321
            NewCoset = GroupSet(list(item))
3322
            q.append(NewCoset)
3323
        return GroupSet(q) 
3324
    # We can save a little time by starting with the original subgroup H   
3325
    q = [H_coset]
3326
    flat = H_coset._List
3327
    for g in G_coset:
3328
        if not(g in flat):
3329
            newCoset = GroupSet([g]) * H_coset
3330
            q.append(newCoset)  
3331
            flat = Flatten(q)
3332
    return GroupSet(q)       
3333
 
3334
def RtCoset(G_list, H_list):   
3335
    # both G and H must be lists or GroupSets for this to work
3336
    # if they are lists, convert to GroupSets (which will remove wrappers if present)
3337
    # We can assume that H is a subgroup of G
3338
    # Note:  The cosets in the output MUST be sorted in order to allow for MultTable to work
3339
    if isinstance(G_list, list):
3340
        G_coset = GroupSet(G_list)
3341
    else:
3342
        G_coset = G_list
3343
    if isinstance(H_list, list):
3344
        H_coset = GroupSet(H_list)
3345
    else:
3346
        H_coset = H_list        
3347
    if not (isinstance(G_coset, GroupSet) and isinstance(H_coset, GroupSet)):
3348
        return "Error: both parameters must be lists or subgroups" 
3349
    H_coset.Sort()
3350
    if isinstance(G_coset._List[0], GapInstance):   # May delegate the task to Gap
3351
        # If G is a list of GapElements, then H will be too.  We can delegate the task to GAP
3352
        Ggap = gap(G_coset._List)
3353
        Ggap.OriginalGroup()    
3354
        Hgap = gap(H_coset._List)
3355
        OutGap = Hgap.RtCoset()
3356
        ## OutGap will be a gap list of lists.  We need to convert to a GroupSet of GroupSets
3357
        q = []
3358
        OutList = list(OutGap)
3359
        for item in OutList:
3360
            NewCoset = GroupSet(list(item))
3361
            q.append(NewCoset)
3362
        return GroupSet(q) 
3363
    # We can save a little time by starting with the original subgroup H   
3364
    q = [H_coset]
3365
    flat = H_coset._List
3366
    for g in G_coset:
3367
        if not(g in flat):
3368
            newCoset = H_coset * GroupSet([g])
3369
            q.append(newCoset)  
3370
            flat = Flatten(q)
3371
    return GroupSet(q)  
3372
 
3373
def Coset(R_list, H_list):   
3374
    # both R and H must be lists or GroupSets for this to work
3375
    # if they are lists, convert to GroupSets (which will remove wrappers if present)
3376
    # We can assume that H is a subring of R
3377
    # Note:  The cosets in the output MUST be sorted in order to allow for MultTable to work
3378
    if isinstance(R_list, list):
3379
        R_coset = GroupSet(R_list)
3380
    else:
3381
        R_coset = R_list
3382
    if isinstance(H_list, list):
3383
        H_coset = GroupSet(H_list)
3384
    else:
3385
        H_coset = H_list        
3386
    if not (isinstance(R_coset, GroupSet) and isinstance(H_coset, GroupSet)):
3387
        return "Error: both parameters must be lists or subrings" 
3388
    H_coset.Sort() 
3389
    if isinstance(R_coset._List[0], GapInstance):   # May delegate the task to Gap
3390
        # If R is a list of GapElements, then H will be too.  We can delegate the task to GAP
3391
        Rgap = gap(R_coset._List)
3392
        Rgap.OriginalGroup()    
3393
        Hgap = gap(H_coset._List)
3394
        OutGap = Hgap.RingCoset()
3395
        ## OutGap will be a gap list of lists.  We need to convert to a GroupSet of GroupSets
3396
        q = []
3397
        OutList = list(OutGap)
3398
        for item in OutList:
3399
            NewCoset = GroupSet(list(item))
3400
            q.append(NewCoset)
3401
        return GroupSet(q) 
3402
    # We can save a little time by starting with the original subgroup H   
3403
    q = [H_coset]
3404
    flat = H_coset._List
3405
    for g in G_coset:
3406
        if not(g in flat):
3407
            newCoset = GroupSet([g]) + H_coset
3408
            q.append(newCoset)  
3409
            flat = Flatten(q)
3410
    return GroupSet(q)        
3411
        
3412
def RingHomo(domain, target):
3413
    funct = Homomorph(domain, target)
3414
    funct.IsRing = True
3415
    return funct        
3416
                
3417
def HomoDef(Fun, input, output):
3418
    if isinstance(Fun, FieldHomo): # this case, there is nothing to check
3419
        Fun.Set(input, output)
3420
    elif isinstance(Fun, Homomorph):    
3421
        ## we have to remove the wrappers from input and output
3422
        ## Also, if input or output is a list, convert to a GroupSet - either domain or target is a quotient group.
3423
        inputx = input
3424
        if isinstance(input, GapElement):
3425
            inputx = input._x
3426
        if isinstance(input, list):
3427
            inputx = GroupSet(input)   
3428
        outputx = output
3429
        if isinstance(output, GapElement):
3430
            outputx = output._x   
3431
        if isinstance(output, list):
3432
            outputx = GroupSet(output)     
3433
        elif not(inputx in Fun.Domain):
3434
            print "Second argument must be in the domain group." 
3435
        elif not(outputx in Fun.Target):
3436
            print "Last argument must be in the target group."
3437
        else:    
3438
            Fun.Set(inputx, outputx)       
3439
    else:
3440
        print "First argument must be a homomorphism."
3441
        
3442
def FinishHomo(Fun):
3443
    return Fun.Finish()   
3444

3445
def CheckHomo(Fun):
3446
    return Fun.Finish()     
3447
    
3448
def Image(Fun, L):
3449
    # L might be a GroupSet, a list, or a single element.  
3450
    # If it is a GroupSet, we can safely assume that we want the corresponding list.
3451
    # If it is a list, then we must strip the GapElement wrappers.
3452
    if isinstance(L, list):
3453
        LL = []
3454
        for item in L:
3455
            if isinstance(item, GapElement):
3456
                LL.append(item._x)
3457
            else:
3458
                LL.append(item)
3459
    elif isinstance(L, GroupSet):
3460
        LL = L._List
3461
    elif isinstance(L, GapElement):
3462
        LL = [L._x]    
3463
    else:    
3464
        return Fun(L)
3465
    R = []
3466
    for g1 in LL:
3467
        z = Fun.Out[Fun.In.index(g1)]
3468
        if not z in R:
3469
            R.append(z)
3470
    return GroupSet(sorted(R))
3471

3472
def HomoInv(Fun, L):
3473
    # L might be a GroupSet, a list, or a single element.  
3474
    # If it is a GroupSet, we can safely assume that we want the corresponding list.
3475
    # If it is a list, then we must strip the GapElement wrappers.
3476
    if isinstance(L, list):
3477
        LL = []
3478
        for item in L:
3479
            if isinstance(item, GapElement):
3480
                LL.append(item._x)
3481
            else:
3482
                LL.append(item)
3483
    elif isinstance(L, GroupSet):
3484
        LL = L._List
3485
    elif isinstance(L, GapElement):
3486
        LL = [L._x]    
3487
    else:    
3488
        LL = [L] 
3489
    return GroupSet(Fun.Inv(LL))
3490

3491
def HomoRange(Fun):
3492
    return Fun.GetRange()
3493
    
3494
def Kernel(Fun):
3495
    R = Fun.Target
3496
    if Fun.IsRing:
3497
        Q = [R[0] + (-R[0])]
3498
    else:    
3499
        Q = [R[0]^0]                   
3500
    return GroupSet(Fun.Inv(Q))        
3501
           
3502
def CycleToPerm(Cyc):
3503
    if isinstance(Cyc, Cycle):
3504
        return Cyc.ToPerm()
3505
    else:
3506
        return NotImplemented    
3507

3508
def PermToCycle(P_Perm):
3509
    if isinstance(P_Perm, Perm):
3510
        return P_Perm.ToCycle()
3511
    else:
3512
        return NotImplemented  
3513
  
3514
def PermToInt(P_Perm):
3515
    if isinstance(P_Perm, Perm):
3516
        return P_Perm.ToInt()
3517
    else:
3518
        return NotImplemented  
3519
      
3520
def Signature(C_Cycle):
3521
    if isinstance(C_Cycle, Cycle):
3522
        count = 1
3523
        for item in C_Cycle._C:
3524
            if len(item) %2 == 0:
3525
                count = - count
3526
        return count
3527
    elif isinstance(C_Cycle, Perm):
3528
        return Signature(C_Cycle.ToCycle()) # Probably a faster way to do this, but this will work.
3529
    else:
3530
        return NotImplemented
3531
        
3532
def Unfactorial(n):
3533
    # returns the smallest number for which i! >= n
3534
    index = 0
3535
    Temp = 1
3536
    while Temp < n:
3537
        index = index + 1
3538
        Temp = Temp * index
3539
    return index
3540
   
3541
def NthPerm(n_int):
3542
   global Temp
3543
   global digit
3544
   if n <= 0:
3545
       return "Argument must be a positive integer"
3546
   u = Unfactorial(n_int)
3547
   Temp = range(u,0,-1) # goes from [u, u-1, u-2, ... , 1]
3548
   S = []
3549
   for k in range(u, 0, -1):
3550
       digit = (((n_int-1) % factorial(k)) // factorial(k-1))
3551
       S.append(Temp[digit])
3552
       del Temp[digit]
3553
   S.reverse()
3554
   return Perm(S)
3555
    
3556
        
3557
def DirectProduct(HList, KList):
3558
    if isinstance(HList, (list, GroupSet)) and isinstance(KList, (list, GroupSet)):
3559
        L  = []
3560
        for i in HList:
3561
            for j in KList:
3562
                L.append(Twople(i,j))
3563
        return GroupSet(L)
3564
    else:
3565
        return NotImplemented   
3566
        
3567
def GroupCenter(Glist):
3568
    # This finds the elements of Glist that commute with all of the elements of Glist
3569
    # If Glist is a group or subgroup, the result we automatically be a subgroup.
3570
    # We will keep order of the elements the same as the order in Glist
3571
    # Glist might either be a GroupSet or a list.  If it is a list, we may have to unwrap the elements.
3572
    if isinstance(Glist, GroupSet):
3573
        G = Glist
3574
    else:
3575
        G = GroupSet(Glist)
3576
    g = []
3577
    for i in G._List:
3578
        include = true
3579
        for j in G._List:
3580
            if (i*j) != (j*i):
3581
                include = false # i is not in the center
3582
                break
3583
        if include:
3584
            g.append(i)
3585
    return GroupSet(g)
3586
   
3587
def Normalizer(Glist, Hlist):
3588
    # This finds the elements of Glist for which g.h.g^-1 is in Hlist for all h in Hlist.
3589
    # If Glist is a group or subgroup, the result we automatically be a subgroup.
3590
    # We will keep order of the elements the same as the order in Glist
3591
    # Glist might either be a GroupSet or a list.  If it is a list, we may have to unwrap the elements.
3592
    if isinstance(Glist, GroupSet):
3593
        G = Glist
3594
    else:
3595
        G = GroupSet(Glist)
3596
    if isinstance(Hlist, GroupSet):
3597
        H = Hlist
3598
    elif isinstance(Hlist, list):  
3599
        H = GroupSet(Hlist)
3600
    else:
3601
        H = GroupSet([Hlist]) # in case H is a single element  
3602
    answer = []
3603
    for g in G._List:
3604
        include = true
3605
        for h in H._List:
3606
            if not((g * h) * g^(-1) in H._List):
3607
                include = false # i is not in the normalizer
3608
                break
3609
        if include:
3610
            answer.append(g)
3611
    return GroupSet(answer)
3612

3613
def NormalClosure(Glist, Slist):
3614
    # This finds the smallest normal subgroup of Glist that contains the elements of Slist.  
3615
    if isinstance(Glist, GroupSet):
3616
        G = Glist
3617
    else:
3618
        G = GroupSet(Glist)
3619
    if isinstance(Slist, GroupSet):
3620
        S = Slist
3621
    elif isinstance(Slist, list):  
3622
        S = GroupSet(Slist)
3623
    else:
3624
        S = GroupSet([Slist]) # in case S is a single element 
3625
    # if the elements are GAP elements, deligate the take to GAP 
3626
    if isinstance(G._List[0], GapInstance):
3627
        Ggap = gap(G._List)
3628
        Ggap.OriginalGroup()  
3629
        Sgap = gap(S._List)
3630
        Ngap = Sgap.NormClosure()
3631
        return GroupSet(list(Ngap)) 
3632
    # find all conjugates of the elements in S
3633
    gen = []
3634
    n = len(G._List) # size of original group
3635
    for g in G._List:
3636
        for h in S._List:
3637
            u = (g * h) * g^(-1)
3638
            if not(u in gen):
3639
                gen.append(u)
3640
    # gen is now the set of generators for the group
3641
    grp = copy(gen) 
3642
    for g1 in grp:
3643
        for g2 in gen:
3644
            u = g1 * g2                  
3645
            if not (u in grp):
3646
                grp.append(u)
3647
                if len(grp)*2 > n:
3648
                    return G
3649
    return GroupSet(sorted(grp)) # TODO: determine sorting algorithm       
3650

3651
def Ideal(Rlist, Slist):
3652
    # This finds the smallest ideal of Rlist that contains the elements of Slist.  
3653
    if isinstance(Rlist, GroupSet):
3654
        R = Rlist
3655
    else:
3656
        R = GroupSet(Rlist)
3657
    if isinstance(Slist, GroupSet):
3658
        S = Slist
3659
    elif isinstance(Slist, list):  
3660
        S = GroupSet(Slist)
3661
    else:
3662
        S = GroupSet([Slist]) # in case S is a single element 
3663
    # if the elements are GAP elements, deligate the take to GAP 
3664
    if isinstance(R._List[0], GapInstance):
3665
        Rgap = gap(R._List)
3666
        Rgap.OriginalGroup()  
3667
        Sgap = gap(S._List)
3668
        Igap = Sgap.IdealClosure()
3669
        return GroupSet(list(Igap)) 
3670
    # first find the union of S, R.S, S.R, and R.S.R
3671
    gen = (R * S) * R
3672
    for item in S:
3673
        if not(item in gen):
3674
            T.append(item)
3675
    T = R * S
3676
    for item in T:
3677
        if not(item in gen):
3678
            T.append(item)
3679
    T = S * R
3680
    for item in T:
3681
        if not(item in gen):
3682
            T.append(item)      
3683
    n = len(R._List) # size of original group
3684
    # gen is now the set of additive generators for the ideal
3685
    id = copy(gen) 
3686
    for g1 in grp:
3687
        for g2 in id:
3688
            u = g1 + g2                  
3689
            if not (u in grp):
3690
                grp.append(u)
3691
                if len(grp)*2 > n:
3692
                    return R
3693
    return GroupSet(sorted(grp)) # TODO: determine sorting algorithm  
3694
    
3695
def ConjugacyClasses(Glist):
3696
    # This finds the conjugacy classes of the group in Glist.
3697
    if isinstance(Glist, GroupSet):
3698
        G = Glist
3699
    else:
3700
        G = GroupSet(Glist)    
3701
    # If the elements are Gap elements, delegate the task to Gap.
3702
    if isinstance(G._List[0], GapInstance):
3703
        Ggap = gap(G._List)
3704
        LL = Ggap.ConjClasses() # will produce a list of lists, in Gap
3705
        L = list(LL)
3706
        CC = []
3707
        for item in L:
3708
            CC.append(GroupSet(list(item)))
3709
        return CC
3710
    CC = []
3711
    FlatC = []    
3712
    for g in G._List:
3713
        if not(g in FlatC):
3714
            Class = []
3715
            for h in G._List:
3716
                u = (h * g) * h^(-1)
3717
                if not(u in Class):
3718
                   Class.append(u)
3719
            Class.sort()
3720
            CC.append(GroupSet(Class))
3721
            FlatC = Flatten(CC)
3722
    return GroupSet(CC)               
3723
                
3724
# This command assumes that both HList and KList are subgroups of some larger group.  
3725
# If one of the groups is a subgroup of the other, Commutator is faster.                
3726
def MutualCommutator(HList, KList):
3727
    if isinstance(HList, GroupSet):
3728
        H = HList
3729
    else:
3730
        H = GroupSet(HList)
3731
    if isinstance(KList, GroupSet):
3732
        K = KList
3733
    else:
3734
        K = GroupSet(KList)
3735
    #If the elements are in GAP, delegate the job to GAP 
3736
    if isinstance(H._List[0], GapInstance) and isinstance(K._List[0], GapInstance):
3737
        Hgap = gap(H._List)
3738
        Hgap.OriginalGroup() 
3739
        Kgap = gap(K._List)
3740
        Agap = Kgap.CommSubgroup()
3741
        return GroupSet(list(Agap))
3742
    gen = []
3743
    for h in H._List:
3744
        for k in K._List:
3745
            u = (h^(-1) * k^(-1)) * (h * k) 
3746
            if not(u in gen):
3747
                gen.append(u)
3748
    grp = copy(gen)  
3749
    for g1 in grp:
3750
        for g2 in gen:
3751
            z = g1 * g2                  
3752
            if not (z in grp):
3753
                grp.append(z)
3754
    return GroupSet(sorted(grp))  
3755
 
3756
# This command assumes that KList is a subset of GList.  
3757
# If K generates the group G, then this gives the derived group G'               
3758
def Commutator(GList, KList):
3759
    if isinstance(GList, GroupSet):
3760
        G = GList
3761
    else:
3762
        G = GroupSet(GList)
3763
    if isinstance(KList, GroupSet):
3764
        K = KList
3765
    else:
3766
        K = GroupSet(KList)
3767
    #If the elements are in GAP, delegate the job to GAP 
3768
    if isinstance(G._List[0], GapInstance):
3769
        Ggap = gap(G._List)
3770
        Ggap.OriginalGroup() 
3771
        Kgap = gap(K._List)
3772
        Agap = Kgap.CommSubgroup()
3773
        return GroupSet(list(Agap))
3774
    n = len(G._List)    
3775
    gen = []
3776
    for h in G._List:
3777
        for k in K._List:
3778
            u = (h^(-1) * k^(-1)) * (h * k) 
3779
            if not(u in gen):
3780
                gen.append(u)         
3781
    grp = copy(gen)  
3782
    for g1 in grp:
3783
        for g2 in gen:
3784
            z = g1 * g2                  
3785
            if not (z in grp):
3786
                grp.append(z)
3787
                if len(grp)*2 > n:
3788
                    # over half the elements, return original list
3789
                    return G
3790
    return GroupSet(sorted(grp))  
3791
 
3792
def ExpressAsWord(genList, target):
3793
    # This routine interfaces with GAP, to find the shortest way of expressing a target permutation in terms of a list of permutations.
3794
    # genList must be a list of strings, each of which represents the variable that a permutation has been defined.
3795
    # target is the permutation that must be reached.
3796
    #
3797
    # There is a twist.  GAP interprets permutations multiplied from left to right, assuming (fog)(x) = g(f(x)).
3798
    # If we multiply permutations from right to left (assuming (fog)(x) = f(g(x)) ), then we have to invert all of the input and output
3799
    # permutations to compensate.  
3800
    #
3801
    # First we evaluate the strings in genList to find the permutations.
3802
    permList = []
3803
    for item in genList:
3804
        perm = eval(item)
3805
        if isinstance(perm, Perm):
3806
            perm = PermToCycle(perm)
3807
        ## Take the inverse of the permutation to compensate GAP.
3808
        perm = perm^-1    
3809
        permList.append(perm)    
3810
    tmpstr = "Perms := Group(" + str(permList[0])
3811
    for i in range(1, len(genList)):
3812
        tmpstr = tmpstr + ", " + str(permList[i])
3813
    tmpstr = tmpstr + ")"
3814
    gap.eval(tmpstr)
3815
    #print tmpstr
3816
    tmpstr = 'phi := EpimorphismFromFreeGroup(Perms:names:=["' + genList[0] + '"'
3817
    for i in range(1, len(genList)):
3818
        tmpstr = tmpstr + ', "' + genList[i] + '"'
3819
    tmpstr = tmpstr + "])"
3820
    gap.eval(tmpstr)
3821
    #print tmpstr
3822
    if isinstance(target, Perm):
3823
        tar = PermToCycle(target)
3824
    else:
3825
        tar = target
3826
    ## Take the inverse of the target to compensate GAP.
3827
    tar = tar^-1    
3828
    tmpstr = "PreImagesRepresentative( phi," + str(tar) + ")"
3829
    outstr = gap.eval(tmpstr)
3830
    #print tmpstr
3831
    return outstr    
3832

3833
def AddRingVar(str):
3834
    # str is a string which will be the name of the variable.
3835
    tmpstr = str + ' := Indeterminate(CurrentRing, "' + str + '")'
3836
    gap.eval(tmpstr)
3837
    tmpstr = "global " + str + "; " + str + ' = GapElement(gap("' + str + '"))'
3838
    exec(tmpstr)
3839
    return None
3840

3841
GenList = []
3842
DegreeList = []
3843

3844
def InitDomain(arg1, *args):
3845
    global LastInit
3846
    global UnivFieldVar
3847
    global BaseCharacteristic
3848
    global CurrentField
3849
    global GenList
3850
    global DegreeList
3851
    GenList = []
3852
    DegreeList = []
3853
    # if there is one argument, it must be an integer, 0 or positive prime.
3854
    # if there are two arguments, the second one must be a string, giving a name of the universal variable.
3855
    if not(isinstance(arg1, (int, long, IntType))):
3856
        return "First argument must be an integer."
3857
    if (arg1 < 0):
3858
        return "First argument must be positive."
3859
    if (arg1 > 0) and not(arg1 in Primes()):
3860
        return str(arg1) + " is not prime."
3861
    if arg1 == 0:
3862
        CurrentField = QQ
3863
    else:
3864
        CurrentField = GF(arg1)
3865
    BaseCharacteristic = arg1    
3866
    LastInit = "InitDomain"
3867
    if len(args) == 0:
3868
        UnivFieldVar = ""
3869
        return None # we are done
3870
    UnivFieldVar = args[0]
3871
    tmpstr = "global " + UnivFieldVar + '; TempDomain = CurrentField["' + UnivFieldVar + '"]; ' + UnivFieldVar + " = TempDomain._first_ngens(1)[0]"
3872
    exec(tmpstr)
3873
    return None
3874
       
3875
def AddFieldVar(str):
3876
    if (BaseCharacteristic > 0) and (len(GenList) > 1):
3877
        if DegreeList[0] > 0:
3878
            return "Sage is unable to define a variable on an extension of an extension of a finite field."
3879
    global LastFieldVar
3880
    LastFieldVar = str
3881
    tmpstr = "global " + str + '; TempDomain = CurrentField["' + str + '"]; ' + str + " = TempDomain._first_ngens(1)[0]"
3882
    exec(tmpstr)   
3883
    return None
3884
    
3885
def RationalFunctions(str):
3886
    global GenList
3887
    global DegreeList
3888
    tmpstr = "global " + str + "; global CurrentField; CurrentField = FunctionField(CurrentField, names = ('" + str + "',)); "
3889
    tmpstr = tmpstr + str + " = CurrentField._first_ngens(1)[0]"
3890
    exec(tmpstr)
3891
    GenList.append(str)
3892
    DegreeList.append(0)
3893
    
3894
def ListField():
3895
    if BaseCharacteristic == 0:
3896
        return "Field is infinite!"
3897
    if len(GenList) == 2:  # in this case, list will not work.  Use a less efficient method.
3898
        return Ring(GenList[0], GenList[1])    
3899
    if len(GenList) == 1 and DegreeList[0] == 0: 
3900
        return "Field is infinite!"   
3901
    return GroupSet(sorted(list(CurrentField)))
3902
    
3903
def Compon(exp, L_List):
3904
    inp = copy(exp)
3905
    out = []
3906
    for item in L_List:
3907
        a = inp.coefficient(item)
3908
        inp = expand(inp - a*item)
3909
        out.append(a)
3910
    out.append(inp)
3911
    return out    
3912
   
3913
def Vectorize(exp):
3914
    try:
3915
        L = vector(exp) 
3916
    except: # not a vector - we are at the lowest level.
3917
        return exp   
3918
    output = []
3919
    for item in L:
3920
        if isinstance(item, (Rational, IntModType, int, long, IntType)):
3921
            output.append(item)
3922
        else:    
3923
            output.append(Vectorize(item))
3924
    return Flatten(output)    
3925
                                               
3926
def InterpolationPolynomial( pointsList, varname):
3927
    Rpolyring = PolynomialRing(CurrentField, varname)
3928
    return Rpolyring.lagrange_polynomial(pointsList)                            
3929
            
3930
def Cyclotomic(n, str):
3931
    TempR = QQ[str]
3932
    return TempR.cyclotomic_polynomial(n)
3933
            
3934
def ToBasis(*args):
3935
    if (len(args) == 0 or len(args) > 2):
3936
        print "One or two arguments are required"
3937
        return None
3938
    if len(args) == 1:
3939
        B = Basis(args[0])
3940
    else:
3941
        B = Basis(args[0], args[1])
3942
    if not (B._worked):
3943
        print "Error: linearly dependent."
3944
        return false
3945
    else:
3946
        print "Sucessful mapping constructed."
3947
        return B    
3948
 
3949
def Coefficients(Bas, value):
3950
    if isinstance(Bas, Basis):
3951
        return Bas.Coeff(value)
3952
    print "First argument must be a basis set by ToBasis."
3953
    return false            
3954
            
3955
def SimpleExtension(arg1, arg2):
3956
    # finds a single element so that Q(w) = Q(arg1, arg2)
3957
    # Note that Sage has several ways of computing this w, but
3958
    # we want to be consistent with the textbook, and return
3959
    # arg1 + n*arg2 for the smallest positive n that works.
3960
    if BaseCharacteristic > 0:
3961
        print "Algorithm does not work on finite fields."
3962
        return false
3963
    w = arg1 + arg2
3964
    v = Vectorize(w)
3965
    n = len(v)
3966
    cont = true
3967
    while cont:
3968
        v = Vectorize(w)
3969
        n = len(v)
3970
        B = [w^i for i in range(n)]
3971
        NB = Basis(B)
3972
        if NB._worked:
3973
            cont = false
3974
        else:
3975
            w = w + arg2 # try again with a different w
3976
    return w    
3977
    
3978
                 
3979
            
3980
print("Initialization Done")
3981

3982
Product

Resources

Company
