程序员用的数学

程序员用的数学
作者：韩守庆著
摘要/内容：
第1章第三方数学软件：V5.5测试（20131204）
第2章 OpenGL和D3D数学基础：线性代数、空间解析几何、局部微分几何、泛函分析、计算机图形学等3D数学基础知识
第3章算法分析
第4章 Eigen库
第5章 MATLAB
第6章 Python图论
第7章 GAP
第8章 MIRACL大数运算库
第9章
第10章
第11章
第12章
第13章
第14章
附录1：数学英语

chap1
功能:求符号函数对某一变量的导数（Symbolic Expression Derivation）
f="x^8"
f = "x^8"
SymExpressionDerivation(f)
ans = "8*x^7"
【
用Mathematica进行求导运算
在Mathematica系统中，用D[f,x]表示f(x)对x的一阶导数，用D[f,{x,n}]表示f(x)对x的n阶导数.
In[1]: =D[x^8,x]
Out[1]=8x^7
In[2]:=D[x^8*Sin[x],x]
Out[2]=x^8Cos[x] +8x^7Sin[x]
】

功能:一元连续函数的积分函数.
格式:Integral(f,a,b)
说明:f为符号函数句柄,其未知数变量默认为x,积分区间为[a,b],当a或b为inf或-inf时,认为是无穷远
f="x/(4+x^2)"
f = "x/(4+x^2)"
Integral(f,0,1)
ans =
[ 0.11157177565710 ]
【
Mathematica数值积分函数为NIntegrate，其调用格式为：
NIntegrate[f,{x,xmin,xmax}]
解
In[1]：=NIntegrate[x/(4+x^2),{x,0,1}]
Out[1]=0.111572
】
功能:数据拟合
格式:polyfit(x,y,n,erro);polyfit(x,y,n);polyfit(x,y)
进行三次多项式拟合
a=[1 2 3 4]
a =
[ 1.00000000000000 2.00000000000000 3.00000000000000 4.00000000000000 ]
b=[16 25 46 85]
b =
[ 16.0000000000000 25.0000000000000 46.0000000000000 85.0000000000000 ]
polyfit(a,b,3)
ans=
"1.00000190734863*x^3-7.62939453125E-06*x^2+2*x+13.0000305175781"
NormResiduals=2.21203107924332E-05
【
Mathematica求插值多项式的函数为InterpolatingPolynomial ，其调用格式为：
InterpolatingPolynomial[数据表,变量]
已知自变量x=1, 2 , 3, 4时，因变量y=16, 25, 46, 85，求一个三次插值多项式逼近该函数，并求y(2.5).
解 In[1]：=data={{1,16},{2,25},{3,46},{4,85}}
In[2]：=f=InterpolatingPolynomial[data,x]
In[3]：=f=Expand[%]
In[4]：=f/.x→2.5
Out[1]={{1, 16}, {2, 25}, {3, 46}, {4, 85}}
Out[2]=16 + (–1 + x) (9 + (–2 + x) (3 + x))
Out[3]=13 + 2 x + x3
Out[4]=33.625
】

功能:采用Steffensen算法求解一元非线性方程根
f="sin(x)"
f = "sin(x)"
x=SolveRootSteffensen(f)
x =
[ 3.14159265358979 1.2246063538E-16 ]
【
Mathematica求方程sinx=0在x=3附近的一个根.
解
In[1]：=FindRoot[Sin[x]= =0,{x,3}]
Out [1]={x →3.14159}
】

功能:求矩阵的秩
格式:rank(a);rank(a,n)
说明:a为矩阵变量,n为误差纠正值[1-18],建议n=10
a=[1 2 3 4
5 6 7 8
9 10 11 12
13 14 15 16
-17 18 19 20]
a =
[ 1.00000000000000 2.00000000000000 3.00000000000000 4.00000000000000
5.00000000000000 6.00000000000000 7.00000000000000 8.00000000000000
9.00000000000000 10.0000000000000 11.0000000000000 12.0000000000000
13.0000000000000 14.0000000000000 15.0000000000000 16.0000000000000
-17.0000000000000 18.0000000000000 19.0000000000000 20.0000000000000 ]
rank(a)
ans =
[ 3.00000000000000 ]
a=[1 5 9 13
17 2 6 10
14 18 3 7
11 15 19 4
8 12 16 20]
a =
[ 1.00000000000000 5.00000000000000 9.00000000000000 13.0000000000000
17.0000000000000 2.00000000000000 6.00000000000000 10.0000000000000
14.0000000000000 18.0000000000000 3.00000000000000 7.00000000000000
11.0000000000000 15.0000000000000 19.0000000000000 4.00000000000000
8.00000000000000 12.0000000000000 16.0000000000000 20.0000000000000 ]
rank(a)
ans =
[ 4.00000000000000 ]
b=[1 5 9 13 17
2 6 10 14 18
3 7 11 15 19
4 8 12 16 20]
b =
[ 1.00000000000000 5.00000000000000 9.00000000000000 13.0000000000000 17.0000000000000
2.00000000000000 6.00000000000000 10.0000000000000 14.0000000000000 18.0000000000000
3.00000000000000 7.00000000000000 11.0000000000000 15.0000000000000 19.0000000000000
4.00000000000000 8.00000000000000 12.0000000000000 16.0000000000000 20.0000000000000 ]
rank(b)
ans =
[ 2.00000000000000 ]
Java通过JNI计算矩阵的秩：
C:\>javac Foo.java
C:\>java Foo
ret45=3
ret54=4
public native static int brank(double[] a,int m,int n);
double a[]=new double[]{1,5,9,13,17,2,6,10,14,18,3,7,11,15,19,4,8,12,16,20};
int ret45=foo.brank(a,4,5);
System.out.println("ret45="+ret45);
int ret54=foo.brank(a,5,4);
System.out.println("ret54="+ret54);
JNIEXPORT jint JNICALL Java_Foo_brank(JNIEnv *env, jclass, jdoubleArray a, jint m, jint n)
{
jdouble* pa = (env)->GetDoubleArrayElements(a, 0 );
int ret=brank(pa,m,n);
(env)->ReleaseDoubleArrayElements(a, pa, 0 );
return ret;
}

Γ(x+yi)=？

Γ(1.0+0.1i)=0.9902-0.05702i

Γ(1.1+0.1i)=0.9959+0.0422i

Γ(1.2+0.1i)=0.912-0.02614i

Γ(1.3+0.1i)=0.8923-0.01492i

Γ(1.4+0.1i)=0.8827-0.005273i

Γ(1.5+0.1i)=0.8821+0.00334i

Γ(1.6+0.1i)=0.8896+0.01132i

Γ(1.7+0.1i)=0.9048+0.01896i

Γ(1.8+0.1i)=0.9276+0.02652i

Γ(1.9+0.1i)=0.9579+0.0342i

功能:伽马函数求值

格式:Gamma(x)

说明:函数成功返回对应的伽马函数值.但是注意,x不能为负整数.

Gamma(1.0+0.1i)

ans = 0.990206629588383 - 0.0568238087537157i

Gamma(1.1+0.1i)

ans = 0.943739327403351 - 0.0397232864251407i

Gamma(1.2+0.1i)

ans = 0.912006443791646 - 0.0261438920264667i

Gamma(1.3+0.1i)

ans = 0.892279466286237 - 0.0149202316065137i

Gamma(1.4+0.1i)

ans = 0.882717571841338 - 0.00527341904680261i

Gamma(1.5+0.1i)

ans = 0.882093197098701 + 0.00334036310197866i

Gamma(1.6+0.1i)

ans = 0.889620571044311 + 0.0113181814096069i

Gamma(1.7+0.1i)

ans = 0.904846700776421 + 0.018964119204293i

Gamma(1.8+0.1i)

ans = 0.927581717215409 + 0.0265233851521887i

Gamma(1.9+0.1i)

ans = 0.957855012678436 + 0.0342048829172188i

功能:求四元变量的自然对数
格式:QuaternionLog(A)
a=[1 2 3 4]
a =
[ 1.00000000000000 2.00000000000000 3.00000000000000 4.00000000000000 ]
QuaternionLog(a)
ans =
[ 1.70059869083108 0.51519029266408 0.77278543899612 1.03038058532817 ]
a=[3 2 0 0]
a =
[ 3.00000000000000 2.00000000000000 0.00000000000000 0.00000000000000 ]
QuaternionLog(a)
ans =
[ 1.28247467873077 0.58800260354756 0.00000000000000 0.00000000000000 ]
用VB_D3DXQuaternionLn计算：
Ln{1,2,3,4}={1,2,3,0}
Ln{3,2,0,0}={4.712389,3.141593,0,0}
复数3+2i的自然对数为：1.282475,0.588003

功能:求四元变量的指数次方
格式:QuaternionExp(A)
a=[1 2 3 4]
a =
[ 1.00000000000000 2.00000000000000 3.00000000000000 4.00000000000000 ]
QuaternionExp(a)
ans =
[ 1.69392272368330 -0.78955962454155 -1.18433943681234 -1.57911924908312 ]
a=[1 2.1 0 0]
a =
[ 1.00000000000000 2.10000000000000 0.00000000000000 0.00000000000000 ]
QuaternionExp(a)
ans =
[-1.37231399230213 2.34644633551727 0.00000000000000 0.00000000000000 ]
用VB_D3DXQuaternionExp计算：
Exp{1,2,3,4}={-.1509214,-.3018427,-.4527641,-.825299}
Exp{1,2.1,0,0}={.3130659,.6574383,0,-.6853938}
exp(1,2.1)= -1.37231+2.34645i

功能:求2个数的算术-几何平均数(Arithmetic - Geometric Mean)
MeanArithmeticGeometric(1.414214,1)
ans = 1.19814043607983
MeanArithmeticGeometric(1,1.414214)
ans = 1.19814043607983
MeanArithmeticGeometric(1,0)
ans = 2.77555756156289E-17
MeanArithmeticGeometric(1,0.707107)
ans = 0.847213204591294
【
//AGM(A0=1.414214+0.000000i;G0=1.000000+0.000000i)=1.198157+0.000000i或AGM(G0=1.000000+0.000000i;A0=1.414214+0.000000i)=1.198157+0.000000i,注意：A1=1.2071067811865475244008443621048，G1=1.1892071150027210667174999705605
//AGM(A0=1.000000+0.000000i;G0=0.000000+0.000000i)=0.000122+0.000000i，即AGM(1,0)=0=>K(1)=+∞
//AGM(A0=1.000000+0.000000i;G0=0.707107+0.000000i)=0.847225+0.000000i，即AGM(1,sqrt(2)/2)≈0.847=>K(sqrt(2)/2)≈1.854
】

双纽线周率~ω是Gauss（1799年）由研究算术几何平均数时发现的量，作为自守形式/椭圆函数论发端的事件是应该被纪念的。

高斯发现的是AGM(sqrt(2),1)=π/~ω这个等式。更一般地Gauss得到了AGM(a,b)的公式。

矩阵运算

Beta(0.25,0.5)

ans =

[ 5.24411510858474 ]

算式解析

5.24411510858474 /2

ans = 2.62205755429237

算式解析
功能:求Jacobi椭圆函数
格式:JacobiSn(u,k)
u:为任意的复数变量
k:必须是介于0与1之间的实数
说明:一般我们用Sn(u,k)来表示本函数.本函数相当于求第一类椭圆积分EllipticF函数的反函数.
JacobiSn(0.5,0.866025)
ans = 0.464426226873146
即sn(u=0.500000+0.000000i;k=0.866025+0.000000i)=0.466428+0.000000i

