PARI/GP编程示例

标签：constant casename GP 示例 PARI vector numbers printf print

素材主要来自我在RosettaCode上贡献的PARI/GP代码

default(realprecision, 110)控制实数计算的有效位数
default(parisize, "32M");增加PARI/GP申请的栈内存大小
顺序结构、选择结构、循环结构
for-loop中设置步长
vector实现map操作
格式化打印
switch-case
for-each one list
应该没有布尔类型

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]);
}

顺序结构、选择结构、循环结构

/* 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中设置步长

这里,和;的位置要熟背

for (n = 1, 11, n += 2; print(n))      \\ 会打印出3,6,9,12

vector实现map操作

\\ 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)))
);
}

应该没有布尔类型

直接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));
}

标签：constant,casename,GP,示例,PARI,vector,numbers,printf,print
From： https://blog.51cto.com/u_15247503/9451558