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