Prime Pride

I was born in 1949, a prime year ending in 49 and the only such year in the thousand year period from 1549 to 2549, both of which mark prime years. Here is a list of prime years ending in 49 over the centuries:

149 349 449 1049 1249 1549 1949 2549 2749 3049 3449 4049 4349 4549 4649 5449 5749 5849 6449 6949 7349 7549 7649 7949 8849 9049 9349 9649 9749 9949

Of the 100 possible candidates that end in 49 over the above ten thousand year period, 30 of them are prime. As a percentage, this is 30%.

In my case, I was born on April 3rd 1949 and concatenating the digits 4, 3 and 49 (as was commonly done during the second half of the twentieth century) produces 4349 which is also prime. If the date is written as 3rd April 1949, then concatenation of the digits 3, 4 and 1949 produces 3449 which is again prime. Composite numbers arise from concatenating 3 and 4 or 4 and 3 with 1949.

Furthermore, I was born at a prime hour: 7 o'clock in the morning and this number remains the same with either a 12 or 24 hour clock. At that time, the Sun was in the 13th degree of Aries (12º47'13" to be precise) and thus in a prime degree. The zodiac is divided into 360 parts, measured from the First Point of Aries. The first degree begins at this point and the 360th degree ends at this same point.

In addition, my birth name (xxxx xxxxxxx xxxxxx) consists of 17 letters. I changed my first name later in life but even so my new name still contained four letters and so the total was still 17.

I'm cherry-picking I guess because there have been plenty of composite numbers in my life. I lived for the first two years of my life at my grandparents' house numbered 69 and another twenty years or more at my parents' house numbered 21. Nevertheless, it's interesting to discover the prime numbers associated with my name and time/date of birth.

Unitary Divisors

My diurnal age today is 25245 and the OEIS A127666 mentions this number as belonging to the sequence of odd infinitary abundant numbers. Unfortunately, I had no idea what was meant by the term infinitary but I was determined to find out. This led me into deeper waters very quickly and I realised that it might be best to start in the paddle pool first by investigating the term unitary.

To quote from Wikipedia:

A natural number a is a unitary divisor (or Hall divisor) of a number b if a is a divisor of b and if a and b/a are coprime, having no common factor other than 1. Thus, 5 is a unitary divisor of 60, because 5 and 60/5 =12 have only 1 as a common factor, while 6 is a divisor but not a unitary divisor of 60, as 6 and 60/10 have a common factor other than 1, namely 2. 1 is a unitary divisor of every natural number. 

Equivalently, a given divisor a of b is a unitary divisor if and only if every prime factor of a has the same multiplicity in a as it has in b. 

The sum of unitary divisors function is denoted by the lowercase Greek letter sigma thus: $\sigma  ^*(n)$. The sum of the $k \, th$ powers of the unitary divisors is denoted by $\sigma_k ^*(n)$: $$\sigma_k^*(n) = \sum_{d\mid n \atop \gcd(d,n/d)=1} d^k$$ If the proper unitary divisors of a given number add up to that number, then that number is called a unitary perfect number.

Now 60 turns out to be a unitary perfect number, one of very few in fact that are known. The list runs:

6, 60, 90, 87360, 146361946186458562560000 

"Some perfect numbers are not unitary perfect numbers, and some unitary perfect numbers are not regular perfect numbers ... It is not known whether or not there are infinitely many unitary perfect numbers, or indeed whether there are any further examples beyond the five already known" from Wikipedia.

The program code I used to determine the unitary divisors of a number in SageMath is (permalink):

#FIND UNITARY DIVISORS AND TEST FOR UNITARY PERFECT NUMBERS

number = 87360

total = 0

D = divisors(number)

print("Unitary Divisors are:")

for i in range(0, len(D)):

    if gcd(D[i], number/D[i]) == 1:

        print(D[i], end=" ")

        total = total+D[i]

print()

print("Sum of unitary divisors less",number,"is",total-number)

if total-number == number:

    print("Therefore",number,"is a unitary perfect number")

else:

    print("Therefore",number,"is not a unitary perfect number")

Unitary Divisors are:

1 3 5 7 13 15 21 35 39 64 65 91 105 192 195 273 320 448 455 832 960 1344 1365 2240 2496 4160 5824 6720 12480 17472 29120 87360 

Sum of unitary divisors less 87360 is 87360

