This site is supported by donations to The OEIS Foundation.

User:Stanislav Sykora

From OeisWiki

Jump to: navigation, search

A physicist, born in 1941 in Czechoslovachia, and presently resident in Italy. Educated in Physics (1963, CVUT in Prague) and Chemical Physics (1967, Czechoslovak Academy of Sciences), worked in various places around the world, mostly as a specialist in Magnetic Resonance. Longtime interest in Mathematics (both applied and a hobby). More details are on my website [1], including a detailed CV (

PARI/GP scripts

In this Section I want to share some of my PARI/GP scripts. Notes:

1. All scripts were tested using the latest PARI version. However, if you find any problem, please let me know.
2. The free PARI and GP software is available at
3. All names of functions defined here, unlike those in PARI library, start with a capital letter.
4. Function names starting with D, like Dname(x), denote derivatives, like that of name(x).
5. Function names starting with AD, like ADname(x), denote anti-derivatives, like that of name(x).
6. Anti-derivatives are unique up to a constant; they are intended for the evaluation of definite integrals.

=== Generic UTILITIES


This file defines two global variables Eps_ and MaxIter_, used in a number of the author's PARI/GP scripts. Eps_ is used primarily to test near-equality of REAL values and MaxIter_ is used to stop runaway iterations. When loaded, it sets the current realprecision to 50 digits, and initializes Eps_ and MaxIter_ accordingly. To change the settings, use the SetEbDefaults function.
SetEbDefaults(dgs=-1,eps=0.0,imx=-1) - Sets desired default precision and the global variables Eps_ and MaxIter_
Floor(z) - Extends the PARI floor function to all numeric type arguments (including t_COMPLEX)
Round(z) - Implements absolute-value rounding of t_REAL and t_COMPLEX arguments
IsInt(z) - Tests whether z matches an integer (within the +-Eps_ tollerance interval)
IsIntLt0(z) - Tests whether z matches a negative integer (within the +-Eps_ tollerance interval)
IsIntLe0(z) - Tests whether z matches a nonpositive integer (within the +-Eps_ tollerance interval)
IsIntGt0(z) - Tests whether z matches a positive integer (within the +-Eps_ tollerance interval)
IsIntGe0(z) - Tests whether z matches a nonnegative integer (within the +-Eps_ tollerance interval)
IsEq(z1,z2) - Tests whether z1 and z2 match each other (within the +-Eps_ tollerance interval)



Provides support for various types of formatted output into files to be used elsewhere (e.g., as input to Excell. Matlab, or as OEIS entries, etc).
Presently contains only output (dumping) functions.
ExpFormat(z,d=0,head="",sep=" ",tail="") - Returns a string with the numeric argument formatted as mantissa-exponent with d digits.
Real2ExpFormat(x,d=0) - Returns a string with the t_REAL argument x formatted in mantissa-exponent with d digits.
TabFuncArr(f,x,d=0,s=" ",file="",head="") - Tabulates the values of a function for an array of arguments.
TabFuncLinX(f,xini,xstep,xend,d=0,s=" ",file="",head="") - Tabulates the values of a function computed for a linear grid of arguments.
Real2OEIS(file,cnst,dgts=1000) - Writes into a file all one needs to submit to OEIS a real constant.
IntSeq2OEIS(file,vect,offset,dumpall=0) - Writes into a text file all one needs to submit to OEIS an integer sequence.



Utilities to test and handle vectors of various types
- Testing vectors:
IsIncreasingVector(v,strict=0) - tests the monotonicity of v (increasing or non-decreasing).
IsDecreasingVector(v,strict=0) - tests the monotonicity of v (decreasing or non-increasing.
- Vector comparisons:
VectorsDifferAt(v1,v2) - returns first index k where v1[k]!=v2[k] (or 0 if none).
VectorsMatchAt(v1,v2) - returns first index k where v1[k]==v2[k] (or 0 if none).
VectorLtVectorAt(v1,v2) - returns first index k where v1[k]< v2[k] (or 0 if none).
VectorLeVectorAt(v1,v2) - returns first index k where v1[k]<=v2[k] (or 0 if none).
VectorGtVectorAt(v1,v2) - returns first index k where v1[k]> v2[k] (or 0 if none).
VectorGeVectorAt(v1,v2) - returns first index k where v1[k]>=v2[k] (or 0 if none).
VectorLtExprAt(v1,expr) - returns first index k where v1[k]< expr(k) (or 0 if none).
VectorLeExprAt(v1,expr) - returns first index k where v1[k]<=expr(k) (or 0 if none).
VectorGtExprAt(v1,expr) - returns first index k where v1[k]> expr(k) (or 0 if none).
VectorGeExprAt(v1,expr) - returns first index k where v1[k]>=expr(k) (or 0 if none).
- Operating on vectors:
ReverseVector(v) - returns vector with reversed order of elements 1..#v (or i..j).
RunningSum(v) - returns the running sum of v.
RunningSumInverse(v) - inverse of the running sum mapping.
RunningProduct(v) - returns the running product of v.
RunningProductInverse(v) - inverse of the running product mapping.
- Lists of natural numbers 1<=a(n)<=nmax satisfying a condition:
ListNCond0(nmax,condition,anmax=0) - lists natural numbers satisfying a simple condition(n) with no parameters. Stops when n reaches anmax (if anmax>0) or when nmax entries were found.
ListVFCond0(cond,file,fmonit,monit=1000,nini=1,idxini=1,nmax=0,idxmax=0) - as above, but writes the results into a file and permits to start at any initial value. Optionally allows to monitor the progress.
DensityNCond0(nmax,condition,monit=0) - density of numbers satisfying the condition among n = 1..nmax; with optional progress monitoring.
ListNCond1(nmax,condition,p,anmax=0) - lists natural numbers satisfying a condition(n,m) with a settable parameter p. Stops when n reaches anmax (if anmax>0) or when nmax entries were found.
ListVFCond1(cond,p,file,fmonit,monit=1000,nini=1,idxini=1,nmax=0,idxmax=0) - as above, but writes the results into a file and permits to start at any initial value; with optional progress monitoring.
DensityNCond1(nmax,condition,m,monit=0) - density of numbers satisfying the condition among n = 1..nmax; with optional progress monitoring.
- Generating functions management:
Egf2CoefSeq(func,nc=20,file="",c0=0.12345) - Generates a vector of the real e.g.f. coefficients for an analytical function func(z). Displays (monitors) the progress and optionally writes the results into an optional file.



Scripts, mostly very simple, related to primes and factorization:
- Boolean prime-related functions:
IsMoebiusZero(n) - Equivalent to "is divisible by a square".
IsMoebiusNonZero(n) - Equivalent to "is not divisible by a square".
- Integer valued prime-related functions:
LargestPrimeFactor(n) - known also as gpf(n).
SumOfPrimeFactors(n) - with multiplicity.
SumOfPrimeFactorsN(n) - without multiplicity.
ProductOfPrimeFactorsN(n) - without multiplicity (else the result would be always n).
PreviousPrime(n) - as precprime(n-1), but for n<=2 returns n.
FollowingPrime(n) - as nextprime(n+1), but for n<2 returns n.
NumberOfNonCoprimes(n) - number of m, 1<m<n, which are not coprime to n.
IntervalPrimesProduct(m,n) - product of all primes in the interval (m,n].
Lambda(n) - central binomial(n,(n\2)), purged of all prime factors exceeding (n+1)\2.
A233511(n) - replace largest prime factor > 2 (if any) with PreceedingPrime.
A233570(n) - replace smallest prime factor (if any) with FollowingPrime.
AltSum1DivPrimePwr(x,eps) - helps to evaluate the series C(x) = -Sum[k]((-1)^k/p(k)^x), where p(k) is the k-th prime.
- These use DigitBasedFunctions.txt:
IsReversedPrime(n,b) - returns 1 if the number obtained from n by reversing its digits in base b is prime.
IsReversiblePrime(n,b) - returns 1 if n is prime AND IsReversedPrime(n,b) returns 1.



Generators of sequences which have something to do with primes and factorizations:
PrimeByMplusD(nmax,m,d) - lists primes p for which m*p+d is also a prime (draft)
CodedPrimesSequence(nmax) - numbers which are either a non-twin prime (odd or 2), or are surrounded by twin primes (even>2).



These functions are expansion-base dependent (the default is usually 10):
DigitsToNumber(d,b) - converts a vector of digits d into an integer, assuming base b.
ReverseDigits(n,b) - in base b, reverses the digits of n.
IsReversedGE(n,b) - tests if, in base b, ReverseDigits(n) >= n.
IsPalindrome(n,b) - tests if, in base b, ReverseDigits(n) == n.
NumberOfDigits(n,b) - number of digits of n in base b.
SumOfDigits(n,b) - sum of digits of n in base b.
ProductOfDigits(n,b) - product of digits of n in base b.
DigitsSumProduct(n,b) - sum of digits of n, multiplied by their product, in base b.
DigitalRoot(n,b) - limit of the iterated sum of digits of n, in base b.
RTruncIsCoprime(n,b) - true when n is coprime to the number obtained by truncating its rightmost digit in base b.
RTruncIsNotCoprime(n,b) - logical negation of RTruncIsCoprime(n,b).
- These refer to binary expansion digits:
Parity(n) - returns 1 when the sum of binary digits is odd, 0 when it is even.
GrayCode(n) - the binary Gray code corresponding to n.
GrayDecode(g) - inverse of the above; the number n corresponding to the Gray code g.



Sequences for endomorphism E(n)->n\b:
GT_Trunc1(nmax,property,b) - Preservation of a property(n) under right-truncation in base b.
GT_Trunc2(nmax,property,b) - Preservation of a property(n,b) under right-truncation in base b.
GT_Trunc3(bmin,bmax,property,outfile) - Preservation of a property(n,b) under right-truncation in bases bmin..bmax; multiple data output.



Generators of various integer sequences (see the file for pertinent sequence numbers):
LiMinusNpDivq(nmax,p,q)(n) - Related to polylogarithm Li(order,x): a(n) = Li(-n,-p/q)*(p+q)^(n+1)/q.



Functions for various optimization tasks (searches for minimum, maximum, etc):
MinSearch(ax,bx,cx,f,monit=0,itermax=10000) - searches for a minimum of a function with one real argument.



Functions related to the Lambert W-function:
LambertW0(x) - the upper branch.
LambertW1(x) - the lower branch.
DLambertW0(x) - derivative of LambertW0.
DLambertW1(x) - derivative of LambertW1.
ADLambertW0(x) - anti-derivative of LambertW0.
ADLambertW1(x) - anti-derivative of LambertW1.


Functions related to the centered normal (Gaussian) errors distribution with unit variance:
NormalErrorDf(x) - normal error density function.
NormalErrorCdf(x) - normal error cumulative distribution function, also called Phi(x).
NormalErrorPercentile(p) - normal error percentile function, also called probit(p).
RngNormal(mean,stddev)(z) - generator of random values with a normal distribution.
Erf(z) - the basic error function (PARI library contains only the complementary erfc(z).
Erfi(z) - imaginary error function.
Fadeevaw(z) - Fadeeva w(z) function.
DErf - first derivative of erf(z).
Derfc(z) - first derivative of erfc(z).
ADErf(z) - anti-derivative of erf(z).
ADerfc(z) - anti-derivative of erfc(z).
Dawson(z) - Dawson integral F(z).
DDawson(z) - first derivative of the Dawson integral.
D2Dawson(z) - second derivative of the Dawson integral.


Brute-force computation of the infinite power tower of x, i.e., T(x) = lim[n->inf] x^^n.


Before using this, please consult Blazys Expansions and Continued Fractions.
Particularly, pay attention to the required precision of input values.
Bx(x,nmax) - returns nmax integers of the Blazys expansion of an irrational x.
Ebx(x,nmax) - the same but using the extended expansion.
Bf(v) - inverse of Bx; evaluates the Blazys continued fraction.
Ebf(v) - the same as above but using the extended version.
BxPeriodicSqrtN(nmax) - lists numbers n for which Bx(sqrt(n)) is periodic.
BxAperiodicSqrtN(nmax) - lists numbers n for which Bx(sqrt(n)) is a-periodic.


Functions related to the Exponential and Logarithmic Integrals.
(complex-valued, accepting any real or complex arguments).
Ei(z) - exponential integral function.
li(z) - logarithmic integral function.
Si(z) - sine integral function.
Ci(z) - cosine integral function.
Cin(z) - entire cosine integral function.
Personal tools