User talk:BenBE: Difference between revisions

→‎Language file: example script
m (→‎Language file: rm spaces)
(→‎Language file: example script)
Line 19:
 
I haven't yet figured out how to handle the escapes but will look at the PHP file and see if that helps.
 
I don't see how to upload a file to SourceForge without opening a new issue, so for now I'll leave it here.
 
<div style="height:30ex;overflow:scroll"><lang php><?php
Line 205 ⟶ 207:
 
[[User:CRGreathouse|CRGreathouse]] 01:32, 8 July 2011 (UTC)
 
Oh, and you had asked for sample source code:
 
<div style="height:30ex;overflow:scroll"><lang parigp>fibmod(n:int, m:int)={
my(f,t);
if (m < 6,
if (m < 0, m = -m);
if (m == 0, error("division by 0"));
if (m == 1, return(0));
if (m == 2, return(n%3>0));
if (m == 3, return(fibonacci(n%8)%3));
if (m == 4, return(fibonacci(n%6)%4));
if (m == 5, return(fibonacci(n%20)%5));
);
if (n < 1000,
return(fibonacci(n)%m);
);
f = factor(m);
t = Mod(fibonacci(n%Pisano(f[1,1], f[1,2])),f[1,1]^f[1,2]);
for(i=2,#f[,1],
t = chinese(t,Mod(fibonacci(n%Pisano(f[i,1], f[i,2])),f[i,1]^f[i,2]));
);
lift(t)
};
addhelp(fibmod,"fibmod(n,m): Returns the nth Fibonacci number mod m. Same as finonacci(n)%m, but faster for large n.");
 
 
\\ Helper function for fibmod; returns a multiple of the period of Fibonacci numbers mod p^e.
\\ Assumes p is prime and e > 0.
Pisano(p:int,e:small)={
my(t);
if(p==5, return(4 * 5^e));
t = p%5;
if (t == 2 || t == 3,
t = 2 * (p + 1)
,
t = p - 1
);
t * p^(e-1)
};
 
