MATH36032 Project 1
Student ID:10299420
February 25, 2021
1 Estimating the Euler-Mascheroni Constant
The Euler-Mascheroni Constant (E-M Constant) γ defined as
γ = lim
n→∞
(
−ln(n) +
n∑
k=1
1
k
)
≈ 0.577215664901533
is a mathematical constant frequently arising in number theory and special integrals. Despite
its appearance in many seemingly unrelated integrals and functions, the E-M constant remains
a mystery and is not as well-known as other mathematical constants like π. Basic questions
like whether it is irrational or not are still open. In many situations, it is more convenient to
approximate complicated constants like this by a rational approximation.
We will explore algorithms which find the best rational approximation p/q, for p, q ∈ N such that
p + q ≤ N for any given N ∈ N. Here, the error ε is measured with respect to the L1 norm,
giving ε = |p
q
− γ|. It will also return p/q in its simplest form, which corresponds to the algorithm
implementing a simple check for whether gcd(p, q) = 1.
The simplest algorithm would iterate over all possible values of q for every value of p. However,
this is computationally expensive for large N as the number of pairs checked would equal
N(N−1)
2
.
A way to visualise the search range would be an N×N matrix A, with entries (aij) = i/j. The
algorithm searches across each row from left to right.
A =
1
1
1
2
. . . 1
N−2
1
N−1
2
1
2
2
. . . 1
N−1 0
…
… . .
.
0
…
N−2
1
N−2
2
. .
. …
…
N−1
1
0 . . . 0 0
The best estimates lie on a narrow diagonal band across the matrix, so an effective method is
to check all pairs (p, q) satisfying p + q = 2, . . . , N , searching diagonally across the matrix. We
introduce upper and lower bounds u and l respectively, corresponding to the best estimates above
and below the E-M constant found so far, so that the algorithm will only check values of p/q that
lie between the two bounds. Initial values are u = 1 for N = 1 and l = 1
2
for N = 2.
1
punctuate equation
Any references?
You haven’t explained what the E-M constant is?
In practice, this would impose limits on the range of p and q.
Let p+ q ≤ N and p/q ≤ u. Substituting for q gives an upper bound for q.
p(u+ 1)
u
≤ N ⇔ p ≤ N
u
u+ 1
Let p be fixed. Then l ≤ p/q ≤ u gives the bounds for q
p
u
≤ q ≤
p
l
This upper bound can be combined with p+ q ≤ N to give
q ≤ min{
p
l
,N − p}
As the algorithm progresses, the bounds u and l can be improved as follows:
For fixed p, p/q decreases as q increases. The limits on p and q shown above ensure we can update
the upper bound u as p/q decreases. If p/q ≤ γ, then we update the lower bound l and can skip
to the next iteration, since any larger values of q will make p/q ≤ l.
Implementing the MATLAB code (see A.1) for N = 2021 gives the output:
>> [p,q] = AppEm (2021)
p =
228
q =
395
2
Use same fonts as in equations
See appendix A1
2 Searching for MyLucky Numbers
An integer n is said to be a MyLucky number if
i) It has only odd prime factors that are distinct.
ii) The number itself is not prime.
iii) For any of its prime factors p, p− 1 divides n− 1.
An example is 10585 = 5× 29× 73. All its primes are odd and distinct, and
10584
4
= 2646,
10584
28
= 378,
10584
72
= 147
We propose the following algorithm to find the smallest MyLucky number greater than or equal
to any given N ∈ N. Since a MyLucky number only has odd factors, it is necessarily odd, so the
algorithm needs to only check odd numbers greater than or equal to N. The bulk of computation
involves the prime factorisation of each odd number. This can be implemented through the
following pseudo-code:
Input: an odd integer N
Output: all prime factors of N
let p be an empty list
for i = 3 to sqrt(N) step 2 do
while a mod n = 0
N = N/i
append i to p
if N=1 do
return p
if N != 1 do
append N to p
return p
In the case of MATLAB, the built-in function factor() can be used, which returns an array
containing the prime factors. Once it has been decomposed into its prime factors, it is an easy
task to check whether the three conditions hold. For condition (ii), it is sufficient to check if the
number of factors is greater than 1.
Below we display the output of the MATLAB code (see A.2) for N = 2021.
>> n = MyLuckynum (2021)
n =
2465
3
3 Finding Beautiful Square Numbers
An integer n is said to be a beautiful square number if its square consists of all nine non-zero digits
exactly once. An example is 118262 = 139854276. We will discuss methods to find the beautiful
square number n, such that n2 is closest to any given N ∈ N.
There are 9! = 362880 permutations of 123456789. If we narrow it to possible square numbers
(ending in 1,4,5,6,9), we still have 5 × 8! = 201600. It is much faster to search through possible
values of n, since
card({n ∈ N | 123456789 ≤ n2 ≤ 987654321}) =
⌈√
123456789
⌉
−
⌊√
987654321
⌋
= 20314
When checking whether n2 consists of all nine non-zero digits, in many programming languages
(including MATLAB) it is often convenient to first convert the number to a string, and then
compare them. In the MATLAB, strings are stored as ASCII values, with ASCII values for
numbers ranging from ‘0’ = 48 to ‘9’ = 57. The MATLAB code (see A.3) utilises this to perform
calculations directly on the ASCII values to retrieve the original digits of n2 and compare them.
Once all beautiful square numbers have been found, it is simple to check which beautiful square
number n minimises |N − n2|. A possible improvement would be to store a list of containing all
of beautiful square numbers, rather to begin the search all over every time the algorithm is run.
Running the MATLAB code (see A.3) for N = 360322021:
>> n = Beautisqnum (360322021)
n =
19023
4
Appendix A MATLAB code
A.1 Euler-Mascheroni Constant
function [p, q] = AppEm(N)
% AppEm: Approximate the E-M constant by the rational number p/q,
% among all pairs of positive integers (p,q) such that p + q <= N emconst = double(eulergamma); error = inf; % Minimum error attained bd = [0.5 ,1]; % Initial bounds num = 1; % p/q > E-M constant past this point
while num <= N*bd(2)/(bd(2) +1) % Limits p/q to between bounds and ensures p + q <= N for dnom = ceil(num/bd(2)):min(floor(num/bd(1)), N-num) % Only checks fractions in simplest form if gcd(num ,dnom) > 1
continue
% Updates values when a better approximation is found
elseif abs(num/dnom – emconst) <= error error = abs(num/dnom - emconst); p = num; q = dnom; end % Update upperbound as p/q decreases if num/dnom > emconst
bd(2) = num/dnom;
% If p/q < E-M constant , update lowerbound and increment p else bd(1) = num/dnom; break end end num = num + 1; end end 5 A.2 MyLucky Numbers function n = MyLuckynum(N) % MyLuckynum: Find the smallest MyLucky number that is greater than % or equal to N n = N; % Searches only odd numbers if mod(N,2)==0 n = n+1; end while true pf = factor(n); if length(pf) > 1 && all(diff(pf)) && all(mod(n-1,pf -1) ==0)
return
end
n = n+2;
end
end
A.3 Beautiful Square Numbers
function n = Beautisqnum(N)
% Beautisqnum: Return n such that n^2 is a perfect square number
% with all nine distinct digits from 1 to 9,
% which is closest to the input N, another 9-digit number
% (can be assumed to be between 10^8 and 10^9 – 1)
nums = 11112:31426;
% Checks if all 9 digits are unique
isbeaut = @(x) (numel(unique(num2str(x^2) – 48)) == 9);
beauty = nums(arrayfun(isbeaut , nums));
% Find closest beautiful square
[~,index] = min(abs(beauty .^2 – N));
n = beauty(index);
end
6
Code is very compact. Good to explain the line with arrayfun.