JacobiSn(0.5,0.866025)
ans = 0.464426226873146
JacobiSn(1+2i,0)
计算出现非数字的情况,请注意检查表达式是否出现除0的情形
JacobiSn(1+2i,0.5)
ans = 2.11350173673824 + 0.16995353310106i
【
//sn(u=0.500000+0.000000i;k=0.866025+0.000000i)=0.466428+0.000000i
//sn(u=1.000000+2.000000i;k=0.000000+0.000000i)=3.165779+1.959601i
】

功能:计算第一类完全椭圆积分(Complete Elliptic Integral)
格式:EllipticK(z)
说明:函数等价于Int(1/Sqrt(1-z^2*Sin(t)^2),0,pi/2),表示对1/sqrt(1-z^2*Sin(t)^2)函数中的变量t从0到pi/2进行定积分.
EllipticK(0.707)
ans = 1.85394676925858
即K(e=0.707)=1.853921
EllipticK(0.5)
ans = 1.6857503548126
即K'(e=0.707)=K(e=0.5)=1.685750

K(e=0.707)=1.853921
K'(e=0.707)=K(e=0.5)=1.685750

EllipticK(0.707)
ans = 1.85394676925858
EllipticK(0.5)
ans = 1.6857503548126
【
//K(k=0.707000+0.000000i)=1.853921+0.000000i
//K(k=0.500000+0.000000i)=1.685750+0.000000i
】

功能:计算第二类椭圆积分(Complete Elliptic Integral)
格式:EllipticE(z)、EllipticE(s,z)
说明:
1、第二类完全椭圆积分EllipticE(z)函数等价于Int(Sqrt(1-z^2*Sin(t)^2),0,pi/2),表示对Sqrt(1-z^2*Sin(t)^2)函数中的变量t从0到pi/2进行定积分.
2、第二类不完全椭圆积分函数EllipticE(s,z)=EllipticRF(c-1,c-z^2,c)-z^2/3*EllipticRD(c-1,c-z^2,c),其中c=1/sin(s)^2
EllipticE(0)
ans = 1.5707963267949
EllipticE(1)
ans = 0.785398163397448
EllipticE(0.707107)
ans = 1.3506437252615
【
椭圆周长积分表达式（第二类椭圆积分）
椭圆的周长是4aE(c/a)，这里的函数E是第二类完全椭圆积分。
E(0)=pi/2（1/4圆弧弧长）

E(1)=1（1/4椭圆弧长的下确界）
E(sqrt(2)/2)=pi^(3/2)Γ(1/4)^(-2)+Γ(1/4)^2/(8sqrt(2))（b=c时的1/4椭圆弧长）
】

功能:Simpson积分
格式:IntegralSimpson(f,x,a,b,n);IntegralSimpson(f,x,a,b);IntegralSimpson(f,a,b)
说明:f为符号函数句柄,x为未知数名称,积分区间[a,b],n为复化积分时的等分数默认为500
f="log(1.0+x)/(1.0+x^2)"
f = "log(1.0+x)/(1.0+x^2)"
IntegralSimpson(f,x,0.0,1.0)
ans =
[ 0.27219826128798 ]

数学变换与滤波中只有
功能:一维傅里叶逼近(Fourier approximation)
格式:
fourapp(f,n,a,b,x)
fourapp(f,n,x)
fourapp(f,n)
功能:(离散)傅里叶变换
格式:Dft(a,n);Dft(a)
功能:(离散)逆傅里叶变换
格式:IDft(a,n);IDft(a)
功能:快速傅里叶变换
格式:FFt(a,n);FFT(a)
功能:快速傅里叶变换逆变换
格式:IFFt(a,n);IFFT(a)

kfour 傅里叶级数逼近
kkfft 快速傅里叶变换
kkfwt 快速沃什变换
kkspt 快速三次平滑
klman 离散随机系统的卡尔曼滤波
kkabg α-β-γ滤波

功能:对方阵进行LU分解.
格式:lu(a)//a为方阵变量
说明:函数自动生成名称为l的下三角矩阵,和名称为u的上三角矩阵,其中u的对角线全为1.
a=[1.0 2.0 -1.0
-2.0 -5.0 3.0
-1.0 -3.0 0]
a =
[ 1.00000000000000 2.00000000000000 -1.00000000000000
-2.00000000000000 -5.00000000000000 3.00000000000000
-1.00000000000000 -3.00000000000000 0.00000000000000 ]
lu(a)
l =
[ 1.00000000000000 0.00000000000000 0.00000000000000
-2.00000000000000 -1.00000000000000 0.00000000000000
-1.00000000000000 -1.00000000000000 -2.00000000000000 ]
u =
[ 1.00000000000000 2.00000000000000 -1.00000000000000
0.00000000000000 1.00000000000000 -1.00000000000000
0.00000000000000 0.00000000000000 1.00000000000000 ]
a=[7 8 9
14 46 51
28 82 163]
a =
[ 7.00000000000000 8.00000000000000 9.00000000000000
14.0000000000000 46.0000000000000 51.0000000000000
28.0000000000000 82.0000000000000 163.000000000000 ]
lu(a)
l =
[ 7.00000000000000 0.00000000000000 0.00000000000000
14.0000000000000 30.0000000000000 0.00000000000000
28.0000000000000 50.0000000000000 72.0000000000000 ]
u =
[ 1.00000000000000 1.14285714285714 1.28571428571429
0.00000000000000 1.00000000000000 1.10000000000000
0.00000000000000 0.00000000000000 1.00000000000000 ]

功能:求方阵行列式.此方法精度高,但运算的方阵的阶数最好不要超过5.
格式:Det(a)
说明:a为方阵

功能:采用LU分解算法求方阵行列式.适合阶数较大的矩阵.
格式:Det2(a)
说明:a为方阵

功能:采用QR分解算法求行列式,算出的结果可能正负号相反,本函数适合大阶数运算,但一般针对非奇异方阵.
格式:Det4(a)
说明:a为方阵

a=[1.0 2.0 -1.0
-2.0 -5.0 3.0
-1.0 -3.0 0]
a =
[ 1.00000000000000 2.00000000000000 -1.00000000000000
-2.00000000000000 -5.00000000000000 3.00000000000000
-1.00000000000000 -3.00000000000000 0.00000000000000 ]
b=Det(a)
b =
[ 2.00000000000000 ]

b=Det2(a)
b =
[ 2.00000000000000 ]

b=Det4(a)
b =
[ 2.00000000000000 ]

功能:矩阵相乘
格式:mul(a,b)
a=[1 0 0
-2 1 0
-1 1 1]
a =
[ 1.00000000000000 0.00000000000000 0.00000000000000
-2.00000000000000 1.00000000000000 0.00000000000000
-1.00000000000000 1.00000000000000 1.00000000000000 ]
b=[1 2 -1
0 -1 1
0 0 -2]
b =
[ 1.00000000000000 2.00000000000000 -1.00000000000000
0.00000000000000 -1.00000000000000 1.00000000000000
0.00000000000000 0.00000000000000 -2.00000000000000 ]
c=mul(a,b)
c =
[ 1.00000000000000 2.00000000000000 -1.00000000000000
-2.00000000000000 -5.00000000000000 3.00000000000000
-1.00000000000000 -3.00000000000000 0.00000000000000 ]
【
矩阵相乘等数值计算　2009-02-27 23:55
Now print resource matrix b[3][3]=
1 0 0
-2 1 0
-1 1 1
Now print resource matrix c[3][3]=
1 2 -1
0 -1 1
0 0 -2
Now printm multiply results matrix a[3][3]=B*C:
1 2 -1
-2 -5 3
-1 -3 0
】

a=[3 0
1 -1]
a =
[ 3.00000000000000 0.00000000000000
1.00000000000000 -1.00000000000000 ]
eig(a)
目前此函数只支持求方阵的特征值。

a=[-4.26923076923077 -4.57692307692308
5.19230769230769 6.26923076923077]
a =
[-4.26923076923077 -4.57692307692308
5.19230769230769 6.26923076923077 ]
eig(a)
目前此函数只支持求方阵的特征值。

【
原矩阵为{{3,0},{1,-1}},特征值为λ_1=3,λ_2=-1,特征向量为V_1={{4},{1}},V_2={{0},{1}} 。

原矩阵为{{-4.26923076923077,-4.57692307692308},{5.19230769230769,6.26923076923077}},特征值为λ_1=3,λ_2=-1,特征向量为V_1={{-3.26923076923077},{5.19230769230769}},V_2={{-7.26923076923077},{5.19230769230769}} 。

原矩阵为{{-3.92307692307692,-7.30769230769231},{2.76923076923077,5.92307692307692}},特征值为λ_1=2.99999999999999,λ_2=-.999999999999989,特征向量为V_1={{-2.92307692307693},{2.76923076923077}},V_2={{-6.92307692307691},{2.76923076923077}} 。

A(x,y)=(3x,x-y)的特征值为λ_1=3,λ_2=-1，在基(I_11=1,I_21=0),(I_12=0,I_22=1)下的矩阵为A。
A在基(E_11,E_21),(E_12,E_22)下的矩阵为E^(-1)AE={{-51/13,-95/13},{36/13,77/13}}，其中E={{E_11,E_12},{E_21,E_22}}={{2,5},{5,6}}
A在基(H_11,H_21),(H_12,H_22)下的矩阵为H^(-1)AH={{17/3,20/3},{-8/3,-11/3}}，其中H={{H_11,H_12},{H_21,H_22}}={{8,8},{4,7}}
从基E到基H的过渡矩阵E^(-1)H，其逆矩阵为H^(-1)E。
H^(-1)E(E^(-1)AE)E^(-1)H=H^(-1)AH
E^(-1)H={{-2.153846,-1},{2.4615,2}}={{-28/13,-1},{32/13,2}}
H^(-1)E={{-13/12,-13/24},{4/3,7/6}}
H^(-1)AH={{17/3,20/3},{-8/3,-11/3}}

原矩阵为{{5.66666666666667,6.66666666666667},{-2.66666666666667,-3.66666666666667}},特征值为λ_1=3,λ_2=-1,特征向量为V_1={{6.66666666666667},{-2.66666666666667}},V_2={{2.66666666666667},{-2.66666666666667}} 。
】

该软件目前还不支持求空间三点确定的三角形面积。
x=[1 2 3
2 5 3
6 1 8]
x =
[ 1.00000000000000 2.00000000000000 3.00000000000000
2.00000000000000 5.00000000000000 3.00000000000000
6.00000000000000 1.00000000000000 8.00000000000000 ]
PgArea(x)
PgArea函数接收的参数必须是行数为2的矩阵。

功能:求一元四次以内多项式的精确解
格式:Polyroot(f)
说明:f为符号句柄或者多项式的系数
f="x^3+x+1"
f=
"x^3+x+1"
polyroot(f)
方程:x^3+x^1+1=0 的根如下
x(1)= 0.34116390191401 - 1.16154139999725i
x(2)= 0.34116390191401 + 1.16154139999725i
x(3)= -0.682327803828019

ans =
[ 0.34116390191401 -1.1615413999972
0.34116390191401 1.16154139999725
-0.6823278038280 0.00000000000000 ]

