无心人 发表于 2023-11-22 16:15:04

PARI/GP的一些笔记

在我没写完前,所有回复帖子一律删除
1、通用设置
1)提高默认堆栈容量
执行目录下有个gprc.txt
加入parisizemax = 100000000
可明显减少对大数进行复杂运算时候报堆栈错误的可能

2) 默认实数的精度
交互环境下输入\p 1000 可提高默认的实数精度,
也可以在gprc.txt里增加 realprecision = 1000
实际的精度要比输入值略大,以符合机器字的整数倍数
推荐gprc.txt默认精度不要太大,否则一次显示不全,也没必要
需要时候临时用\p 10000这种调节

3)启用计时器
交互环境下,输入#回车即可启用计时器,每次执行命令结果都会附加执行时间
再输入一次#即可关闭
##显示最后一次的执行时间

4)交互环境下,执行结果会存入%开始后跟数字的变量里
GP>x=1
%1=1
GP>y=2
%2=2
这些变量可以在输入中直接使用,
特别的单独的%代表最后一次的值

5)在线帮助
输入?会显示帮助列表,? n显示对应主题帮助
? function 会显示function函数的帮助
?? function 是更详细的function帮助
??? function 是所有和function有关联的函数清单
?\ 是显示快捷键的列表,比如#,##等
?. 是显示相关类型的成员函数


2、常见的常量
Pi 圆周率
oo 正无穷
-oo 负无穷
I 虚数单位i
Euler 欧拉常数
Catalan 卡塔兰常数

3、类型
(1)t_INT    : Z
无限精度整数,0x开头可以定有16进制数,0b开头定义二进制数
例如:12345678, -234567, 0x1E2365, 0xe3, 0b10100111
当输入整数时候,不带小数点,否则类型是实数
GP > type(1.)
%1 = "t_REAL"
GP > type(1)
%1 = "t_INT"

(2)t_REAL   : R
实数分三部分
符号s:+1, 0, -1
底数部分m:一个多精度整数,精度由系统变量realprecision控制
指数部分e:根据机器字长决定的有符号整数,范围[-2^B, 2^B], 字长32位时候B=31,64时候B=63
表示为实际值 s*m*2^e
例如:1.0, 23.56, 3.14159, 2.3E-18, 23e29

(3)t_INTMOD : Z/nZ
Mod(n, m)定义n模m,其中n, m必须是整数,可以为负数,m不为0
Mod(n, -m)=Mod(n, m)
Mod(-n, m)=Mod(m-n, m)
通常的+-*均可用于Mod(n, m)
但是/,如果m不是素数可能没有结果而保存,同样倒数也可能报错

gp > x=Mod(3,5)
%1 = Mod(3, 5)
gp > y=Mod(2,5)
%2 = Mod(2, 5)
gp > x*y
%3 = Mod(1, 5)
gp > x/y
%4 = Mod(4, 5)
gp > x+y
%5 = Mod(0, 5)

gp > x=Mod(4,8)
%6 = Mod(4, 8)
gp > y=Mod(5,8)
%7 = Mod(5, 8)
gp > x*y
%8 = Mod(4, 8)
gp > x+y
%9 = Mod(1, 8)
gp > x/y
%10 = Mod(4, 8)
gp > y/x
    ***   at top-level: y/x
    ***                  ^--
    *** _/_: impossible inverse in Fl_inv: Mod(4, 8).

(4)t_FRAC   : Q
有理数中分母不是1的形式
分数定义为 n/m,其中n, m都是整数,且m不为0
若m整除n,则n/m为t_INT类型
分数运算总会化简到最简形式,这会用到GCD算法,
通常,把分数运算转换到整数运算会提高运行速度

gp > 2/3
%1 = 2/3
gp > 2^(2/3)
%2 = 1.5874010519681994747517056392723082603914933278998530098082857618252165056242191732735442132622209570229347616813220179034976599
gp > 27^(2/3)
%3 = 9
gp > 3/4+5/7
%4 = 41/28
gp > 1/2-1/3
%5 = 1/6

(5)t_FFELT: Fp
有限域类型

(6)t_COMPLEX: C
复数可以表示成 x + y * I,x和y可以是t_INT, t_REAL, t_FRAC, t_INTMOD, t_PADIC

gp > x=Mod(2,65537)+Mod(1,65537)*I
%1 = Mod(2, 65537) + Mod(1, 65537)*I
gp > x^65537
%2 = Mod(2, 65537) + Mod(1, 65537)*I
gp > x^65535
%3 = Mod(52430, 65537) + Mod(39322, 65537)*I

