This site is supported by donations to The OEIS Foundation.

User:R. J. Mathar/Pari

From OeisWiki

Jump to: navigation, search
/**
* \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)) );
        ) ;
}

Personal tools