chap2
1.空间平面对应于一个三元一次方程Ax+By+Cz+D=0。反之，任意一个三元一次方程也对应于空间中的一个平面。如果平面α的方程是Ax+By+Cz+D=0，其含义是平面α上任意动点(x,y,z)都是Ax+By+Cz+D=0的解。而Ax+By+Cz+D=0的每一组解也对应于α上某一点。按：2004.4.21求三点确定的平面的例子：①求解或求证：由A1(1,1,1),A2(1,1,-1),A3(1,-1,1)三点确定的平面α为x=1(z,y∈R)。解：设平面α:k1x+k2y+k3y=k，则线性方程组为:{{A1_x,A2_x,A3_x},{A1_y,A2_y,A3_y},{A1_z,A2_z,A3_z}}|{{k},{k},{k}}，得到k2=k3=0,k1=k，即x=1。②求由B1(2,-1,1),B2(1,2,-1),B3(0,5,-1)三点确定的平面β。解1：~B1B2=(-1,3,-2)，~B1B3=(-2,6,-2)，~B2B3=(-1,3,0)，B1与~B1B2、~B1B3的正交解空间为：取(A,B,C)∈{c(3,1,0)}，任取(x,y,z)∈β,则3(x-2)+1(y+1)+0(z-1)=0，即平面β:3x+y=5(z∈R)。解2：设平面α:k1x+k2y+k3y=k，则线性方程组为:{{A1_x,A2_x,A3_x},{A1_y,A2_y,A3_y},{A1_z,A2_z,A3_z}}|{{k},{k},{k}}，得到k1=6k/10,k2=2k/10,k3=0k/10，即由{A1_x,A1_y,A1_z},{A2_x,A2_y,A2_z},{A3_x,A3_y,A3_z}确定的平面为3x+y=5。注:由{A1_x,A2_x,A3_x},{A1_y,A2_y,A3_y},{A1_z,A2_z,A3_z}确定的平面为9x-8y+7z=10。按：2004.4.14线性代数的欧氏几何意义的重新发现：三向量{A1_x,A1_y,A1_z},{A2_x,A2_y,A2_z},{A3_x,A3_y,A3_z}共面<=>|{{A1_x,A1_y,A1_z},{A2_x,A2_y,A2_z},{A3_x,A3_y,A3_z}}|=0<=>{A1_x,A1_y,A1_z},{A2_x,A2_y,A2_z},{A3_x,A3_y,A3_z}线性相关；二向量{A1_x,A1_y,A1_z},{A2_x,A2_y,A2_z}共线<=>|{{i,j,k},{A1_x,A1_y,A1_z},{A2_x,A2_y,A2_z}}|=0<=>b_i∝a_i,特别地,|{{i,j,k},{A1_x,A1_y,0},{A2_x,A2_y,0}}|=|{A1_x,A1_y},{A2_x,A2_y}}|·k=0<=>|{A1_x,A1_y},{A2_x,A2_y}}|=0。如,{3,1},{6,2}}、{{1,2},{2,4},{3,6}}均共线(平行),{{1,2},{2,4},{3,6}}不共线,{{1,2,1},{2,4,2},{3,6,4}}线性相关=>共面[按照这里的上下文,这里用矩阵表示几个向量,且用1个列向量表示1个向量]；{{1,2,3},{2,4,6},{1,2,4}}线性相关=>共面[按照这里的上下文,这里用矩阵表示几个向量,且用1个行向量表示1个向量]；({A1_x,A1_y,A1_z},{A2_x,A2_y,A2_z}}=0∩({A1_x,A1_y,A1_z},{A3_x,A3_y,A3_z}}=0=>{A2_x,A2_y,A2_z},{A3_x,A3_y,A3_z}共线。
2.欧氏几何公理的代数描述：
平面α:Ax+By+Cz+D=0上两点(x1,y1,z1),(x2,y2,z2)确定的直线为l:{(x,y(x),z(x))}
对于任意x∈l,Ax+By(x)+Cz(x)+D=Ax+D+[B(y2-y1)+C(z2-Z1)](x-x1)/(x2-x1)+By1+Cz1=Ax+D-A(x-x1)-Ax1-D=0，x∈α。
即公理1:如果一条直线上有两点在一个平面内,则该直线上所有点都在该平面内,也称直线在平面内。
公理2的解析描述：平面α:A1x+B1y+C1z+D1=0∩平面β:A2x+B2y+C2z+D2=0若有解则必有无穷解。
公理3的解析描述：当(x2-x1,y2-y1,z2-z1)与(x3-x2,y3-y2,z3-z2)线性无关即|{{i,j,k},{x2-x1,y2-y1,z2-z1},{x3-x2,y3-y2,z3-z2}}|≠0时，通过此三点的平面为…
3.R^3中曲面的切平面存在定理（2010.9.27）以及用二元实函显函数表示的曲面的切面计算公式（2004.5.13）：设R^3中曲面Σ的方程为F(x,y,z)=0。若函数F(x,y,z)在点X_0(x_0,y_0,z_0)∈Σ处可微，且函数F(x,y,z)在点X_0处的各一阶偏导数不全为0，则曲面Σ在X_0有切平面存在，其方程为F'_x(X_0)(x-x_0)+F'_y(X_0)(y-y_0)+F'_z(X_0)(z-z_0)=0。按：曲面z=f(x,y)在(x_0,y_0,f(x_0,y_0))处的切面为z=f'_x(x_0,y_0)(x-x_0)+f'_y(x_0,y_0)(y-y_0)+f(x_0,y_0)。例如，z=sqrt(1-x^2-y^2)在(0.5,0.5,sqrt(2)/2)处的切面为z=-(sqrt(2)/2)x-(sqrt(2)/2)y+sqrt(2)。
考虑Σ为平面的情形，R^3中的平面α:Ax+By+Cz+D=0在点X_0(x_0,y_0,z_0)∈α(<=>Ax_0+By_0+Cz_0+D=0)处存在切平面，且切平面为A(x_x_0)+B(y-y_0)+C(z-z_0)=Ax+By+Cz+D=0，即为平面α自身。
4.在切平面存在的情况下又可以根据如下定理来求曲面的法线方程和单位法向量m：
设R^3中曲面Σ的方程为F(x,y,z)=0。若函数F(x,y,z)在点X_0(x_0,y_0,z_0)∈Σ处可微，且函数F(x,y,z)在点X_0处的各一阶偏导数不全为0，则曲面Σ在X_0的法线方程为(x-x_0)/F'_x(x_0,y_0,z_0)=(y-y_0)/F'_y(x_0,y_0,z_0)=(z-z_0)/F'_z(x_0,y_0,z_0)，~n=±(δ(y,z)/δ(u,v),δ(z,x)/δ(u,v),δ(x,y)/δ(u,v))|X_0。
考虑Σ为平面的情形（Σ=α），R^3中的平面α:Ax+By+Cz+D=0在点X_0(x_0,y_0,z_0)∈α处的法线方程为(x-x_0)/A=(y-y_0)/B=(z-z_0)/C。
2011.5.20:取△X_3X_2X_1{<}α，利用三维叉积和克莱姆法则，则α的定向用它的法向量N=(N_x,N_y,N_z)=(X_2X_1×X_2X_3)=(X_1-X_2)×(X_3-X_2)=(y1z3+y2z1+y3z2-y1z2-y2z3-y3z1,x1z2+x2z3+x3z1-x1z3-x2z1-x3z2,x1y3+x2y1+x3y2-x1y2-x2y3-x3y1)=(X_2-X_1)×(X_2-X_3)=(△/D)(A,B,C)或法向量(A,B,C)或法向量(-A,-B-C)表示，其中N_x=y1z3+y2z1+y3z2-y1z2-y2z3-y3z1，N_y=x1z2+x2z3+x3z1-x1z3-x2z1-x3z2，N_z=x1y3+x2y1+x3y2-x1y2-x2y3-x3y1，△=|{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3}}是AX_i+BY_i+CZ_i+D=0的系数行列式，A,B,C可以由x_i,y_i,z_i,D表示出来，A=△_A/△=(D/△)N_x，B=△_B/△=(D/△)N_y，C=△_C/△=(D/△)N_z。平面α的单位法向量即m=(A,B,C)/|(A,B,C)|。
若取法向量为N=(A,B,C)，X_0∈α，则D=-N·X_0。则平面方程可以用向量的形式表示：N·X_0=-D。
如果Ax+By+Cz+D>0<=>N·X>-D，则点X=(x,y,z)在平面α的外侧（正面，N或m指向同一侧）；如果Ax+By+Cz+D<0<=>N·X<-D，则点X=(x,y,z)在平面α的内侧（反面，与N或m的方向相反）。例如，对于数学中的右手系和平面α:F(x,y,z)=z=0，i,j∈α,m=k=i×j=(0,0,1),因m·(x,y,2)>-D=0，m·(x,y,-3)<-D=0，所以点X=(x,y,2)在平面α的外侧，X=(x,y,-3)在平面α的内侧。对于D3D中的左手系(i,j,-k),m=-k=-i×j=(0,0,1)(-k),结果不变。
又如，由切点X_0=(1/2,1/2,sqrt(2)/2)和切平面方程F(x,y,z)=(x-1/2)+(y-1/2)+sqrt(2)(z-sqrt(2)/2)=0<=>x+y+sqrt(2)z-2=0就可以求出法线方程x-1/2=y-1/2=(1/sqrt(2))(z-sqrt(2)/2)，取单位主法向量m=(1/2,1/2,sqrt(2)/2)时就规定了一个定向。
5.面面角的计算：①η_1:9x+4y+7z=0与η_2:9x+8y+4z=0的面面角为<(9,4,7),(9,8,4)>=arccos(141/(sqrt(146)sqrt(161)))=0.403569281（合23.12281654°）。②面x+y-z=0与面z=0的面面角等于向量(1,1,-1)与向量(0,0,-1)夹arccos(1/(sqrt(3)-1))。
6.面面距的计算（4.23引理）：平行平面η_1:k_1x+k_2y+k_3z=-D_1与η_2:k_1x+k_2y+k_3z=-D_2的面面距为d(η_1,η_2)=|D_2-D_1|/sqrt((k_1)^2+(k_2)^2+(k_3)^2)。
7.R^n(n=2,3)中其他几何量的计算（>=4.14）：
长度|α-β|称为向量α和β的距离[其实就是欧氏空间R^n的弱抽象--线性度量空间R^n中的非齐次坐标表示的两个点的距离]，记为d(α,β)。
点面距的计算：①点(-2,2,0)到平面18x+26y+11z-42=0的距离为26/sqrt(1121)。②已知截距为a,b,c,则平面:x/a+y/b+z/c=1(或k1x+k2y+k3z=k)，原点到其距离d=|abc|/sqrt(a^2·b^2+a^2·c^2+b^2·c^2)(或d=k/sqrt(k1^2+k2^2+k3^2)) 。③已知x+2y+3z=12,求证x^2+2y^2+3z^2≥24(=|x=2,y=2,z=2)<=>……x+sqrt(2)y+sqrt(3)z=12……x^2+y^2+z^2≥24。利用点面距即可证得：原点到平面的距离为12/sqrt(1+2+3)=2sqrt(6)(=|(x,y,z)=(2,2sqrt(2),2sqrt(3)))。四月最值不等式：x+sqrt(e)y+sqrt(f)z=g=>x^2+y^2+z^2≥g^2/(1+e+f)(=|(x,y,z)=g(1,sqrt(e),sqrt(f))/(1+e+f))<=>x+ey+fz=g=>x^2+ey^2+fz^2≥g^2/(1+e+f)(=|(x,y,z)=g(1,1,1)/(1+e+f))。
四面体体积的计算：(1,1/2,1),(-1/3,1,2),(2/3,2,-2),(-2,2,0)四点确定的四面体体积为：(1/6)|{{1,1/2,1,1},{-1/3,1,2,1},{2/3,2,-2,1},{-2,2,0,1}}|=-13/18，(1/6)|{{-4/3,1/2,1},{-1/3,3/2,-3},{-3,3/2,-1}}|=13/18。
三角形面积的计算：(1,1/2,1),(-1/3,1,2),(2/3,2,-2)三点确定的三角形面积为：(1/2)‖{{e_1,e_2,e_3},{-4/3,1/2,1},{-1/3,3/2,-3}}‖=sqrt(1121)/12。
又如：A=(1,0,1),B=(1,2,-1),O=(0,0,0)确定的平面AOB为z=x-y，S_△AOB=(1/2)‖{{e_1,e_2,e_3},{1,0,1},{1,2,-1}}‖=(1/2)‖(-2,2,2)‖=sqrt(3)。
A=(1,0,1),B=(1,2,-1),O=(0,0,0),C=(-1,2,0)四点确定的四面体体积为：V_四面体=(1/3)S_△AOB*CH=(1/3)sqrt(3)*sqrt(3)=1，其中高CH=sqrt(3)。
A=(1,0,1),B=(1,2,-1),O=(0,0,0),C=(-1,2,0)四点确定的平行六面体体积为：V_平行六面体=(OA,OB,OC)=|A,B,C|=|{{1,1,-1},{0,2,2},{1,-1,0}}|=6。