(7)t_PADIC: Qp   
若p为素数,p_adic数定义为
\sum_{i=k}^{oo} {a_i * p^i}
gp > t= 2*7^(-1) + 3 + 4*7 + 2*7^2 + O(7^3)
%1 = 2*7^-1 + 3 + 4*7 + 2*7^2 + O(7^3)
gp > type(t)
%2 = "t_PADIC"
注意,系统并不检查p是不是素数,所以可以
gp > t = 2 * (1/10 + O(10^5));
gp > type(t)
%4 = "t_PADIC"
但是参与运算时候可能发生错误

(8)t_QUAD   : Q
非平方非零整数 d, d%4=0 或者d%4=1
w = quadgen(d, 'w)
d%4=0时候,w=sqrt(d)/2
d%4=1时候,w=(1+sqrt(d))/2
二次域数表达为 x + y * w,其中x, y是t_INT, t_INTMOD,t_FRAC, t_POL,t_POLMOD类型
令z = x + y*w

z.disc返回d
z.mod返回z的最小多项式T
z.pol返回t_POL类型的x+y*w
real(z)返回x, imag(z)返回y
gp > w=quadgen(-3,'w)
%1 = w
gp > z=1+2*w
%2 = 1 + 2*w
gp > z.disc
%3 = -3
gp > z.mod
%4 = w^2 - w + 1
gp > z.pol
%5 = 2*w + 1
gp > real(z)
%6 = 1
gp > imag(z)
%7 = 2

(9)t_POL    : T      
多项式,系数可以是t_INT, t_INTMOD, t_FAC, t_COMPLEX, t_QUAD

gp > w=quadgen(-3,'w)
%1 = w
gp > z=1+2*w
%2 = 1 + 2*w
gp > p=1+x+1/2*x^2+12*x^3+Mod(1,7)*x^4+z*x^5
%3 = (1 + 2*w)*x^5 + Mod(1, 7)*x^4 + Mod(5, 7)*x^3 + Mod(4, 7)*x^2 + Mod(1, 7)*x + Mod(1, 7)
gp > p1=1.2+x^2
%4 = x^2 + 1.20.....

多项式可以用常规方法输入
p = x^4 + 2*x^3 + 3*x^2 + 4*x + 5
也可以用向量形式,默认变量是x,从最高次开始
gp > P = Pol()
%1 = x^4 + 2*x^3 + 3*x^2 + 4*x + 5
也可以指定变量
gp > P = Pol(, 'y)
%1 = y^4 + 2*y^3 + 3*y^2 + 4*y + 5
也可以从最低次开始
gp > P = Polrev(, 'y)
%1 = 5*y^4 + 4*y^3 + 3*y^2 + 2*y + 1

(10)t_POLMOD : T/(P)         
MOD(x, y)中x, y都是多项式
若z=Mod(x, y),
z.pol返回化简后的x
z.mod返回y
gp > zz=Mod(x^7+1, x^2+x+1)
%1 = Mod(x + 1, x^2 + x + 1)
gp > zz.pol
%2 = x + 1
gp > zz.mod
%3 = x^2 + x + 1

(11)t_SER    : T((X))
类似于p-adic的形式的多项式
gp > f = x^(-3)+2*x^(-2)+1+x^3+O(x^7)
%1 = x^-3 + 2*x^-2 + 1 + x^3 + O(x^7)
gp > type(f)
%2 = "t_SER"
非有理函数会被转化成这种形式
gp > sin(x)
%3 = x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + 1/362880*x^9 - 1/39916800*x^11 + 1/6227020800*x^13 - 1/1307674368000*x^15 + O(x^17)
gp > type(sin(x))
%4 = "t_SER"

(12)t_RFRAC: T(X)
有理函数中分母不是1的形式
gp > type((1+x^2)/x)
%1 = "t_RFRAC"
gp > type((1+x^2))
%2 = "t_POL"

(13)t_QFB    : 二元二次型
通过 Qfb(a, b, c) 定义的 p = a*x^2 + b*x*y + c*y^2 二次型

(14)t_VEC    : T^n(行向量)      
(15)t_COL    : T^n(列向量)
t_VEC和t_COL都是向量,区别是VEC是横向的,COL是纵向的
输出时候,COL有~后缀
gp > v1=Vec()
%1 =                                                             ^-
gp > v2=Vec()
%2 =
gp > c1=Col()
%3 = ~
gp > c2=Col()
%4 = ~
gp > m=matrix(3,3)
%5 =






gp > m=v1
%6 =
gp > m=v2
%7 =
gp > m[,1]=c1
%8 = ~
gp > m[,2]=c2
%9 = ~
gp > m
%10 =





两者都可以通过1开始的下标访问和赋值
v1, v2, c1, c2等
也可以用
v这种范围语法
v[^2]表示第二项移除后的结果
-1表示最后一个元素
因此有
V[^-1]移除最后一个元素
V[-2..-1]表示最后两个元素
也可用vector函数
vector(n, x, seq)
vector(5, i, i+1)=
Col(vector(5, i, i^2)) = ~
或者用

例如
, isprime(x)] =
[ p | p <- primes(10), isprime(p+2) ] =
[ x+y | x <- ; y <- ] =
[ | x <- , isprime(x); y <- , y % 3 == 1 ] = [, ]

(16)t_MAT    : M_{m*n}(T)(矩阵)
矩阵赋值
matrix(m,{n = m},{X},{Y},{expr = 0})
gp > m = matrix(2,2)
%1 =



gp > matrix(5, 5, x, y, x+y)
%2 =









也可以
gp > m = Mat()
%1 =






m, m...表示矩阵的行向量
m[,1], m[,2]...表示矩阵的列向量
m表示前两行
m[^1,]表示移除第一行后其他行
m表示第1行,2-3列
m表示第1列,2-2行
m[^1, 2..3]表示移除第1行后剩余行的2-3列

(17)t_LIST   : 列表            
t_List类似t_VEC, t_COL,可以通过List()从他们转成t_LIST
List() = List([])
List(Vec()) = List()
List(Col()) = List()
List([ x+y | x <- ; y <- ]) = List()
其他例子不再赘述
t_LIST可以通过listput追加元素
gp > l = List()
%1 = List()
gp > listput(~l, 10)
%2 = 10
gp > listput(~l, 20)
%3 = 20
gp > l
%4 = List()
这里的~l表示对l的引用,不需要复制到函数

(18)t_STR    : 字符串
""括起来的字符串,"用\"表示
例如 "123"         

(19)t_VECSMALL: 小整数向量
以当前机器字长限制的小整数的向量,32位下绝对值小于w^31,64位下绝对值小于2^63
gp > x = Vec()
%1 =
gp > y = Vecsmall()
%2 = Vecsmall()
gp > type(x)
%3 = "t_VEC"
gp > type(y)
%4 = "t_VECSMALL"
gp > x=2^65
%5 = 36893488147419103232
gp > x
%6 =
gp > y=2^63
***   at top-level: y=2^63
***               ^---------
***   incorrect type in t_VECSMALL assignment (t_INT).

(20)t_CLOSURE: 函数
所有函数

(21)t_ERROR: 错误文本   
每当发生错误时都会创建这种类型的对象:它包含有关错误和错误上下文的一些信息。通常,会立即打印适当的错误,计算中止,GP进入“中断循环”
gp > 1/0; 1 + 1
***   at top-level: 1/0;1+1
***                  ^------
*** _/_: division by a noninvertible object
***   Break loop: type 'break' to go back to the GP prompt

可以用iferr捕获错误,并打印错误信息
gp > i = Mod(6,12); iferr(1/i, E, print(E)); 1 + 1
error("impossible inverse modulo: Mod(6, 12).")
%1 = 2
gp > i = Mod(6,12); iferr(1/i, E, print(component(E,1)));
Mod(6, 12)

(22)t_INFINITY: -oo 和 +oo
正负无穷      

无心人 发表于 2023-11-22 16:44:59

4、函数定义和常见函数
1)定义

gp > f(x)=x^2+1
%1 = (x)->x^2+1
gp > f(3)
%2 = 10

gp > f(x,y)=x+y
%3 = (x,y)->x+y
gp > f(i,5)
%4 = i + 5

或者以下方式也可以
f=x->x^2+1
f=(x)->x^2+1
f=(x,y)->x+y

也可以定义多行函数,每行用;结束即可,函数值是最后一行的结果

g(x) = y = x^2; z = y % 9;
g(7)
4

如果在某些需要seq参数的函数中,这些行需要用{}括起来


无心人 发表于 2023-11-22 17:05:22

2) if
if(a, seq1)
如果 a != 0 执行seq1
if(a, seq1, seq2)
如果 a != 0 执行seq1,否则执行seq2
if(a1, seq1,
   a2, seq2.
   .....
    defaultseq)
这个执行多个判定条件,类似连续的if判断
所有条件都不满足的话,执行最后的默认seq

无心人 发表于 2023-11-22 17:47:04

3) for
for(x=a, b, seq)
a, b是实数,从a 到 b 执行 seq, 步长是1
如果 a > b 则不执行
如果b 是 +∞(b = oo ),则循环是无限循环

GP > for(x=2.5, 3, print(x))
2.5

GP > for(x=2.5, 4, print(x))
2.5
3.5

可以用break跳出循环
GP> for(i=1,100,{print(i); if(i>=10,break())})
1
2
3
4
5
6
7
8
9
10
可以用next跳过当前循环进入下一个
gp > for(i=1,10,if(isprime(i), next, print(i)))
1
4
6
8
9
10

其他for函数:
forcomposite(n = a, b, seq)
对 区间的合数执行 seq
gp > forcomposite(a=2, 10,print(a))
4
6
8
9
10

forprime(n = a, b, seq)
对 区间的素数执行 seq
gp > forprime(n=2,10, print(n))
2
3
5
7

forprimestep(p = a, b, q, seq)
这里a, b都是unsigned long类型(32或者64位), p是素数, a <= p <= b
若q是t_INT,p满足p=a mod q
若q是t_INTMOD, q = Mod(c, N), 则Mod(p, N)=Mod(c, N)
gp > forprimestep(p = 4, 30, 5, print(p))
   19
   29
gp > forprimestep(p = 4, 30, Mod(1,5), print(p))
   11

fordiv(n, d, seq)
对n的因子d执行 seq
gp > fordiv(100, d, print(d))
1
2
4
5
10
20
25
50
100

fordivfactored(n, d, seq)
对n进行分解,然后对n的所有因子和因子的分解执行seq

gp > fordivfactored(24, d, print(d))

)]
)]
)]
]
)]
]
]
gp > fordivfactored(24, d, print(eulerphi(d)))
1
1
2
2
2
4
4
8

