Sage は Unix 系OS、MacOS、Windows OS いずれにも対応している数式処理システムである。Maple, Mathematica などと同様と考えて良い。決定的な違いは、オープンソースつまり、システムの中身が公開されていて、GPL というライセンス形態で、開発されており、無償で提供されている点である。責任を持つ企業などがないから不安だと考えるひともいるかも知れないが、同時に、たとえば Maple や Mathematica は利用しているアルゴリズムなどが公開されておらず、問題があった場合には、開発元が修正するまでは待つしかない。アルゴリズムが公開されていないということは、学術的に信頼ができるかどうか、peer review の原則にたつ、研究の評価にもある枷がかかることになる。他にも Magma などは、代数計算において優れたシステムであるが、これをシステムの一部にとりこんで、自ら開発する事はできない。新しいクラスなどを定義することができないためである。そう考えたときに、Sage のようなオープンソースのシステムの価値は高い。
これまで、Sage 以外にも、オープンソースのシステムはいくつも存在したが、研究のために開発されたものが多く、汎用性が高いものは少なく、特に、初心者に使いやすいものも少なかった。その意味では、Maple や Mathematica は優れていると言わざるを得なかったであろう。しかし、教育目的を考えた時、これら市販のシステムの値段の高さはいかんともしがたい。日本以外の100を越える国々（日本以外ほとんどすべてと言ってもよい）では、Student License が 100ドル未満の値段で提供されているが、日本は例外的にこの二倍から三倍の価格となっている。それでも安くなった方で、以前はもっと高額であった。わたしは、本社にこのことを訴えたが返事はすべて代理店経由でくるため、ほとんど改善は得られなかった。日本の数学が一般人に受け入れられない理由はいくつもあるだろうが、そのうちの一つだと個人的には認識している。また、学生用以外のものの値段は、よほどの研究費をもっている研究機関でないと、すくなくとも自由には利用できない。たとえ、安価になっても、学生が家で、そして、卒業してからも一生のものとして、色々なカタチで数学を身近なものとするために、利用することをかんがえると、たとえ機能があるていど限定されていても、無償のものの価値は人間の財産として重要である。
Sage の良い点は、既存の オープンソースのシステムをも利用できるように設計されている点である。Maxima がもとだと言われているが、群論関係では有名な GAP も利用できるようになっている。他のオープンソースのシステムを利用してきた研究者にもなじみやすいのではないだろうか。
Sage には、四つのインターフェイスが用意されている。Terminal モード、Notebook モード、SageMathCell、SageMathCloud である。Notebook モードはオンラインでも、インストールしても利用できる。
SageMathCell は、電卓としても（実際に、SageMathCell Server を用いた、SmartPhone のアプリも存在する）簡単なデモ用にも、Interact モードを利用して教育用にも有効である。SageMathCell では、Sage, GAP, GP, HTML, Maxima, Octave, Python, R が利用可能である。Octave は MatLab のコードがほぼ使えるので、線形代数などを使った数値計算にも有効である。
Sage によって数学が、専門家だけのものではなく、また数学を使う人が、専門家だけではなく、すべてのひとにひろがることを願っている。計算力があることは望ましいし、そのことは、研究推進に大きな力になるであろう。しかし、計算力がなくても、すこし数学を理解することで、Sage の利用が可能となり、数学が多くの人に利用されるようになると信じている。スポーツがプロだけのものではなく、肉体的な訓練を十分しなくても楽しめ、かつ人生を豊かにするように、数学もすべてのひとが楽しむことができるようになることを期待している。数学は人間すべてにとってすばらしいもので、利用し、楽しむ人それぞれにとって、意味はすこしずつ違っていても、価値のあるものだと信じるからである。
[INPUT]
def birthday(n): # same as birthday 2
return 1-RR(binomial(365,n)*factorial(n)/(365^n))
[INPUT]
def birthday2(n):
p=1
for i in range(0,n):
p = p*(365-i)/365
return 1-RR(p)
[INPUT]
def birthday3(n):
p=1
for i in range(0,n):
p = p*(366-i)/366
return 1-RR(p)
[INPUT]
n=var('n'); f = 1-RR((365-31)/365)^n
[INPUT]
var('a, b, c')
solve(x^2+b*x+c==0,x)
[INPUT]
var('a,b')
f=expand((a-b)^3*(a+3*b))
f; factor(f)
[INPUT]
x,y,z=var('x,y,z')
eq=(x-y-z==0, 5*x+13*y==8, -13*y+14*z==3)
sol=solve(eq, x,y,z); sol
[INPUT] # Problem 1.6.8 # H, O, Ca, C in this order # H3O + CaCO3 -> H2O + Ca + CO2 A8 = Matrix([[3,0,-2,0,0,0],[1,3,-1,0,-2,0],[0,1,0,-1,0,0],[0,1,0,0,-1,0]]) print(A8) print A8.echelon_form()
[INPUT] # Pb, N, Cr, Mn, O in this order # PbN6 + CrMn2O8 -> Pb3O4 + Cr2O3 + MnO2 + NO # treat as integer matrix A10 = Matrix([[1, 0, -3, 0, 0, 0, 0], [6,0,0,0,0,-1, 0],[0,1,0,-1,0,0,0],[0,1,0,0,-1,0,0],[0,8,-4,-3,-2,-1,0]]) print(A10) print A10.echelon_form()
[INPUT] # Pb, N, Cr, Mn, O in this order # PbN6 + CrMn2O8 -> Pb3O4 + Cr2O3 + MnO2 + NO # treat as a matrix with rational number entries 有理数（分数）行列として扱う A10 = Matrix(QQ, [[1, 0, -3, 0, 0, 0, 0], [6,0,0,0,0,-1, 0],[0,1,0,-1,0,0,0],[0,1,0,0,-1,0,0],[0,8,-4,-3,-2,-1,0]]) print(A10) print A10.echelon_form()
[INPUT]
var('x1,x2,x3,x4,x5,x6')
eq = [x1 - 3*x3 == 0, 6*x1 - x6 ==0, x2 - x4 == 0, x2 - x5 == 0, 8*x2 - 4*x3 - 3*x3 -2 *x5 - x6 == 0]
solve(eq, [x1,x2,x3,x4,x5,x6])
[INPUT] U=matrix([[1,1,1,1],[1,-1,0,0],[1,0,-1,0],[1,0,0,-1]]); print(U); U.det()
[INPUT]
var('a')
U1=matrix([[1,a,a,a],[1,-a,0,0],[1,0,-a,0],[1,0,0,-a]]); show(U1); U1.det()
T = matrix([[-1, -1, 2],[1, 0, -1],[0, 1, 1]]); T
T1=T[:,0]; T2=T[:,1]; T3=T[:,2]; T1; T2; T3
[INPUT] # A definition of a matrix A = Matrix([[1, 0, 2],[2, -1, 3], [1, 1, 8]])
A.det(); A.trace()
factor(A.charpoly())
A.transpose(); A.adjoint()
A*A.adjoint()
[INPUT]
var('x1,x2,x3')
T = matrix([2*x1+3*x2+x3,x1-x2-x3,-x1+3*x2+5*x3])
[INPUT] A1=matrix(3,3); A1[0] = T(x1=1,x2=0,x3=0); A1[1] = T(x1=0,x2=1,x3=0); A1[2] = T(x1=0,x2=0,x3=1) A11 = A1.transpose(); A11
[INPUT] A1T=A11.transpose(); A1T[0].cross_product(A1T[1])
[INPUT] A11.det().abs()
[INPUT]
var('a,b')
A2=matrix([[a,b,b,b],[b,a,b,b],[b,b,a,b],[b,b,b,a]])
B2=matrix([[b,b,b,b],[b,b,b,b],[b,b,b,b],[b,b,b,b]])
A2; B2
[INPUT]
print(A2.determinant())
print("Determinant of A")
factor(A2.determinant())
[INPUT] print(A2.characteristic_polynomial()) A2.eigenvalues()
[INPUT] B2.eigenvectors_right()
[INPUT] D2,T2 = B2.eigenmatrix_right(); print(D2) print(T2)
[INPUT] T2^-1*A2*T2
[INPUT] A3 = matrix([[3,-1,2,0],[1,2,-2,5],[-1,3,1,1],[2,3,1,-2]]) A3
[INPUT] print(A3.det()) factor(A3.det())
[INPUT] b3 = matrix([1,2,3,4]).transpose() print(A3[:,0]) B31 = copy(A3); B32 = copy(A3) print(B31) print(B32) B31.column(1,b3) print(B31) A3
[INPUT] # x は変数として定義されているので定義しなくてもよい。 f1a=cos(arcsec(x))/(x*sqrt(x^2-1)) integral(f1a,x,2/sqrt(3),2)
[INPUT]
var('y')
f1b=1/(y*sqrt(1+(log(y)^2)))
integral(f1b,y,1,exp(1))
[INPUT]
var('t')
f1c=(t^4-4*t^3+2*t^2-3*t+1)/(t^2+1)^3
if1c=integral(f1c,t); if1c; if1c.derivative(); if1c.derivative().factor()
[INPUT] f1d=16*x^3*(log(x))^2 integral(f1d,x)
[INPUT]
var('v')
f1e=sqrt(1-v^2)/v^2
integral(f1e,v)
[INPUT] f1g=1/(1+sin(x)+cos(x)) integral(f1g,x)
[INPUT]
var('n')
f6=(n/(n+1))^n
f6.limit(n=oo)
[INPUT]
var('y')
f9=y^2/(cos(y)-cosh(y))
f9.limit(y=0)
[INPUT] f7=exp(sin(x)) taylor(f7,x,0,2)
[INPUT] f8=arcsinh(x) taylor(f8,x,0,10)
[INPUT] #下の値が1より大きければ収束、小さければ発散 f10=n*((2*n)*(2*n+1)/(2*n-1)^2-1) f10.limit(n=oo)
[INPUT] # Any line after # is a comment. # R.<x> is a polynomial ring with rational coefficients # P.<z> is a power series ring, which is called Laurent Series Ring, with coefficients in R.R.<x>=PolynomialRing(QQ,1) P.<z> = LaurentSeriesRing(R) f = 1-2*x*z + z^2 g = 1/f # Display f and find Coefficients of f print(f) f.coefficients()
[INPUT]
# find the first 6 coefficients starting from 0th coefficient
for i in range(6):
g.coefficients()[i]
[INPUT] # the first 6 coefficients of g = 1/f, simpler way g.coefficients()[:6]
[INPUT] # using Landau's symbol g+O(z^6)
[INPUT]
x=var('x'); f = x^3 - 3*x^2 - 3*x -1; [f(n) for n in range(8)]
[INPUT]
x=var('x')
f = x^5 + 2*x - 5
print(f.subs(x=1))
print(f.subs(x=2))
find_root(f,1,2)
f.derivative()
for i in range(0,9):
factor(2^(2^i)+1)
def fe(n):
v = []
for i in range(0,n):
v.append(factor(2^(2^i)+1))
return v
m = fe(9); m
[INPUT] def me(n): v = [] for i in prime_range(2,n): if is_prime(2^i-1): v.append(i) return v
def fibonatti(n):
f=[1,1]
for i in range(2,n):
f.append(f[i-1]+f[i-2])
return f
fibonatti(10)
[INPUT]
var('n,m,k')
def surj(m,n): # number of onto mappings from an m-set to an n-set
s=0
for k in range(0,n):
s = s+(-1)^k*binomial(n,k)*(n-k)^m
return s
[INPUT]
var('y')
solve(y==3*x/(x-2),x)
単純に x について解くことで逆関数を求めます。x だけは最初から変数として登録してありますから、定義する必要はありませんが、他の変数（この例では、y）は定義しないと、使うことができません。
[INPUT] g(x) = 3*x/(x-2) h(x) = 2*x/(x-3) p1=plot(g,(x,-10,15), ymax=10, ymin=-1, color='red', detect_poles = 'show') p2=plot(h,(x,-10,15), ymax=10, ymin=-1, color='blue', detect_poles = 'show') p3=plot(x,(x,-10,15), color='black') (p1+p2+p3).show(aspect_ratio=1)上の例で求めた、逆関数が h(x) となっています。それも定義して、グラフにしたものです。plot の中のオプションは、大体想像できると思います。まずは、plot ?? と入力して、どのようなものが使えるか、また例を見てみるのが良いでしょう。plot は条件を表す点（この場合はグラフ）を決めるものです。それを表示するのは、show()。一つのグラフとして表すので、p1+p2+p3 としています。p3 は対称軸を表しています。
[INPUT]
f(x) = x/(x^2+1)
plot(f(x),(x,-10,10)).show(figsize=[8,4])
f.derivative().show()
sol=solve(f.derivative(),x)
print(sol)
for i in range(len(sol)):
print("at x = {0}, f(x) = {1}".format(sol[i].rhs(), f(sol[i].rhs())))
minval,xmin=find_minimum_on_interval(f,-5,5)
print("Min at x={0}, Min value = {1}").format(xmin, minval)
関数の極大極小を求めるものです。関数を定義し、グラフを表示し、導関数を計算し表示、導関数が 0になる、つまり臨界点（critical point）を計算し、その値を、表示させています。solve の答えは、リストで、[x==a, x==b] と得られますから、それぞれの右辺を代入した結果を表示させています。最後は、極小の値を、数値で求めています。print の中身は、{0}, {1} は表示させるデータを意味し、最初の例では、sol[i].rhs() だったり、f(sol[i].rhs()) だったりということです。range(2) は 0, 1 を表します。そして、sol[0] が一つ目、sol[1] が二つ目と言うことになります。
[INPUT] print(xgcd(51, 288)) print(xgcd(357,629)) print(xgcd(180, 252)) print(xgcd(4278, 71929)) print(factor(4278)) print(factor(71929))
[INPUT] U18=[g for g in Integers(18) if g.is_unit()]; len(U18);U18;[g.multiplicative_order() for g in Integers(18) if g.is_unit()]
[INPUT] S9=SymmetricGroup(9); A9=AlternatingGroup(9); Set([x.order() for x in S9]); Set([x.order() for x in A9])N=8 or 9 may be a practical limit. See magma section.
[INPUT] S6=SymmetricGroup(6); Q= Set([x for x in S6 if x.order() == 4]); I= Set([x for x in S6 if x.order() == 2]); Q.cardinality(); I.cardinality()
[INPUT] G=SymmetricGroup(6) r=G((1,3,4,5,6)); s=G((1,3,2)) K=G.subgroup([r,s]) K; K.order()
g=G('(1,6,5,2,3)')
g.word_problem(K.gens())
毎回同じ方法（置換）で順番を入れ換えるカードを切る機械がある。A から K までのハートのトランプをこの機械で二回切ると 10, 9, Q, 8, K, 3, 4, A, 5, J, 6, 2, 7 が得られた。一回切ったあとはどうなっているか。b をその置換とすると、c = b・b = (1,8,4,7,13,5,9,2,12,3,6,11,10) だから、b の 13乗が恒等写像であることが、わかる。そこで、b の逆元を a とすると、b の 12乗 すなわち、c の 6乗が a であることが分かるので、その逆元をもとめれば b が分かる。
[INPUT] c = Permutation([8,12,6,7,9,11,13,4,2,1,10,3,5]) b=power(c,6).inverse() b; a.cycle_string()
[INPUT]
a=Permutation('(1,2,8,12,4,3,7,6,13,11,5,10,9)')
b=a.inverse()
s=range(1,14)
b.action(s); (b*b).action(s)
[INPUT]
G721=PermutationGroup([[(1,2)],[(2,3,4)]])
print("Group in question is: \n{0}".format(G721))
print("\nOrder of the group is {0}".format(G721.order()))
print("Is isomorphic to Sym(4)? {0}.".format(G721.\
is_isomorphic(SymmetricGroup(4))))
[INPUT]
G735=PermutationGroup([[(1)], [(1,2),(3,4)], [(1,2,3,4),(5,6)],\
[(1,3),(2,4)], [(1,4,3,2),(5,6)], [(5,6),(1,3)], \
[(1,4),(2,3)], [(2,4),(5,6)]])
print("Group in question is: \n{0}".format(G735))
print("\nOrder of the group is {0}".format(G735.order()))
print("Is Abelian? {0}".format(G735.is_abelian()))
print("\nStabilizer of 1 is: \n{0}".format(G735.stabilizer(1)))
print("Order = {0}".format(G735.stabilizer(1).order()))
print("Orbit of 1 is the set {0}".format(G735.orbit(1)))
print("\nStabilizer of 3 is: \n{0}".format(G735.stabilizer(3)))
print("Order = {0}".format(G735.stabilizer(3).order()))
print("Orbit of 3 is the set {0}".format(G735.orbit(3)))
print("\nStabilizer of 5 is: \n{0}".format(G735.stabilizer(5)))
print("Order = {0}".format(G735.stabilizer(5).order()))
print("Orbit of 5 is the set {0}".format(G735.orbit(5)))
[INPUT]
# The following program takes time to run. Please wait patiently.
G19=AbelianGroup(3,[12,20,10],names = list("abc"))
sub = G19.subgroups()
sub20=[g for g in sub if order(g)==20]
len(sub20)
for i in range(5):
print(sub20[i])
[INPUT]
# The following program takes time to run. Please wait patiently.
G20=AbelianGroup(3,[20,5,60],names = list("stu"))
sub1 = G20.subgroups()
print("Number of subgroups = {0}".format(len(sub1)))
sub1n=[g for g in sub1 if g.is_cyclic()==False]
print("Number of noncyclic subgroups = {0}".format(len(sub1n)))
sub14=[g for g in sub1n if g.invariants()==[2,2]]
print("Number of subgroups isomorphic to Z2xZ2 = {0}".format(len(sub14)))
print(sub14[0])
print(list(sub14[0]))
[INPUT] R=Integers(220) print(R.unit_group_order()) print(factor(euler_phi(200))) print(R.unit_gens()) r220=[a.multiplicative_order() for a in R.unit_gens()] AbelianGroup(r220).elementary_divisors()
[INPUT]
# nth cyclotomic polynomial
for i in range(1,20):
print "n=",i,": ", cyclotomic_polynomial(i)
[OUTPUT]
n= 1 : x - 1
n= 2 : x + 1
n= 3 : x^2 + x + 1
n= 4 : x^2 + 1
n= 5 : x^4 + x^3 + x^2 + x + 1
n= 6 : x^2 - x + 1
n= 7 : x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
n= 8 : x^4 + 1
n= 9 : x^6 + x^3 + 1
n= 10 : x^4 - x^3 + x^2 - x + 1
n= 11 : x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
n= 12 : x^4 - x^2 + 1
n= 13 : x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4
+ x^3 + x^2 + x + 1
n= 14 : x^6 - x^5 + x^4 - x^3 + x^2 - x + 1
n= 15 : x^8 - x^7 + x^5 - x^4 + x^3 - x + 1
n= 16 : x^8 + 1
n= 17 : x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^8
+ x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
n= 18 : x^6 - x^3 + 1
n= 19 : x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10
+ x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
[INPUT]
# coefficient c[m] of x^m of the nth cyclotomic polynomial that is not 1 or -1
for n in range(1,150):
f = cyclotomic_polynomial(n)
c = f.coefficients()
for m in range(len(c)):
if c[m].abs() > 1:
print n, m, c[m]
[OUTPUT]
105 5 -2
105 27 -2
[INPUT] g=graphs.PetersenGraph() pg=g.plot() pg.show()
[INPUT]
# Information of J(n,d) taken from Terwilliger's Subconstituent Algebra paper 1,2,3
n,d,i,j,nu,mu,diam = var('n,d,i,j,nu,mu,diam')
# Information taken from page 192 in subconst 3
# r = -n+d-1; s = -n-2; t = -n*(n-1)/(d*(n-d))
def johnson_bj(n,d,j):
return (d-j)*(n-d-j)
def johnson_cj(n,d,j):
return j^2
def johnson_aj(n,d,j):
return d*(n-d) - j^2 - (d-j)*(n-d-j)
def johnson_array(n,d): # d has to be a positive integer
c = [johnson_cj(n,d,j) for j in range(d+1)]
a = [johnson_aj(n,d,j) for j in range(d+1)]
b = [johnson_bj(n,d,j) for j in range(d+1)]
return matrix(3,d+1,[c,a,b])
def johnson_evj(n,d,j):
s = -n-2
return d*(n-d) + j*(j+1+s)
def johnson_devj(n,d,j):
t = -n*(n-1)/(d*(n-d))
return n-1 + t*j
# Information taken from page 372 in subconst 1
def johnson_bsj(n,d,j):
r = -n+d-1; s = -n-2; t = -n*(n-1)/(d*(n-d))
if d==j:
return 0
else:
return t*(j-d)*(j+1+s)*(j+1+r)/((2*j+1+s)*(2*j+2+s))
def johnson_csj(n,d,j):
r = -n+d-1; s = -n-2; t = -n*(n-1)/(d*(n-d))
if d==j:
return -t*d*(d+s-r)/(2*d+s)
else:
return -t*j*(j+s+d+1)*(j+s-r)/((2*j+1+s)*(2*j+s))
def johnson_asj(n,d,j):
r = -n+d-1; s = -n-2; t = -n*(n-1)/(d*(n-d))
if d==j:
return n-1 - johnson_csj(n,d,j)
else:
return n-1 - johnson_csj(n,d,j) - johnson_bsj(n,d,j)
def johnson_dualarray(n,d): # d has to be a positive integer
c = [johnson_csj(n,d,j) for j in range(d+1)]
a = [johnson_asj(n,d,j) for j in range(d+1)]
b = [johnson_bsj(n,d,j) for j in range(d+1)]
return matrix(3,d+1,[c,a,b])
def johnson_diam(n,d,nu,mu): # n, d, nu, mu have to be a nonnegative integers
d1 = d-2*nu
d2 = min(d-mu,n-d-2*nu)
if (d1 < 0) or (d1 < d - 2*nu) or (d1 > d - mu):
if (d2 < 0) or (d2 < d - 2*nu) or (d2 > d - mu):
return []
else:
return [d2]
else:
if (d2 < 0) or (d2 < d - 2*nu) or (d2 > d - mu) or (d1 == d2):
return [d1]
else:
return [d1, d2]
def johnson_mod_parameters(n,d):
return [(nu,mu,diam)
for nu in range(0,d+1)
for mu in range(nu,d+1)
for diam in johnson_diam(n,d,nu,mu)]
# Information taken from page 198 in subconst 3
def johnson_mod_bj(n,d,nu,mu,diam,j):
return (diam-j)*(n-diam-2*nu-mu-j)
def johnson_mod_cj(n,d,nu,mu,diam,j):
return j*(j+2*nu-mu)
def johnson_mod_aj(n,d,nu,mu,diam,j):
return d*(n-d)+mu*(mu+diam-n-1)+diam*(diam-n+2*nu)+j*(n-4*nu-2*j)
def johnson_mod_array(n,d,nu,mu,diam): # d has to be a positive integer
c = [johnson_mod_cj(n,d,nu,mu,diam,j) for j in range(diam+1)]
a = [johnson_mod_aj(n,d,nu,mu,diam,j) for j in range(diam+1)]
b = [johnson_mod_bj(n,d,nu,mu,diam,j) for j in range(diam+1)]
return matrix(3,diam+1,[c,a,b])
def a2m(a,i,j):
if j==i+1:
return a[0][i+1]
elif j==i:
return a[1][i]
elif j==i-1:
return a[2][i-1]
else:
return 0
def array2matrix(a):
return matrix(len(a[0]),len(a[0]),[a2m(a,i,j) for i in range(len(a[0]))
for j in range(len(a[0]))])
[INPUT]
# distance-regular graphs with classical parameters (d,b,alpha, beta, b)
m, n, d, b, alpha, beta, h, i, j, q= var('m, n, d, b, alpha, beta, h, i, j, q')
def bn(m,b): # (b^m-1)/(b-1)
if m<0:
return 0
else:
return sum(b^i, i, 0, m-1)
def a2m(a,i,j):
if j==i+1:
return a[0][i+1]
elif j==i:
return a[1][i]
elif j==i-1:
return a[2][i-1]
else:
return 0
def array2matrix(a):
return matrix(len(a[0]),len(a[0]),[a2m(a,i,j) for i in range(len(a[0]))
for j in range(len(a[0]))])
def classical_bj(d, b, alpha, beta, j):
return (bn(d,b)-bn(j,b))*(beta-alpha*bn(j,b))
def classical_cj(d, b, alpha, beta, j):
return bn(j,b)*(1+alpha*bn(j-1,b))
def classical_aj(d, b, alpha, beta, j):
return bn(j,b)*(beta-1+alpha*(bn(d,b)-bn(j,b)-bn(j-1,b)))
def classical_array(d, b, alpha, beta): # d has to be a positive integer
c = [classical_cj(d, b, alpha, beta, j) for j in range(d+1)]
a = [classical_aj(d, b, alpha, beta, j) for j in range(d+1)]
b = [classical_bj(d, b, alpha, beta, j) for j in range(d+1)]
return matrix(3,d+1,[c,a,b])
def classical_evj(d, b, alpha, beta, j):
return b^(-j)*classical_bj(d,b,alpha,beta,j) - bn(j,b)
G:=SymmetricGroup(6); r:=(1,3,4,5,6); s:=(1,3,2); K:=Subgroup(G,[r,s]); Size(K);
t:=s*r; Factorization(K,t); Factorization(K,(1,2)*(3,4,5,6)); # Operation は左の元から Factorization(K,(1,2,3,4,5));