'定义向量u=(u_x,u_y,u_z)和向量v=(v_x,v_y,v_z)的叉积运算为：u×v=(u_yv_z-u_zv_y,u_zv_x-u_xv_z,u_xv_y-u_yv_x)
function cross(u1,u2,u3,v1,v2,v3)
cross1=u2*v3-u3v2
cross2=u3*v1-u1v3
cross3=u1*v2-u2v1
cross="(" & cross1 & "," & cross2 & "," & cross3 & ")"
end function
'i×j = k,j×k = i,k×i = j
msgbox "i×j=" & cross(1,0,0,0,1,0)
msgbox "j×k=" & cross(0,1,0,0,0,1)
msgbox "k×i=" & cross(0,0,1,1,0,0)

function Arccos(X)
Arccos=Atn(-X/Sqr(-X*X+1))+2*Atn(1)
end function
'向量u=(u_x,u_y,u_z)和向量v=(v_x,v_y,v_z)的夹角θ=arccos[(u,v)/(‖u‖·‖v‖)],0<=θ<=pi,
function angle(u1,u2,u3,v1,v2,v3)
X=(u1*v1+u2*v2+u3*v3)/(sqr(u1*u1+u2*u2+u3*u3)*sqr(v1*v1+v2*v2+v3*v3))
angle=Arccos(X)
end function
msgbox angle(1,1,-1,0,0,-1)
msgbox Arccos(sqr(3)/3)

chap3
2014-6-21
上界O——大O记法
下界Ω——大Ω记法
在讨论排序时，不好的排序算法是O(N^2)，或称为二次的。
第7章排序
我们对内部排序的考查将指出：
存在几种容易的算法以O(N^2)排序，如插入排序。
有一种算法叫希尔排序，它编程简单，以o(N^2)（小o记法）运行。
还有一些稍微复杂的O(NlogN)的排序算法。
任何通用的排序算法均需要Ω(NlogN)次比较。
排序还是一种能够对其进行精确的算法的范例。
基于比较的排序：<和>操作符、赋值运算。
在STL中，排序是通过使用函数模板sort来完成的。
sort的参数反映了一个容器（的范围）的头尾标志，以及一个可选的比较器。
迭代器必须只是随机访问。sort算法不能保证相等的项保持它们原始的次序（如果这很重要的话，可以使用stable_sort来代替sort）。
sort使用例子：
非降序排列
非升序排列
前半部分按非降序排列
插入排序为O(N^2)，平均情形也是Θ(N^2)（大Θ记法）。
定理7.1：N个互异元素的数组的平均逆序数是N(N-1)/4。
定理7.2：通过交换相邻元素进行排序的任何算法平均需要Ω(N^2)时间。
快速排序的平均运行时间是O(NlogN)，它的最坏情形的性能为O(N^2)。

chap4
首先到Eigen官网上下载Eigen源码包，下载后解压完直接放到自己平时软件所在的目录下（例如E:\EigenRoot，并在VS2010的附加包含目录下指定该路径即可），不需要安装。Eigen下载地址为
http://eigen.tuxfamily.org/index.php?title=Main_Page
http://eigen.tuxfamily.org/dox/unsupported/group__MatrixFunctions__Module.html#title9

使用Eigen的第一个例子：
#include <iostream>
#include <vector>
#include <Eigen/Eigen>
using namespace Eigen;
using namespace std;
int main()
{
/*
v1 =
5
6
result:
20 25
24 30
*/
Eigen::Vector2d v1, v2; //Eigen中的变量
v1 << 5, 6; //默认的向量为列向量
cout << "v1 = " << endl << v1 << endl;
v2 << 4, 5 ;
Matrix2d result = v1*v2.transpose();
cout << "result: " << endl << result << endl;
system("pause");
return 0;
}

第二个例子，计算矩阵的基本初等函数：
#include <iostream>
#if 0
#include <vector>
#include <Eigen/Eigen>
#endif
#include <unsupported/Eigen/MatrixFunctions>
using namespace Eigen;
using namespace std;
int main()
{
/*
v1 =
5
6
result:
20 25
24 30
*/
#if 0
Eigen::Vector2d v1, v2; //Eigen中的变量
v1 << 5, 6; //默认的向量为列向量
cout << "v1 = " << endl << v1 << endl;
v2 << 4, 5 ;
Matrix2d result = v1*v2.transpose();
cout << "result: " << endl << result << endl;
#endif
/*
http://eigen.tuxfamily.org/dox/unsupported/group__MatrixFunctions__Module.html
The matrix A is:
0 -0.785398 0
0.785398 0 0
0 0 0
The matrix exponential of A is:
0.707107 -0.707107 0
0.707107 0.707107 0
0 0 1
对上面的矩阵A求正弦的话会出现如下断言失败
Assertion failed: "Taylor series does not converge" && 0, file e:\eigenroot\unsupported\eigen\src\matrixfunctions\matrixfunctionatomic.h, line 89
也就是泰勒级数sinA=A - A^3/3! + A^5/5! -A^7/7!.....对某些矩阵A不收敛，这一点跟实数域和复数域的情况不一样。
cos(A) =
1.32461 0 0
0 1.32461 0
0 0 1
*/
#if 1
{
const double pi = std::acos(-1.0);
MatrixXd A(3,3);
A << 0, -pi/4, 0,
pi/4, 0, 0,
0, 0, 0;
std::cout << "The matrix A is:\n" << A << "\n\n";
std::cout << "The matrix exponential of A is:\n" << A.exp() << "\n\n";
//A = MatrixXd::Random(3,3);
//MatrixXd sinA = A.sin();
//std::cout << "sin(A) = \n" << sinA << "\n\n";
MatrixXd cosA = A.cos();
std::cout << "cos(A) = \n" << cosA << "\n\n";
// The matrix functions satisfy sin^2(A) + cos^2(A) = I,like the scalar functions.
//std::cout << "sin^2(A) + cos^2(A) = \n" << sinA*sinA + cosA*cosA << "\n\n";
}
#endif
system("pause");
return 0;
}

第三个例子：
/*
验证：
{{-3,1,-1},{-7,5,-1},{-6,6,-2}}^2={{8,-4,4},{-8,12,4},{-12,12,4}}
{{-3,1,-1},{-7,5,-1},{-6,6,-2}}^3={{-20,12,-12},{-84,76,-12},{-72,72,-8}}
矩阵代数是不可交换的，然而是结合的，至今研究过的大部分系统都是结合的。约定+只用于可交换的系统，而×不必受此限制。
m1:
-3 1 -1
-7 5 -1
-6 6 -2
m2:
8 -4 4
-8 12 4
-12 12 4
m3:
-20 12 -12
-84 76 -12
-72 72 -8
*/
#if 1
MatrixXf m1(3,3);
m1<<-3,1,-1,-7,5,-1,-6,6,-2;
cout << "m1: " << endl << m1 << endl;
MatrixXf m2 = m1*m1;
cout << "m2: " << endl << m2 << endl;
MatrixXf m3 = m1*m1*m1;
cout << "m3: " << endl << m3 << endl;
#endif

第四个例子：
由矩阵的幂和无穷级数导出的等式：
sin{{4,6,0},{-3,-5,0},{-3,-6,1}}={{-2,0,-1},{1,0,1},{0,1,1}}{{sin1,0,0},{0,sin1,0},{0,0,-sin2}}{{-1,-1,0},{-1,-2,1},{1,2,0}}
cos{{2,0,0},{1,2,-1},{1,0,1}}发散，真的吗？
20131220注意：
{{-2,0,-1},{1,0,1},{0,1,1}}的逆矩阵是{{-1,-1,0},{-1,-2,1},{1,2,0}}，而不是{{-1,-1,1},{-1,-2,2},{0,1,0}}。
A:
4 6 0
-3 -5 0
-3 -6 1
B:
-2 0 -1
1 0 1
0 1 1
C:
0.841471 0 0
0 0.841471 0
0 0 -0.909297
D:
-1 -1 0
-1 -2 1
1 2 0
sin(A) =
2.59224 3.50154 -1.35828e-008
-1.75077 -2.66007 2.92391e-007
-1.75077 -3.50154 0.841471
BCD:
2.59224 3.50154 0
-1.75077 -2.66007 0
-1.75077 -3.50154 0.841471
#if 1
MatrixXf A(3,3);
A<<4,6,0,-3,-5,0,-3,-6,1;
cout << "A: " << endl << A << endl;
MatrixXf B(3,3);
B<<-2,0,-1,1,0,1,0,1,1;
cout << "B: " << endl << B << endl;
MatrixXf C(3,3);
C<<sinf(1),0,0,0,sinf(1),0,0,0,-sinf(2);
cout << "C: " << endl << C << endl;
#if 0
MatrixXf D = B.inverse();
#else
MatrixXf D(3,3);
D<<-1,-1,0,-1,-2,1,1,2,0;
#endif
cout << "D: " << endl << D << endl;
MatrixXf sinA = A.sin();
std::cout << "sin(A) = \n" << sinA << "\n\n";
MatrixXf BCD = B*C*D;
cout << "BCD: " << endl << BCD << endl;
#endif

A:
2 0 0
1 2 -1
1 0 1
cosA:
-0.416147 -9.20451e-008 7.45058e-008
-0.956449 -0.416147 0.956449
-0.956449 -1.97194e-007 0.540302
#if 1
MatrixXf A(3,3);
A<<2,0,0,1,2,-1,1,0,1;
cout << "A: " << endl << A << endl;
MatrixXf cosA = A.cos();
cout << "cosA: " << endl << cosA << endl;
#endif

chap5
20150420添加：Matlab画3维空间中的2阶、3阶贝塞尔曲线

函数文件bezier2.m：
function bezier2(p0,p1,p2)
t=0:0.001:1;
x=(p2(1)-2*p1(1)+p0(1))*t.^2+2*(p1(1)-p0(1))*t+p0(1);
y=(p2(2)-2*p1(2)+p0(2))*t.^2+2*(p1(2)-p0(2))*t+p0(2);
z=(p2(3)-2*p1(3)+p0(3))*t.^2+2*(p1(3)-p0(3))*t+p0(3);
plot3(x,y,z,'r');

执行：
>> bezier2([1 3 0],[4 18 0],[7 6 0])
>> bezier2([1 3 0],[4 18 0],[7 6 2])
函数文件bezier3.m：
function bezier3(p0,p1,p2,p3)
t=0:0.001:1;
x=(1-t).^3*p0(1)+3*t.*(1-t).^2*p1(1)+3*t.^2.*(1-t)*p2(1)+t.^3*p3(1);
y=(1-t).^3*p0(2)+3*t.*(1-t).^2*p1(2)+3*t.^2.*(1-t)*p2(2)+t.^3*p3(2);
z=(1-t).^3*p0(3)+3*t.*(1-t).^2*p1(3)+3*t.^2.*(1-t)*p2(3)+t.^3*p3(3);
plot3(x,y,z,'r')
执行：
>> bezier3([0 3 0],[5 20 0],[7 2 0],[9 1 0])
>> bezier3([0 3 0],[5 20 0],[7 2 0],[9 1 2]) %这个三阶贝塞尔曲线应该是一条挠曲线