Therefore 87360 is a unitary perfect number

I started off this post by investigating the term unitary and ended up focusing on unitary divisors rather than unitary numbers. Let's address that deficiency now. If the sum of the unitary divisors of a number is greater than the number, the number can be described as a unitary abundant number. Such numbers can be odd or even, with the former being much rarer. If the sum of the unitary divisors is less than the number, the number can be described as a unitary deficient number. As we've seen, if the sum of the unitary divisors is equal to the number, the number can be described as a unitary perfect number.

Of course, I've not examined in this post what are meant by infinitary divisors. That will have to wait for a future post.

Necklaces and Bracelets

At my recent diurnal age of 25136 (I'm 25138 days old as I write this), I noted that:

25236 is also a member of OEIS A141783: number of bracelets (turn over necklaces) with n beads: 1 blue, 12 green, and r = n-13 red (here n=20 and so r=7). This can be restated as "the number of bracelets with 1 blue, 12 green and 7 red beads is 25236". 

I had trouble understanding why the number of bracelets was 25236 and not 25194. According to my initial calculations, the number of necklaces (arrangements of beads that can't be lifted out of the plane and reversed) should be 19!/(12!*7!) = 50388. Here the numerator is 19! because 20!/20 adjusts for the rotational symmetry of the twenty beads and the denominator takes into account the fact that the beads in each of the groups, 12 green and 7 red, are indistinguishable from one another. Because bracelets can be lifted out of the plane and turned over, the number of bracelets should be 50338/2 = 25194.

As with most things, it's better to start at the beginning. The sequence runs: 1, 7, 49, 231, 924, 3108, 9324, 25236, ... and begins with n=13 at which point there is only one possible configuration as there are no red beads. When n=14, there is one red bead, one blue bead and twelve green beads. I would reckon the number of possible necklace configurations to be: 13!/12!=13 but in fact there is one more as can be seen in Diagram 1 and 2:

Diagram 1: blue bead on left, red bead on right

Diagram 2: red bead on left, blue bead on right.

As the two diagrams show, the first configuration with a blue bead on the left and a red bead on the right can be rotated to produce a new configuration where the red bead is on the left and the blue bed is on the right. This means that there are 13+1=14 possible necklaces and thus 14/2=7 possible bracelets.

When there are 20 beads, the situation is of course more complicated but the principle is the same. There will be the 50388 necklaces as calculated above but there are also the extra necklaces created by rotating the blue and red beads by 180° for the symmetric configurations like the one shown in Diagram 3. In every symmetric configuration, there will be nine beads (six green and three red) on either side of the red-blue bead axis. There are $\binom{9}{3} = 84$ ways this can be done and so the total number of necklaces is 50388+84=50472, corresponding to 50472/2=25236 bracelets.

Diagram 3: a symmetrical position with green and red

identically placed on either side of the red-blue axis

Polygonal Number Generating Function and Formula

Today my diurnal age is 25233, a number that happens to be a polygonal number, specifically a 32-gonal number. A WolframAlpha article states that the generating function for the n-gonal numbers is given by: $$G_n(x)= \frac{x \, [(n-3)x+1]}{(1-x)^3}$$ This means that for the 32-gonal numbers, the formula becomes: $$G_{32}(x)= \frac{x \, (33x+1)}{(1-x)^3}$$ In WolframAlpha, the coefficients can then be identified using the series command:

In Sage, the Taylor series is generated using the code:

g(x)=x*(33*x+1)/(1-x)^3
g.taylor(x,0,15).coefficients()

When run, this code produces the following output:

In the case of the 36-gonal numbers, the formula for the n-th term is given by: $$a(n)=n(17n-16)$$ The general formula for the n-th polygonal number is given by: $$a(n)=(P-2) \frac{n(n+1)}{2}-(P-3)n$$ where $P$ is the number of vertices of the polygon.

In the case of $P=36$ (our 36-gon), the formula becomes: $$34\frac{n(n+1)}{2}-33n=n(17n-16)$$ Here's a video describing how the formula is derived:

Number of Divisors

I came across this formula for calculating the number of divisors of a given number: $$\sigma_0(n)=\prod_{i=1}^r(a^i+1)$$ where $r$ is the number of distinct prime factors and $a^i$ are their coefficients. For example, $$\sigma_0(24)=\prod_{i=1}^2(a^i+1)=(3+1)(1+1)=8$$ where $r$ = 2 because there are two distinct prime factors (2 and 3) with indices of 3 and 1 respectively.

There are indeed eight factors of $24: 1, 2, 4, 8, 3, 6, 12, 24$ so this is a quick and easy way to calculate the number of factors based on prime factors of a number and the degrees to which each of them are raised. It's similar to the formula for calculating the number of ways in which a number can be written as the sum of two squares viz.:

The fundamental theorem on sums of two squares is:

Let $n=2^kp_1^{a_1} \dots p_r^{a_r}q_1^{b_1} \dots q_s^{b_s}$ , where the $p_i$ are distinct primes with $p_i \equiv 1 \pmod{4}$ and the $q_i$ are distinct primes with $q_j \equiv 3 \pmod{4}$. Then $n$ is the sum of two squares if and only if all the $b_j$ are even. In that case, the number of distinct solutions is $$\lceil \frac{\scriptstyle{1}}{\scriptstyle{2}} (a_1+1)(a_2+1) \dots (a_r+1) \rceil$$ where $\lceil x \rceil$ is the ceiling function, the smallest integer greater than or equal to $x$.

Of course by this rule, $24$ cannot be written as a sum of two squares because it has a $4k+3$ factor ($3$), raised to an odd power ($1$). However, $25=5^2$ has no $4k+3$ factors and its $4k+1$ factor ($5$) is raised to the $2nd$ power. This means that it can be written as the sum of two squares in two distinct ways because $\lceil \frac{1}{2}(2+1) \rceil =2$. The two ways are $5^2+0^2$ and $4^2+3^2$.

Getting back to the number of divisors however, the formula quoted is actually part of a more general formula (as explained in Wikipedia):

The divisor function is multiplicative, but not completely multiplicative.  The consequence of this is that, if we write $$n = \prod_{i=1}^r p_i^{a_i}$$ where $r=\omega (n)$ is the number of distinct prime factors of $n$, $p_i$ is the $i$th prime factor, and $a_i$ is the maximum power of $p_i$ by which $n$ is divisible, then we have $$\sigma_x(n) = \prod_{i=1}^{r} \frac{p_{i}^{(a_{i}+1)x}-1}{p_{i}^x-1}$$ which is equivalent to the useful formula: $$\sigma_x(n) = \prod_{i=1}^r \sum_{j=0}^{a_i} p_i^{j x} = \prod_{i=1}^r (1 + p_i^x + p_i^{2x} + \cdots + p_i^{a_i x}).$$ It follows (by setting $x = 0$) that $d(n)$ is: $$\sigma_0(n)=\prod_{i=1}^r (a_i+1).$$

A Look at Sums of Multiples of Squares

I certainly have a deeper understanding now of what numbers can be expressed as a sum of two squares and in how many ways. See my earlier posts:

• Sums of Square Revisited
• Fundamental Theorem on Sums of Two Squares
• Numbers that are Sums of Two Squares
• Sums of Squares
• Finding Significance in Insignificant Numbers
• Sum of Two Squares
• Sum of Squares of Integers and Catalan Numbers

There's certainly some overlap in these posts but mostly they've dealt with numbers of the form $x^2+y^2$ rather than $ax^2+by^2$ where $a$ and $b$ are integers. It's the latter type of expression that I want to explore in this post. My interest was stimulated by today's diurnal age number, 25229, that happens to be a 4k+1 prime.

As a 4k+1 prime, it can be expressed as a sum of two squares in one way only, viz. $98^2+125^2$. However, it turns out that it can be expressed a sum of two squares multiplied by varying coefficients in many different ways. I created some Sage code to investigate:

number=25229 #enter any number
x,y = var('x'),var('y')
a,b = var('a'), var('b')
for a in range(100):
    for b in range(100):
        for x in range(sqrt(number)):
            for y in range(sqrt(number)):
                if a * x^2 + b * y^2 == number:
                    if a < b:
                        print(a, b, x, y)

Running this code on the SageMathCell Server generated the following results for a <100 and b < 100. The first number represents a, the second b, the third x and the fourth y. For example, the first entry (1, 4, 125, 49) represents $125^2+4 \times 49^2$.

(1, 4, 125, 49)
(1, 5, 27, 70)
(1, 7, 41, 58)
(1, 13, 101, 34)
(1, 20, 27, 35)
(1, 25, 98, 25)
(1, 28, 41, 29)
(1, 29, 75, 26)
(1, 49, 125, 14)
(1, 52, 101, 17)
(1, 65, 42, 19)
(1, 85, 152, 5)
(1, 91, 127, 10)
(1, 97, 94, 13)
(2, 3, 61, 77)
(2, 11, 105, 17)
(2, 21, 110, 7)
(2, 51, 37, 21)
(3, 26, 9, 31)
(4, 25, 49, 25)
(4, 65, 21, 19)
(4, 85, 76, 5)
(4, 97, 47, 13)
(5, 6, 71, 2)
(5, 9, 70, 9)
(5, 19, 11, 36)
(5, 24, 71, 1)
(5, 46, 61, 12)
(5, 59, 45, 16)
(5, 69, 8, 19)
(5, 76, 11, 18)
(5, 81, 70, 3)
(5, 89, 69, 4)
(6, 53, 56, 11)
(7, 22, 31, 29)
(7, 29, 60, 1)
(7, 46, 53, 11)
(8, 21, 55, 7)
(9, 20, 9, 35)
(9, 29, 25, 26)
(9, 65, 14, 19)
(11, 18, 17, 35)
(11, 50, 17, 21)
(11, 98, 17, 15)
(13, 29, 36, 17) 

It would be interesting to see the full range of possibilities for a, b, x and y. However, for larger values of a and b, the SageMathCell Server seems to time out. It is easy to see how the case of (1, 4, 125, 49) comes about because: $$125^2+(2 \times 49)^2 =125^2+(2 \times 49)^2 = 125^2+ 4 \times 49^2$$ However, in the case of (2, 3, 61, 77) where $(\sqrt{2} \times 61)^2+ (\sqrt{3} \times 77)^2 = 25229$, it's far from clear how the result comes about.

It's probably better to start investigating with a smaller 4k+1 prime such as 61. Now $61=5^2+6^2$ and there is only one possibility with coefficients and that is $4 \times 2^2+5 \times 3^2$. Similarly with $73=3^2+8^2$, there is only one possibility and that is $5 \times 3^2+7 \times 2^2$. With $97=4^2+9^2$, there are two possibilities: $4 \times 2^2 + 9 \times 3^2$ and $5 \times 3^2+13 \times 2^2$.

It can be noted that 4k+3 primes such as 43, while not able to be represented as a sum of two squares, can be represented as a sum of squares with coefficients: $3 \times 3^2+4 \times 2^2$. Essentially, the elements of the set of all multiples of the square numbers (let's call the set S) are being combined with themselves by addition to form a set of new numbers (let's call this set N). This set consists of $1^2 =1$ and the multiples of 1 are the natural numbers. This of course is trivial because it means for any number n can be written as $(n-m) \times 1^2 + m \times 1^2$ where $m<n$.

This result means that it is necessary to exclude 1 from our investigation and begin with $2^2=4$. So all multiples of 4 are permissible, as are all multiples of 9, 25, 36, 49, 64, 81 and so on. In this scheme, the smallest possible number would be $2^2+2 \times 2^2=12$ with $2^2+2^2=8$ being ignored since we know the rules governing the requirements for a number to be a sum of two squares. This means that coefficients must not be perfect squares and thus $2^2+4 \times 2^2 =2^2+4^2 = 20$ is not acceptable. Even coefficients that add to perfect squares when the squared terms are the same must be excluded.

Thus $2^2+3 \times 2^2 = 2^2 \times (1+3) = 2^2 \times 4 = 4^2 =16$ falls victim. From the outset, it's clear that certain numbers (apart from the numbers 1 to 11) cannot be represented by any combination of multiples of squares: 14, 15, 17, 18, 19, 20, 21 and 23. However, it seems clear that every number above 23 can be represented either as:

• a perfect square e.g. $25=5^2$
• a multiple of a perfect square e.g. $50=2 \times 25$
• a sum of two different squares e.g. $25=4^2+3^2$
• a sum of the multiples of two squares e.g. $25=4 \times 2^2 +3^2$. 

I've not proved this by any means but it seems highly likely, given the increase in possibilities as the numbers get larger. So the main takeaway from all this is that $ax^2+by^2=23$ is the last ellipse in the family of ellipses $ax^2+by^2=N$, where N is a natural number and a>1 and b>1, that does not have integer solutions for a, b, x and y.