PARI/GP编程示例
目录
- 推荐一些语法手册
default(realprecision, 110)控制实数计算的有效位数default(parisize, "32M");增加PARI/GP申请的栈内存大小- 顺序结构、选择结构、循环结构
- for-loop中设置步长
vector或者apply实现map操作- 格式化打印
- switch-case
- for-each one list
- 应该没有布尔类型
- 字符串不能直接取某位置的字符，需要
Vecsmall转成数组，然后才能取某位置的元素 - subString和isSubstring
- 哈希表
- 因子分解的后处理
- 借助
list写法可以同时给多个变量进行赋值 - 字符串数组拼接成字符串
- 获取某个数域的某个property
- 反转向量后再切片
- 遇到factorial()函数返回结果不是整数怎么办？
bitand以及hammingweight
推荐一些语法手册
https://math.mit.edu/~brubaker/PARI/PARIrefcard.pdf
https://math.mit.edu/~brubaker/PARI/PARIusers.pdf
default(realprecision, 110)控制实数计算的有效位数
\\ Set the precision to, say, 110 digits
default(realprecision, 110);
\\ Function to display Apéry's constant
my_show_apery_constant(expr, caption) = printf("%s %.101g\n\n", caption, expr);
\\ Apéry's constant via PARI/GP's zeta function
my_show_apery_constant(zeta(3), "Apéry's constant via PARI/GP's Zeta:\n");
\\ Apéry's constant via reciprocal cubes
my_show_apery_constant(sum(k = 1, 1000, 1/k^3), "Apéry's constant via reciprocal cubes:\n");
\\ Apéry's constant via Markov's summation
my_show_apery_constant((5/2) * sum(k = 1, 158, (-1)^(k - 1) * (k!)^2 / ((2*k)! * k^3)), "Apéry's constant via Markov's summation:\n");
\\ Apéry's constant via Wedeniwski's summation
my_show_apery_constant(1/24 * sum(k = 0, 19, (-1)^k * ((2*k + 1)!)^3 * ((2*k)!)^3 * (k!)^3 * (126392*k^5 + 412708*k^4 + 531578*k^3 + 336367*k^2 + 104000*k + 12463) / (((3*k + 2)!) * ((4*k + 3)!)^3)), "Apéry's constant via Wedeniwski's summation:\n");
default(parisize, "32M");增加PARI/GP申请的栈内存大小
\\ Increase the stack size if necessary
default(parisize, "32M"); \\ Increase to 32MB, adjust if necessary
Riordan(N) = {
my(a = vector(N));
a[1] = 1; a[2] = 0; a[3] = 1;
for (n = 3, N-1,
a[n+1] = ((n - 1) * (2 * a[n] + 3 * a[n-1]) \ (n + 1)); \\ Integer division
);
return(a);
}
rios = Riordan(10000);
\\ Now print the first 32 elements in the desired format
for (i = 1, 32,{
print1(rios[i]," ")
});
print("")
\\ Print the number of digits for the 1000th and 10000th Riordan numbers
print("The 1,000th Riordan has ", #digits(rios[1000]), " digits.");
print("The 10,000th Riordan has ", #digits(rios[10000]), " digits.");
KlamerRado(N) = {
my(kr = vector(100 * N), ret = [], idx = 1);
kr[1] = 1;
while (idx <= #kr / 3,
if (kr[idx],
if (2 * idx + 1 <= #kr, kr[2 * idx + 1] = 1);
if (3 * idx + 1 <= #kr, kr[3 * idx + 1] = 1);
);
idx++;
);
for (n = 1, #kr,
if (kr[n], ret = concat(ret, n));
);
ret
}
default(parisize, "1024M");
{
kr1m = KlamerRado(1000000);
print("First 100 Klarner-Rado numbers:");
for (i = 1, 100,
print1(kr1m[i], " ");
);
print();
print("The 1,000th Klarner-Rado number is ", kr1m[1000]);
print("The 10,000th Klarner-Rado number is ", kr1m[10000]);
print("The 100,000th Klarner-Rado number is ", kr1m[100000]);
print("The 1000,000th Klarner-Rado number is ", kr1m[1000000]);
}
顺序结构、选择结构、循环结构
isDuffinian(n) = (!isprime(n)) && (gcd(n, sigma(n)) == 1);
testDuffinians()=
{
print("First 50 Duffinian numbers:");
count = 0; n = 2;
while(count < 50,
if (isDuffinian(n),
print1(n, " ");
count++;
);
n++;
);
print("\n\nFifteen Duffinian triplets:");
count = 0; n = 2;
while (count < 15,
if (isDuffinian(n) && isDuffinian(n + 1) && isDuffinian(n + 2),
print(n, " ", n + 1, " ", n + 2);
count++;
);
n++;
);
}
testDuffinians();
/* Define the Cullen and Woodall number functions */
cullen(n) = n * 2^n + 1;
woodall(n) = n * 2^n - 1;
{
/* Generate the first 20 Cullen and Woodall numbers */
print(vector(20, n, cullen(n)));
print(vector(20, n, woodall(n)));
/* Find the first 5 Cullen prime numbers */
cps = [];
for(i = 1, +oo,
if(isprime(cullen(i)),
cps = concat(cps, i);
if(#cps >= 5, break);
);
);
print(cps);
/* Find the first 12 Woodall prime numbers */
wps = [];
for(i = 1, +oo,
if(isprime(woodall(i)),
wps = concat(wps, i);
if(#wps >= 12, break);
);
);
print(wps);
}
for-loop中设置步长
这里,用来分隔函数里的各个入参。
;用来"显式"地表示;前面是一个语句。
\\ 会打印出1,2,3,4,5
for (n = 1, 5,
;
print(n);
)
\\ 左闭右闭
left=10;
right=1;
step=-3;
forstep(x=left,right,step,
print(x)
)
vector或者apply实现map操作
apply(x->x^2, [1,2,3,4])
\\ Define the Jacobsthal function
Jacobsthal(n) = (2^n - (-1)^n) / 3;
\\ Define the JacobsthalLucas function
JacobsthalLucas(n) = 2^n + (-1)^n;
\\ Define the JacobsthalOblong function
JacobsthalOblong(n) = Jacobsthal(n) * Jacobsthal(n + 1);
{
\\ Generate and print Jacobsthal numbers for 0 through 29
print(vector(30, n, Jacobsthal(n-1)));
\\ Generate and print JacobsthalLucas numbers for 0 through 29
print(vector(30, n, JacobsthalLucas(n-1)));
\\ Generate and print JacobsthalOblong numbers for 0 through 19
print(vector(20, n, JacobsthalOblong(n-1)));
\\ Find the first 20 prime numbers in the Jacobsthal sequence
myprimes = [];
i = 0;
while(#myprimes < 40,
if(isprime(Jacobsthal(i)), myprimes = concat(myprimes, [i, Jacobsthal(i)]));
i++;
);
for (i = 1, #myprimes\2, print(myprimes[2*i-1] " " myprimes[2*i]); );
}
格式化打印
printf("floor: %d, field width 3: %3d, with sign: %+3d\n", Pi, 1, 2);
printf("%.5g %.5g %.5g\n",123,123/456,123456789);
printf("%-2.5s:%2.5s:%2.5s\n", "P", "PARI", "PARIGP");
x = 23; y=-1/x; printf("x=%+06.2f y=%+0*.*f\n", x, 6, 2, y);
for (i = 2, 5, printf("%05d\n", 10^i))
for (i = 2, 5, printf("%05d\n", 10^i))
printf("%4d", [1,2,3]);
printf("%5.2f", mathilbert(3));
print("fsasfafaT " 2314 " afdgj") \\ 直接拼接
switch-case
这里,和;的位置要熟背
\\ 最外层的花括号不能省略
{
casename = "place_holder";
r = 2;
if (r == 2,
casename = "triangular";
,
r == 3,
casename = "tetrahedral";
,
r == 4,
casename = "pentatopic";
,
r == 12,
casename = "12-simplex";
);
print(casename);
}
for-each one list
V = [[1,2,3,4], [5,6,7,8], [9,10,11,12]];
foreach(V, x, my([a,b,c,d] = x); print([a,b,c,d, a+b^2+c^3+d^4]))
我语法没学好，我还是习惯下面这么写：
/* Polytopic number generation function */
polytopic(r, range) = {
vector(#range, i, binomial(range[i] + r - 1, r))
}
/* Triangular root function */
triangularRoot(x) = {
(sqrt(8*x + 1) - 1)/2
}
/* Tetrahedral root function */
tetrahedralRoot(x) = {
N = (3*x + sqrt(9*x^2 - 1/27))^(1/3) + (3*x - sqrt(9*x^2 - 1/27))^(1/3) - 1;
precision(N, 18)
}
/* Pentatopic root function */
pentatopicRoot(x) = {
(sqrt(5 + 4*sqrt(24*x + 1)) - 3)/2
}
{
/* Displaying polytopic numbers */
r_sel=[2,3,4,12];
for(i = 1, #r_sel,
r=r_sel[i];
casename="place_holder";
if(r == 2, casename="triangular"; ,
r == 3, casename="tetrahedral"; ,
r == 4, casename="pentatopic"; ,
r == 12, casename="12-simplex";
);
printf("\nFirst 30 %s numbers:\n %s\n" , Str(casename) , Str(polytopic(r, [0..29])) )
);
/* Displaying roots of specific numbers */
nums = [7140, 21408696, 26728085384, 14545501785001];
for(i = 1, #nums,
n = nums[i];
printf("\nRoots of %s:\n", Str(n));
printf(" triangular-root: %s\n", Str(triangularRoot(n)));
printf(" tetrahedral-root: %s\n", Str(tetrahedralRoot(n)));
printf(" pentatopic-root: %s\n", Str(pentatopicRoot(n)))
);
}
应该没有布尔类型
官方文档里说
Any nonzero value is interpreted as true and any zero as false (this includes empty vectors or matrices).
我这里直接用0,1了。
chowla(n) = {
if (n < 1, error("Chowla function argument must be positive"));
if (n < 4, return(0));
my(divs = divisors(n));
sum(i=1, #divs, divs[i]) - n - 1;
}
\\ Function to count Chowla numbers
countchowlas(n, asperfect = 1, verbose = 1) = {
my(count = 0, chow, i);
for (i = 2, n,
chow = chowla(i);
if ( (asperfect && (chow == i - 1)) || ((!asperfect) && (chow == 0)),
count++;
if (verbose, print("The number " i " is " if (asperfect, "perfect.", "prime.")));
);
);
count;
}
\\ Main execution block
{
print("The first 37 chowla numbers are:");
for (i = 1, 37, printf("Chowla(%s) is %s\n", Str(i), Str(chowla(i)) ) );
m=100;
while(m<=10000000, print("The count of the primes up to " m " is " countchowlas(m, 0, 0)); m=m*10);
print("The count of perfect numbers up to 35,000,000 is " countchowlas(35000000, 1, 1));
}
字符串不能直接取某位置的字符，需要Vecsmall转成数组，然后才能取某位置的元素
对应的逆操作是str = Str(Vec(guessvector));
colors = ["grey", "yellow", "green"];
wordle(answer, guess) =
{
n = #guess;
if (#answer != n, error("The words must be of the same length."));
answervector = Vecsmall(answer);
guessvector = Vecsmall(guess); \\ Convert guess to a vector of ASCII values
result = vector(n, i, 0);
for (i = 1, n,
if (guessvector[i] == answervector[i],
answervector[i] = 0;
result[i] = 2;
);
);
for (i = 1, n,
c = guessvector[i];
for (j = 1, n,
if (answervector[j] == c,
answervector[j] = 0;
result[i] = 1;
break;
);
);
);
result;
}
{
testpairs = [["ALLOW", "LOLLY"], ["BULLY", "LOLLY"], ["ROBIN", "ALERT"], ["ROBIN", "SONIC"], ["ROBIN", "ROBIN"]];
for (i = 1, #testpairs,
pair = testpairs[i];
res = wordle(pair[1], pair[2]);
res2 = vector(#res, j, colors[res[j] + 1]);
print(pair[1] " v " pair[2] " => " res " => " res2);
)
}
subString和isSubstring
/* Returns a substring of str starting at s with length n */
ssubstr(str, s = 1, n = 0) = {
my(vt = Vecsmall(str), ve, vr, vtn = #str, n1);
if (vtn == 0, return(""));
if (s < 1 || s > vtn, return(str));
n1 = vtn - s + 1; if (n == 0, n = n1); if (n > n1, n = n1);
ve = vector(n, z, z - 1 + s); vr = vecextract(vt, ve); return(Strchr(vr));
}
/* Checks if subStr is a substring of mainStr */
isSubstring(mainStr, subStr) = {
mainLen = #Vecsmall(mainStr);
subLen = #Vecsmall(subStr);
for (startPos = 1, mainLen - subLen + 1,
if (ssubstr(mainStr, startPos, subLen) == subStr,
return(1); /* True: subStr found in mainStr */
)
);
return(0); /* False: subStr not found */
}
/* Determines if a number's factors, all > 9, are substrings of its decimal representation */
contains_its_prime_factors_all_over_9(n) = {
if (n < 10 || isprime(n), return(0)); /* Skip if n < 10 or n is prime */
strn = Str(n); /* Convert n to string */
pfacs = factor(n)[, 1]; /* Get unique prime factors of n */
for (i = 1, #pfacs,
if (pfacs[i] <= 9, return(0)); /* Skip factors ≤ 9 */
if (!isSubstring(strn, Str(pfacs[i])), return(0)); /* Check if factor is a substring */
);
return(1); /* All checks passed */
}
/* Main loop to find and print numbers meeting the criteria */
{
found = 0; /* Counter for numbers found */
for (n = 0, 30 * 10^6, /* Iterate from 0 to 30 million */
if (contains_its_prime_factors_all_over_9(n),
found += 1; /* Increment counter if n meets criteria */
print1(n, " "); /* Print n followed by a space */
if (found % 10 == 0, print("")); /* Newline every 10 numbers */
if (found == 20, break); /* Stop after finding 20 numbers */
);
);
}
哈希表
mymap = Map();
mapisdefined(mymap,key)
mapget(mymap, key)
mapput(mymap, key, value)
就这些
\\ Initialize an associative array equivalent
pisanos = Map();
\\ Function to calculate the Pisano period for a given prime p
pisano(p) = {
local(lastn, n, i);
if (p < 2, return(1));
if (mapisdefined(pisanos, p),
return(mapget(pisanos, p));
);
lastn = 0; n = 1;
for (i = 1, p^2,
[lastn, n] = [n, Mod(lastn + n, p)];
if (lastn == 0 && n == 1,
mapput(pisanos, p, i);
return(i);
);
);
return(1);
}
\\ Function to calculate Pisano period for a prime p and a power k
pisanoprime(p, k) = {
my(i = pisano(p));
if (!isprime(p), error("p must be prime"));
p^(k-1) * i;
}
\\ Function to calculate Pisano period for a composite number n
pisanotask(n) = {
my(factors = factor(n));
\\print("factors=" factors);
\\apply(x -> print("x=" x), factors);
lcm( vector(#factors~, i, pisanoprime(factors[i, 1], factors[i, 2])) );
}
{
\\ Print Pisano periods for prime numbers up to 15 with k=2
for (i = 1, 15,
if (isprime(i),
print("pisanoPrime[" i ", 2] = " pisanoprime(i, 2))
);
);
\\ Print Pisano periods for prime numbers up to 180 with k=1
for (i = 1, 180,
if (isprime(i),
print("pisanoPrime[" i ", 1] = " pisanoprime(i, 1))
);
);
\\ Print Pisano periods for numbers 2 to 180
print("\nPisano[n] for n from 2 to 180:");
print(vector(179, i, pisano(i+1)));
\\ Print Pisano periods using pisanotask for numbers 2 to 180
print("\nPisano[n] using pisanoPrime for n from 2 to 180:");
print(vector(179, i, pisanotask(i+1)));
}
因子分解的后处理
factor得到的是一个Matrix，这里所谓的后处理就是在根据这个Matrix进行计算。
my(factors = factor(n));
\\print("factors=" factors);
\\apply(x -> print("x=" x), factors);
lcm( vector(#factors~, i, pisanoprime(factors[i, 1], factors[i, 2])) );
借助list写法可以同时给多个变量进行赋值
[lastn, n] = [n, Mod(lastn + n, p)]
字符串数组拼接成字符串
homeprimechain(n) = {
my(chain = [], concat_result, factors, factor_str);
while (!isprime(n),
chain = concat(chain, n); /* Append n to the chain */
factors = factor(n);
/* Correctly build the concatenated string of factors */
factor_str = Str(concat(Vec(vector(#factors~, i,
concat(Vec(concat(vector(factors[i, 2], j, Str(factors[i, 1])
)))
)))));
concat_result = factor_str;
\\print("concat_result=" concat_result);
n = eval(concat_result); /* Convert string back to number */
);
concat(chain, n); /* Append the final prime to the chain */
}
printHPiter(n, numperline = 4) = {
my(chain = homeprimechain(n), len = #chain, i);
for (i = 1, len,
if (i < len,
print1("HP" , chain[i] , " (" , len - i , ") = " , if(i % numperline == 0, "\n", "")),
print(chain[i], " is prime.\n\n");
);
);
}
{
/* Iterate over a set of numbers */
forstep(i = 2, 20, 1, print("Home Prime chain for ", i, ": "); printHPiter(i););
printHPiter(65);
}
获取某个数域的某个property
\\ https://www.lmfdb.org/NumberField/3.1.45387.1
K = bnfinit(y^3 - 41, 1)
print(K.pol)
print("Degree: " poldegree(K.pol))
print("Signature: " K.sign)
print("Discriminant: " K.disc)
print("Root discriminant: " abs(K.disc)^(1/poldegree(K.pol)))
print("Ramified primes: " factor(abs(K.disc))[,1]~)
print("Integral basis: " K.zk)
print("Class group and class number: " K.clgp)
print("Fundamental unit: " K.fu)
print("Torsion generator: " K.tu[2])
print("Regulator: " K.reg)
# self-contained Pari/GP code snippet to compute the analytic class number formula
K = bnfinit(x^3 - 41, 1);
print([polcoeff (lfunrootres (lfuncreate (K))[1][1][2], -1), 2^K.r1 * (2*Pi)^K.r2 * K.reg * K.no / (K.tu[1] * sqrt (abs (K.disc)))])
print("Galois group: " polgalois(K.pol))
L = nfsubfields(K); L[2..length(b)]
print("Intermediate fields" L)
\\ to obtain a list of $[e_i,f_i]$ for the factorization of the ideal $p\mathcal{O}_K$ for $p=7$ in Pari:
p = 7; pfac = idealprimedec(K, p);
print("Frobenius cycle types: " vector(length(pfac), j, [pfac[j][3], pfac[j][4]]));
反转向量后再切片
Vecrev(b)[j..j+i-1]
遇到factorial()函数返回结果不是整数怎么办？
虽然不知道为什么，但是打表先存一下好像可以解决该问题。
default("parisizemax", "1024M");
\\ Define the function wilsonprimes with a default limit of 11000
wilsonprimes(limit) = {
\\ Set the default limit if not specified
my(limit = if(limit, limit, 11000));
\\ Precompute factorial values up to the limit to save time
my(facts = vector(limit, i, i!));
\\ Sign variable for adjustment in the formula
my(sgn = 1);
print(" n: Wilson primes\n--------------------");
\\ Loop over the specified range (1 to 11 in the original code)
for(n = 1, 11,
print1(Str(" ", n, ": "));
sgn = -sgn; \\ Toggle the sign
\\ Loop over all primes up to the limit
forprime(p = 2, limit,
\\ Check the Wilson prime condition modified for PARI/GP
index=1;
if(n<2,index=1,index=n-1);
if(p > n && Mod(facts[index] * facts[p - n] - sgn, p^2) == 0,
print1(Str(p, " "));
)
);
print1("\n");
);
}
\\ Execute the function with the default limit
wilsonprimes();
bitand以及hammingweight
f(n)=(-1)^hammingweight(bitand(n,n*2));
/* Loop from 0 to 20 and print the result of applying f */
{
for(n = 0, 20, print(n " , " f(n)))
}