rad(n:int)={
if(n < 0, n = -n);
if (issquarefree(n), return(n));
vecprod(factor(n)[,1])
};
addhelp(rad, "rad(n): Radical of n, the largest squarefree number dividing n. Sloane's A007947.");
 
 
isprimepower(n:int)={
ispower(n,,&n);
isprime(n)
};
addhelp(isprimepower, "isprimepower(n): Is n a prime power? Characteristic sequence of Sloane's A000961.");
 
 
isPowerful(n:int)={
if (bitand(n,3)==2,return(0));
if (n%3 == 0 && n%9 > 0, return(0));
if (n%5 == 0 && n%25 > 0, return(0));
vecmin(factor(n)[,2])>1
};
addhelp(isPowerful, "isPowerful(n): Is n powerful (min exponent 2)? Sloane's A112526; characteristic function of Sloane's A001694.");
 
 
prp(n:int,b:int=2)={
Mod(b,n)^(n-1)==1
};
addhelp(prp, "prp(n,b=2): Is n a b-probable prime?");
 
 
sprp(n:int,b:int=2)={
my(d = n, s = 0);
until(bitand(d,1), d >>= 1; s++);
d = Mod(b, n)^d;
if (d == 1, return(1));
for(i=1,s-1,
if (d == -1, return(1));
d = d^2;
);
d == -1
};
addhelp(sprp, "sprp(n,b=2): Is n a b-strong probable prime?");
 
 
sopfr(n:int)={
my(f=factor(n));
sum(i=1,#f[,1],f[i,1]*f[i,2])
};
addhelp(sopfr, "sopfr(n): Sum of prime factors of n (with multiplicity). Sloane's A001414.");
 
 
primorial(n)={
my(t1=1,t2=1,t3=1,t4=1);
forprime(p=2,n>>3,t1=t1*p);
forprime(p=(n>>3)+1,n>>2,t2=t2*p);
t1=t1*t2;t2=1;
forprime(p=(n>>2)+1,(3*n)>>3,t2=t2*p);
forprime(p=((3*n)>>3)+1,n>>1,t3=t3*p);
t1=t1*(t2*t3);t2=1;t3=1;
forprime(p=(n>>1)+1,(5*n)>>3,t2=t2*p);
forprime(p=((5*n)>>3)+1,(3*n)>>2,t3=t3*p);
t2=t2*t3;t3=1;
forprime(p=((3*n)>>2)+1,(7*n)>>3,t3=t3*p);
forprime(p=((7*n)>>3)+1,n,t4=t4*p);
t1*(t2*(t3*t4))
};
addhelp(primorial, "Returns the product of each prime less than or equal to n. Sloane's A034386.");
 
 
gpf(n:int)={
my(f=factor(n)[,1]);
if (n <= 1,
if (n >= -1, return (1)); \\ Convention
n = -n
);
f[#f]
};
addhelp(gpf, "The greatest prime factor of a number. Sloane's A006530.");
 
 
lpf(n:int)={
if (n <= 1,
if (n >= -1, return (1)); \\ Convention
n = -n
);
forprime(p=2,2e4,
if(n%p==0,return(p));
);
factor(n)[1,1]
};
addhelp(lpf,"The least prime factor of a number. Sloane's A020639.");
 
isTriangular(n:int)={
issquare(n<<3+1)
};
addhelp(isTriangular, "isTriangular(n): Is n a triangular number? Sloane's A010054.");
 
 
isFibonacci(n:int)={
my(k=n^2);
k+=((k + 1) << 2);
issquare(k) || (n > 0 && issquare(k-8))
};
addhelp(isFibonacci, "isFibonacci(n): Is n a Fibonacci number? Sloane's A010056; characteristic function of Sloane's A000045.");
 
 
largestSquareFactor(n:int)={
my(f=factor(n));
prod(i=1,#f[,1],f[i,1]^(f[i,2]>>1))^2
};
addhelp(largestSquareFactor, "largestSquareFactor(n): Largest square dividing n. Sloane's A008833.");
 
 
hamming(n:int)={
my(v=binary(n));
sum(i=1,#v,v[i])
};
addhelp(hamming, "hamming(n): Hamming weight of n (considered as a binary number). Sloane's A000120.");
 
 
istwo(n:int)={
my(f);
if (n < 3,
return (n >= 0);
);
f = factor(oddres(n));
for (i=1,#f[,1],
if(bitand(f[i,2],1)==1 && bitand(f[i, 1], 3)==3, return(0))
);
1
};
addhelp(istwo, "Is the number a sum of two squares? Characteristic function of Sloane's A001481.");
 
 
ways2(n:int)={
\\sum(k=0,sqrtint(n>>1),issquare(n-k^2))
my(f=factor(oddres(n)),t=1);
for(i=1,#f[,1],
if (f[i,1]%4==3,
if (f[i,2]%2, return(0))
,
t *= f[i,2] + 1;
)
);
(t + 1) >> 1
};
addhelp(ways2, "Number of ways that n can be represented as a sum of two squares. Sloane's A000161.");
 
 
isthree(n:int)={
my(tmp=valuation(n, 2));
bitand(tmp, 1) || bitand(n >> tmp, 7) != 7
};
addhelp(isthree, "isthree(n): Is the number the sum of three squares? Sloane's A071374; characteristic function of Sloane's A000378.");
 
 
ways3(n:int)={
my(t);
sum(k=0,sqrtint(n\3),
t=n-k^2;
sum(m=k,sqrtint(t>>1),
issquare(t-m^2)
)
)
};
addhelp(ways3, "Number of ways that n can be represented as a sum of three squares. Sloane's A000164.");
 
 
msb(n:int)={
my(k:int=0);
n=floor(n);
while(n>15,
n>>=4;
k+=4
);
if(n>3,
n>>=2;
k+=2
);
min(n,2)<<k
};
addhelp(msb, "msb(n): Most significant bit of n: returns the greatest power of 2 <= the number. Sloane's A053644.");
 
 
Faulhaber(e:small,a='x)={
substpol(sum(i=0,e,binomial(e+1,i)*bernfrac(i)*'x^(e+1-i))/(e+1) + 'x^e, 'x, a)
};
addhelp(Faulhaber, "Faulhaber(e,{a='x}): Returns the polynomial for the sum 1^e + 2^e + ... + x^e, evaluated at a.");</lang></div>