This site is supported by donations to The OEIS Foundation.
User:R. J. Mathar/Pari
From OeisWiki
/**
* \mainpage
* @since 2010-09-08
* @author Richard J. Mathar
*
* <a href="http://oeis.org/">OEIS</a> functions in the
* \cite{PARI} interpreter language PARI/gp.
*
*/
/** GENERIC AUXILIARY FUNCTIONS
*/
/** Write the vector into a b-file.
* The elements of the vector are appended to the file.
* The file is created of not existing.
* @param[in] v The vector of values
* @param[in] f The file name to be created.
* @param[in] offs The index of the first element of v.
* @since 2010-09-08
* @todo f should be replaced by the A-number and the %I and %H lines be printed.
* @todo Do not append to an existing file.
* @author R. J. Mathar
*/
vecToBfile(v,f,offs)={
for(i=1,length(v),
write(f,i+offs-1" "v[i]) ;
)
}
/** n-th term of a linear recurrences with constant coefficients.
* @param[in] vi The vector of the initial values.
* @param[in] vc The vector of the coefficients o fthe recurrence.
* @param[in] n The index of the value returned, 1-based.
* @return a(n) defined by a(i)=vi[i] for n<=length(vi), otherwise a(n)=sum_{c=1..leng(vc)} vc[c]*a(n-c)
* @since 2010-09-08
* @author R. J. Mathar
*/
linHRec(vi,vc,n)={
local(vshf,j,s,lenvi) ;
j=length(vi) ;
lenvi=j ;
if(n <= j,
return( vi[n] )
) ;
if( length(vc) > lenvi,
error("LinHRec reccurence ", length(vc), " deep but only ",j," initial elements") ;
) ;
vshf=vector(j,i,vi[i]) ;
while(j <n,
s=sum(i=1,length(vc),vc[i]*vshf[lenvi+1-i]) ;
/* move up one element */
for(i=1,lenvi-1,
vshf[i]=vshf[i+1]
);
vshf[lenvi] = s ;
j++;
) ;
vshf[lenvi] ;
}
/** The first n terms of a linear recurrences with constant coefficients.
* @param[in] vi The vector of the initial values.
* @param[in] vc The vector of the coefficients of the recurrence.
* @param[in] n The length of the vector returned.
* @return a[1..n] defined by a(i)=vi[i] for n<=length(vi), otherwise a(n)=sum_{c=1..leng(vc)} vc[c]*a(n-c)
* @since 2010-09-08
* @author R. J. Mathar
*/
linHRec_vec(vi,vc,n)={
local(vshf,j,s,lenvi) ;
j=length(vi) ;
lenvi=j ;
if( length(vc) > lenvi,
error("LinHRec reccurence ", length(vc), " deep but only ",j," initial elements") ;
) ;
vshf=vector(n) ;
for(i=1, min(lenvi,n),
vshf[i] = vi[i]
) ;
while(j <n,
s=sum(i=1,length(vc),vc[i]*vshf[j+1-i]) ;
vshf[j+1] = s ;
j++;
) ;
vshf ;
}
/** A list of distinct factorizations of n.
* @param n The argument >=1 .
* @param maxf The maximum factor of n in the factorizations.
* @return A list of the form [list1,list2,list3,..] where all numbers in listi are <= maxf,
* where the product of all numbers in listi is n, and where the numbers in each listi are non-increasing.
* Example: n=100 and maxf=10 returns [[[5, 5, 2, 2], [5, 5, 4], [10, 5, 2], [10, 10]].
* This is mainly an auxiliary function to A005179, and the multiplicative analog to partitions.
* @since 2008-05-26
* @author R. J. Mathar
*/
prodR(n, maxf)={
local(dfs, a=[], r, tmp ) ;
dfs=divisors(n) ;
for(i=2, length(dfs),
if( dfs[i]<=maxf,
if(dfs[i]==n,
a=concat(a, [[n]]),
r=prodR(n/dfs[i], min(dfs[i], maxf)) ;
for(j=1, length(r),
tmp=concat(dfs[i], r[j]) ;
a=concat(a, [tmp]) ;
) ;
) ;
) ;
) ;
return(a) ;
}
/** NUMBERED SEQUENCES IN A-NUMBER ORDER
*/
/** Tau(n), number of divisors of n.
* @param[in] n The argument the divisors of which are counted.
* @return tau(n).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000005(n)={
/* sigma(n,0) */
numdiv(n)
}
addhelp(A000005,"tau(n), number of divisors of n") ;
/** Number of partitions of n into distinct parts.
* @param[in] n The argument n>=0.
* @return q(n).
* @since 2000
* @author N. J. A. Sloane
*/
A000009(n)={
polcoeff( prod(k=1,n, 1+x^k, 1+x*O(x^n)), n)
}
addhelp(A000009,"q(n), number of partitions of n into distinct parts") ;
/** Vector with number of partitions of 0..n-1 into distinct parts.
* @param[in] n The lenght of the vector.
* @return q(0..n-1).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000009_vec(n)={
local( gf= prod(k=1,n, 1+x^k, 1+x*O(x^n)), a=vector(n)) ;
for(i=0,n-1,
a[i+1]=polcoeff(gf,i)
) ;
a
}
addhelp(A000009_vec,"q(0..n-1), vector of partitions of 0..n-1 into distinct parts") ;
/** Euler totient function phi(n).
* @param[in] n The upper limit of primes counted.
* @return phi(n), count of numbers <=n and relative prime to n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000010(n)={
eulerphi(n)
}
addhelp(A000010,"phi(n), Euler's totient function.") ;
/** All 1's sequence.
* @param[in] n The (dummy) argument >=0.
* @return 1
* @since 2010-09-08
* @author R. J. Mathar
*/
A000012(n)={
1
}
/** The natural numbes.
* @param[in] n The argument >=1.
* @return n
* @since 2010-09-08
* @author R. J. Mathar
*/
A000027(n)={
if( n<1,
error("A000027 argument ", n, " <=0.") ;
) ;
n
}
/** The n-th Lucas number.
* @param[in] n The index into the sequences
* @return L(n), which L(0)=2, L(1)=1 and L(n)=L(n-1)+L(n-2).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000032(n)={
/* recurrence a(0,1) = (2,1) and a(n)=a(n-1)+a(n-2), offset 0 */
linHRec([2,1],[1,1],n+1)
}
addhelp(A000032,"Lucas(n)") ;
/** The first n Lucas numbers, L(0..n-1).
* @param[in] n The length of the vecto returned
* @return L(0..n-1), where L(0)=2, L(1)=1 and L(n)=L(n-1)+L(n-2).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000032_vec(n)={
linHRec_vec([2,1],[1,1],n)
}
addhelp(A000032_vec,"Vector Lucas(0..n-1)") ;
/** Test the argument for being prime.
* @param[in] n The value to be tested.
* @return true if n is prime.
* @since 2010-09-08
* @author R. J. Mathar
*/
isA000040(n)={
isprime(n)
}
/** The prime numbers.
* @param[in] n The index into the primes
* @return prime(n)
* @since 2010-09-08
* @author R. J. Mathar
*/
A000040(n)={
prime(n)
}
addhelp(A000040,"prime(n)") ;
/** The smallest prime larger than n.
* @param[in] n The number smaller than the prime returned.
* @return The prime p > n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000040_next(n)={
nextprime(n+1)
}
addhelp(A000040_next,"Smallest prime larger than n.") ;
/** Vector of the first n prime numbers.
* @param[in] n The number of vector components.
* @return prime(1..n)
* @since 2010-09-08
* @author R. J. Mathar
*/
A000040_vec(n)={
vector(n,i,A000040(i)) ;
}
addhelp(A000040_vec,"The vector of the first n primes") ;
/** The number of unrestricted partitions of n.
* @param[in] n The index into the P(n).
* @return P(n). This is 0 if n=0 or 1, 2 if n=2.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000041(n)={
numbpart(n)
}
addhelp(A000041,"Number of partitions of n.") ;
/** Exponents of Mersenne primes.
* @param[in] n The index into the list.
* @return The prime exponent: the n-th element in [2,3,4,7,13,..]
* This is simply taken from a finite lookup table.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000043(n)={
local(a=
[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937,
21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917]
) ;
if( n>length(a),
error("A000043 argument ", n, " too large.") ,
if(n <1,
error("A000043 argument ", n, " <1e.") ,
a[n]
)
) ;
}
addhelp(A000043,"Prime exponent of n-th Mersenne prime.") ;
/** The n-th fibonacci number.
* @param[in] n The index into the F(n).
* @return F(n). This is 0 if n=0, 1 if n=1 or 2.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000045(n)={
fibonacci(n)
}
addhelp(A000045,"fibonacci(n)") ;
/** The first n Fibonacci numbers.
* @param[in] n The length of the vector returned.
* @return Fibonacci(0..n-1).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000045_vec(n)={
vector(n,i,fibonacci(i-1))
}
/** The n-th power of 2.
* @param[in] n The exponent. >=0.
* @return 2^n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000079(n)={
2^n
}
addhelp(A000079,"2^n") ;
/** Catalan(n).
* @param[in] n The index into the Catalan numbers, >=0.
* @return C(n)
* @since 2010-09-08
* @author R. J. Mathar
*/
A000108(n)={
binomial(2*n,n)/(n+1)
}
addhelp(A000108,"Catalan(n)") ;
/** The first n Catalan numbers, starting C(0).
* @param[in] n The length of the vector returned.
* @return C(0..n-1)
* @since 2010-09-08
* @author R. J. Mathar
*/
A000108_vec(n)={
vector(n,i, A000108(i-1))
}
addhelp(A000108_vec,"The vector Catalan(0..n-1)") ;
/** n!, the factorial of n.
* @param[in] n The argument.
* @return 1*2*3*..*n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000142(n)={
n!
}
addhelp(A000142,"Factorial(n)");
/** Vector of the first n Bell numbers.
* @param[in] n The length of the vector to be created.
* @return Bell(0..n-1).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000110_vec(n)={
local(a=vector(n),s) ;
a[1] = 1 ;
a[2] = 1 ;
a[3] = 2 ;
if( n <= 3,
return(vector(n,i,a[i])) ;
) ;
for(nloc=3,n-1,
s = sum(k=1,nloc, a[nloc-k+1]/(nloc-k)!/(k-1)! );
a[nloc+1]= (nloc-1)!*s
) ;
a ;
}
addhelp(A000110_vec,"Vector of Bell(0..n-1)");
/** Bell(n). The number of ways of placing n labeled balls into n non-labelled boxes.
* @param[in] n The argument >=0.
* @return Bell(n).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000110(n)={
local v;
if( n < 0,
error("A002110 argument ", n, " <0.") ,
v = A000110_vec(n+1) ;
v[n+1] ;
) ;
}
addhelp(A000110,"Bell(n)");
/** Hamming weight of n. Number of 1's in the binary expansion of n.
* @param[in] n The argument.
* @return The number of bits set in the binary expansion of n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000120(n)={
local(a=0) ;
while(n>0,
n>>=1 ;
a += bittest(n,0) ;
);
a ;
}
addhelp(A000120,"Binary weight of n. Number of bits set in n.");
/** Pell(n).
* @param[in] n The index into the sequence.
* @return The element a(n) in a(0)=0, a(1)=1, a(n)=2*a(n-1)+a(n-2).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000129(n)={
linHRec([0,1],[2,1],n+1)
}
addhelp(A000129,"Pell(n).") ;
/** The vector of the first n Pell numbers.
* @param[in] n The length of the vector to be returned.
* @return The Pell(0..n-1), starting 0,1,2,5,12,29.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000129_vec(n)={
linHRec_vec([0,1],[2,1],n)
}
addhelp(A000129_vec,"The first n Pell numbers.") ;
/** sigma(n), sum of divisors of n.
* @param[in] n The n of which the divisors are taken.
* @return sigma_1(n).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000203(n)={
sigma(n,1)
}
addhelp(A000203,"sigma(n), sum of divisors of n") ;
/** The n-th triangular number.
* @param[in] n The index into the list, n>=0.
* @return n*(n+1)/2.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000217(n)={
n*(n+1)/2
}
addhelp(A000217,"n-th Triangular number") ;
/** The vector of the first n triangular number.
* @param[in] n The length of the vector returned.
* @return The first n elements of [0,1,3,6,10,...]
* @since 2010-09-08
* @author R. J. Mathar
*/
A000217_vec(n)={
vector(n,i, A000217(i-1))
}
addhelp(A000217_vec,"Vector of Triangular numbers 0 to n-1") ;
/** 2^n-1.
* @param[in] n The exponent n.
* @return 2^n minus 1.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000225(n)={
if( n < 0,
error("A000225 argument ", n, " <0.") ;
) ;
2^n-1
}
addhelp(A000225,"2^n-1.") ;
/** The n-th power of 3.
* @param[in] n The exponent. >=0.
* @return 3^n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000244(n)={
3^n
}
addhelp(A000244,"3^n.") ;
/** n squared.
* @param[in] n The basis.
* @return n^2.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000290(n)={
n^2
}
addhelp(A000290,"The square of n.") ;
/** n*(n+1)*(n+2)/6. n-th tetrahedral number.
* @param[in] n The argument, >=0.
* @return The polynomial evaluated at n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000292(n)={
if( n < 0,
error("A000292 argument ", n, " <0.") ;
) ;
n*(n+1)*(n+2)/6
}
addhelp(A000292,"N-th tetrahedral number.") ;
/** Vector of the first n th tetrahedral numbers.
* @param[in] n The length of the vector constructed.
* @return The vector starting [0,1,4,10,20,...]
* @since 2010-09-08
* @author R. J. Mathar
*/
A000292_vec(n)={
vector(n,i,A0000292(i-1)) ;
}
addhelp(A000292_vec,"Vector with the n smallest tetrahedral numbers.") ;
/** Test whether the argument is a tetrahedral number.
* @param[in] n The argument.
* @return True if n is of the form k(k+1)(k+2)/6 with k>=0.
* @since 2010-09-08
* @author R. J. Mathar
*/
isA000292(n)={
local(k) ;
if( n < 0,
return(0),
if(n<=1,
return(1)
)
) ;
/* start from an estimate T(k) >= k^3/6. Then check k downwards.
*/
k=floor((6*n)^(1/3)) ;
while( A000292(k)> n,
k--
) ;
A000292(k) ==n
}
addhelp(isA000292,"True if n is a tetrahedral number.") ;
/** n cubed.
* @param[in] n The basis.
* @return n^3.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000578(n)={
n^3
}
addhelp(A000578,"n^3.") ;
/** TEst whether the argument is a positive 3rd power.
* @param[in] n The prospective perfect cube.
* @return True if n is of the form k^3, k>=0.
* @since 2010-09-08
* @author R. J. Mathar
*/
isA000578(n)={
if(n<0,
return(0),
ispower(n,3)
)
}
addhelp(isA000578,"True if n is a positive perfect cube.") ;
/** Prime-pi(n). The number of primes <=n.
* @param[in] n The limit of the prime count.
* @return pi(n).
* @since 2010-09-08
* @author R. J. Mathar
*/
A000720(n)={
primepi(n)
}
addhelp(A000720,"pi(n), the number of primes <=n.") ;
/** Test whether the argument is a prime power p^k, k>=0.
* @param n The argument >=1.
* @return True if n is of the form p^k, p a prime, k>=0.
* @since 2010-09-08
* @author R. J. Mathar
*/
isA000961(n)={
if( n < 0,
return(0),
/* true if either 1 or having one distinct prime factor
*/
if(n==1,
return(1),
omega(n)==1
) ;
) ;
}
addhelp(isA000961,"True if prime power p^k, k>=0.") ;
/** The next prime power larger than n.
* @param n The argument and lower limit.
* @return The smallest p^k>n, p a prime and k>=0.
* @since 2010-09-08
* @author R. J. Mathar
*/
A000961_next(n)={
local(a=n+1) ;
if( n < 1,
return(1),
while( !isA000961(a),
a++
) ;
a
) ;
}
addhelp(A000961_next,"Smallest prime power >n.") ;
/** Vector of the n smallest prime powers.
* @param n The length of the vector returned.
* @return The first elements of [1,2,3,4,5,7,8,9,11,13,..]
* @since 2010-09-08
* @author R. J. Mathar
*/
A000961_vec(n)={
local(a=vector(n), pk=1) ;
for(i=1,n,
a[i] = pk ;
pk = A000961_next(pk) ;
) ;
a
}
addhelp(A000961_vec,"The vector of the n first prime powersn.") ;
/** Central binomial coefficient binomia(2n,n).
* @param n The argument >=0.
* @return binomia(2*n,n)
* @since 2010-09-08
* @author R. J. Mathar
*/
A000984(n)={
if( n < 0,
error("A000984 argument ", n, " <0.") ;
) ;
binomial(2*n,n)
}
addhelp(A000984,"binomia(2*n,n), central binomial coefficient.") ;
/** n-th Motzkin number. Nonintersection chords between n labeled points on a circle.
* @param[in] n The index into the sequence.
* @return The element a(n) in a(0)=a(1)=1, a(2)=2 etc
* @since 2010-09-08
* @author R. J. Mathar
*/
A001006(n)={
sum(k=0,n,binomial(n,2*k)*A000108(k) )
}
addhelp(A001006,"Motzkin(n).") ;
/** Jacobsthal(n).
* @param[in] n The index into the sequence.
* @return The element a(n) in a(0)=0, a(1)=1, a(n)=a(n-1)+2*a(n-2).
* @since 2010-09-08
* @author R. J. Mathar
*/
A001045(n)={
linHRec([0,1],[1,2],n+1)
}
addhelp(A001045,"Jacobsthal(n).") ;
/** The vector of the first n Jacobsthal numbers.
* @param[in] n The lenght of the vector to be returned.
* @return The J(0..n-1).
* @since 2010-09-08
* @author R. J. Mathar
*/
A001045_vec(n)={
linHRec_vec([0,1],[1,2],n)
}
addhelp(A001045_vec,"The first n Jacobsthal numbers.") ;
/** Central binomial coefficient binomial(n,floor(n/2)).
* @param n The argument >=0.
* @return binomial(n,[n/2])
* @since 2010-09-08
* @author R. J. Mathar
*/
A001405(n)={
if( n < 0,
error("A001405 argument ", n, " <0.") ;
) ;
binomial(n,n\/2)
}
addhelp(A001405,"binomial(n,[n/2]), central binomial coefficient.") ;
/** The first n Central binomial coefficients binomial(k,floor(k/2)).
* @param n The length of the vector returned.
* @return The vector starting 1,1,2,3,6,10,....
* @since 2010-09-08
* @author R. J. Mathar
*/
A001405_vec(n)={
if( n < 0,
error("A001405_vec argument ", n, " <0.") ;
) ;
vector(n,i,A001405(i-1))
}
addhelp(A001405_vec,"Vector of binomial(k,[k/2]), 0<=k<n.") ;
/** The 1*3*5*...(2n-1) double factorial.
* @param[in] n The index into the odd double factorials, >=0.
* @return (2n-1)!!
* @since 2010-09-08
* @author R. J. Mathar
*/
A001147(n)={
if( n < 0,
error("A001147 argument ", n, " <0.") ;
) ;
prod(i=1,n,2*i-1)
}
addhelp(A001147,"double factorial of 2n-1.") ;
/** omega(n). The number of distinct prime divisors of n.
* @param[in] n The number to be factorized.
* @return Omega(n).
* @since 2010-09-08
* @author R. J. Mathar
*/
A001221(n)={
omega(n)
}
addhelp(A001221,"omega(n), number of distinct prime divisors.") ;
/** The vector of the first n omega(.).
* @param[in] n The length of the vector returned.
* @return omega(1..n).
* @since 2010-09-08
* @author R. J. Mathar
*/
A001221_vec(n)={
vector(n,i,omega(i))
}
addhelp(A001221_vec,"The vector omega(1..n) with counts of distinct prime divisors.") ;
/** Bigomega(n). The number of prime divisors, counted with multiplicity.
* @param[in] n The number to be factoized.
* @return big-Omega(n).
* @since 2010-09-08
* @author R. J. Mathar
*/
A001222(n)={
bigomega(n)
}
addhelp(A001222,"Big-omega(n), number of prime divisors with multiplicity.") ;
/** Vector of the first n big-omega(.).
* @param[in] n The length of the vector to be returned.
* @return big-Omega(1..n).
* @since 2010-09-08
* @author R. J. Mathar
*/
A001222_vec(n)={
vector(n,i, bigomega(i))
}
/** The square of Lucas(n).
* @param[in] n The index into the L(n).
* @return L(n)^2. This is 4 if n=0, 1 if n=1, 9 if n=2...
* @since 2010-09-08
* @author R. J. Mathar
*/
A001254(n)={
(A000032(n))^2
}
addhelp(A001254,"Lucas(n) squared") ;
/** Test for being a semiprime.
* @param[in] n The number to be tested.
* @return true if a semiprime, else false.
* @since 2010-09-08
* @author R. J. Mathar
*/
isA001358(n)={
if(n<3,
0,
bigomega(n)==2
) ;
}
addhelp(A001358,"Semiprime(n), starting with 4 at n=1.") ;
/** The next semiprime after the argument.
* @param[in] n The currently known semiprime.
* @return The smallest semiprime > n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A001358_next(n)={
local(a=n+1) ;
while( !isA001358(a),
a++
) ;
a
}
addhelp(A001358_next,"Smallest semiprime larger than n.") ;
/** The n smallest semiprimes.
* @param[in] n The length of the vector to be returned.
* @return A vector starting [4,6,9,10,..]
* @since 2010-09-08
* @author R. J. Mathar
*/
A001358_vec(n)={
local(a=vector(n)) ;
if(n>0,
a[1]=4;
for(j=2,n,
a[j] = A001358_next(a[j-1]) ;
) ;
);
a
}
addhelp(A001358_vec,"Vector of the first n semiprimes") ;
/** Product of the first n primes.
* @param n the number of factors, >=0.
* @return Primorial(n). This is 1 if n==0.
* @since 2010-09-08
* @author R. J. Mathar
*/
A002110(n)={
if( n < 0,
error("A002110 argument ", n, " <0.") ;
) ;
prod(i=1,n,prime(i))
}
/** R(n), repunits in base 10.
* @param n the index >=0.
* @return The number represented by n 1's in base 10.
* @since 2010-09-08
* @author R. J. Mathar
*/
A002275(n)={
if( n < 0,
error("A002275 argument ", n, " <0.") ;
) ;
(10^n-1)/9
}
/** Test n against being oblong.
* @param[in] n The argument.
* @return True if n=k*(k+1) for some k>=0.
* @since 2010-09-08
* @author R. J. Mathar
*/
isA002378(n)={
if(n<0,
return(0),
issquare(1+4*n)
)
}
addhelp(isA002378,"Test argument against the oblong form k*(k+1).") ;
/** n*(n+1). Oblong number.
* @param[in] n The argument.
* @return n*(n+1).
* @since 2010-09-08
* @author R. J. Mathar
*/
A002378(n)={
n*(n+1)
}
addhelp(A002378,"Oblong number n*(n+1).") ;
/** Test if the argument is a composite number.
* @param n The argument to be tested for primeness.
* @return True if a number of the form x*y wiht x,y>1.
* @since 2010-09-08
* @author R. J. Mathar
*/
isA002808(n)={
if ( n <4,
0,
! isprime(n)
);
}
addhelp(isA002808,"Test wether n is composite.") ;
/** The smallest composite number >n.
* @param n The lower limit set.
* @return The smallest composite >n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A002808_next(n)={
if ( n <4,
4,
if( isprime(n+1),
n+2,
n+1
);
);
}
addhelp(A002808_next,"Smallest composite >n.") ;
/** The n-th composite number.
* @param n The index into the list, >=1.
* @return The n-th composite number.
* @since 2010-09-08
* @author R. J. Mathar
*/
A002808(n)={
local (a=3) ;
if( n < 0,
error("A002808 argument ", n, " <0.") ;
) ;
/* This is the very slow implementation which should
* be replaced by a primepi() call.
*/
for(i=1,n,
a = A002808_next(a)
);
a
}
addhelp(A002808,"The n-th composite number.") ;
/** Test the argument against being square-free.
* @param n The argument to be checked >=1.
* @return True if n is squarefree.
* @since 2010-09-08
* @author R. J. Mathar
*/
isA005117(n)={
if( n < 1,
error("isA005117 argument ", n, " <1.") ;
) ;
issquarefree(n)
}
addhelp(isA005117,"True if n is squarefree >=1.") ;
/** The smallest squarefree number >n.
* @param n The argument and lower limit.
* @return The smallest number larger than n which is squarefree.
* @since 2010-09-08
* @author R. J. Mathar
*/
A005117_next(n)={
local(a=n+1) ;
if( n < 1,
return(1)
) ;
while( !issquarefree(a),
a++
) ;
a
}
addhelp(A005117_next,"The smallest squarefree number >n.") ;
/** The n-th squarefree number.
* @param n The argument >=1 .
* @return The n-th element in [1,2,3,5,6,7,10,...]
* @since 2010-09-08
* @author R. J. Mathar
*/
A005117(n)={
local(a=0) ;
if( n < 1,
error("A005117 argument ", n, " <1.") ;
) ;
for(i=1,n,
a = A005117_next(a)
);
a
}
addhelp(A005117,"The n-th squarefree number.") ;
/** The smallest number with exactly n divisors.
* @param n The argument >=1 .
* @return The n-th element in [1,2,4,6,16,12,64,...]
* @since 2008-05-26
* @author R. J. Mathar
*/
A005179(n)={
local(pf=prodR(n, n), a=1, b) ;
for(i=1, length(pf),
b=prod(j=1, length(pf[i]), prime(j)^(pf[i][j]-1)) ;
if(b<a || i==1,
a=b ) ;
) ;
return(a) ;
}
addhelp(A005179,"The smallest number with n divisors.") ;
/** 2n+1.
* @param n the index >=0.
* @return the associated odd number.
* @since 2010-09-08
* @author R. J. Mathar
*/
A005408(n)={
if( n < 0,
error("A005408 argument ", n, " <0.") ;
) ;
2*n+1
}
addhelp(A005408,"Odd number 2n+1.") ;
/** 2*n.
* @param n The argument >=0 .
* @return 2 times n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A005843(n)={
if( n < 0,
error("A005843 argument ", n, " <0.") ;
) ;
2*n
}
addhelp(A005843,"2*n.") ;
/** Largest prime dividing n.
* @param n The number to be factored.
* @return The largest of the prime factors, 1 if n=1.
* @since 2010-09-08
* @author R. J. Mathar
*/
A006530(n)={
if( n < 1,
error("A006530 argument ", n, " <=0.") ;
) ;
if( n==1,
return(1)
) ;
vecmax(factor(n)[,1]) ;
}
addhelp(A006530,"Largest prime factor of n.") ;
/** Vector of largest prime factors for for 1..n.
* @param n The length of the vector returned.
* @return The vector [1,2,3,2,5,3...]
* @since 2010-09-08
* @author R. J. Mathar
*/
A006530_vec(n)={
vector(n,i,A006530(i))
}
addhelp(A006530_vec,"Vector of largest prime factors of 1..n.") ;
/** Binomial(n,m). Pascal's triangle.
* @param[in] n The numerator in the expression.
* @param[in] m The denominator in the expression.
* @return (n choose m).
* @since 2010-09-08
* @author R. J. Mathar
*/
A007318(n,m)=
{
binomial(n,m)
}
addhelp(A007318,"binomial(n,m)") ;
/** Stirling2(n,k).
* @param[in] n
* @param[in] k
* @return Stirling number of the second kind.
* @since 2010-09-08
* @author R. J. Mathar
*/
A008277(n,k)=
{
sum(i=0,k, (-1)^(k-i)*binomial(k,i)*i^n) /k!
}
addhelp(A008277,"Stirling2(n,m).") ;
/** Vector of Stirling2(n,m) for 1<=m<=n.
* @param[in] n
* @return Stirling numbers of the second kind.
* @since 2010-09-08
* @author R. J. Mathar
*/
A008277_row(n)=
{
vector(n,i, A008277(n,i))
}
addhelp(A008277_row,"Stirling2(n,1..n).") ;
/** Mu(n), Moebius function.
* @param[in] n The argument.
* @return mu(n), 0 if n is not squarefree, otherwise (-1) to the power of the number of distinct prime divisors of n.
* @since 2010-09-08
* @author R. J. Mathar
*/
A008683(n)={
moebius(n)
}
addhelp(A008684,"Moebius(n).") ;
/** Vector of the Moebius function of 1..n.
* @param[in] n The length of the vector returned.
* @return mu(1..n).
* @since 2010-09-08
* @author R. J. Mathar
*/
A008683_vec(n)={
vector(n,i,moebius(i))
}
/** Unitary sigma.
* @param[in] n The argument
* @return usigma(n), sum of the divisors d where gcd(d,n/d)=1.
* @since 2000
* @author Rick Shepherd
*/
A034448(n)={
sumdiv(n, d, if(gcd(d, n/d)==1, d))
}
/** (-1)sigma(n)
* @param[in] n The parameter >=1.
* @return (-1)^omega(n) * sum_{d|n} d*(-1)^omega(d).
* @author R. J. Mathar
* @since 2006-10-12
*/
A049060(n)={
local(i,resul,rmax,p) ;
if(n==1,
return(1),
i=factor(n) ;
rmax=matsize(i)[1] ;
resul=1 ;
for(r=1,rmax,
p=sum(j=1,i[r,2],i[r,1]^j) ;
resul *= p-1 ;
) ;
resul ;
) ;
}
addhelp(A049060,"(-1)sigma(n)") ;
/** Coefficient [x^m] of Chebyshev's S(n,x).
* @param[in] n The lower index parameter n >=0.
* @param[in] m The exponent of x.
* @return The coefficient of x^m of S(n,x).
* @author R. J. Mathar
* @since 2010-09-08
*/
A049310(n,m)={
if(n<0 || m<0 ||m>n || ((n+m)% 2 ),
return(0)
);
(-1)^((n+m)/2+m)*binomial((n+m)/2,m)
}
addhelp(A049310,"Coefficient [x^m] Chebyshev-S(n,x).") ;
/** Coefficients Chebyshev's S(n,x) for x^0 up to x^n.
* @param[in] n The lower index parameter n >=0.
* @return The vector of coefficients, exponents of x increasing.
* @author R. J. Mathar
* @since 2010-09-08
*/
A049310_row(n)={
vector(n+1,i, A049310(n,i-1))
}
addhelp(A049310_row,"Coefficients x^0, x^1 etc of Chebyshev-S(n,x).") ;
/** Increase the preceding number, of index n-1, by the smallest number which
* has the same number of divisors as the preceding number.
* @param[in] n The argument.
* @since 2010-10-02
* @author R. J. Mathar
*/
A175300(n)={
local(a=1,aprev) ;
if(n==1,
return(1),
aprev=A175300(n-1) ;
return (aprev+A005179(sigma(aprev,0)) );
) ;
}
