| Title: | Number-Theoretic Functions |
|---|---|
| Description: | Provides number-theoretic functions for factorization, prime numbers, twin primes, primitive roots, modular logarithm and inverses, extended GCD, Farey series and continued fractions. Includes Legendre and Jacobi symbols, some divisor functions, Euler's Phi function, etc. |
| Authors: | Hans Werner Borchers |
| Maintainer: | Hans W. Borchers <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.8-5 |
| Built: | 2024-09-06 06:49:59 UTC |
| Source: | CRAN |
Provides number-theoretic functions for factorization, prime numbers, twin primes, primitive roots, modular logarithm and inverses, extended GCD, Farey series and continued fractions. Includes Legendre and Jacobi symbols, some divisor functions, Euler's Phi function, etc.
The DESCRIPTION file:
| Package: | numbers |
| Type: | Package |
| Title: | Number-Theoretic Functions |
| Version: | 0.8-5 |
| Date: | 2022-11-22 |
| Author: | Hans Werner Borchers |
| Maintainer: | Hans W. Borchers <[email protected]> |
| Depends: | R (>= 4.1.0) |
| Suggests: | gmp (>= 0.5-1) |
| Description: | Provides number-theoretic functions for factorization, prime numbers, twin primes, primitive roots, modular logarithm and inverses, extended GCD, Farey series and continued fractions. Includes Legendre and Jacobi symbols, some divisor functions, Euler's Phi function, etc. |
| License: | GPL (>= 3) |
| NeedsCompilation: | no |
| Packaged: | 2022-11-22 18:59:22 UTC; hwb |
| Repository: | CRAN |
| Date/Publication: | 2022-11-23 08:10:02 UTC |
Index of help topics:
GCD GCD and LCM Integer Functions Primes Prime Numbers Sigma Divisor Functions agm Arithmetic-geometric Mean bell Bell Numbers bernoulli_numbers Bernoulli Numbers carmichael Carmichael Numbers catalan Catalan Numbers cf2num Generalized Continous Fractions chinese Chinese Remainder Theorem collatz Collatz Sequences contfrac Continued Fractions coprime Coprimality div Integer Division divisors List of Divisors dropletPi Droplet Algorithm for pi and e egyptian_complete Egyptian Fractions - Complete Search egyptian_methods Egyptian Fractions - Specialized Methods eulersPhi Eulers's Phi Function extGCD Extended Euclidean Algorithm fibonacci Fibonacci and Lucas Series hermiteNF Hermite Normal Form iNthroot Integer N-th Root isIntpower Powers of Integers isNatural Natural Number isPrime isPrime Property isPrimroot Primitive Root Test legendre_sym Legendre and Jacobi Symbol mersenne Mersenne Numbers miller_rabin Miller-Rabin Test mod Modulo Operator modinv Modular Inverse and Square Root modlin Modular Linear Equation Solver modlog Modular (or: Discrete) Logarithm modpower Power Function modulo m moebius Moebius Function necklace Necklace and Bracelet Functions nextPrime Next Prime numbers-package Number-Theoretic Functions omega Number of Prime Factors ordpn Order in Faculty pascal_triangle Pascal Triangle periodicCF Periodic continued fraction previousPrime Previous Prime primeFactors Prime Factors primroot Primitive Root pythagorean_triples Pythagorean Triples quadratic_residues Quadratic Residues ratFarey Farey Approximation and Series rem Integer Remainder solvePellsEq Solve Pell's Equation stern_brocot_seq Stern-Brocot Sequence twinPrimes Twin Primes zeck Zeckendorf Representation
Although R does not have a true integer data type, integers can be represented exactly up to 2^53-1 . The numbers package attempts to provided basic number-theoretic functions that will work correcty and relatively fast up to this level.
Hans Werner Borchers
Maintainer: Hans W. Borchers <[email protected]>
Hardy, G. H., and E. M. Wright (1980). An Introduction to the Theory of Numbers. 5th Edition, Oxford University Press.
Riesel, H. (1994). Prime Numbers and Computer Methods for Factorization. Second Edition, Birkhaeuser Boston.
Crandall, R., and C. Pomerance (2005). Prime Numbers: A Computational Perspective. Springer Science+Business.
Shoup, V. (2009). A Computational Introduction to Number Theory and Algebra. Second Edition, Cambridge University Press.
Arndt, J. (2010). Matters Computational: Ideas, Algorithms, Source Code. 2011 Edition, Springer-Verlag, Berlin Heidelberg.
Forster, O. (2014). Algorithmische Zahlentheorie. 2. Auflage, Springer Spektrum Wiesbaden.
The arithmetic-geometric mean of real or complex numbers.
agm(a, b)agm(a, b)
a, b
|
real or complex numbers. |
The arithmetic-geometric mean is defined as the common limit of the two
sequences and .
Returnes one value, the algebraic-geometric mean.
The calculation of the AGM is continued until the two values of a and
b are identical (in machine accuracy).
Borwein, J. M., and P. B. Borwein (1998). Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity. Second, reprinted Edition, A Wiley-interscience publication.
Arithmetic, geometric, and harmonic mean.
## Gauss constant 1 / agm(1, sqrt(2)) # 0.834626841674073 ## Calculate the (elliptic) integral 2/pi \int_0^1 dt / sqrt(1 - t^4) f <- function(t) 1 / sqrt(1-t^4) 2 / pi * integrate(f, 0, 1)$value 1 / agm(1, sqrt(2)) ## Calculate pi with quadratic convergence (modified AGM) # See algorithm 2.1 in Borwein and Borwein y <- sqrt(sqrt(2)) x <- (y+1/y)/2 p <- 2+sqrt(2) for (i in 1:6){ cat(format(p, digits=16), "\n") p <- p * (1+x) / (1+y) s <- sqrt(x) y <- (y*s + 1/s) / (1+y) x <- (s+1/s)/2 } ## Not run: ## Calculate pi with arbitrary precision using the Rmpfr package require("Rmpfr") vpa <- function(., d = 32) mpfr(., precBits = 4*d) # Function to compute \pi to d decimal digits accuracy, based on the # algebraic-geometric mean, correct digits are doubled in each step. agm_pi <- function(d) { a <- vpa(1, d) b <- 1/sqrt(vpa(2, d)) s <- 1/vpa(4, d) p <- 1 n <- ceiling(log2(d)); for (k in 1:n) { c <- (a+b)/2 b <- sqrt(a*b) s <- s - p * (c-a)^2 p <- 2 * p a <- c } return(a^2/s) } d <- 64 pia <- agm_pi(d) print(pia, digits = d) # 3.141592653589793238462643383279502884197169399375105820974944592 # 3.1415926535897932384626433832795028841971693993751058209749445923 exact ## End(Not run)## Gauss constant 1 / agm(1, sqrt(2)) # 0.834626841674073 ## Calculate the (elliptic) integral 2/pi \int_0^1 dt / sqrt(1 - t^4) f <- function(t) 1 / sqrt(1-t^4) 2 / pi * integrate(f, 0, 1)$value 1 / agm(1, sqrt(2)) ## Calculate pi with quadratic convergence (modified AGM) # See algorithm 2.1 in Borwein and Borwein y <- sqrt(sqrt(2)) x <- (y+1/y)/2 p <- 2+sqrt(2) for (i in 1:6){ cat(format(p, digits=16), "\n") p <- p * (1+x) / (1+y) s <- sqrt(x) y <- (y*s + 1/s) / (1+y) x <- (s+1/s)/2 } ## Not run: ## Calculate pi with arbitrary precision using the Rmpfr package require("Rmpfr") vpa <- function(., d = 32) mpfr(., precBits = 4*d) # Function to compute \pi to d decimal digits accuracy, based on the # algebraic-geometric mean, correct digits are doubled in each step. agm_pi <- function(d) { a <- vpa(1, d) b <- 1/sqrt(vpa(2, d)) s <- 1/vpa(4, d) p <- 1 n <- ceiling(log2(d)); for (k in 1:n) { c <- (a+b)/2 b <- sqrt(a*b) s <- s - p * (c-a)^2 p <- 2 * p a <- c } return(a^2/s) } d <- 64 pia <- agm_pi(d) print(pia, digits = d) # 3.141592653589793238462643383279502884197169399375105820974944592 # 3.1415926535897932384626433832795028841971693993751058209749445923 exact ## End(Not run)
Generate Bell numbers.
bell(n)bell(n)
n |
integer, asking for the n-th Bell number. |
Bell numbers, commonly denoted as , are defined as the number of
partitions of a set of n elements. They can easily be calculated
recursively.
Bell numbers also appear as moments of probability distributions, for example
B_n is the n-th momentum of the Poisson distribution with mean 1.
A single integer, as long as n<=22.
sapply(0:10, bell) # 1 1 2 5 15 52 203 877 4140 21147 115975sapply(0:10, bell) # 1 1 2 5 15 52 203 877 4140 21147 115975
Generate the Bernoulli numbers.
bernoulli_numbers(n, big = FALSE)bernoulli_numbers(n, big = FALSE)
n |
integer; starting from 0. |
big |
logical; shall double or GMP big numbers be returned? |
Generate the n+1 Bernoulli numbers B_0,B_1, ...,B_n,
i.e. from 0 to n. We assume B1 = +1/2.
With big=FALSE double integers up to 2^53-1 will be used,
with big=TRUE GMP big rationals (through the 'gmp' package).
B_25 is the highest such number that can be expressed as an
integer in double float.
Returns a matrix with two columns, the first the numerator, the second the denominator of the Bernoulli number.
M. Kaneko. The Akiyama-Tanigawa algorithm for Bernoulli numbers. Journal of Integer Sequences, Vol. 3, 2000.
D. Harvey. A multimodular algorithm for computing Bernoulli numbers. Mathematics of Computation, Vol. 79(272), pp. 2361-2370, Oct. 2010. arXiv 0807.1347v2, Oct. 2018.
bernoulli_numbers(3); bernoulli_numbers(3, big=TRUE) ## Big Integer ('bigz') 4 x 2 matrix: ## [,1] [,2] [,1] [,2] ## [1,] 1 1 [1,] 1 1 ## [1,] 1 2 [2,] 1 2 ## [2,] 1 6 [3,] 1 6 ## [3,] 0 1 [4,] 0 1 ## Not run: bernoulli_numbers(24)[25,] ## [1] -236364091 2730 bernoulli_numbers(30, big=TRUE)[31,] ## Big Integer ('bigz') 1 x 2 matrix: ## [,1] [,2] ## [1,] 8615841276005 14322 ## End(Not run)bernoulli_numbers(3); bernoulli_numbers(3, big=TRUE) ## Big Integer ('bigz') 4 x 2 matrix: ## [,1] [,2] [,1] [,2] ## [1,] 1 1 [1,] 1 1 ## [1,] 1 2 [2,] 1 2 ## [2,] 1 6 [3,] 1 6 ## [3,] 0 1 [4,] 0 1 ## Not run: bernoulli_numbers(24)[25,] ## [1] -236364091 2730 bernoulli_numbers(30, big=TRUE)[31,] ## Big Integer ('bigz') 1 x 2 matrix: ## [,1] [,2] ## [1,] 8615841276005 14322 ## End(Not run)
Checks whether a number is a Carmichael number.
carmichael(n)carmichael(n)
n |
natural number |
A natural number n is a Carmichael number if it is a Fermat
pseudoprime for every a, that is a^(n-1) = 1 mod n, but
is composite, not prime.
Here the Korselt criterion is used to tell whether a number n is
a Carmichael number.
Returns TRUE or FALSE
There are infinitely many Carmichael numbers, specifically there should
be at least n^(2/7) Carmichael numbers up to n (for n large enough).
R. Crandall and C. Pomerance. Prime Numbers - A Computational Perspective. Second Edition, Springer Science+Business Media, New York 2005.
carmichael(561) # TRUE ## Not run: for (n in 1:100000) if (carmichael(n)) cat(n, '\n') ## 561 2821 15841 52633 ## 1105 6601 29341 62745 ## 1729 8911 41041 63973 ## 2465 10585 46657 75361 ## End(Not run)carmichael(561) # TRUE ## Not run: for (n in 1:100000) if (carmichael(n)) cat(n, '\n') ## 561 2821 15841 52633 ## 1105 6601 29341 62745 ## 1729 8911 41041 63973 ## 2465 10585 46657 75361 ## End(Not run)
Generate Catalan numbers.
catalan(n)catalan(n)
n |
integer, asking for the n-th Catalan number. |
Catalan numbers, commonly denoted as , are defined as
and occur regularly in all kinds of enumeration problems.
A single integer, as long as n<=30.
C <- numeric(10) for (i in 1:10) C[i] <- catalan(i) C[5] #=> 42C <- numeric(10) for (i in 1:10) C[i] <- catalan(i) C[5] #=> 42
Evaluate a generalized continuous fraction as an alternating sum.
cf2num(a, b = 1, a0 = 0, finite = FALSE)cf2num(a, b = 1, a0 = 0, finite = FALSE)
a |
numeric vector of length greater than 2. |
b |
numeric vector of length 1 or the same length as a. |
a0 |
absolute term, integer part of the continuous fraction. |
finite |
logical; shall Algorithm 1 be applied. |
Calculates the numerical value of (simple or generalized) continued fractions of the form
by converting it into an alternating sum and then applying the accelleration Algorithm 1 of Cohen et al. (2000).
The argument is by default set to ,
that is the continued fraction is treated in its simple form.
With finite=TRUE the accelleration is turned off.
Returns a numerical value, an approximation of the continued fraction.
This function is not vectorized.
H. Cohen, F. R. Villegas, and Don Zagier (2000). Experimental Mathematics, Vol. 9, No. 1, pp. 3-12. <www.emis.de/journals/EM>
## Examples from Wolfram Mathworld print(cf2num(1:25), digits=16) # 0.6977746579640077, eps() a = 2*(1:25) + 1; b = 2*(1:25); a0 = 1 # 1/(sqrt(exp(1))-1) cf2num(a, b, a0) # 1.541494082536798 a <- b <- 1:25 # 1/(exp(1)-1) cf2num(a, b) # 0.5819767068693286 a <- rep(1, 100); b <- 1:100; a0 <- 1 # 1.5251352761609812 cf2num(a, b, a0, finite = FALSE) # 1.525135276161128 cf2num(a, b, a0, finite = TRUE) # 1.525135259240266## Examples from Wolfram Mathworld print(cf2num(1:25), digits=16) # 0.6977746579640077, eps() a = 2*(1:25) + 1; b = 2*(1:25); a0 = 1 # 1/(sqrt(exp(1))-1) cf2num(a, b, a0) # 1.541494082536798 a <- b <- 1:25 # 1/(exp(1)-1) cf2num(a, b) # 0.5819767068693286 a <- rep(1, 100); b <- 1:100; a0 <- 1 # 1.5251352761609812 cf2num(a, b, a0, finite = FALSE) # 1.525135276161128 cf2num(a, b, a0, finite = TRUE) # 1.525135259240266
Executes the Chinese Remainder Theorem (CRT).
chinese(a, m)chinese(a, m)
a |
sequence of integers, of the same length as |
m |
sequence of natural numbers, relatively prime to each other. |
The Chinese Remainder Theorem says that given integers and
natural numbers , relatively prime (i.e., coprime) to each other,
there exists a unique solution such that the following
system of linear modular equations is satisfied:
More generally, a solution exists if the following condition is satisfied:
This version of the CRT is not yet implemented.
Returns th (unique) solution of the system of modular equalities as an
integer between 0 and M=prod(m).
m <- c(3, 4, 5) a <- c(2, 3, 1) chinese(a, m) #=> 11 # ... would be sufficient # m <- c(50, 210, 154) # a <- c(44, 34, 132) # x = 4444m <- c(3, 4, 5) a <- c(2, 3, 1) chinese(a, m) #=> 11 # ... would be sufficient # m <- c(50, 210, 154) # a <- c(44, 34, 132) # x = 4444
Generates Collatz sequences with n -> k*n+l for n odd.
collatz(n, k = 3, l = 1, short = FALSE, check = TRUE)collatz(n, k = 3, l = 1, short = FALSE, check = TRUE)
n |
integer to start the Collatz sequence with. |
k, l
|
parameters for computing |
short |
logical, abbreviate stps with |
check |
logical, check for nontrivial cycles. |
Function n, k, l generates iterative sequences starting with
n and calculating the next number as n/2 if n is
even and k*n+l if n is odd. It stops automatically when
1 is reached.
The default parameters k=3, l=1 generate the classical Collatz
sequence. The Collatz conjecture says that every such sequences will end
in the trivial cycle ...,4,2,1. For other parameters this does not
necessarily happen.
k and l are not allowed to be both even or both odd – to make
k*n+l even for n odd. Option short=TRUE calculates
(k*n+l)/2 when n is odd (as k*n+l is even in this case),
shortening the sequence a bit.
With option check=TRUE will check for nontrivial cycles, stopping
with the first integer that repeats in the sequence. The check is disabled
for the default parameters in the light of the Collatz conjecture.
Returns the integer sequence generated from the iterative rule.
Sends out a message if a nontrivial cycle was found (i.e. the sequence is not ending with 1 and end in an infinite cycle). Throws an error if an integer overflow is detected.
The Collatz or 3n+1-conjecture has been experimentally verified
for all start numbers n up to 10^20 at least.
See the Wikipedia entry on the 'Collatz Conjecture'.
collatz(7) # n -> 3n+1 ## [1] 7 22 11 34 17 52 26 13 40 20 10 5 16 8 4 2 1 collatz(9, short = TRUE) ## [1] 9 14 7 11 17 26 13 20 10 5 8 4 2 1 collatz(7, l = -1) # n -> 3n-1 ## Found a non-trivial cycle for n = 7 ! ## [1] 7 20 10 5 14 7 ## Not run: collatz(5, k = 7, l = 1) # n -> 7n+1 ## [1] 5 36 18 9 64 32 16 8 4 2 1 collatz(5, k = 7, l = -1) # n -> 7n-1 ## Info: 5 --> 1.26995e+16 too big after 280 steps. ## Error in collatz(5, k = 7, l = -1) : ## Integer overflow, i.e. greater than 2^53-1 ## End(Not run)collatz(7) # n -> 3n+1 ## [1] 7 22 11 34 17 52 26 13 40 20 10 5 16 8 4 2 1 collatz(9, short = TRUE) ## [1] 9 14 7 11 17 26 13 20 10 5 8 4 2 1 collatz(7, l = -1) # n -> 3n-1 ## Found a non-trivial cycle for n = 7 ! ## [1] 7 20 10 5 14 7 ## Not run: collatz(5, k = 7, l = 1) # n -> 7n+1 ## [1] 5 36 18 9 64 32 16 8 4 2 1 collatz(5, k = 7, l = -1) # n -> 7n-1 ## Info: 5 --> 1.26995e+16 too big after 280 steps. ## Error in collatz(5, k = 7, l = -1) : ## Integer overflow, i.e. greater than 2^53-1 ## End(Not run)
Evaluate a continued fraction or generate one.
contfrac(x, tol = 1e-12)contfrac(x, tol = 1e-12)
x |
a numeric scalar or vector. |
tol |
tolerance; default |
If x is a scalar its continued fraction will be generated up to
the accuracy prescribed in tol. If it is of length greater 1, the
function assumes this to be a continued fraction and computes its value
and convergents.
The continued fraction is assumed to be
finite and neither periodic nor infinite. For implementation uses the
representation of continued fractions through 2-by-2 matrices
(i.e. Wallis' recursion formula from 1644).
If x is a scalar, it will return a list with components cf
the continued fraction as a vector, rat the rational approximation,
and prec the difference between the value and this approximation.
If x is a vector, the continued fraction, then it will return a list
with components f the numerical value, p and q the
convergents, and prec an estimated precision.
This function is not vectorized.
Hardy, G. H., and E. M. Wright (1979). An Introduction to the Theory of Numbers. Fifth Edition, Oxford University Press, New York.
contfrac(pi) contfrac(c(3, 7, 15, 1)) # rational Approx: 355/113 contfrac(0.555) # 0 1 1 4 22 contfrac(c(1, rep(2, 25))) # 1.414213562373095, sqrt(2)contfrac(pi) contfrac(c(3, 7, 15, 1)) # rational Approx: 355/113 contfrac(0.555) # 0 1 1 4 22 contfrac(c(1, rep(2, 25))) # 1.414213562373095, sqrt(2)
Determine whether two numbers are coprime, i.e. do not have a common prime divisor.
coprime(n,m)coprime(n,m)
n, m
|
integer scalars |
Two numbers are coprime iff their greatest common divisor is 1.
Logical, being TRUE if the numbers are coprime.
coprime(46368, 75025) # Fibonacci numbers are relatively prime to each other coprime(1001, 1334)coprime(46368, 75025) # Fibonacci numbers are relatively prime to each other coprime(1001, 1334)
Integer division.
div(n, m)div(n, m)
n |
numeric vector (preferably of integers) |
m |
integer vector (positive, zero, or negative) |
div(n, m) is integer division, that is discards the fractional part,
with the same effect as n %/% m.
It can be defined as floor(n/m).
A numeric (integer) value or vector/matrix.
div(c(-5:5), 5) div(c(-5:5), -5) div(c(1, -1), 0) #=> Inf -Inf div(0,c(0, 1)) #=> NaN 0div(c(-5:5), 5) div(c(-5:5), -5) div(c(1, -1), 0) #=> Inf -Inf div(0,c(0, 1)) #=> NaN 0
Generates a list of divisors of an integer number.
divisors(n)divisors(n)
n |
integer whose divisors will be generated. |
The list of all divisors of an integer n will be calculated
and returned in ascending order, including 1 and the number itself.
For n>=1000 the list of prime factors of n will be
used, for smaller n a total search is applied.
Returns a vector integers.
divisors(1) # 1 divisors(2) # 1 2 divisors(2^5) # 1 2 4 8 16 32 divisors(1000) # 1 2 4 5 8 10 ... 100 125 200 250 500 1000 divisors(1001) # 1 7 11 13 77 91 143 1001divisors(1) # 1 divisors(2) # 1 2 divisors(2^5) # 1 2 4 8 16 32 divisors(1000) # 1 2 4 5 8 10 ... 100 125 200 250 500 1000 divisors(1001) # 1 7 11 13 77 91 143 1001
Generates digits for pi resp. the Euler number e.
dropletPi(n) dropletE(n)dropletPi(n) dropletE(n)
n |
number of digits after the decimal point; should not exceed 1000 much as otherwise it will be very slow. |
Based on a formula discovered by S. Rabinowitz and S. Wagon.
The droplet algorithm for pi uses the Euler transform of the alternating Leibniz series and the so-called “radix conversion".
String containing “3.1415926..." resp. “2.718281828..." with
n digits after the decimal point (i.e., internal decimal places).
Borwein, J., and K. Devlin (2009). The Computer as Crucible: An Introduction to Experimental Mathematics. A K Peters, Ltd.
Arndt, J., and Ch. Haenel (2000). Pi – Algorithmen, Computer, Arithmetik. Springer-Verlag, Berlin Heidelberg.
## Example dropletE(20) # [1] "2.71828182845904523536" print(exp(1), digits=20) # [1] 2.7182818284590450908 dropletPi(20) # [1] "3.14159265358979323846" print(pi, digits=20) # [1] 3.141592653589793116 ## Not run: E <- dropletE(1000) table(strsplit(substring(E, 3, 1002), "")) # 0 1 2 3 4 5 6 7 8 9 # 100 96 97 109 100 85 99 99 103 112 Pi <- dropletPi(1000) table(strsplit(substring(Pi, 3, 1002), "")) # 0 1 2 3 4 5 6 7 8 9 # 93 116 103 102 93 97 94 95 101 106 ## End(Not run)## Example dropletE(20) # [1] "2.71828182845904523536" print(exp(1), digits=20) # [1] 2.7182818284590450908 dropletPi(20) # [1] "3.14159265358979323846" print(pi, digits=20) # [1] 3.141592653589793116 ## Not run: E <- dropletE(1000) table(strsplit(substring(E, 3, 1002), "")) # 0 1 2 3 4 5 6 7 8 9 # 100 96 97 109 100 85 99 99 103 112 Pi <- dropletPi(1000) table(strsplit(substring(Pi, 3, 1002), "")) # 0 1 2 3 4 5 6 7 8 9 # 93 116 103 102 93 97 94 95 101 106 ## End(Not run)
Generate all Egyptian fractions of length 2 and 3.
egyptian_complete(a, b, show = TRUE)egyptian_complete(a, b, show = TRUE)
a, b
|
integers, a != 1, a < b and a, b relatively prime. |
show |
logical; shall solutions found be printed? |
For a rational number 0 < a/b < 1, generates all Egyptian fractions
of length 2 and three, that is finds integers x1, x2, x3 such that
a/b = 1/x1 + 1/x2a/b = 1/x1 + 1/x2 + 1/x3.
All solutions found will be printed to the console if show=TRUE;
returns invisibly the number of solutions found.
https://www.ics.uci.edu/~eppstein/numth/egypt/
egyptian_complete(6, 7) # 1/2 + 1/3 + 1/42 egyptian_complete(8, 11) # no solution with 2 or 3 fractions # TODO # 2/9 = 1/9 + 1/10 + 1/90 # is not recognized, as similar cases, # because 1/n is not considered in m/n.egyptian_complete(6, 7) # 1/2 + 1/3 + 1/42 egyptian_complete(8, 11) # no solution with 2 or 3 fractions # TODO # 2/9 = 1/9 + 1/10 + 1/90 # is not recognized, as similar cases, # because 1/n is not considered in m/n.
Generate Egyptian fractions with specialized methods.
egyptian_methods(a, b)egyptian_methods(a, b)
a, b
|
integers, a != 1, a < b and a, b relatively prime. |
For a rational number 0 < a/b < 1, generates Egyptian fractions
that is finds integers x1, x2, ..., xk such that
a/b = 1/x1 + 1/x2 + ... + 1/xk
using the following methods:
‘greedy’
Fibonacci-Sylvester
Golomb (same as with Farey sequences)
continued fractions (not yet implemented)
No return value, all solutions found will be printed to the console.
https://www.ics.uci.edu/~eppstein/numth/egypt/
egyptian_methods(8, 11) # 8/11 = 1/2 + 1/5 + 1/37 + 1/4070 (Fibonacci-Sylvester) # 8/11 = 1/2 + 1/6 + 1/21 + 1/77 (Golomb-Farey) # Other solutions # 8/11 = 1/2 + 1/8 + 1/11 + 1/88 # 8/11 = 1/2 + 1/12 + 1/22 + 1/121egyptian_methods(8, 11) # 8/11 = 1/2 + 1/5 + 1/37 + 1/4070 (Fibonacci-Sylvester) # 8/11 = 1/2 + 1/6 + 1/21 + 1/77 (Golomb-Farey) # Other solutions # 8/11 = 1/2 + 1/8 + 1/11 + 1/88 # 8/11 = 1/2 + 1/12 + 1/22 + 1/121
Euler's Phi function (aka Euler's ‘totient’ function).
eulersPhi(n)eulersPhi(n)
n |
Positive integer. |
The phi function is defined to be the number of positive integers
less than or equal to n that are coprime to n, i.e.
have no common factors other than 1.
Natural number, the number of coprime integers <= n.
Works well up to 10^9.
eulersPhi(9973) == 9973 - 1 # for prime numbers eulersPhi(3^10) == 3^9 * (3 - 1) # for prime powers eulersPhi(12*35) == eulersPhi(12) * eulersPhi(35) # TRUE if coprime ## Not run: x <- 1:100; y <- sapply(x, eulersPhi) plot(1:100, y, type="l", col="blue", xlab="n", ylab="phi(n)", main="Euler's totient function") points(1:100, y, col="blue", pch=20) grid() ## End(Not run)eulersPhi(9973) == 9973 - 1 # for prime numbers eulersPhi(3^10) == 3^9 * (3 - 1) # for prime powers eulersPhi(12*35) == eulersPhi(12) * eulersPhi(35) # TRUE if coprime ## Not run: x <- 1:100; y <- sapply(x, eulersPhi) plot(1:100, y, type="l", col="blue", xlab="n", ylab="phi(n)", main="Euler's totient function") points(1:100, y, col="blue", pch=20) grid() ## End(Not run)
The extended Euclidean algorithm computes the greatest common divisor and solves Bezout's identity.
extGCD(a, b)extGCD(a, b)
a, b
|
integer scalars |
The extended Euclidean algorithm not only computes the greatest common
divisor of and , but also two numbers and
such that .
This algorithm provides an easy approach to computing the modular inverse.
a numeric vector of length three, c(d, n, m), where d is the
greatest common divisor of a and b, and n and m
are integers such that d = n*a + m*b.
There is also a shorter, more elegant recursive version for the extended Euclidean algorithm. For R the procedure suggested by Blankinship appeared more appropriate.
Blankinship, W. A. “A New Version of the Euclidean Algorithm." Amer. Math. Monthly 70, 742-745, 1963.
extGCD(12, 10) extGCD(46368, 75025) # Fibonacci numbers are relatively prime to each otherextGCD(12, 10) extGCD(46368, 75025) # Fibonacci numbers are relatively prime to each other
Rational approximation of real numbers through Farey fractions.
ratFarey(x, n, upper = TRUE) farey_seq(n)ratFarey(x, n, upper = TRUE) farey_seq(n)
x |
real number. |
n |
integer, highest allowed denominator in a rational approximation. |
upper |
logical; shall the Farey fraction be grater than |
Rational approximation of real numbers through Farey fractions,
i.e. find for x the nearest fraction in the Farey series of
rational numbers with denominator not larger than n.
farey_seq(n) generates the full Farey sequence of rational
numbers with denominators not larger than n. Returns the
fractions as 'big rational' class in 'gmp'.
Returns a vector with two natural numbers, nominator and denominator.
farey_seq is very slow even for n > 40, due to the
handling of rational numbers as 'big rationals'.
Hardy, G. H., and E. M. Wright (1979). An Introduction to the Theory of Numbers. Fifth Edition, Oxford University Press, New York.
contFrac
ratFarey(pi, 100) # 22/7 0.0013 ratFarey(pi, 100, upper = FALSE) # 311/99 0.0002 ratFarey(-pi, 100) # -22/7 ratFarey(pi - 3, 70) # pi ~= 3 + (3/8)^2 ratFarey(pi, 1000) # 355/113 ratFarey(pi, 10000, upper = FALSE) # id. ratFarey(pi, 1e5, upper = FALSE) # 312689/99532 - pi ~= 3e-11 ratFarey(4/5, 5) # 4/5 ratFarey(4/5, 4) # 1/1 ratFarey(4/5, 4, upper = FALSE) # 3/4ratFarey(pi, 100) # 22/7 0.0013 ratFarey(pi, 100, upper = FALSE) # 311/99 0.0002 ratFarey(-pi, 100) # -22/7 ratFarey(pi - 3, 70) # pi ~= 3 + (3/8)^2 ratFarey(pi, 1000) # 355/113 ratFarey(pi, 10000, upper = FALSE) # id. ratFarey(pi, 1e5, upper = FALSE) # 312689/99532 - pi ~= 3e-11 ratFarey(4/5, 5) # 4/5 ratFarey(4/5, 4) # 1/1 ratFarey(4/5, 4, upper = FALSE) # 3/4
Generates single Fibonacci numbers or a Fibonacci sequence; or generates a Lucas series based on the Fibonacci series.
fibonacci(n, sequence = FALSE) lucas(n)fibonacci(n, sequence = FALSE) lucas(n)
n |
an integer. |
sequence |
logical; default: FALSE. |
Generates the n-th Fibonacci number, or the whole Fibonacci sequence
from the first to the n-th number; starts with (1, 1, 2, 3, ...).
Generates only single Lucas numbers. The Lucas series can be extenden to
the left and starts as (... -4, 3, -1, 2, 1, 3, 4, ...).
The recursive version is too slow for values n>=30. Therefore, an
iterative approach is used. For numbers n > 78 Fibonacci numbers
cannot be represented exactly in R as integers (>2^53-1).
A single integer, or a vector of integers.
fibonacci(0) # 0 fibonacci(2) # 1 fibonacci(2, sequence = TRUE) # 1 1 fibonacci(78) # 8944394323791464 < 9*10^15 lucas(0) # 2 lucas(2) # 3 lucas(76) # 7639424778862807 # Golden ratio F <- fibonacci(25, sequence = TRUE) # ... 46368 75025 f25 <- F[25]/F[24] # 1.618034 phi <- (sqrt(5) + 1)/2 abs(f25 - phi) # 2.080072e-10 # Fibonacci numbers w/o iteration fibo <- function(n) { phi <- (sqrt(5) + 1)/2 fib <- (phi^n - (1-phi)^n) / (2*phi - 1) round(fib) } fibo(30:33) # 832040 1346269 2178309 3524578 for (i in -8:8) cat(lucas(i), " ") # 47 -29 18 -11 7 -4 3 -1 2 1 3 4 7 11 18 29 47 # Lucas numbers w/o iteration luca <- function(n) { phi <- (sqrt(5) + 1)/2 luc <- phi^n + (1-phi)^n round(luc) } luca(0:10) # [1] 2 1 3 4 7 11 18 29 47 76 123 # Lucas primes # for (j in 0:76) { # l <- lucas(j) # if (isPrime(l)) cat(j, "\t", l, "\n") # } # 0 2 # 2 3 # ... # 71 688846502588399fibonacci(0) # 0 fibonacci(2) # 1 fibonacci(2, sequence = TRUE) # 1 1 fibonacci(78) # 8944394323791464 < 9*10^15 lucas(0) # 2 lucas(2) # 3 lucas(76) # 7639424778862807 # Golden ratio F <- fibonacci(25, sequence = TRUE) # ... 46368 75025 f25 <- F[25]/F[24] # 1.618034 phi <- (sqrt(5) + 1)/2 abs(f25 - phi) # 2.080072e-10 # Fibonacci numbers w/o iteration fibo <- function(n) { phi <- (sqrt(5) + 1)/2 fib <- (phi^n - (1-phi)^n) / (2*phi - 1) round(fib) } fibo(30:33) # 832040 1346269 2178309 3524578 for (i in -8:8) cat(lucas(i), " ") # 47 -29 18 -11 7 -4 3 -1 2 1 3 4 7 11 18 29 47 # Lucas numbers w/o iteration luca <- function(n) { phi <- (sqrt(5) + 1)/2 luc <- phi^n + (1-phi)^n round(luc) } luca(0:10) # [1] 2 1 3 4 7 11 18 29 47 76 123 # Lucas primes # for (j in 0:76) { # l <- lucas(j) # if (isPrime(l)) cat(j, "\t", l, "\n") # } # 0 2 # 2 3 # ... # 71 688846502588399
Greatest common divisor and least common multiple
GCD(n, m) LCM(n, m) mGCD(x) mLCM(x)GCD(n, m) LCM(n, m) mGCD(x) mLCM(x)
n, m
|
integer scalars. |
x |
a vector of integers. |
Computation based on the Euclidean algorithm without using the extended version.
mGCD (the multiple GCD) computes the greatest common divisor for
all numbers in the integer vector x together.
A numeric (integer) value.
The following relation is always true:
n * m = GCD(n, m) * LCM(n, m)
GCD(12, 10) GCD(46368, 75025) # Fibonacci numbers are relatively prime to each other LCM(12, 10) LCM(46368, 75025) # = 46368 * 75025 mGCD(c(2, 3, 5, 7) * 11) mGCD(c(2*3, 3*5, 5*7)) mLCM(c(2, 3, 5, 7) * 11) mLCM(c(2*3, 3*5, 5*7))GCD(12, 10) GCD(46368, 75025) # Fibonacci numbers are relatively prime to each other LCM(12, 10) LCM(46368, 75025) # = 46368 * 75025 mGCD(c(2, 3, 5, 7) * 11) mGCD(c(2*3, 3*5, 5*7)) mLCM(c(2, 3, 5, 7) * 11) mLCM(c(2*3, 3*5, 5*7))
Hermite normal form over integers (in column-reduced form).
hermiteNF(A)hermiteNF(A)
A |
integer matrix. |
An mxn-matrix of rank r with integer entries is said to be
in Hermite normal form if:
(i) the first r columns are nonzero, the other columns are all zero;
(ii) The first r diagonal elements are nonzero and d[i-1] divides d[i]
for i = 2,...,r .
(iii) All entries to the left of nonzero diagonal elements are non-negative
and strictly less than the corresponding diagonal entry.
The lower-triangular Hermite normal form of A is obtained by the following three types of column operations:
(i) exchange two columns
(ii) multiply a column by -1
(iii) Add an integral multiple of a column to another column
U is the unitary matrix such that AU = H, generated by these operations.
List with two matrices, the Hermite normal form H and the unitary
matrix U.
Another normal form often used in this context is the Smith normal form.
Cohen, H. (1993). A Course in Computational Algebraic Number Theory. Graduate Texts in Mathematics, Vol. 138, Springer-Verlag, Berlin, New York.
n <- 4; m <- 5 A = matrix(c( 9, 6, 0, -8, 0, -5, -8, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, -5, 0), n, m, byrow = TRUE) Hnf <- hermiteNF(A); Hnf # $H = 1 0 0 0 0 # 1 2 0 0 0 # 28 36 84 0 0 # -35 -45 -105 0 0 # $U = 11 14 32 0 0 # -7 -9 -20 0 0 # 0 0 0 1 0 # 7 9 21 0 0 # 0 0 0 0 1 r <- 3 # r = rank(H) H <- Hnf$H; U <- Hnf$U all(H == A %*% U) #=> TRUE ## Example: Compute integer solution of A x = b # H = A * U, thus H * U^-1 * x = b, or H * y = b b <- as.matrix(c(-11, -21, 16, -20)) y <- numeric(m) y[1] <- b[1] / H[1, 1] for (i in 2:r) y[i] <- (b[i] - sum(H[i, 1:(i-1)] * y[1:(i-1)])) / H[i, i] # special solution: xs <- U %*% y # 1 2 0 4 0 # and the general solution is xs + U * c(0, 0, 0, a, b), or # in other words the basis are the m-r vectors c(0,...,0, 1, ...). # If the special solution is not integer, there are no integer solutions.n <- 4; m <- 5 A = matrix(c( 9, 6, 0, -8, 0, -5, -8, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, -5, 0), n, m, byrow = TRUE) Hnf <- hermiteNF(A); Hnf # $H = 1 0 0 0 0 # 1 2 0 0 0 # 28 36 84 0 0 # -35 -45 -105 0 0 # $U = 11 14 32 0 0 # -7 -9 -20 0 0 # 0 0 0 1 0 # 7 9 21 0 0 # 0 0 0 0 1 r <- 3 # r = rank(H) H <- Hnf$H; U <- Hnf$U all(H == A %*% U) #=> TRUE ## Example: Compute integer solution of A x = b # H = A * U, thus H * U^-1 * x = b, or H * y = b b <- as.matrix(c(-11, -21, 16, -20)) y <- numeric(m) y[1] <- b[1] / H[1, 1] for (i in 2:r) y[i] <- (b[i] - sum(H[i, 1:(i-1)] * y[1:(i-1)])) / H[i, i] # special solution: xs <- U %*% y # 1 2 0 4 0 # and the general solution is xs + U * c(0, 0, 0, a, b), or # in other words the basis are the m-r vectors c(0,...,0, 1, ...). # If the special solution is not integer, there are no integer solutions.
Determine the integer n-th root of .
iNthroot(p, n)iNthroot(p, n)
p |
any positive number. |
n |
a natural number. |
Calculates the highest natural number below the n-th root of
p in a more integer based way than simply floor(p^{1/n}).
An integer.
iNthroot(0.5, 6) # 0 iNthroot(1, 6) # 1 iNthroot(5^6, 6) # 5 iNthroot(5^6-1, 6) # 4 ## Not run: # Define a function that tests whether isNthpower <- function(p, n) { q <- iNthroot(p, n) if (q^n == p) { return(TRUE) } else { return(FALSE) } } ## End(Not run)iNthroot(0.5, 6) # 0 iNthroot(1, 6) # 1 iNthroot(5^6, 6) # 5 iNthroot(5^6-1, 6) # 4 ## Not run: # Define a function that tests whether isNthpower <- function(p, n) { q <- iNthroot(p, n) if (q^n == p) { return(TRUE) } else { return(FALSE) } } ## End(Not run)
Determine whether p is the power of an integer.
isIntpower(p) isSquare(p) isSquarefree(p)isIntpower(p) isSquare(p) isSquarefree(p)
p |
any integer number. |
isIntpower(p) determines whether p is the power of an integer
and returns a tupel (n, m) such that p=n^m where m is
as small as possible. E.g., if p is prime it returns c(p,1).
isSquare(p) determines whether p is the square of an integer;
and isSquarefree(p) determines if p contains a square number
as a divisor.
A 2-vector of integers.
isIntpower(1) # 1 1 isIntpower(15) # 15 1 isIntpower(17) # 17 1 isIntpower(64) # 8 2 isIntpower(36) # 6 2 isIntpower(100) # 10 2 ## Not run: for (p in 5^7:7^5) { pp <- isIntpower(p) if (pp[2] != 1) cat(p, ":\t", pp, "\n") } ## End(Not run)isIntpower(1) # 1 1 isIntpower(15) # 15 1 isIntpower(17) # 17 1 isIntpower(64) # 8 2 isIntpower(36) # 6 2 isIntpower(100) # 10 2 ## Not run: for (p in 5^7:7^5) { pp <- isIntpower(p) if (pp[2] != 1) cat(p, ":\t", pp, "\n") } ## End(Not run)
Natural number type.
isNatural(n)isNatural(n)
n |
any numeric number. |
Returns TRUE for natural (or: whole) numbers between 1 and 2^53-1.
Boolean
IsNatural <- Vectorize(isNatural) IsNatural(c(-1, 0, 1, 5.1, 10, 2^53-1, 2^53, Inf)) # isNatural(NA) ?IsNatural <- Vectorize(isNatural) IsNatural(c(-1, 0, 1, 5.1, 10, 2^53-1, 2^53, Inf)) # isNatural(NA) ?
Vectorized version, returning for a vector or matrix of positive integers a vector of the same size containing 1 for the elements that are prime and 0 otherwise.
isPrime(x)isPrime(x)
x |
vector or matrix of nonnegative integers |
Given an array of positive integers returns an array of the same size of 0 and 1, where the i indicates a prime number in the same position.
array of elements 0, 1 with 1 indicating prime numbers
x <- matrix(1:10, nrow=10, ncol=10, byrow=TRUE) x * isPrime(x) # Find first prime number octett: octett <- c(0, 2, 6, 8, 30, 32, 36, 38) - 19 while (TRUE) { octett <- octett + 210 if (all(isPrime(octett))) { cat(octett, "\n", sep=" ") break } }x <- matrix(1:10, nrow=10, ncol=10, byrow=TRUE) x * isPrime(x) # Find first prime number octett: octett <- c(0, 2, 6, 8, 30, 32, 36, 38) - 19 while (TRUE) { octett <- octett + 210 if (all(isPrime(octett))) { cat(octett, "\n", sep=" ") break } }
Determine whether g generates the multiplicative group
modulo p.
isPrimroot(g, p)isPrimroot(g, p)
g |
integer greater 2 (and smaller than p). |
p |
prime number. |
Test is done by determining the order of g modulo p.
Returns TRUE or FALSE.
isPrimroot(2, 7) isPrimroot(2, 71) isPrimroot(7, 71)isPrimroot(2, 7) isPrimroot(2, 71) isPrimroot(7, 71)
Legendre and Jacobi Symbol for quadratic residues.
legendre_sym(a, p) jacobi_sym(a, n)legendre_sym(a, p) jacobi_sym(a, n)
a, n
|
integers. |
p |
prime number. |
The Legendre Symbol (a/p), where p must be a prime number,
denotes whether a is a quadratic residue modulo p or not.
The Jacobi symbol (a/p) is the product of (a/p) of all prime
factors p on n.
Returns 0, 1, or -1 if p divides a, a is a quadratic
residue, or if not.
Lsym <- Vectorize(legendre_sym, 'a') # all quadratic residues of p = 17 qr17 <- which(Lsym(1:16, 17) == 1) # 1 2 4 8 9 13 15 16 sort(unique((1:16)^2 %% 17)) # the same ## Not run: # how about large numbers? p <- 1198112137 # isPrime(p) TRUE x <- 4652356 a <- mod(x^2, p) # 520595831 legendre_sym(a, p) # 1 legendre_sym(a+1, p) # -1 ## End(Not run) jacobi_sym(11, 12) # -1Lsym <- Vectorize(legendre_sym, 'a') # all quadratic residues of p = 17 qr17 <- which(Lsym(1:16, 17) == 1) # 1 2 4 8 9 13 15 16 sort(unique((1:16)^2 %% 17)) # the same ## Not run: # how about large numbers? p <- 1198112137 # isPrime(p) TRUE x <- 4652356 a <- mod(x^2, p) # 520595831 legendre_sym(a, p) # 1 legendre_sym(a+1, p) # -1 ## End(Not run) jacobi_sym(11, 12) # -1
Determines whether is a Mersenne number, that is such that
is prime.
mersenne(p)mersenne(p)
p |
prime number, not very large. |
Applies the Lucas-Lehmer test on p. Because intermediate numbers will
soon get very large, uses ‘gmp’ from the beginning.
Returns TRUE or FALSE, indicating whether p is a Mersenne number or not.
https://mathworld.wolfram.com/Lucas-LehmerTest.html
mersenne(2) ## Not run: P <- Primes(32) M <- c() for (p in P) if (mersenne(p)) M <- c(M, p) # Next Mersenne numpers with primes are 521 and 607 (below 1200) M # 2 3 5 7 13 17 19 31 61 89 107 gmp::as.bigz(2)^M - 1 # 3 7 31 127 8191 131071 ... ## End(Not run)mersenne(2) ## Not run: P <- Primes(32) M <- c() for (p in P) if (mersenne(p)) M <- c(M, p) # Next Mersenne numpers with primes are 521 and 607 (below 1200) M # 2 3 5 7 13 17 19 31 61 89 107 gmp::as.bigz(2)^M - 1 # 3 7 31 127 8191 131071 ... ## End(Not run)
Probabilistic Miller-Rabin primality test.
miller_rabin(n)miller_rabin(n)
n |
natural number. |
The Miller-Rabin test is an efficient probabilistic primality test based
on strong pseudoprimes. This implementation uses the first seven prime
numbers (if necessary) as test cases. It is thus exact for all numbers
n < 341550071728321.
Returns TRUE or FALSE.
miller_rabin() will only work if package gmp has been loaded
by the user separately.
https://mathworld.wolfram.com/Rabin-MillerStrongPseudoprimeTest.html
miller_rabin(2) ## Not run: miller_rabin(4294967297) #=> FALSE miller_rabin(4294967311) #=> TRUE # Rabin-Miller 10 times faster than nextPrime() N <- n <- 2^32 + 1 system.time(while (!miller_rabin(n)) n <- n + 1) # 0.003 system.time(p <- nextPrime(N)) # 0.029 N <- c(2047, 1373653, 25326001, 3215031751, 2152302898747, 3474749660383, 341550071728321) for (n in N) { p <- nextPrime(n) T <- system.time(r <- miller_rabin(p)) cat(n, p, r, T[3], "\n")} ## End(Not run)miller_rabin(2) ## Not run: miller_rabin(4294967297) #=> FALSE miller_rabin(4294967311) #=> TRUE # Rabin-Miller 10 times faster than nextPrime() N <- n <- 2^32 + 1 system.time(while (!miller_rabin(n)) n <- n + 1) # 0.003 system.time(p <- nextPrime(N)) # 0.029 N <- c(2047, 1373653, 25326001, 3215031751, 2152302898747, 3474749660383, 341550071728321) for (n in N) { p <- nextPrime(n) T <- system.time(r <- miller_rabin(p)) cat(n, p, r, T[3], "\n")} ## End(Not run)
Modulo operator.
mod(n, m) modq(a, b, k)mod(n, m) modq(a, b, k)
n |
numeric vector (preferably of integers) |
m |
integer vector (positive, zero, or negative) |
a, b
|
whole numbers (scalars) |
k |
integer greater than 1 |
mod(n, m) is the modulo operator and returns n mod m.
mod(n, 0) is n, and the result always has the same sign
as m.
modq(a, b, k) is the modulo operator for rational numbers and
returns a/b mod k. b and k must be coprime,
otherwise NA is returned.
a numeric (integer) value or vector/matrix, resp. an integer number
The following relation is fulfilled (for m != 0):
mod(n, m) = n - m * floor(n/m)
mod(c(-5:5), 5) mod(c(-5:5), -5) mod(0, 1) #=> 0 mod(1, 0) #=> 1 modq(5, 66, 5) # 0 (Bernoulli 10) modq(5, 66, 7) # 4 modq(5, 66, 13) # 5 modq(5, 66, 25) # 5 modq(5, 66, 35) # 25 modq(-1, 30, 7) # 3 (Bernoulli 8) modq( 1, -30, 7) # 3 # Warning messages: # modq(5, 66, 77) : Arguments 'b' and 'm' must be coprime. # Error messages # modq(5, 66, 1) : Argument 'm' mustbe a natural number > 1. # modq(5, 66, 1.5) : All arguments of 'modq' must be integers. # modq(5, 66, c(5, 7)) : Function 'modq' is *not* vectorized.mod(c(-5:5), 5) mod(c(-5:5), -5) mod(0, 1) #=> 0 mod(1, 0) #=> 1 modq(5, 66, 5) # 0 (Bernoulli 10) modq(5, 66, 7) # 4 modq(5, 66, 13) # 5 modq(5, 66, 25) # 5 modq(5, 66, 35) # 25 modq(-1, 30, 7) # 3 (Bernoulli 8) modq( 1, -30, 7) # 3 # Warning messages: # modq(5, 66, 77) : Arguments 'b' and 'm' must be coprime. # Error messages # modq(5, 66, 1) : Argument 'm' mustbe a natural number > 1. # modq(5, 66, 1.5) : All arguments of 'modq' must be integers. # modq(5, 66, c(5, 7)) : Function 'modq' is *not* vectorized.
Computes the modular inverse of n modulo m.
modinv(n, m) modsqrt(a, p)modinv(n, m) modsqrt(a, p)
n, m
|
integer scalars. |
a, p
|
integer modulo p, p a prime. |
The modular inverse of n modulo m is the unique natural
number 0 < n0 < m such that n * n0 = 1 mod m. It is a
simple application of the extended GCD algorithm.
The modular square root of a modulo a prime p is a number
x such that x^2 = a mod p. If x is a solution, then
p-x is also a solution module p. The function will always
return the smaller value.
modsqrt implements the Tonelli-Shanks algorithm which also works
for square roots modulo prime powers. The general case is NP-hard.
A natural number smaller m, if n and m are coprime,
else NA. modsqrt will return 0 if there is no solution.
modinv(5, 1001) #=> 801, as 5*801 = 4005 = 1 mod 1001 Modinv <- Vectorize(modinv, "n") ((1:10)*Modinv(1:10, 11)) %% 11 #=> 1 1 1 1 1 1 1 1 1 1 modsqrt( 8, 23) # 10 because 10^2 = 100 = 8 mod 23 modsqrt(10, 17) # 0 because 10 is not a quadratic residue mod 17modinv(5, 1001) #=> 801, as 5*801 = 4005 = 1 mod 1001 Modinv <- Vectorize(modinv, "n") ((1:10)*Modinv(1:10, 11)) %% 11 #=> 1 1 1 1 1 1 1 1 1 1 modsqrt( 8, 23) # 10 because 10^2 = 100 = 8 mod 23 modsqrt(10, 17) # 0 because 10 is not a quadratic residue mod 17
Solves the modular equation a x = b mod n.
modlin(a, b, n)modlin(a, b, n)
a, b, n
|
integer scalars |
Solves the modular equation a x = b mod n. This eqation is solvable
if and only if gcd(a,n)|b. The function uses the extended greatest
common divisor approach.
Returns a vector of integer solutions.
modlin(14, 30, 100) # 95 45 modlin(3, 4, 5) # 3 modlin(3, 5, 6) # [] modlin(3, 6, 9) # 2 5 8modlin(14, 30, 100) # 95 45 modlin(3, 4, 5) # 3 modlin(3, 5, 6) # [] modlin(3, 6, 9) # 2 5 8
Realizes the modular (or discrete) logarithm modulo a prime number
, that is determines the unique exponent such that
, g a primitive root.
modlog(g, x, p)modlog(g, x, p)
g |
a primitive root mod p. |
x |
an integer. |
p |
prime number. |
The method is in principle a complete search, cut short by "Shank's
trick", the giantstep-babystep approach, see Forster
(1996, pp. 65f). g has to be a primitive root modulo p,
otherwise exponentiation is not bijective.
Returns an integer.
Forster, O. (1996). Algorithmische Zahlentheorie. Friedr. Vieweg u. Sohn Verlagsgesellschaft mbH, Wiesbaden.
modlog(11, 998, 1009) # 505 , i.e., 11^505 = 998 mod 1009modlog(11, 998, 1009) # 505 , i.e., 11^505 = 998 mod 1009
Calculates powers and orders modulo m.
modpower(n, k, m) modorder(n, m)modpower(n, k, m) modorder(n, m)
n, k, m
|
Natural numbers, |
modpower calculates n to the power of k modulo
m.
Uses modular exponentiation, as described in the Wikipedia article.
modorder calculates the order of n in the multiplicative
group module m. n and m must be coprime.
Uses brute force, trick to use binary expansion and square is not more
efficient in an R implementation.
Natural number.
This function is not vectorized.
modpower(2, 100, 7) #=> 2 modpower(3, 100, 7) #=> 4 modorder(7, 17) #=> 16, i.e. 7 is a primitive root mod 17 ## Gauss' table of primitive roots modulo prime numbers < 100 proots <- c(2, 2, 3, 2, 2, 6, 5, 10, 10, 10, 2, 2, 10, 17, 5, 5, 6, 28, 10, 10, 26, 10, 10, 5, 12, 62, 5, 29, 11, 50, 30, 10) P <- Primes(100) for (i in seq(along=P)) { cat(P[i], "\t", modorder(proots[i], P[i]), proots[i], "\t", "\n") } ## Not run: ## Lehmann's primality test lehmann_test <- function(n, ntry = 25) { if (!is.numeric(n) || ceiling(n) != floor(n) || n < 0) stop("Argument 'n' must be a natural number") if (n >= 9e7) stop("Argument 'n' should be smaller than 9e7.") if (n < 2) return(FALSE) else if (n == 2) return(TRUE) else if (n > 2 && n %% 2 == 0) return(FALSE) k <- floor(ntry) if (k < 1) k <- 1 if (k > n-2) a <- 2:(n-1) else a <- sample(2:(n-1), k, replace = FALSE) for (i in 1:length(a)) { m <- modpower(a[i], (n-1)/2, n) if (m != 1 && m != n-1) return(FALSE) } return(TRUE) } ## Examples for (i in seq(1001, 1011, by = 2)) if (lehmann_test(i)) cat(i, "\n") # 1009 system.time(lehmann_test(27644437, 50)) # TRUE # user system elapsed # 0.086 0.151 0.235 ## End(Not run)modpower(2, 100, 7) #=> 2 modpower(3, 100, 7) #=> 4 modorder(7, 17) #=> 16, i.e. 7 is a primitive root mod 17 ## Gauss' table of primitive roots modulo prime numbers < 100 proots <- c(2, 2, 3, 2, 2, 6, 5, 10, 10, 10, 2, 2, 10, 17, 5, 5, 6, 28, 10, 10, 26, 10, 10, 5, 12, 62, 5, 29, 11, 50, 30, 10) P <- Primes(100) for (i in seq(along=P)) { cat(P[i], "\t", modorder(proots[i], P[i]), proots[i], "\t", "\n") } ## Not run: ## Lehmann's primality test lehmann_test <- function(n, ntry = 25) { if (!is.numeric(n) || ceiling(n) != floor(n) || n < 0) stop("Argument 'n' must be a natural number") if (n >= 9e7) stop("Argument 'n' should be smaller than 9e7.") if (n < 2) return(FALSE) else if (n == 2) return(TRUE) else if (n > 2 && n %% 2 == 0) return(FALSE) k <- floor(ntry) if (k < 1) k <- 1 if (k > n-2) a <- 2:(n-1) else a <- sample(2:(n-1), k, replace = FALSE) for (i in 1:length(a)) { m <- modpower(a[i], (n-1)/2, n) if (m != 1 && m != n-1) return(FALSE) } return(TRUE) } ## Examples for (i in seq(1001, 1011, by = 2)) if (lehmann_test(i)) cat(i, "\n") # 1009 system.time(lehmann_test(27644437, 50)) # TRUE # user system elapsed # 0.086 0.151 0.235 ## End(Not run)
The classical Moebius and Mertens functions in number theory.
moebius(n) mertens(n)moebius(n) mertens(n)
n |
Positive integer. |
moebius(n) is +1 if n is a square-free positive integer
with an even number of prime factors, or +1 if there are an odd
of prime factors. It is 0 if n is not square-free.
mertens(n) is the aggregating summary function, that sums up all
values of moebius from 1 to n.
For moebius, 0, 1 or -1, depending on the prime
decomposition of n.
For mertens the values will very slowly grow.
Works well up to 10^9, but will become very slow for the Mertens
function.
sapply(1:16, moebius) sapply(1:16, mertens) ## Not run: x <- 1:50; y <- sapply(x, moebius) plot(c(1, 50), c(-3, 3), type="n") grid() points(1:50, y, pch=18, col="blue") x <- 1:100; y <- sapply(x, mertens) plot(c(1, 100), c(-5, 3), type="n") grid() lines(1:100, y, col="red", type="s") ## End(Not run)sapply(1:16, moebius) sapply(1:16, mertens) ## Not run: x <- 1:50; y <- sapply(x, moebius) plot(c(1, 50), c(-3, 3), type="n") grid() points(1:50, y, pch=18, col="blue") x <- 1:100; y <- sapply(x, mertens) plot(c(1, 100), c(-5, 3), type="n") grid() lines(1:100, y, col="red", type="s") ## End(Not run)
Necklace and bracelet problems in combinatorics.
necklace(k, n) bracelet(k, n)necklace(k, n) bracelet(k, n)
k |
The size of the set or alphabet to choose from. |
n |
the length of the necklace or bracelet. |
A necklace is a closed string of length n over a set of size
k (numbers, characters, clors, etc.), where all rotations are
taken as equivalent. A bracelet is a necklace where strings may also
be equivalent under reflections.
Polya's enumeration theorem can be utilized to enumerate all necklaces
or bracelets. The final calculation involves Euler's Phi or totient
function, in this package implemented as eulersPhi.
Returns the number of necklaces resp. bracelets.
https://en.wikipedia.org/wiki/Necklace_(combinatorics)
necklace(2, 5) necklace(3, 6) bracelet(2, 5) bracelet(3, 6)necklace(2, 5) necklace(3, 6) bracelet(2, 5) bracelet(3, 6)
Find the next prime above n.
nextPrime(n)nextPrime(n)
n |
natural number. |
nextPrime finds the next prime number greater than n, while
previousPrime finds the next prime number below n.
In general the next prime will occur in the interval [n+1,n+log(n)].
In double precision arithmetic integers are represented exactly only up to 2^53 - 1, therefore this is the maximal allowed value.
Integer.
p <- nextPrime(1e+6) # 1000003 isPrime(p) # TRUEp <- nextPrime(1e+6) # 1000003 isPrime(p) # TRUE
Number of prime factors resp. sum of all exponents of prime factors in the prime decomposition.
omega(n) Omega(n)omega(n) Omega(n)
n |
Positive integer. |
'omega(n)' returns the number of prime factors of 'n' while 'Omega(n)' returns the sum of their exponents in the prime decomposition. 'omega' and 'Omega' are identical if there are no quadratic factors.
Remark: (-1)^Omega(n) is the Liouville function.
Natural number.
Works well up to 10^9.
omega(2*3*5*7*11*13*17*19) #=> 8 Omega(2 * 3^2 * 5^3 * 7^4) #=> 10omega(2*3*5*7*11*13*17*19) #=> 8 Omega(2 * 3^2 * 5^3 * 7^4) #=> 10
Calculates the order of a prime number p in n!, i.e. the
highest exponent e such that p^e|n!.
ordpn(p, n)ordpn(p, n)
p |
prime number. |
n |
natural number. |
Applies the well-known formula adding terms floor(n/p^k).
Returns the exponent e.
ordpn(2, 100) #=> 97 ordpn(7, 100) #=> 16 ordpn(101, 100) #=> 0 ordpn(997, 1000) #=> 1ordpn(2, 100) #=> 97 ordpn(7, 100) #=> 16 ordpn(101, 100) #=> 0 ordpn(997, 1000) #=> 1
Generates the Pascal triangle in rectangular form.
pascal_triangle(n)pascal_triangle(n)
n |
integer number |
Pascal numbers will be generated with the usual recursion formula and stored in a rectangular scheme.
For n>50 integer overflow would happen, so use the arbitrary
precision version gmp::chooseZ(n, 0:n) instead for calculating
binomial numbers.
Returns the Pascal triangle as an (n+1)x(n+1) rectangle with zeros filled in.
See Wolfram MathWorld or the Wikipedia.
n <- 5; P <- pascal_triangle(n) for (i in 1:(n+1)) { cat(P[i, 1:i], '\n') } ## 1 ## 1 1 ## 1 2 1 ## 1 3 3 1 ## 1 4 6 4 1 ## 1 5 10 10 5 1 ## Not run: P <- pascal_triangle(50) max(P[51, ]) ## [1] 126410606437752 ## End(Not run)n <- 5; P <- pascal_triangle(n) for (i in 1:(n+1)) { cat(P[i, 1:i], '\n') } ## 1 ## 1 1 ## 1 2 1 ## 1 3 3 1 ## 1 4 6 4 1 ## 1 5 10 10 5 1 ## Not run: P <- pascal_triangle(50) max(P[51, ]) ## [1] 126410606437752 ## End(Not run)
Generates a periodic continued fraction.
periodicCF(d)periodicCF(d)
d |
positive integer that is not a square number |
The function computes the periodic continued fraction of the square root of an integer that itself shall not be a square (because otherwise the integer square root will be returned). Note that the continued fraction of an irrational quadratic number is always a periodic continued fraction.
The first term is the biggest integer below sqrt(d) and the rest is
the period of the continued fraction. The period is always exact, there is
no floating point inaccuracy involved (though integer overflow may happen
for very long fractions).
The underlying algorithm is sometimes called "The Fundamental Algorithm for Quadratic Numbers". The function will be utilized especially when solving Pell's equation.
Returns a list with components
cf |
the continued fraction with integer part and first period. |
plen |
the length of the period. |
Integer overflow may happen for very long continued fractions.
Hans Werner Borchers
Mak Trifkovic. Algebraic Theory of Quadratic Numbers. Springer Verlag, Universitext, New York 2013.
periodicCF(2) # sqrt(2) = [1; 2,2,2,...] = [1; (2)] periodicCF(1003) ## $cf ## [1] 31 1 2 31 2 1 62 ## $plen ## [1] 6periodicCF(2) # sqrt(2) = [1; 2,2,2,...] = [1; (2)] periodicCF(1003) ## $cf ## [1] 31 1 2 31 2 1 62 ## $plen ## [1] 6
Find the next prime below n.
previousPrime(n)previousPrime(n)
n |
natural number. |
previousPrime finds the next prime number smaller than n,
while nextPrime finds the next prime number below n.
In general the previousn prime will occur in the interval
[n-1,n-log(n)].
In double precision arithmetic integers are represented exactly only up to 2^53 - 1, therefore this is the maximal allowed value.
Integer.
p <- previousPrime(1e+6) # 999983 isPrime(p) # TRUEp <- previousPrime(1e+6) # 999983 isPrime(p) # TRUE
primeFactors computes a vector containing the prime factors of
n. radical returns the product of those unique prime
factors.
primeFactors(n) radical(n)primeFactors(n) radical(n)
n |
nonnegative integer |
Computes the prime factors of n in ascending order,
each one as often as its multiplicity requires, such that
n == prod(primeFactors(n)).
## radical() is used in the abc-conjecture:
# abc-triple: 1 <= a < b, a, b coprime, c = a + b
# for every e > 0 there are only finitely many abc-triples with
# c > radical(a*b*c)^(1+e)
Vector containing the prime factors of n, resp.
the product of unique prime factors.
divisors, gmp::factorize
primeFactors(1002001) # 7 7 11 11 13 13 primeFactors(65537) # is prime # Euler's calculation primeFactors(2^32 + 1) # 641 6700417 radical(1002001) # 1001 ## Not run: for (i in 1:99) { for (j in (i+1):100) { if (coprime(i, j)) { k = i + j r = radical(i*j*k) q = log(k) / log(r) # 'quality' of the triple if (q > 1) cat(q, ":\t", i, ",", j, ",", k, "\n") } } } ## End(Not run)primeFactors(1002001) # 7 7 11 11 13 13 primeFactors(65537) # is prime # Euler's calculation primeFactors(2^32 + 1) # 641 6700417 radical(1002001) # 1001 ## Not run: for (i in 1:99) { for (j in (i+1):100) { if (coprime(i, j)) { k = i + j r = radical(i*j*k) q = log(k) / log(r) # 'quality' of the triple if (q > 1) cat(q, ":\t", i, ",", j, ",", k, "\n") } } } ## End(Not run)
Eratosthenes resp. Atkin sieve methods to generate a list of prime numbers
less or equal n, resp. between n1 and n2.
Primes(n1, n2 = NULL) atkin_sieve(n)Primes(n1, n2 = NULL) atkin_sieve(n)
n, n1, n2
|
natural numbers with |
The list of prime numbers up to n is generated using the "sieve of
Eratosthenes". This approach is reasonably fast, but may require a lot of
main memory when n is large.
Primes computes first all primes up to sqrt(n2) and then
applies a refined sieve on the numbers from n1 to n2, thereby
drastically reducing the need for storing long arrays of numbers.
The sieve of Atkins is a modified version of the ancient prime number sieve of Eratosthenes. It applies a modulo-sixty arithmetic and requires less memory, but in R is not faster because of a double for-loop.
In double precision arithmetic integers are represented exactly only up to 2^53 - 1, therefore this is the maximal allowed value.
vector of integers representing prime numbers
A. Atkin and D. Bernstein (2004), Prime sieves using quadratic forms. Mathematics of Computation, Vol. 73, pp. 1023-1030.
isPrime, gmp::factorize, pracma::expint1
Primes(1000) Primes(1949, 2019) atkin_sieve(1000) ## Not run: ## Appendix: Logarithmic Integrals and Prime Numbers (C.F.Gauss, 1846) library('gsl') # 'European' form of the logarithmic integral Li <- function(x) expint_Ei(log(x)) - expint_Ei(log(2)) # No. of primes and logarithmic integral for 10^i, i=1..12 i <- 1:12; N <- 10^i # piN <- numeric(12) # for (i in 1:12) piN[i] <- length(primes(10^i)) piN <- c(4, 25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534, 455052511, 4118054813, 37607912018) cbind(i, piN, round(Li(N)), round((Li(N)-piN)/piN, 6)) # i pi(10^i) Li(10^i) rel.err # -------------------------------------- # 1 4 5 0.280109 # 2 25 29 0.163239 # 3 168 177 0.050979 # 4 1229 1245 0.013094 # 5 9592 9629 0.003833 # 6 78498 78627 0.001637 # 7 664579 664917 0.000509 # 8 5761455 5762208 0.000131 # 9 50847534 50849234 0.000033 # 10 455052511 455055614 0.000007 # 11 4118054813 4118066400 0.000003 # 12 37607912018 37607950280 0.000001 # -------------------------------------- ## End(Not run)Primes(1000) Primes(1949, 2019) atkin_sieve(1000) ## Not run: ## Appendix: Logarithmic Integrals and Prime Numbers (C.F.Gauss, 1846) library('gsl') # 'European' form of the logarithmic integral Li <- function(x) expint_Ei(log(x)) - expint_Ei(log(2)) # No. of primes and logarithmic integral for 10^i, i=1..12 i <- 1:12; N <- 10^i # piN <- numeric(12) # for (i in 1:12) piN[i] <- length(primes(10^i)) piN <- c(4, 25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534, 455052511, 4118054813, 37607912018) cbind(i, piN, round(Li(N)), round((Li(N)-piN)/piN, 6)) # i pi(10^i) Li(10^i) rel.err # -------------------------------------- # 1 4 5 0.280109 # 2 25 29 0.163239 # 3 168 177 0.050979 # 4 1229 1245 0.013094 # 5 9592 9629 0.003833 # 6 78498 78627 0.001637 # 7 664579 664917 0.000509 # 8 5761455 5762208 0.000131 # 9 50847534 50849234 0.000033 # 10 455052511 455055614 0.000007 # 11 4118054813 4118066400 0.000003 # 12 37607912018 37607950280 0.000001 # -------------------------------------- ## End(Not run)
Find the smallest primitive root modulo m, or find all primitive roots.
primroot(m, all = FALSE)primroot(m, all = FALSE)
m |
A prime integer. |
all |
boolean; shall all primitive roots module p be found. |
For every prime number there exists a natural number that
generates the field , i.e. are
all different.
The computation here is all brute force. As most primitive roots are relatively small, so it is still reasonable fast.
One trick is to factorize and test only for those prime factors.
In R this is not more efficient as factorization also takes some time.
A natural number if m is prime, else NA.
This function is not vectorized.
Arndt, J. (2010). Matters Computational: Ideas, Algorithms, Source Code. Springer-Verlag, Berlin Heidelberg Dordrecht.
P <- Primes(100) R <- c() for (p in P) { R <- c(R, primroot(p)) } cbind(P, R) # 7 is the biggest prime root here (for p=71)P <- Primes(100) R <- c() for (p in P) { R <- c(R, primroot(p)) } cbind(P, R) # 7 is the biggest prime root here (for p=71)
Generates all primitive Pythagorean triples of integers
such that , where are coprime (have no
common divisor) and .
pythagorean_triples(c1, c2)pythagorean_triples(c1, c2)
c1, c2
|
lower and upper limit of the hypothenuses |
If is a primitive Pythagorean triple, there are integers
with such that
with and being odd.
Returns a matrix, one row for each Pythagorean triple, of the form
(m n a b c).
https://mathworld.wolfram.com/PythagoreanTriple.html
pythagorean_triples(100, 200) ## [,1] [,2] [,3] [,4] [,5] ## [1,] 10 1 99 20 101 ## [2,] 10 3 91 60 109 ## [3,] 8 7 15 112 113 ## [4,] 11 2 117 44 125 ## [5,] 11 4 105 88 137 ## [6,] 9 8 17 144 145 ## [7,] 12 1 143 24 145 ## [8,] 10 7 51 140 149 ## [9,] 11 6 85 132 157 ## [10,] 12 5 119 120 169 ## [11,] 13 2 165 52 173 ## [12,] 10 9 19 180 181 ## [13,] 11 8 57 176 185 ## [14,] 13 4 153 104 185 ## [15,] 12 7 95 168 193 ## [16,] 14 1 195 28 197pythagorean_triples(100, 200) ## [,1] [,2] [,3] [,4] [,5] ## [1,] 10 1 99 20 101 ## [2,] 10 3 91 60 109 ## [3,] 8 7 15 112 113 ## [4,] 11 2 117 44 125 ## [5,] 11 4 105 88 137 ## [6,] 9 8 17 144 145 ## [7,] 12 1 143 24 145 ## [8,] 10 7 51 140 149 ## [9,] 11 6 85 132 157 ## [10,] 12 5 119 120 169 ## [11,] 13 2 165 52 173 ## [12,] 10 9 19 180 181 ## [13,] 11 8 57 176 185 ## [14,] 13 4 153 104 185 ## [15,] 12 7 95 168 193 ## [16,] 14 1 195 28 197
List all quadratic residues of an integer.
quadratic_residues(n)quadratic_residues(n)
n |
integer. |
Squares all numbers between 0 and n/2 and generate a unique list of
all these numbers modulo n.
Vector of integers.
quadratic_residues(17)quadratic_residues(17)
Integer remainder function.
rem(n, m)rem(n, m)
n |
numeric vector (preferably of integers) |
m |
must be a scalar integer (positive, zero, or negative) |
rem(n, m) is the same modulo operator and returns .
mod(n, 0) is NaN, and the result always has the same sign
as n (for n != m and m != 0).
a numeric (integer) value or vector/matrix
rem(c(-5:5), 5) rem(c(-5:5), -5) rem(0, 1) #=> 0 rem(1, 1) #=> 0 (always for n == m) rem(1, 0) # NA (should be NaN) rem(0, 0) #=> NaNrem(c(-5:5), 5) rem(c(-5:5), -5) rem(0, 1) #=> 0 rem(1, 1) #=> 0 (always for n == m) rem(1, 0) # NA (should be NaN) rem(0, 0) #=> NaN
Sum of powers of all divisors of a natural number.
Sigma(n, k = 1, proper = FALSE) tau(n)Sigma(n, k = 1, proper = FALSE) tau(n)
n |
Positive integer. |
k |
Numeric scalar, the exponent to be used. |
proper |
Logical; if |
Total sum of all integer divisors of n to the power of k,
including 1 and n.
For k=0 this is the number of divisors, for k=1
it is the sum of all divisors of n.
tau is Ramanujan's tau function, here computed using
Sigma(., 5) and Sigma(., 11).
A number is called refactorable, if tau(n) divides n,
for example n=12 or n=18.
Natural number, the number or sum of all divisors.
Works well up to 10^9.
https://en.wikipedia.org/wiki/Divisor_function
https://en.wikipedia.org/wiki/Ramanujan_tau_function
sapply(1:16, Sigma, k = 0) sapply(1:16, Sigma, k = 1) sapply(1:16, Sigma, proper = TRUE)sapply(1:16, Sigma, k = 0) sapply(1:16, Sigma, k = 1) sapply(1:16, Sigma, proper = TRUE)
Find the basic, that is minimal, solution for Pell's equation, applying the technique of (periodic) continued fractions.
solvePellsEq(d)solvePellsEq(d)
d |
non-square integer greater 1. |
Solving Pell's equation means to find integer solutions (x,y)
for the Diophantine equation
for a
non-square integer. These solutions are important in number theory and
for the theory of quadratic number fields.
The procedure goes as follows: First find the periodic continued
fraction for , then determine the convergents of this
continued fraction. The last pair of convergents will provide the
solution for Pell's equation.
The solution found is the minimal or fundamental solution. All other solutions can be derived from this one – but the numbers grow up very rapidly.
Returns a list with components
x, y
|
solution (x,y) of Pell's equation. |
plen |
length of the period. |
doubled |
logical: was the period doubled? |
msg |
message either "Success" or "Integer overflow". |
If 'doubled' was TRUE, there exists also a solution for the negative Pell equation
Integer overflow may happen for the convergents, but very rarely.
More often, the terms x^2 or y^2 can overflow the
maximally representable integer 2^53-1 and checking Pell's
equation may end with a value differing from 1, though in
reality the solution is correct.
Hans Werner Borchers
H.W. Lenstra Jr. Solving the Pell Equation. Notices of the AMS, Vol. 49, No. 2, February 2002.
See the "List of fundamental solutions of Pell's equations" in the Wikipedia entry for "Pell's Equation".
s = solvePellsEq(1003) # $x = 9026, $y = 285 9026^2 - 1003*285^2 == 1 # TRUEs = solvePellsEq(1003) # $x = 9026, $y = 285 9026^2 - 1003*285^2 == 1 # TRUE
The function generates the Stern-Brocot sequence up to
length n.
stern_brocot_seq(n)stern_brocot_seq(n)
n |
integer; length of the sequence. |
The Stern-Brocot sequence is a sequence S of natural
numbers beginning with
1, 1, 2, 1, 3, 2, 3, 1, 4, 3, 5, 2, 5, 3, 4, 1, ...
defined with S[1] = S[2] = 1 and the following rules:
S[k] = S[k/2] if k is even S[k] = S[(k-1)/2] + S[(k+1)/2] if k is not even
The Stern-Brocot has the remarkable properties that
(1) Consecutive values in this sequence are coprime;
(2) the list of rationals S[k+1]/S[k] (all in reduced
form) covers all positive rational numbers once and once only.
Returns a sequence of length n of natural numbers.
N. Calkin and H.S. Wilf. Recounting the rationals. The American Mathematical Monthly, Vol. 7(4), 2000.
Graham, Knuth, and Patashnik. Concrete Mathematics - A Foundation for Computer Science. Addison-Wesley, 1989.
( S <- stern_brocot_seq(92) ) # 1, 1, 2, 1, 3, 2, 3, 1, 4, 3, 5, 2, 5, 3, 4, 1, 5, 4, 7, # 3, 8, 5, 7, 2, 7, 5, 8, 3, 7, 4, 5, 1, 6, 5, 9, 4, 11, 7, 10, # 3, 11, 8, 13, 5, 12, 7, 9, 2, 9, 7, 12, 5, 13, 8, 11, 3, 10, 7, 11, # 4, 9, 5, 6, 1, 7, 6, 11, 5, 14, 9, 13, 4, 15, 11, 18, 7, 17, 10, 13, # 3, 14, 11, 19, 8, 21, 13, 18, 5, 17, 12, 19, 7, ... table(S) ## S ## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 17 18 19 21 ## 7 5 9 7 12 3 11 5 5 3 7 3 5 2 1 2 2 2 1 which(S == 1) # 1 2 4 8 16 32 64 ## Not run: # Find the rational number p/q in S # note that 1/2^n appears in position S[c(2^(n-1), 2^(n-1)+1)] occurs <- function(p, q, s){ # Find i such that (p, q) = s[i, i+1] inds <- seq.int(length = length(s)-1) inds <- inds[p == s[inds]] inds[q == s[inds + 1]] } p = 3; q = 7 # 3/7 occurs(p, q, S) # S[28, 29] '%//%' <- function(p, q) gmp::as.bigq(p, q) n <- length(S) S[1:(n-1)] %//% S[2:n] ## Big Rational ('bigq') object of length 91: ## [1] 1 1/2 2 1/3 3/2 2/3 3 1/4 4/3 3/5 ## [11] 5/2 2/5 5/3 3/4 4 1/5 5/4 4/7 7/3 3/8 ... as.double(S[1:(n-1)] %//% S[2:n]) ## [1] 1.000000 0.500000 2.000000 0.333333 1.500000 0.666667 3.000000 ## [8] 0.250000 1.333333 0.600000 2.500000 0.400000 1.666667 0.750000 ... ## End(Not run)( S <- stern_brocot_seq(92) ) # 1, 1, 2, 1, 3, 2, 3, 1, 4, 3, 5, 2, 5, 3, 4, 1, 5, 4, 7, # 3, 8, 5, 7, 2, 7, 5, 8, 3, 7, 4, 5, 1, 6, 5, 9, 4, 11, 7, 10, # 3, 11, 8, 13, 5, 12, 7, 9, 2, 9, 7, 12, 5, 13, 8, 11, 3, 10, 7, 11, # 4, 9, 5, 6, 1, 7, 6, 11, 5, 14, 9, 13, 4, 15, 11, 18, 7, 17, 10, 13, # 3, 14, 11, 19, 8, 21, 13, 18, 5, 17, 12, 19, 7, ... table(S) ## S ## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 17 18 19 21 ## 7 5 9 7 12 3 11 5 5 3 7 3 5 2 1 2 2 2 1 which(S == 1) # 1 2 4 8 16 32 64 ## Not run: # Find the rational number p/q in S # note that 1/2^n appears in position S[c(2^(n-1), 2^(n-1)+1)] occurs <- function(p, q, s){ # Find i such that (p, q) = s[i, i+1] inds <- seq.int(length = length(s)-1) inds <- inds[p == s[inds]] inds[q == s[inds + 1]] } p = 3; q = 7 # 3/7 occurs(p, q, S) # S[28, 29] '%//%' <- function(p, q) gmp::as.bigq(p, q) n <- length(S) S[1:(n-1)] %//% S[2:n] ## Big Rational ('bigq') object of length 91: ## [1] 1 1/2 2 1/3 3/2 2/3 3 1/4 4/3 3/5 ## [11] 5/2 2/5 5/3 3/4 4 1/5 5/4 4/7 7/3 3/8 ... as.double(S[1:(n-1)] %//% S[2:n]) ## [1] 1.000000 0.500000 2.000000 0.333333 1.500000 0.666667 3.000000 ## [8] 0.250000 1.333333 0.600000 2.500000 0.400000 1.666667 0.750000 ... ## End(Not run)
Generate a list of twin primes between n1 and n2.
twinPrimes(n1, n2)twinPrimes(n1, n2)
n1, n2
|
natural numbers with |
twinPrimes uses Primes and uses diff to find all
twin primes in the given interval.
In double precision arithmetic integers are represented exactly only up to 2^53 - 1, therefore this is the maximal allowed value.
Returnes a nx2-matrix, where nis the number of twin primes
found, and each twin tuple fills one row.
twinPrimes(1e6+1, 1e6+1001)twinPrimes(1e6+1, 1e6+1001)
Generates the Zeckendorf representation of an integer as a sum of Fibonacci numbers.
zeck(n)zeck(n)
n |
integer. |
According to Zeckendorfs theorem from 1972, each integer can be uniquely represented as a sum of Fibonacci numbers such that no two of these are consecutive in the Fibonacci sequence.
The computation is simply the greedy algorithm of finding the highest
Fibonacci number below n, subtracting it and iterating.
List with components fibs the Fibonacci numbers that add sum up to
n, and inds their indices in the Fibonacci sequence.
zeck( 10) #=> 2 + 8 = 10 zeck( 100) #=> 3 + 8 + 89 = 100 zeck(1000) #=> 13 + 987 = 1000zeck( 10) #=> 2 + 8 = 10 zeck( 100) #=> 3 + 8 + 89 = 100 zeck(1000) #=> 13 + 987 = 1000