forfactored(n = a, b, seq)
从a到b,对n和n的分解执行seq
gp > forfactored(n=2,10, print(n))
)]
)]
)]
)]
]
)]
)]
)]
]


forsquarefree(N=a, b, seq)
无平方因子的N和N的分解
gp > s=0;forsquarefree(n=1,100, s+=n);s
%1 = 2967
gp > forsquarefree(n=1,10, print(n));

)]
)]
)]
]
)]
]

foreach(v, x, seq)
对v中每个元素x执行seq, v的类型是t_LIST, t_VEC,t_COL, t_MAT
如果类型是t_MAT,x是矩阵的行
gp > foreach(factor(24), d, print(d))
~
~


forpart(X = k,seq,{a = k1},{n = k2})
X是满足x_1+x_2...+x_n = k的vec,如果可选项a=k1存在,则x_i最大是k1,如果可选项n=k2存在,则最多k2项
gp > forpart(X=4, print(X))
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
gp > forpart(v=5,print(Vec(v)), 4, 3)
   
   
   
   
有2-4项,小于5,填充0
gp > forpart(v=5,print(Vec(v)),,)
   
   
   
   
   

forperm(a, p, seq)
若a是t_INT,对1..a的全排列执行seq
若a是t_VECS,对a的全排列执行seq
gp > forperm(3, p, print(p))
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
gp > forperm(, p, print(p))
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()