20150403-20150412补充
用到的函数文件mo.m
function rchamo=mo(r) %求向量r的模
sum=0;
for i=1:3
sum=sum+r(i)^2;
end
rchamo=simple(sum^(1/2));

求曲线在任一点处的曲率与挠率
空间曲线曲率与挠率的计算：
求曲线r=(t,t^2/2,t^3/3)在点t=1上的曲率向量k和曲率|k|。
>> syms t;r=[t t^2/2 t^3/3];r1=diff(r)

r1 =

[ 1, t, t^2]

>> r1mo=mo(r1)

r1mo =

(1+t^2+t^4)^(1/2)

>> T=r1/r1mo

T =

[ 1/(1+t^2+t^4)^(1/2), t/(1+t^2+t^4)^(1/2), t^2/(1+t^2+t^4)^(1/2)]

>> T1=diff(T)

T1 =

[ -1/2/(1+t^2+t^4)^(3/2)*(2*t+4*t^3), 1/(1+t^2+t^4)^(1/2)-1/2*t/(1+t^2+t^4)^(3/2)*(2*t+4*t^3), 2*t/(1+t^2+t^4)^(1/2)-1/2*t^2/(1+t^2+t^4)^(3/2)*(2*t+4*t^3)]

>> k=T1/r1mo

k =

[ -1/2/(1+t^2+t^4)^2*(2*t+4*t^3), (1/(1+t^2+t^4)^(1/2)-1/2*t/(1+t^2+t^4)^(3/2)*(2*t+4*t^3))/(1+t^2+t^4)^(1/2), (2*t/(1+t^2+t^4)^(1/2)-1/2*t^2/(1+t^2+t^4)^(3/2)*(2*t+4*t^3))/(1+t^2+t^4)^(1/2)]

>> t=1;kzhi=eval(k)

kzhi =