forstep(X = a,b,s,seq)
a<=X<=b
若s是整数,则X=a, a+s, a+2s,a+3s......
若s是Mod(c, N),则Mod(X, N) = Mod(c, N)
若s是vec, 比如s =
则X = a, a+1, a+1+2, a+1+2+1, a+1+2+1+2...
gp > forstep(x=5, 10, 2, print(x))
   5
   7
   9
gp > forstep(x=5, 10, Mod(1,3), print(x))
   7
   10
gp > forstep(x=5, 10, , print(x))
   5
   6
   8
   9
gp >forstep(x=5, 20, , print(x))
5
6
8
11
12
14
17
18
20

forsubgroup(H = G,{bound},seq)
群G的阿贝尔子群,执行seq

forsubset(nk, s, seq)
若nk是整数,的所有子集执行seq
若nk=, 的所有k个元素的子集执行seq
gp > forsubset(, s, print(s))
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
gp > forsubset(3, s, print(s))
   Vecsmall([])
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()
   Vecsmall()

forvec(X = v,seq,{flag = 0})
v是n个组成的,
函数生成n个元素的vec,每个元素v_i,都满足a_i<=v_i<=b_i
如果flag=1,vec中元素满足v_i<=v_j, i<j
如果flag=2,vec中元素满足v_i<v_j, i<j