-0.3333 0.0000 0.3333 %曲率向量k=(-1/3,0,1/3)
>> kzhimo=sqrt(kzhi*kzhi')

kzhimo =

0.4714 %曲率|k|=sqrt(2)/3
求曲线r=(3t-t^3,3t^2,3t+t^3)的曲率和挠率。
>> syms t;r=[3*t-t^3 3*t^2 3*t+t^3];r1=diff(r)

r1 =

[ 3-3*t^2, 6*t, 3+3*t^2]

>> r2=diff(r1)

r2 =

[ -6*t, 6, 6*t]

>> r1mo=mo(r1)

r1mo =

3*2^(1/2)*(1+t^2)

>> r1r2=cross(r1,r2)

r1r2 =

[ 18*t^2-18, -6*(3+3*t^2)*t-6*(3-3*t^2)*t, 18+18*t^2]

>> r1r2mo=mo(r1r2)

r1r2mo =

18*2^(1/2)*(1+t^2)

>> k=r1r2mo/r1mo^3

k =

1/3/(1+t^2)^2

>> r3=diff(r2)

r3 =

[ -6, 0, 6]

>> r123=dot(cross(r1,r2),r3)

r123 =

216

>> tau=r123/r1r2mo^2

tau =

1/3/(1+t^2)^2 %注意，沿此曲线，|k|=tau=1/(3*(1+t^2)^2)

用符号积分求解二重积分∫(-1,2)∫(y*y,y+2)xydxdy
>> syms x y;%定义两个符号变量
>> a=int(int(x*y,x,y*y,y+2),y,-1,2);%积分限x：y*y,y+2 ,y:,-1,2
>> b=simple(a);%化简
>> c=vpa(b,4)%得到4位近似解，也可以任意N位解

c =

5.625
a =

45/8
b =

45/8
>> int(int(x*y,x,y*y,y+2),y,-1,2)

ans =

45/8
dblquad计算二重数值积分
>> f=@(x,y)x.*y;
>> dblquad(f,y*y,y+2,-1,2)
出错
>> dblquad(f,pi/4,1,2,4)

ans =

1.1494

多重积分的MATLAB实现
http://wenku.baidu.com/link?url=hk9pQjn2em9epQpLXaKreqU3FwFKI7SB1g5-NE8klWTE2ywsOIY2u8s4pDy4w8LplEg_EFgtuza2dwZEu9_sGZn4bm59UVljxEkDSrEC-6m

将积分分为定积分、二重积分、三重积分、第一型曲线积分、第二型曲线积分、第一型曲面积分和第二型曲面积分等七种积分
例1计算二重积分，积分对象是(x^2+y^2-x)dxdy=(x^2+y^2-x)dx∧dy，积分区域是由直线y=2,y=x,y=2x所围成的闭区域
>> int(int(x^2+y^2-x,x,y/2,y),y,0,2)

ans =

13/6
例9计算三重积分，积分对象是xyzdxdydz=xyzdx∧dy∧dz，积分区域是球面x^2+y^2+z^2=1及三个坐标面所围成的在第一象限内的区域
>> syms x y z;int(int(int(x*y*z,z,0,sqrt(1-x^2-y^2)),y,0,sqrt(1-x^2)),x,0,1)

ans =

1/48
计算定积分，函数y=x^2在(1,2)区间内的投影面积
>> f=@(x)x.^2;quad(f,1,2) %算子（泛函）quad(f,a,b)表示在区间(a,b)内对x求积分

ans =

2.3333
MATLAB积分
http://wenku.baidu.com/link?url=vyGVVVoPFxXtd-CLC7dwdqlX6S140Q0gAhDNnUwU3i6Ilm57aR7PnKrr5Ja6uiymFInaof6h0hzA1vhyiZCJGVx1nYzFqBFwGANTTz_gX7K
triplequad计算三重数值积分
>> f=@(x,y,z)(y*sin(x)+z*cos(x));triplequad(f,0,pi,0,1,-1,1)

ans =

2.0000
实验七+曲线、曲面积分及其应用
http://wenku.baidu.com/link?url=NZJ54kMV82K117znxE62kN5UprabK758Xj0GfkH--NKU_aTeVYrLl-TYKIy0xXUIUQ_vWCaboouputuLfiRE7dPMh1S3FqG8HF48G9ufzlW
用第二型曲线积分求椭圆的面积（利用格林公式计算平面图形的面积）
>> syms x y t dx dy a b
>> x=a*cos(t);y=b*sin(t);
>> dx=diff(x);dy=diff(y);
>> int(x*dy-y*dx,t,0,2*pi)

ans =

2*a*b*pi

>> (1/2)*int(x*dy-y*dx,t,0,2*pi)

ans =

a*b*pi
——椭圆x^2/a^2+y^2/b^2=1的面积为abpi

用第一型曲线积分求椭圆的周长
>> ds=sqrt(dx^2+dy^2)

ds =

(a^2*sin(t)^2+b^2*cos(t)^2)^(1/2)

>> int(ds,t,0,2*pi)

ans =

4*EllipticE((-(-a^2+b^2)/a^2)^(1/2))*(b^2/a^2)^(1/2)*a^2/(b^2)^(1/2)
>> simple(ans)

ans =

4*EllipticE(1/a*(a^2-b^2)^(1/2))*a

雅可比用q来表示模和周期，例如
k=4sqrt(q)[(1+q^2)(1+q^4)(1+q^6)…/((1+q)(1+q^3)(1+q^5)…)]^4。
k=kq(0.207879576350762)=0.985171431009416
>> syms kq1 q;kq1=maple('product(((1+q^(2*n))/(1+q^(2*n-1))),n=1..100)');q=0.207879576350762;kq=eval(kq1);k=4*sqrt(q)*kq^4

k =

0.9852

chap6
2015-11-20 17:07添加：
http://networkx.github.io/documentation/latest/
NetworkX安装：
以管理员身份在命令行下输入
C:\Python27\Scripts>easy_install.exe C:\networkx-1.10.tar.gz
会出现：
……
Installed c:\python27\lib\site-packages\decorator-4.0.4-py2.7.egg
Finished processing dependencies for networkx==1.10
打开C:\Python27\python.exe，输入import networkx as nx、print nx，就有相关显示，如果没有成功，会提示不认识networkx类
Python 2.7.9 (default, Dec 10 2014, 12:24:55) [MSC v.1500 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import networkx as nx
>>> print nx
<module 'networkx' from 'C:\Python27\lib\site-packages\networkx-1.10-py2.7.egg\networkx\__init__.pyc'>

定义：在有向图G=<V，E>中，对任意结点v∈V，以v为始结点的弧的条数，称为结点v的出度，记为d+(v)；以v为终结点的弧的条数，称为v的入度，记作d-(v)；结点v的出度与入度之和，称为结点的度数，记为d(v)，显然d(v)=d+(v)+d-(v)。
>>> DG= nx.DiGraph()
>>> DG.add_node(1)
>>> DG.add_node(2)
>>> DG.add_nodes_from([3,4,5,6])
>>> DG.add_cycle([1,2,3,4])
>>> DG.add_edge(1,3)
>>> DG.add_edges_from([(3,5),(3,6),(6,7)])
>>> nx.degree(DG)
{1: 3, 2: 2, 3: 5, 4: 2, 5: 1, 6: 2, 7: 1}
>>> nx.degree_histogram(DG)
[0, 2, 3, 1, 0, 1]
>>> DG.number_of_edges()
8
>>> len(DG)
7
>>> DG.adj
{1: {2: {}, 3: {}}, 2: {3: {}}, 3: {4: {}, 5: {}, 6: {}}, 4: {1: {}}, 5: {}, 6:{7: {}}, 7: {}}
对于无向图G=<V，E>，结点v∈V的度数等于联结它的边数，也记为d(v)。若v点有环，规定该点度因环而增加2。
>>> G=nx.Graph()
>>> G=nx.Graph(name='G1')
>>> e=[(1,2),(1,3),(2,3),(3,4),(4,5)]
>>> G=nx.Graph(e)
>>> nx.degree(G)
{1: 2, 2: 2, 3: 3, 4: 2, 5: 1}
>>> nx.degree_histogram(G)
[0, 1, 3, 1]
对于无向图G=<V，E>，记Δ(G)或Δ=max{d(v)|v∈V}，δ(G)或δ=min{d(v)|v∈V}，它们分别称为图G的最大度和最小度。
这里图G的最大度是3，最小度是1。
关于无向图中的结点的度，欧拉给出一个定理，这是图论中的第一个定理。
定理：给定无向图G=<V，E>，则∑_v∈V(d(v))=2|E|
定理：在任何无向图中，奇度结点的数目为偶数。
>>> G.number_of_edges()
5
>>> len(G) #返回G中节点数目
5
>>> G.adj
{1: {2: {}, 3: {}}, 2: {1: {}, 3: {}}, 3: {1: {}, 2: {}, 4: {}}, 4: {3: {}, 5: {}}, 5: {4: {}}}
>>> 1 in G
True
>>> 6 in G
False

>>> nx.density(G)
0.5
>>> nx.degree_centrality(G)
{1: 0.5, 2: 0.5, 3: 0.75, 4: 0.5, 5: 0.25}
>>> nx.betweenness_centrality(G)
{1: 0.0, 2: 0.0, 3: 0.6666666666666666, 4: 0.5, 5: 0.0}
>>> nx.triangles(G)
{1: 1, 2: 1, 3: 1, 4: 0, 5: 0}
>>> nx.transitivity(G)
0.5
>>> nx.clustering(G)
{1: 1.0, 2: 1.0, 3: 0.3333333333333333, 4: 0.0, 5: 0.0}

chap7
20151024添加：
GAP内建命令LogFFE( z, r )求离散对数
生成有限域F_2[x]/(x^5+x^3+1)，已知f=x^5+x^3+1在F_2上本原。设a为f的一个根，计算a^2+a+1对a的离散对数。
gap> F2:=GF(2);x:=Indeterminate(F2,"x");f:=x^5+x^3+1;F:=GF(F2,f);a:=PrimitiveRoot(F);Print(LogFFE(a^2+a+1,a),"\n");
GF(2)
x
x^5+x^3+Z(2)^0
GF(2^5)
Z(2^5)
11
a^2+a+1对a的离散对数为11。
在有限域F_2[x]/(x^7+x^6+1)中，已知f=x^7+x^6+1在F_2上本原。设a为f的一个根，求a^3+a+1对a的离散对数。
gap> F2:=GF(2);x:=Indeterminate(F2, "x");f:=x^7+x^6+1;F:=GF(F2,f);a:=PrimitiveRoot(F);Print(LogFFE(a^3+a+1,a),"\n");
GF(2)
x
x^7+x^6+Z(2)^0
GF(2^7)
Z(2^7)
31
有两种有限域的构造方法：
一个是从某个素域通过单扩张得到有限域。定理：一个有限域一定是它所含素域的单扩域。
另一个是通过添加某个多项式的根得到分裂域。定理：设F_p是含有p个元素的素域，q=p^n(n>=1),那么多项式f(x)=x^q-x在F_q上的分裂域是一个含有q个元素的有限域。

求多项式x^2+x+1在F_16中的子域，零元，单位，本原元。
gap> F16:=GF(2,4);x:=Indeterminate(F16,"x");f:=x^2+x+1;F:=GF(F16,f);PrimitiveRoot(F);Zero(f);One(f);Characteristic(f);Order(f);Inverse(f);DegreeFFE(f);
GF(2^4)
x
x^2+x+Z(2)^0
AsField( GF(2^4), GF(2^8) )
Z(2^8)
0*Z(2)
Z(2)^0
2

模n的原根有许多，但模n的最小原根是唯一的. 在GAP中有两个函数可以很方便地计算模n的最小原根：IntFFE(Z(p));和IntFFE(PrimitiveRoot(GF(p)))。
模23的最小原根等于5。
gap> IntFFE(Z(23));
5
gap> IntFFE(PrimitiveRoot(GF(23)));
5

Fermat小定理：若p是素数，a和p互素，则a^p≡a(modp)。
Fermat小定理可用于素性测试，这样就可以测试F5是合数还是素数。
gap> F0:=2^(2^0)+1;F1:=2^(2^1)+1;F2:=2^(2^2)+1;F3:=2^(2^3)+1;F4:=2^(2^4)+1;2^F0 mod F0;2^F1 mod F1;2^F2 mod F2;2^F3 mod F3;2^F4 mod F4;
3
5
17
257
65537
2
2
2
2
2
gap> 3^F0 mod F0;3^F1 mod F1;3^F2 mod F2;3^F3 mod F3;3^F4 mod F4;
0
3
3
3
3
gap> F5:=2^(2^5)+1;
4294967297
gap> PowerMod(2,F0,F0);PowerMod(2,F1,F1);PowerMod(2,F2,F2);PowerMod(2,F3,F3);PowerMod(2,F4,F4);PowerMod(2,F5,F5);
2
2
2
2
2
2
gap> PowerMod(3,F0,F0);PowerMod(3,F1,F1);PowerMod(3,F2,F2);PowerMod(3,F3,F3);PowerMod(3,F4,F4);PowerMod(3,F5,F5);
0
3
3
3
3
497143886
gap> F5 mod 3;
2
所以F_5不是素数。事实上
gap> Factors(F5);
[ 641, 6700417 ]
gap> PowerMod(4,F0,F0);PowerMod(4,F1,F1);PowerMod(4,F2,F2);PowerMod(4,F3,F3);PowerMod(4,F4,F4);PowerMod(4,F5,F5);
1
4
4
4
4
4

第9个费马数是合数
gap> F9:=2^(2^9)+1;
1340780792994259709957402499820584612747936582059239337772356144372176403007354697680187429816690342769003185818648605\
0853753882811946569946433649006084097
gap> PowerMod(2,F9,F9);
2
gap> PowerMod(3,F9,F9);
1334675769651040385428187301046783395533879238453972810857773932975597967623327836963094381826193117863973778281140461\
7970328133474715265143960685411614642
gap> Factors(F9);
[ 2424833, 7455602825647884208337395736200454918783366342657,
741640062627530801524787141901937474059940781097519023905821316144415759504705008092818711693940737 ]

用欧拉的二次不可约多项式x^2+x+41计算x从10到50的数，查看是否全部为素数；
gap> for i in [10..50] do x:=i;y:=x^2+x+41;Print(y);Print(":");Print(IsPrime(y));Print(Factors(y));Print("\n"); od;
151:true[ 151 ]
173:true[ 173 ]
197:true[ 197 ]
223:true[ 223 ]
251:true[ 251 ]
281:true[ 281 ]
313:true[ 313 ]
347:true[ 347 ]
383:true[ 383 ]
421:true[ 421 ]
461:true[ 461 ]
503:true[ 503 ]
547:true[ 547 ]
593:true[ 593 ]
641:true[ 641 ]
691:true[ 691 ]
743:true[ 743 ]
797:true[ 797 ]
853:true[ 853 ]
911:true[ 911 ]
971:true[ 971 ]
1033:true[ 1033 ]
1097:true[ 1097 ]
1163:true[ 1163 ]
1231:true[ 1231 ]
1301:true[ 1301 ]
1373:true[ 1373 ]
1447:true[ 1447 ]
1523:true[ 1523 ]
1601:true[ 1601 ]
1681:false[ 41, 41 ]
1763:false[ 41, 43 ]
1847:true[ 1847 ]
1933:true[ 1933 ]
2021:false[ 43, 47 ]
2111:true[ 2111 ]
2203:true[ 2203 ]
2297:true[ 2297 ]
2393:true[ 2393 ]
2491:false[ 47, 53 ]
2591:true[ 2591 ]

二次互反律中引入勒让德符号：
Lengendre符号L(n/p)=1，当存在整数x满足x的平方模p余n；反之为-1。
gap> IsPrime(593);
true
克罗内克-雅可比符号J(438,593)=-1即同余式x^2≡286(mod563)无解，n=438是模数m2=593的二次非剩余。
设n=0或1(mod4),n非平方数且m>0，则Kronecker符号J(n,m)定义为：
若m整除n，即m|n,则J(n,m)=0,
若n=1(mod8),则J(n,2)=1
若n=5(mod8),则J(n,2)=-1
若m为奇素数且m不整除n，则J(n,m)转化为勒让德符号L(n/m)。
Jacobi符号Jacobi(n,m)=J(n/m)（1837年），n>2且为奇整数,取值范围为{0,1,-1}，是勒让德符号的一种推广，但是根据雅可比符号的值不能判断同余式是否有解；勒让德符号L(n/p)在GAP4中的另一种推广Legendre(n,m)=L(n/m)能准确判断二次同余方程x^2≡n(mod m)一定是否有解。
gap> n:=11;;m:=35;;Legendre(n,m);Jacobi(n,m);#9^2 = 11 mod 35
1//准确的判断是否有解
1//雅可比符号的值
gap> n:=-1;;m:=9;;Legendre(n,m);Jacobi(n,m);#x^2 = -1 mod 3,9均无解
-1//准确的判断是否有解
1//雅可比符号的值
gap> n:=2;;m:=3599;;Legendre(n,m);Jacobi(n,m);#x^2 = 2 mod 3599无解
-1//准确的判断是否有解
1//雅可比符号的值
gap> n:=3;;m:=35;;Legendre(n,m);Jacobi(n,m);#x^2 = 3 mod 35无解
-1//准确的判断是否有解
1//雅可比符号的值
gap> n1:=286;;n2:=438;;m1:=563;;m2:=593;;IsPrime(m1);Legendre(n1,m1);Jacobi(n1,m1);IsPrime(m2);Legendre(n2,m2);Jacobi(n2,m2);
true
-1
-1
true
-1
-1
即使当m=p时，GAP4中的这两个符号还是有细微区别的。
gap> M:=[3,5,7,11,13,17];for i in M do Print(Legendre(-3,i));Print(",");Print(Jacobi(-3,i));Print("\n"); od;
[ 3, 5, 7, 11, 13, 17 ]
1,0
-1,-1
1,1
-1,-1
1,1
-1,-1
gap> n:=5;;M:=[3,5,7,11,13,17];;for i in M do Print(Legendre(n,i));Print(",");Print(Jacobi(n,i));Print("\n"); od;
-1,-1
1,0
-1,-1
1,1
-1,-1
-1,-1
gap> n:=-7;;M:=[3,5,7,11,13,17];;for i in M do Print(Legendre(n,i));Print(",");Print(Jacobi(n,i));Print("\n"); od;
-1,-1
-1,-1
1,0
1,1
-1,-1
-1,-1
gap> n:=-11;;M:=[3,5,7,11,13,17];;for i in M do Print(Legendre(n,i));Print(",");Print(Jacobi(n,i));Print("\n"); od;
1,1
1,1
-1,-1
1,0
-1,-1
-1,-1
模m是2或合数的情形：
gap> n:=13;;M:=[2,4,6,8,9,10,12,14,15,16];;for i in M do Print(Legendre(n,i));Print(",");Print(Jacobi(n,i));Print("\n"); od;
1,-1
1,1
1,-1
-1,-1
1,1
-1,1
1,1
-1,1
-1,-1
-1,1
gap> n:=-3;;M:=[2,4,6,8,9,10,12,14,15,16];;for i in M do Print(Legendre(n,i));Print(",");Print(Jacobi(n,i));Print("\n"); od;
1,-1
1,1
1,0
-1,-1
-1,0
-1,1
1,0
1,-1
-1,0
-1,1
gap> m:=8;;for n in [5] do Print("n=");Print(n);Print(",m=8:");Print(Legendre(n,m));Print(",");Print(Jacobi(n,m));Print("\n"); od;
n=5,m=8:-1,-1
gap> m:=5;;for n in [8] do Print("n=");Print(n);Print(",m=5:");Print(Legendre(n,m));Print(",");Print(Jacobi(n,m));Print("\n"); od;
n=8,m=5:-1,-1
http://math.fau.edu/richman/jacobi.htm
5,45都是奇数，该网站求得的雅可比符号值为：
(8/5) = (3/5) = (2/3) = -(1/3) = -1
(11/45) = (1/11) = 1
Not a residue
请输入一个素数 p=5 和一个与它互素的正整数 m=8 :
勒让德符号 (m=8 / p=5 ) 的值为 -1
gap> m:=45;;for n in [1..60] do Print("n=");Print(n);Print(",m=45:");Print(Legendre(n,m));Print(",");Print(Jacobi(n,m));Print("\n"); od;
n=1,m=45:1,1
n=2,m=45:-1,-1
n=3,m=45:-1,0
n=4,m=45:1,1
n=5,m=45:-1,0
n=6,m=45:-1,0
n=7,m=45:-1,-1
n=8,m=45:-1,-1
n=9,m=45:1,0
n=10,m=45:1,0
n=11,m=45:-1,1
n=12,m=45:-1,0
n=13,m=45:-1,-1
n=14,m=45:-1,1
n=15,m=45:-1,0
n=16,m=45:1,1
n=17,m=45:-1,-1
n=18,m=45:-1,0
n=19,m=45:1,1
n=20,m=45:-1,0
n=21,m=45:-1,0
n=22,m=45:-1,-1
n=23,m=45:-1,-1
n=24,m=45:-1,0
n=25,m=45:1,0
n=26,m=45:-1,1
n=27,m=45:-1,0
n=28,m=45:-1,-1
n=29,m=45:-1,1
n=30,m=45:-1,0
n=31,m=45:1,1
n=32,m=45:-1,-1
n=33,m=45:-1,0
n=34,m=45:1,1
n=35,m=45:-1,0
n=36,m=45:1,0
n=37,m=45:-1,-1
n=38,m=45:-1,-1
n=39,m=45:-1,0
n=40,m=45:1,0
n=41,m=45:-1,1
n=42,m=45:-1,0
n=43,m=45:-1,-1
n=44,m=45:-1,1
n=45,m=45:1,0
n=46,m=45:1,1
n=47,m=45:-1,-1
n=48,m=45:-1,0
n=49,m=45:1,1
n=50,m=45:-1,0
n=51,m=45:-1,0
n=52,m=45:-1,-1
n=53,m=45:-1,-1
n=54,m=45:1,0
n=55,m=45:1,0
n=56,m=45:-1,1
n=57,m=45:-1,0
n=58,m=45:-1,-1
n=59,m=45:-1,1
n=60,m=45:-1,0

2015-10-12 20:32添加：
矩阵A的特征多项式为|xE-A|
gap> A:=[ [ 1, 0,0], [ 1,0, 1 ], [ 0,1, 0] ];y:=CharacteristicPolynomial(A);Factors(y);A1:=TransposedMat(A);y1:=CharacteristicPolynomial(A1);Factors(y1);
[ [ 1, 0, 0 ], [ 1, 0, 1 ], [ 0, 1, 0 ] ]
x_1^3-x_1^2-x_1+1
[ x_1-1, x_1-1, x_1+1 ]
[ [ 1, 1, 0 ], [ 0, 0, 1 ], [ 0, 1, 0 ] ]
x_1^3-x_1^2-x_1+1
[ x_1-1, x_1-1, x_1+1 ]
gap> A100:=A^100;
[ [ 1, 0, 0 ], [ 50, 1, 0 ], [ 50, 0, 1 ] ]

向量内积
gap> theta:=[2,0,-1];;t:=ScalarProduct(theta,theta);
5
gap> theta1:=[2,0,0,0,-1,-1];;t1:=ScalarProduct(theta1,theta1);
6
gap> theta1:=[2,0,0,0,-1,-1];;theta2:=[1,-1,-1,-1,1,1];;theta3:=[1,1,1,1,1,1];;t12:=ScalarProduct(theta1,theta2);t13:=ScalarProduct(theta1,theta3);t23:=ScalarProduct(theta2,theta3);
0
0
0
gap> theta:=[1,E(3),E(3)^2];;t:=ScalarProduct(theta,theta);
0

//线性卷积conv
//1.000000 6.000000 16.000000 26.000000 31.000000 20.000000
gap> a:=[1,2,3,4];;b:=[1,4,5];;ProductCoeffs(a,b);
[ 1, 6, 16, 26, 31, 20 ]

gap> lo:= LieObject( [ [ 1, 0 ], [ 0, 1 ] ] );m:=UnderlyingRingElement(lo);lo*lo;m*m;
LieObject( [ [ 1, 0 ], [ 0, 1 ] ] )
[ [ 1, 0 ], [ 0, 1 ] ]
LieObject( [ [ 0, 0 ], [ 0, 0 ] ] )
[ [ 1, 0 ], [ 0, 1 ] ]
三维叉积对加法的分配律，线性性和雅可比恒等式表明：具有向量加法和叉积的R^3构成了一个Lie代数。
a × (b × c) + b × (c × a) + c × (a × b) = 0
雅可比恒等式就是下列等式：[X,[Y,Z]]+[Y,[Z,X]]+[Z,[X,Y]]=0
Lie代数是满足雅可比恒等式的代数结构的一个主要例子。
注意，满足雅可比恒等式的代数结构不一定满足反交换律。

x^3/6+x^2/2+x+1无有理根
2x^3-5x^2+5x-3的有理根为3/2
x^3-6x^2+15x-14的有理根为2
x^3/6+x^2/2+x+1无重根
x^3-3*x-2有重根
gap> x:=Indeterminate(Rationals);RootsOfUPol(x^3/6+x^2/2+x+1);
x_1
[ ]
gap> x:=Indeterminate(Rationals);RootsOfUPol(2*x^3-5*x^2+5*x-3);
x_1
[ 3/2 ]
gap> x:=Indeterminate(Rationals);RootsOfUPol(x^3-6*x^2+15*x-14);
x_1
[ 2 ]
gap> x:=Indeterminate(Rationals,"x");;y:=x^3/6+x^2/2+x+1;y1:=Derivative(y);Gcd(y,y1);GcdOp(y,y1);RootsOfUPol(y);
1/6*x^3+1/2*x^2+x+1
1/2*x^2+x+1
1
1
[ ]
gap> x:=Indeterminate(Rationals,"x");;y:=x^3-3*x-2;y1:=Derivative(y);Gcd(y,y1);GcdOp(y,y1);RootsOfUPol(y);
x^3-3*x-2
3*x^2-3
x+1
x+1
[ 2, -1, -1 ]

有理数域上的一元多项式环模不可约多项式生成的理想得到的无限域
gap> R:=PolynomialRing(Rationals,1);x:=Indeterminate(Rationals);poly:=x^2+x+1;IsIrreducibleRingElement(R,poly);Degree(poly);Factors(poly);F:=FieldByPolynomial(poly);Size(F);IsField(F);IsPrimeField(F);IsFinite(F);IsRationalsPolynomialRing(F);
Rationals[x_1]
x
x^2+x+1
true
2
[ x^2+x+1 ]
<algebraic extension over the Rationals of degree 2>
infinity
true
false
false
false

问题：GAP求基域为有理数域的可约多项式（例如f(x)=x^4-5*x^2+6=0）的伽罗瓦群？
注意：GAP软件中多项式环的功能尚不完善
R:=PolynomialRing(Rationals,1);IsUniqueFactorizationRing(R);IsIntegralRing(R);IsEuclideanRing(R);IsIrreducibleRingElement(R,x^2-2);IsIrreducibleRingElement(R,x^4-10*x^2+1);IsIrreducibleRingElement(R,50-45*x-6*x^2+x^3);IsIrreducibleRingElement(R,x^4-5*x^2+6);
R:=PolynomialRing(Rationals,2);IsUniqueFactorizationRing(R);IsEuclideanRing(R);
R:=PolynomialRing(Integers,1);IsUniqueFactorizationRing(R);IsEuclideanRing(R);
R:=PolynomialRing(Integers,2);IsUniqueFactorizationRing(R);IsEuclideanRing(R);
R:=PolynomialRing(GaussianIntegers,1);IsUniqueFactorizationRing(R);
R:=PolynomialRing(GaussianIntegers,2);IsUniqueFactorizationRing(R);
R:=PolynomialRing(GaussianRationals,1);IsUniqueFactorizationRing(R);IsIntegralRing(R);IsEuclideanRing(R);IsIrreducibleRingElement(R,x^2-2);
R:=PolynomialRing(GaussianRationals,2);IsUniqueFactorizationRing(R);IsEuclideanRing(R);
R:=UnivariatePolynomialRing(Rationals,"x");IsUniqueFactorizationRing(R);IsIntegralRing(R);IsIrreducibleRingElement(R,x^2-2);
R:=UnivariatePolynomialRing(Integers,"x");IsUniqueFactorizationRing(R);
R:=UnivariatePolynomialRing(GaussianRationals,"x");IsUniqueFactorizationRing(R);IsIntegralRing(R);IsIrreducibleRingElement(R,x^2-2);
R:=UnivariatePolynomialRing(GaussianIntegers,"x");IsUniqueFactorizationRing(R);

求方程x^3-2x-5=0的根。
gap> x:=Indeterminate(Rationals);RootsOfUPol(x^3-2*x-5);
x_1
[ ]
求方程x^3+2x^2+10x-20的根。
gap> x:=Indeterminate(Rationals);RootsOfUPol(x^3+2*x^2+10*x-20);
x_1
[ ]
gap> x:=Indeterminate(Rationals);RootsOfUPol(50-45*x-6*x^2+x^3);
x_1
[ 10, 1, -5 ]
gap> poly:=x^3-2*x-5;IsIrreducibleRingElement(PolynomialRing(Rationals,1),poly);IsSeparablePolynomial(poly);Degree(poly);
x^3-2*x-5
true
true
3
gap> poly:=x^3+2*x^2+10*x-20;IsIrreducibleRingElement(PolynomialRing(Rationals,1),poly);IsSeparablePolynomial(poly);Degree(poly);
x^3+2*x^2+10*x-20
true
true
3
gap> poly:=50-45*x-6*x^2+x^3;IsIrreducibleRingElement(PolynomialRing(Rationals,1),poly);IsSeparablePolynomial(poly);Degree(poly);
x^3-6*x^2-45*x+50
false
true
3
gap> poly:=x^4-5*x^2+6;IsIrreducibleRingElement(PolynomialRing(Rationals,1),poly);IsSeparablePolynomial(poly);Degree(poly);Factors(poly);
x^4-5*x^2+6
false
true
4
[ x^2-3, x^2-2 ]
定义：设F为域，称F[x]中多项式f(x)是F上的可分多项式，如果它在F[x]中的任一不可约因子都无重根（允许f(x)有重根）。否则，称f(x)是F上的不可分多项式。
任一无重根的不可约多项式必是可分多项式，不可分多项式的任一倍式必是不可分多项式。
完全域F上的多项式必是可分多项式。因为特征数为0的域和有限域必是完全域，因此，只有在特征数为p的无限域上才可能存在不可分多项式和不可分元。
gap> poly:=x^2+x+1;IsIrreducibleRingElement(PolynomialRing(GF(2),1),poly);Degree(poly);Factors(poly);
x^2+x+Z(2)^0
true
2
[ x^2+x+Z(2)^0 ]

gap> f:= MinimalPolynomial( Rationals, Sqrt(2) );
x^2-2
sqrt(2)在Q上的最小多项式为x^2-2。
gap> f:= MinimalPolynomial( Rationals, Sqrt(-1) );
x^2+1
gap> f:= MinimalPolynomial( Rationals, Sqrt(2)+Sqrt(3) );
x^4-10*x^2+1
sqrt(2)+sqrt(3)在Q上的最小多项式为f(x)=x^4-10x^2+1

http://www.gap-system.org/Manuals/pkg/MatricesForHomalg/doc/chap3.html
gap> LoadPackage("RingsForHomalg");
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Loading AutoDoc 2014.08.21 (Generate documentation from GAP source code)
by Sebastian Gutsche (http://wwwb.math.rwth-aachen.de/~gutsche/) and
Max Horn (http://www.quendi.de/math).
Homepage: http://wwwb.math.rwth-aachen.de/~gutsche/gap_packages/AutoDoc
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Loading ToolsForHomalg 2014.12.08 (Provides special methods and knowledge propagation tools)
by Mohamed Barakat (http://www.mathematik.uni-kl.de/~barakat/),
Sebastian Gutsche (http://wwwb.math.rwth-aachen.de/~gutsche/), and
Markus Lange-Hegermann (http://wwwb.math.rwth-aachen.de/~markus/).
Homepage: http://wwwb.math.rwth-aachen.de/~gutsche/gap_packages/ToolsForHomalg
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Loading MatricesForHomalg 2015.02.03 (Lazy evaluated matrices with clever operations for the homalg project)
by Mohamed Barakat (http://www.mathematik.uni-kl.de/~barakat/),
Markus Lange-Hegermann (http://wwwb.math.rwth-aachen.de/~markus/),
Martin Leuner (http://wwwb.math.rwth-aachen.de/Mitarbeiter/leuner.php), and
Vinay Wagh (http://www.iitg.ernet.in/vinay.wagh/).
Homepage: http://homalg.math.rwth-aachen.de/index.php/core-packages/matricesforhomalg
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Loading HomalgToCAS 2015.02.03 (A window to the outer world)
by Thomas B?chler (http://wwwb.math.rwth-aachen.de/~thomas/),
Mohamed Barakat (http://www.mathematik.uni-kl.de/~barakat/),
Simon G?rtzen (http://wwwb.math.rwth-aachen.de/~simon/),
Sebastian Gutsche (http://wwwb.math.rwth-aachen.de/~gutsche/), and
Vinay Wagh (http://www.iitg.ernet.in/vinay.wagh/).
Homepage: http://homalg.math.rwth-aachen.de/index.php/core-packages/homalgtocas
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Loading RingsForHomalg 2015.02.10 (Dictionaries of external rings for the homalg project)
by Mohamed Barakat (http://www.mathematik.uni-kl.de/~barakat/),
Simon Goertzen (http://wwwb.math.rwth-aachen.de/goertzen/),
Markus Kirschmer (http://www.math.rwth-aachen.de/~Markus.Kirschmer/),
Markus Lange-Hegermann (http://wwwb.math.rwth-aachen.de/~markus/),
Oleksandr Motsak (http://www.mathematik.uni-kl.de/~motsak/),
Daniel Robertz (http://wwwb.math.rwth-aachen.de/~daniel/),
Hans Sch?nemann (http://www.mathematik.uni-kl.de/~hannes/),
Andreas Steenpa? (), and
Vinay Wagh (http://www.iitg.ernet.in/vinay.wagh/).
Homepage: http://homalg.math.rwth-aachen.de/index.php/core-packages/ringsforhomalg
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
true
gap> ZZ:=HomalgRingOfIntegers();
Z
gap> IsHomalgRing(ZZ);
true
gap> IsHomalgRing(Integers);
false
gap> IsPreHomalgRing(ZZ);
false
gap> IsPreHomalgRing(Integers);
false
gap> Display( ZZ );
<An internal ring>
gap> Z256 := ZZ / 2^8;
Z/( 256 )
gap> Display( Z256 );
<A residue class ring>
gap> Z2 := Z256 / 6;
Z/( 256, 6 )
gap> BasisOfRows( MatrixOfRelations( Z2 ) );
<An unevaluated non-zero 1 x 1 matrix over an internal ring>
gap> Z2;
Z/( 2 )
gap> Display( Z2 );
<A residue class ring>
gap> F2 := HomalgRingOfIntegers()/2;
Z/( 2 )
gap> Z6 := HomalgRingOfIntegers()/6;
Z/( 6 )
gap> Z4:= HomalgRingOfIntegers()/4;
Z/( 4 )

gap> IsZero(ZZ);
false
gap> ContainsAField(ZZ);ContainsAField(F2);
false
true
gap> IsRationalsForHomalg(ZZ);IsRationalsForHomalg(F2);
false
false
gap> IsFieldForHomalg(ZZ);IsFieldForHomalg(F2);IsFieldForHomalg(Z6);IsFieldForHomalg(Z4);
false
true
false
false
gap> IsDivisionRingForHomalg(ZZ);IsDivisionRingForHomalg(F2);IsDivisionRingForHomalg(Z6);
false
true
false
gap> IsResidueClassRingOfTheIntegers(ZZ);IsResidueClassRingOfTheIntegers(F2);IsResidueClassRingOfTheIntegers(Z6);IsResidueClassRingOfTheIntegers(Z4);
true
true
true
true
gap> IsBezoutRing(ZZ);
true
gap> IsIntegrallyClosedDomain(ZZ);
true
gap> IsUniqueFactorizationDomain(ZZ);IsUniqueFactorizationDomain(F2);
true
true
gap> IsKaplanskyHermite(ZZ);
true
gap> IsDedekindDomain(ZZ);IsDedekindDomain(F2);IsDedekindDomain(Z6);
true
true
false=>主理想环PIR和主理想整环PID是不同的
每个主理想整环PID都是戴德金整环。
代数数论（algebraic number theory）中的最基本事实是:
一个代数数域K中的代数整数全体（称为K的整数环）/满足理想唯一分解条件的整环是一个戴德金整环（Dedekind Domain）[是整数环Z的弱抽象，是整环的强抽象]。
每个理想可以唯一表示为素理想的乘积。
理想全体是素理想生成的乘法半群，其分式理想全体是一个乘法群，称为K的理想群。
正是这种“理想的唯一素分解”可部分弥补“复数一般不能唯一素因子分解”的不足，在历史上使代数数论发展起来。
理想的唯一分解定理：代数整数环中每一个非平凡理想都可以唯一分解为有限多个素理想的乘积。
gap> IsIntegralDomain(ZZ);IsIntegralDomain(F2);IsIntegralDomain(Z6);
true
true
false
gap> IsNoetherian(ZZ);IsNoetherian(F2);IsNoetherian(Z6);
true
true
true=>诺特环不一定要求是整环
交换代数中的希尔伯特基定理（1888年）：若R为诺特环，则R[x]亦然。
[一个系数环上的任何多个变量的多项式所成的环，当这系数的环有一个单位元素和一组有限基时，这多项式环本身也有一组有限基。]
gap> IsLeftNoetherian(ZZ);IsLeftNoetherian(Z6);
true
true
gap> IsRightNoetherian(ZZ);IsRightNoetherian(Z6);
true
true
gap> IsArtinian(ZZ);IsArtinian(F2);
false
整数环是非阿廷环的诺特环
true
gap> IsLeftArtinian(ZZ);
false
gap> IsRightArtinian(ZZ);
false
gap> IsPrincipalIdealRing(ZZ);IsPrincipalIdealRing(F2);IsPrincipalIdealRing(Z6);
true
true
true
证明：整数环Z的每个理想都是主理想。
证明：模m剩余类环的每个理想都是主理想。=>主理想环PIR不一定要求是整环
gap> IsLeftPrincipalIdealRing(ZZ);
true
gap> IsRightPrincipalIdealRing(ZZ);
true
gap> IsRegular(ZZ);IsRegular(F2);
true
true
冯·诺伊曼引入的正则环(regular ring)->戴德金环

同调代数http://www.docin.com/p-67295980.html
gap> LoadPackage("RingsForHomalg");
true
gap> ZZ:=HomalgRingOfIntegers();
Z
gap> M := HomalgMatrix( "[ 2, 3, 4, 5, 6, 7 ]", 2, 3, ZZ );
<A 2 x 3 matrix over an internal ring>
gap> LoadPackage("Modules");
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Loading Modules 2015.02.10 (A homalg based package for the Abelian category of finitely presented modules over computable rings)
by Thomas B?chler (http://wwwb.math.rwth-aachen.de/~thomas/),
Mohamed Barakat (http://www.mathematik.uni-kl.de/~barakat/),
Florian Diebold (),
Sebastian Gutsche (http://wwwb.math.rwth-aachen.de/~gutsche/), and
Markus Lange-Hegermann (http://wwwb.math.rwth-aachen.de/~markus/).
Homepage: http://homalg.math.rwth-aachen.de/index.php/core-packages/modules
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
true
gap> M := LeftPresentation( M );
<A non-torsion left module presented by 2 relations for 3 generators>
gap> Display( M );
[ [ 2, 3, 4 ],
[ 5, 6, 7 ] ]

Cokernel of the map

Z^(1x2) --> Z^(1x3),

currently represented by the above matrix
gap> N := TorsionObject( M );
<A cyclic torsion left module presented by yet unknown relations for a cyclic generator>
gap> iota := TorsionObjectEmb( M );
<A monomorphism of left modules>
gap> psi := Ext( 1, iota, N );
<A homomorphism of right modules>
gap> ByASmallerPresentation( psi );
<A non-zero homomorphism of right modules>
gap> Display( psi );
[ [ 2 ] ]

the map is currently represented by the above 1 x 1 matrix
gap> extNN := Range( psi );
<A non-zero cyclic torsion right module on a cyclic generator satisfying 1 relation>
gap> IsIdenticalObj( extNN, Ext( 1, N, N ) );
true
gap> extMN := Source( psi );
<A non-zero cyclic torsion right module on a cyclic generator satisfying 1 relation>
gap> IsIdenticalObj( extMN, Ext( 1, M, N ) );
true
gap> Display( extNN );
Z/< 3 >
gap> Display( extMN );
Z/< 3 >

附录1：数学英语
多复变函数论theory of analytic functions of several variables
Ricci张量或Einstein张量
m阶逆变，n阶协变的张量
Riemann-Christoffel张量
代数闭包（algebraic closure）
最小多项式或极小多项式（minimal polynomial）
代数整数（algebraic integer）
代数数域(algebraic number field)
本原多项式（primitive polynomial）
数域P上的多项式函数（polynomial function）
n元多项式环P[x_1,x_2,...,x_n]
最大公因式（highest common factor）
公因式（common factor）
辅助方程式(Auxiliary Equation)

主元素(Identity Element)
逆元素(Inverse Element)
结合律(Associative Law)
置换群(Substitution Group)
伽罗华函数(Galois Function)
伽罗华分解式(Galois Resolvent)
立方倍积(The Duplication of the Cube)
元数(Order，即是元素的个数)
不变约群(Invarient Subgroup)
变形(Transform)
极大不变真约群(Maximal Invarient Proper Subgroup)
组合因数(Composition Factors)
可解群(Solvable Group)
循环群(Cyclic Group)
正置换群(Regular Substitution Group)
循环正置换群(Regular Cyclic Group)
元素(Element)
运算(Operation)
置换(Substitution)
群(Group)

巴比伦人(Babylonians)
可约的(reducible,就是说可以分解因数)
不可约的(irreducible)
数域(Field)
实数域(Field of Real Numbers)
复数域(Field of Complex Numbers)
可以用根式解(Solvable by Radicals)
约群(Subgroup)
不动置换(Identity Substitution)

群论(Theory of Groups)

调和harmony

几何geometric

算术arithmetic

平方square

连续函数（Continuous function）：开集的原像是开集的函数。该定义和epsilon-delta语言给出的定义是等价的。
连通集合（Connected set）
(Path connected)：集合中任意两点都存在连续路径相连。
中值定理(Intermediate Value Theorem)
极值定理(Extreme Value Theory)
一致收敛定理(Uniform Convergence Theorem)
基本群(Fundamental Group)
紧集（Compact set）
有界数列必然存在收敛子列
实数空间中有界闭集是紧的
紧集的任意开覆盖存在有限子覆盖
紧集中的数列必存在收敛子列
Compactness

数学英语
多复变函数论theory of analytic functions of several variables
Ricci张量或Einstein张量
m阶逆变，n阶协变的张量
Riemann-Christoffel张量
代数闭包（algebraic closure）
最小多项式或极小多项式（minimal polynomial）
代数整数（algebraic integer）
代数数域(algebraic number field)
本原多项式（primitive polynomial）
数域P上的多项式函数（polynomial function）
n元多项式环P[x_1,x_2,...,x_n]
最大公因式（highest common factor）
公因式（common factor）
辅助方程式(Auxiliary Equation)

群论(Theory of Groups)

调和harmony

几何geometric

算术arithmetic

平方square

三次曲线cubic curve
三角级数trigonometric series
一一映射bijection，bijective mapping
一致收敛序列uniformly convergent sequence
二元二次型binary quadratic form
二次曲面quadric，quadratic surface
二重积分double integral
子域subfield
子集subset
子群subgroup
平面一般线性群general linear group of the plane GL(2,R)
无穷集infinite set
开集open set
正规子群normal subgroup
本原解primitive solution
平坦环面flat torus
右陪集right coset
左陪集left coset
全同构holohedral isomorphism
自同构automorphism
同胚homeomorphism
同态homomorphism
同构isomorphism
四色问题four colour problem
四元数quaternion
可解群soluble group
可数集countable set
多面体polyhedron
多项式polynomial
交换群commutative group
交换环commutative ring
仿射群affine group
仿射变换affine transformation
满射surjection
圆锥曲线conic
p进赋值p-adic valuation
p进绝对值p-adic absolute
p进距离p-adic distance