gp > forvec (X=[,[-1,1]], print(X));
   
   
   
   
   
   
gp > forvec (X=[,[-1,1]], print(X), 1);
   
   
   
   ? forvec (X=[,[-1,1]], print(X), 2)
   

gp >forvec (X=[,[-1,1], ], print(X));



























gp >forvec (X=[,[-1,1], ], print(X), 1);









gp >forvec (X=[,[-1,1], ], print(X), 2);




无心人 发表于 2023-11-23 08:39:52

4)while
while(a, seq)
当a不为0时候,执行seq,直到a是0

gp > x=1;while(x<=10, {print(x);x=x+1});
1
2
3
4
5
6
7
8
9
10

注意,
如果a永远为1,将无限循环seq
如果a一开始为0,则不执行seq

无心人 发表于 2023-11-23 08:47:02

5) until
until(a, seq)
执行seq直到a为1
需要注意的是,seq会至少执行一次

gp > x=1;until(x>0, {print(x);x=x+1});
1
gp > x=1;until(x>10, {print(x);x=x+1});
1
2
3
4
5
6
7
8
9
10

until 和 while 是相反的条件

无心人 发表于 2023-11-27 14:35:37

6)输入输出

(1)input(): 获取命令行输入
gp > x = input()
1
%1 = 1

(2)print({str}*),print1({str}*), printp({str}*): 输出
gp > x = 1
%1 = 1
gp > print(x)
1
gp > m =
%83 =




gp > print(m) //以raw方式显示,等价于print1

gp > printp(m) //prettymatrix方式显示





(3) printsep(sep,{str}*)
用sep字符间隔的输出
printsep1(sep,{str}*)
和printsep类似,但是输出成raw形式

gp > printsep(":", 1,2,3,4)
1:2:3:4
gp > printsep(".", 1,2,3,4)
1.2.3.4
gp > printsep(",", 1,2,3,4)
1,2,3,4
gp > printsep(" ", 1,2,3,4)
1 2 3 4
gp > printsep1(":", , )
:

(4) printf(fmt, {x}*)
根据fmt字符串格式化输出,类似C的同名函数
不带%的原样输出
%的根据转换规则输出
即%开始,
    后跟一个或多个可选标注(#, 0, -, +, " "),
    一个可选的输出最小输出字符串长度,一个可选的小数精度,这两个用 . 分割,
    可选的转换类型字符 (l, d, i, o, u, x, X, p, e, E, f, F,g, G, s)
#表示16进制输出x格式加前缀"0x", X格式加前缀"0X",8进制o格式加前缀"o"
0表示数字类型的,规定最小输出字符串长度后,不足部分补0,
" "表示,规定最小输出字符串长度后,不足部分补空格
-表示输出会对齐在左边,默认对齐右边,- 和 0同时出现,忽略0
      即当输出长度大于实际数字宽度时候,- 会导致在输出位置先输出数字,后出现补齐的空格
+表示 + 和 - 符号总是会输出,+ 和" "同时出现,会忽略" "
m.n 宽度限定,表示最少输出m个字符,小数固定n位,m可省略
类型部分
l: 有符号机器字十进制整数
d, i:有符号十进制整数
u:无符号十进制整数
(当用u输出大的负整数时候,显示的是补码表示)
o: 无符号八进制整数
x, X:无符号十六进制表示,其中A-F字母,x格式用小写a-f表示,X格式用大写A-F表示
p: GEN类型
e, E:数字表达成小数 d * 10^n 形式,其中1<= d < 10, n是有符号整数, e格式输出成d en, E格式输出成d En
f, F:当数字绝对值小于1时候,总是输出成小数 d.dddd 的形式,
      当数字绝对值大于等于1时候,如果小数形式精度可以输出到小数点,用小数形式,
      否则输出和e, E一致(这里要依赖于realprecision的设定)
                例如realprecision = 19, 输出 123456789012345678.9可以用小数形式,但1234567890123456789就得用指数形式
g, G: 当输出数字绝对值小于1,如果绝对值小于1E-4,用小数形式输出,否则等同于e, E
      当输出数字绝对值大于等于1,按照 f, F形式处理
s: 字符串

gp > printf("floor: %d, field width 3: %3d, with sign: %+3d\n", Pi, 1, 2);
   floor: 3, field width 3:   1, with sign:+2

gp > printf("%.5g %.5g %.5g\n",123,123/456,123456789);
   123.00 0.26974 1.2346 e8

gp >x = 23; y=-1/x; printf("x=%+06.2f y=%+0*.*f\n", x, 6, 2, y);
   x=+23.00 y=-00.04

gp > for (i = 2, 5, printf("%05d\n", 10^i))
   00100
   01000
   10000
   100000\\ don't truncate fields whose length is larger than the minimum width
gp > printf("%.2f|%06.2f|", Pi,Pi)
   3.14|3.14|

gp > printf("%4d", );
   [   1,   2,   3]
gp >printf("%5.2f", mathilbert(3));
   [ 1.000.500.33]

   [ 0.500.330.25]

   [ 0.330.250.20]

(5) strprintf(fmt, {x}*)
和printf类似,但是会输出到一个字符串里

(6) printtex({x}*), strtex({x}*): 输出成tex形式
printtex在标准输出设备上输出,strtex输出到字符串
gp >v =
   %1
   ? strtex(v)
   %2 = "\\pmatrix{ 1&2&3\\cr}\n"

无心人 发表于 2023-11-27 16:37:55

7) strjoin 和 strsplit
strjoin(v,{p = ""}):
把v中元素用p连接起成字符串,如果p未指定,则v直接连接成字符串
gp > strjoin(, "-")
%1 = "1-2-3-4-5"
gp > strjoin(, "*")
%2 = "1*2*3*4*5"
gp > strjoin(, ".")
%3 = "1.2.3.4.5"
gp > strjoin()
%4 = "12345"
gp > strjoin(["hello", "world"], ",")
%5 = "hello,world"


strsplit(s,{p = ""}):
和strjoin相反的操作,把字符串 s 通过 p,分割成 vec, 如果 p 未指定,则分割成单个字符序列
gp > strsplit("1-2-3-4-5", "-")
%1 = ["1", "2", "3", "4", "5"]
gp > strsplit("1-2-3-4-5")
%1 = ["1", "-", "2", "-", "3", "-", "4", "-", "5"]
如果开始,结束位置有p字符,将分割出0长度字符""
gp > strsplit("-1-2-3-4-5-", "-")
%1 = ["", "1", "2", "3", "4", "5", ""]

无心人 发表于 2023-11-27 16:40:00

8) apply, select

apply(f, A)
如果A是数字或者字符串,返回f(A)
如果A是t_POL或者t_SER,返回对应类型,但是对应系数a_i系数变成f(a_i)
如果A是t_VEC, t_VECSMALL, t_LIST, t_MAT,对应位置的值x变成f(x)
gp > f=x->x^2+1
%1 = (x)->x^2+1
gp > apply(f, 2)
%2 = 5
gp > apply(f, )
%3 =
gp > apply(f, y^2+2*y+1)
%4 = 2*y^2 + 5*y + 2
gp > apply(f, )
%5 =
[ 25]



select(f, A, {flag = 0})、
A是t_VEC 或者 t_LIST
如果flag=0,返回 A 中满足 f(x) 不是0 的元素组成的新 t_VEC 或者 t_LIST
如果flag=1,返回 A 中满足 f(x) 不是0 的元素在原变量里的位置索引组成的 t_VECSMALL
gp > select(isprime, vector(100, i, i))
%1 =
gp > select(isprime, vector(100, i, i^2+1))
%2 =
gp > select(isprime, vector(100, i, i^2+1), 1)
%3 = Vecsmall()
gp > select(isprime, List())
%4 = List()

无心人 发表于 2023-11-27 16:56:22

9)fold, prod, sum

fold(f, A)
f是二元函数f(x, y),A是至少有两个元素的t_VEC
返回f(...f(f(A,A),A)...,A[#A])
gp > f=(x,y)->x+y
%1 = (x,y)->x+y
gp > fold(f, )
%2 = 15

prod(X=a, b, expr, {x=1})
遍历a, b然后计算expr, 再把expr结果乘起来
如果最后的参数x默认不是1,再乘以x
gp > prod(i=1, 10, i)
%1 = 3628800
gp > prod(i=1, 10, i, 10)
%2 = 36288000

sum(X = a,b,expr,{x = 0}):
遍历a, b然后计算expr, 再把expr结果加起来
如果最后的参数x默认不是0,再加x
gp > sum(i=1, 10, i)
%1 = 55
gp > sum(i=1, 10, i, 10)
%2 = 65

页: [1] 2
查看完整版本: PARI/GP的一些笔记