Duffinian Numbers

Duffinian numbers are defined as numbers that are composite and relatively prime to the sum of their divisors. The first such number is 4 because it is composite and the sum of its divisors (1, 2 and 4) is 7 and relatively prime to 4. The next composite number is 6 whose sum of divisors (1, 2, 3 and 6) is 12. However, since 12 and 6 are not relatively prime it is not Duffinian. The next composite number is 8 with divisors 1, 2, 4 and 8 whose sum is 15 and thus relatively prime to 8. Thus the first two Duffinian numbers are 4 and 8.

The first Duffinian numbers up to 247 are:

4, 8, 9, 16, 21, 25, 27, 32, 35, 36, 39, 49, 50, 55, 57, 63, 64, 65, 75, 77, 81, 85, 93, 98, 100, 111, 115, 119, 121, 125, 128, 129, 133, 143, 144, 155, 161, 169, 171, 175, 183, 185, 187, 189, 201, 203, 205, 209, 215, 217, 219, 221, 225, 235, 237, 242, 243, 245, 247

Since there are 59 numbers listed above, their density is 59/247 or nearly 24%. Thus they are more common than the prime numbers that occur with a frequency of around 18% in the same interval. There is an infinity of such numbers. As explained in Numbers Aplenty:

Any number of the form  \( p^k \)  where  \( p \)  is prime and  \( k \ge 2 \)  is Duffinian, since  \( \sigma(p^k)=1+p+p^2+\cdots+p^k \)  cannot be divisible by \(p \).

Of course, every prime \( p \) is relatively prime to the sum of its divisors \( p+1 \) but by definition the prime numbers are excluded. However, OEIS A014567 lists numbers \(n \) such that \( n \) and \( \sigma(n) \) are relatively prime, without \( n \) needing to be composite. Such numbers are sometimes called sigma primes and of course include all the primes. OEIS A003624 lists the Duffinian numbers.

Numbers Aplenty goes on to say that:

P. Heichelheim proved that there exists a run of five consecutive Duffinian numbers starting at 202605639573839041, and that there cannot exist a longer such run. 

Rose Mary Zbiek has proved that every even Duffinian number is either a square or twice a square. 

The smallest 3 × 3 magic square whose entries are consecutive Duffinian numbers is:

It's a simple enough matter in SAGE to determine whether or not a number is Duffinian. Here is the code (with "gcd" standing for "greatest common divisor"):

number=25387
if gcd(number,sigma(number,1))==1:
    print number, "is Duffinian"
else:
    print number, "is not Duffinian"

The input 25387 produces the output: 25387 is Duffinian

Ulam Numbers

I've been encountering these numbers for a while now, ever since making use of the website NUMBERS APLENTY which categorises numbers in a wide variety of ways and one of those categories is that of the Ulam number. Yesterday for example I was 25384 days old and 25384 is an Ulam number defined by the site as:

The Ulam sequence is defined by \(U_1=1, U_2=2 \) and, for \( k>2, U_k \) is the smallest integer that can be written in exactly one way as \(U_i+U_j \) with \( i<j<k \). The members of the Ulam sequence are called Ulam numbers. 

For example, after the first 4 terms which are trivially 1, 2, 3 and 4, the value of \(U_5 \) cannot be 5, since \(5=U_1+U_4=U_2+U_3 \), but it is 6, which can be obtained only as  \( U_2+U_4\)  (not as \(U_3+U_3 \) because the terms added must be distinct).

The problem with the Online Encyclopaedia of Integer Sequences (OEIS) is that it will not flag large numbers like 25384 as Ulam numbers because the numbers in the Ulam sequence (OEIS A002858) are only listed up to 339.

1, 2, 3, 4, 6, 8, 11, 13, 16, 18, 26, 28, 36, 38, 47, 48, 53, 57, 62, 69, 72, 77, 82, 87, 97, 99, 102, 106, 114, 126, 131, 138, 145, 148, 155, 175, 177, 180, 182, 189, 197, 206, 209, 219, 221, 236, 238, 241, 243, 253, 258, 260, 273, 282, 309, 316, 319, 324, 339

There is a list of the first 10,000 Ulam numbers in the links but these numbers do not show up in a search of the OEIS. That's why sites like NUMBERS APLENTY are so useful as a supplement to the OEIS. The site always provides a link to the relevant sequence in the OEIS as well as a Mathworld and Wikipedia link.

So what's special about the Ulam numbers? To begin with they are not distributed uniformly. As the NUMBERS APLENTY website states:

Ulam numbers are not distributed uniformly, but their density has peaks at an average distance of 21.6 (see Figure 1):

A plot of the Ulam numbers up to 1082, arranged line by line in a square 108×108. 

The plot shows the peaks in their density that occur with a frequency roughly equal to 108/5=21.6.

Wikipedia claims that the Ulam numbers have a density of approximately 0.07398 or 7.4%. It also states that:

There are infinitely many Ulam numbers. For, after the first n numbers in the sequence have already been determined, it is always possible to extend the sequence by one more element: \( U_{n − 1} + U_n \) is uniquely represented as a sum of two of the first n numbers, and there may be other smaller numbers that are also uniquely represented in this way, so the next element can be chosen as the smallest of these uniquely representable numbers.

Wikipedia goes on to say that:

The idea can be generalised as (u, v)-Ulam numbers by selecting different starting values (u, v). A sequence of (u, v)-Ulam numbers is regular if the sequence of differences between consecutive numbers in the sequence is eventually periodic. When v is an odd number greater than three, the (2, v)-Ulam numbers are regular. When v is congruent to 1 (mod 4) and at least five, the (4, v)-Ulam numbers are again regular. However, the Ulam numbers themselves do not appear to be regular. 

A sequence of numbers is said to be s-additive if each number in the sequence, after the initial 2s terms of the sequence, has exactly s representations as a sum of two previous numbers. Thus, the Ulam numbers and the (u, v)-Ulam numbers are 1-additive sequences. 

If a sequence is formed by appending the largest number with a unique representation as a sum of two earlier numbers, instead of appending the smallest uniquely representable number, then the resulting sequence is the sequence of Fibonacci numbers. 

This last point about the similarity with the Fibonacci numbers is interesting. Of course the Fibonacci numbers get very big very quickly while the Ulam numbers, because the smallest sum of any two previous numbers is sought, remain relatively smaller.

 Mathworld has a table of some of the other Ulam sequences:

I was struggling to get my SAGE code to generate these sequences of terms but came across some Python 2 code (source) that does the job, although I'm struggling to understand how it works:

ulams=[ ]
ulams.append(1)
ulams.append(2)
limit=100
sums=[0 for i in range(2*limit)]
newUlam=2
sumIndex=1
while(newUlam <limit):
  for i in ulams:
if (i<newUlam):
sums[i+newUlam] += 1
while(sums[sumIndex]!=1):
sumIndex += 1
newUlam = sumIndex
sumIndex += 1
ulams.append(newUlam) for i in ulams:
print i

With the limit set to 100, the following output is produced: 1 2 3 4 6 8 11 13 16 18 26 28 36 38 47 48 53 57 62 69 72 77 82 87 97 99 102

The code can be easily modified to generate other sequences with different starting terms e.g. 2 and 3:

ulams=[ ]
ulams.append(2)
ulams.append(3)
limit=100
sums=[0 for i in range(2*limit)]
newUlam=3
sumIndex=1
while(newUlam <limit):
  for i in ulams:
if (i<newUlam):
sums[i+newUlam] += 1
while(sums[sumIndex]!=1):
sumIndex += 1
newUlam = sumIndex
sumIndex += 1
ulams.append(newUlam) for i in ulams:
print i

This code produces the output: 2 3 5 7 8 9 13 14 18 19 24 25 29 30 35 36 40 41 46 51 56 63 68 72 73 78 79 83 84 89 94 115

ADDENDUM: I finally developed my own code which I do understand and it works! Here are the terms below 200 in the Ulam sequence beginning 2,3:

ulam=[i for i in range(1,200)]
ulam[0], ulam[1]=2,3
k=2
while k<len(ulam):
    count=0
    for x in range(k):
        for y in range(k):
            if ulam[x] > ulam[y]:
                if ulam[x] + ulam[y] == ulam[k]:
                    count+=1
    if count!=1:
        ulam.remove(ulam[k])
        k-=1
    k+=1
print ulam  

This code produces the output: 2, 3, 5, 7, 8, 9, 13, 14, 18, 19, 24, 25, 29, 30, 35, 36, 40, 41, 46, 51, 56, 63, 68, 72, 73, 78, 79, 83, 84, 89, 94, 115, 117, 126, 153, 160, 165, 169, 170, 175, 176, 181, 186, 191

Here is a link to an interesting article that looks at patterns in these 1-additive sequences.

Apéry's Constant

Today I turned 25355 days old and this number turns up in the Engel expansion of \( \zeta(3) \). Firstly however, let's remind ourselves that \( \zeta \) is the Riemann zeta function and can be expressed as:$$ \zeta(3)=\sum_1^{\infty} \frac{1}{n^3} $$ $$ \text{or  } \zeta(3)=\lim_{n \rightarrow \infty} \left( \frac{1}{1^3}+\frac{1}{2^3}+ \cdots + \frac{1}{n^3} \right) $$This works out to around 1.202056903159594285399738161511449990764986292 and is known as Apèry's constant. I've written about this constant before in a post titled The Basel Problem and Beyond on May 7th 2017. According to Wikipedia:

This constant arises naturally in a number of physical problems, including in the second-and third-order terms of the electron's gyromagnetic ratio using quantum electrodynamics. It also arises in the analysis of random minimum spanning trees and in conjunction with the gamma function when solving certain integrals involving exponential functions in a quotient which appear occasionally in physics, for instance when evaluating the two-dimensional case of the Debye model and the Stefan–Boltzmann law.

This doesn't mean much to me but there you have it. As for the Engel expansion, I made a blog post about it in 2016. The expansion consists of the terms that go to make up the denominators of the fractions that when added will approximate \( \zeta(3) \). OEIS A053980 lists the first few terms as: 1, 5, 98, 127, 923, 5474, 16490, 25355

A053980

Engel expansion of zeta(3) = 1.20206...

This means that \( \zeta(3) \) can be approximated as:$$ \zeta(3) \approx \frac{1}{1}+\frac{1}{1 \times 5} +\frac{1}{1 \times 5 \times 98}+\frac{1}{1 \times 5 \times 98 \times 127} +\frac{1}{1 \times 5 \times 98 \times 127} + \cdots $$Of course, once the zeta function is touched upon, one can find oneself in very deep water very quickly so I'm not going to say too much more except to include a screenshot (Figure 1) from Wolfram MathWorld showing the Engel expansions for some of the other constants:

Figure 1

The SAGE code for generating these sequences is fairly straightforward. Figure 2 shows what's involved for Apèry's constant:

Figure 2: permalink

Of course, simply replacing u in the above code by say \(e\) or \(\pi \) will produce the associated Engel expansion.

on May 12th 2021

Fermat Pseudoprimes

Today I turned 25351 days old and the number 25351 turns out to have some quite interesting properties connected with pseudoprimes. To quote from Wikipedia:

In number theory, the Fermat pseudoprimes make up the most important class of pseudoprimes that come from Fermat's little theorem. Fermat's little theorem states that if p is prime and a is coprime to p, then \(a^{p−1} − 1 \) is divisible by p. 

For an integer a > 1, if a composite integer x divides \( a^{x−1} − 1, \) then x is called a Fermat pseudoprime to base a. In other words, a composite integer is a Fermat pseudoprime to base a if it successfully passes the Fermat primality test for the base a. It follows that if x is a Fermat pseudoprime to base a, then x is coprime to a. 

The smallest base-2 Fermat pseudoprime is 341. It is not a prime, since it equals 11·31, but it satisfies Fermat's little theorem: \( 2^{340} ≡ 1 \) (mod 341) and thus passes the Fermat primality test for the base 2. Pseudoprimes to base 2 are sometimes called Poulet numbers, after the Belgian mathematician Paul Poulet, Sarrus numbers, or Fermatians (sequence A001567 in the OEIS). 

A Fermat pseudoprime is often called a pseudoprime, with the modifier Fermat being understood. An integer x that is a Fermat pseudoprime for all values of a that are coprime to x is called a Carmichael number.

Now 25351 is a member of OEIS A005936: Pseudoprimes to base 5. This is because \( 5^{25350} ≡ 1 \) (mod 25351). However, it turns out that 25351 is not simply a Fermat pseudoprime. It is also a strong pseudoprime defined by Wikipedia as:

In number theory, a probable prime is a number that passes a primality test. A strong probable prime is a number that passes a strong version of a primality test. All primes pass these tests, but a small fraction of composites also pass, making them "false primes". Unlike the Fermat pseudoprimes, for which there exist numbers that are pseudoprimes to all coprime bases (the Carmichael numbers), there are no composites that are strong pseudoprimes to all bases.

Formally, an odd composite number \(n = d \times 2^s + 1 \) with d also odd is called a strong (Fermat) pseudoprime to a base a when one of the following conditions holds:$$ a^{d} \equiv 1\! \mod n \text{ or}$$$$a^{d\cdot 2^{r}}\equiv -1 \!\mod n\quad {\mbox{ for some }}0\leq r < s $$ If a number n satisfies one of the above conditions and we don't yet know whether it is prime, it is more precise to refer to it as a strong probable prime to base a. But if we know that n is not prime, then one may use the term strong pseudoprime. The definition is trivially met if a ≡ ±1 mod n so these trivial bases are often excluded. Guy mistakenly gives a definition with only the first condition, which is not satisfied by all primes.

Now 25351 = 2 * 12675 + 1 satisfies the \(n = 12675\times 2^s + 1 \) form and 

\( 5^{12675} ≡ 1 \) (mod 25351). Hence 25351 is a strong pseudoprime. This condition may seem somewhat arbitrary but there is good explanation of why it's so at WolframAlpha.

Things get more complicated still when we learn that 25351 is also a Euler pseudoprime but I'm not going to delve any deeper at the moment. That's probably enough except that I'll finish up by saying a little about the Carmichael numbers that satisfy Fermat's little theorem to all coprime bases. The first of these numbers are 561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841 and 29341. Carmichael numbers have at least three positive prime factors e.g. 561 = 3 * 11 * 17.

ADDENDUM: entered on March 11th 2020

I'm revisiting pseudoprimes because it's a topic about which I always seem to get confused. I've mentioned them previously in the following posts:

• Fermat Pseudoprimes on the 30th August 2018
• Carmichael Numbers on the 30th January 2019  
• The Original Taxi Cab Number in a New Light on the 21st December 2019

What brought to mind again was the number 25909, the number of days old that I was yesterday. One of its properties is that it's a member of OEIS A020267: strong pseudoprimes to base 41. Now it's easy enough to show that it's a pseudoprime to base 41 because$$41^{25908}-1 \equiv 0 \bmod{25909}$$However, to show 25909 is a strong pseudoprime, it needs to fulfil one of two possible conditions. Formally, an odd composite number \(n = d \times 2^s + 1\) with d also odd is called a strong (Fermat) pseudoprime to a base, when one of the following conditions holds:$$ a^{d} \equiv 1\! \mod n \text{ or}$$$$a^{d\cdot 2^{r}}\equiv -1 \!\mod n\quad {\mbox{ for some }}0\leq r \leq s $$So in the case of 25909 we have: $$25909=25908+1=2^2 \times 3 \times 17 \times 127+1=6477 \times 2^2+1$$The value of \(d\) here is 6477 so let's test it out using the first condition. We find that:$$41^{6477} \equiv 8805 \bmod{25909}$$Clearly it does not pass the first test, so let's try the second.$$41^{6477*2} \equiv -1 \bmod 25909$$It passes and so 25909 can claim to be a strong pseudoprime. Note that \( 0 \leq r < s\) and \(r\) must be 0 or 1 and the latter satifies.

Quadratic Reciprocity

I just read a most interesting article in Quanta Magazine about the new Fields medallist Peter Scholze. His discoveries are quite beyond my comprehension but in the article mention is made of quadratic reciprocity. I've come across the concept repeatedly but have always avoided taking the time to understand it. However, in this article it was explained simply and clearly. 

Reciprocity laws are generalisations of the 200-year-old quadratic reciprocity law, a cornerstone of number theory and one of Scholze’s personal favourite theorems. The law states that given two prime numbers p and q, in most cases p is a perfect square on a clock with q hours exactly when q is a perfect square on a clock with p hours. For example, five is a perfect square on a clock with 11 hours, since \(5 = 16 = 4^2 \), and 11 is a perfect square on a clock with five hours, since \(11 = 1 = 1^2 \).

Here is a more detailed explanation of quadratic reciprocity (source):

It's easy to find integer solutions of the equation \( y^2 = 2x^2 + 7k \), but there are no integer solutions at all of the seemingly similar equation \( y^2 = 3x^2 + 7k \).  This is the sort of thing that experienced number theorists can tell at a glance, simply by noting that the equations written modulo 7 are \( y^2 = 2x^2 \) and \( y^2 = 3x^2 \) respectively, and since division is unique in the field of integers (mod p) we have \( (y/x)^2 \equiv n \pmod{p} \), which implies that n must be the square of some element of the field of integers modulo p.  But the integers mod 7 are just 0, 1, 2, 3, 4, 5, 6, whose squares (mod 7) are 0, 1, 4, 2, 2, 4, 1. Thus the equation \( y^2 = nx^2 + 7k \) can have integer solutions only if n is congruent to 0, 1, 2, or 4 (mod 7).  These are called the quadratic residues (mod 7), and the remaining numbers 3, 5, 6 are called the non-quadratic residues.  

It's often extremely important when dealing with problems in number theory to know whether a certain prime p is a square (i.e., a quadratic residue) modulo some other particular prime q. Legendre defined a symbol to represent this information, which we will denote as

This is a remarkable fact, and not at all self-evident.  Legendre succeeded in proving some special cases of this, and also gave what he thought was a complete proof, but his argument relied on the premise that every arithmetic progression contains infinitely many primes. This is in fact true, as subsequently shown by Dirichlet, but at the time of Legendre's proof it wasn't known, so Gauss pointed out that Legendre's proof was incomplete. Gauss went even further, and gave several (valid) proofs of this remarkable theorem, which he called the Fundamental Theorem of number theory.

Thus [p\q] equals either +1 or -1 depending on whether p is or isn't a square modulo q. From examining many individual cases, Euler had previously noticed a striking relationship between the [p\q] and [q\p].  Specifically, he noticed that [p\q] = [q\p] except when p and q are both of the form 4k-1, in which case [p\q] = -[q\p]. 

Let's take a specific example. Yesterday I turned 25339 days old. The previous prime is 25321. It can seen that \( 25339 \equiv 18 \pmod{25321} \) and \( 25321 \equiv 25321 \pmod{253339} \), neither residue being square. Therefore \([25339\25321] = [25321\25339] = -1\). Note that 25321 is a 4k+1 prime while 25339 is a 4k-1 prime, which is in accord with the theorem.

I've download a PDF titled Quadratic Residues and Non-Residues Selected Topics by Steve Wright from the Department of Mathematics and Statistics Oakland University, Rochester, Michigan and dated 21st October 2106. It's freely available here. I'll have a read over this and see if I can deepen my understanding of this topic.

Here are some SAGE commands that can be used when dealing with quadratic residues (double click to enlarge):

Source: http://doc.sagemath.org/html/en/constructions/number_theory.html

ADDENDUM:

When trying to tackle the above Quadratic Residues and Non-Residues Selected Topics, I quickly realised that my foundations in modular arithmetic were a little shaky and needed to be reinforced a little. I began watching a series of videos by Polar Pi on YouTube, beginning with this video on a linear congruence equation with one solution:

The next video in this series relates to linear congruence equations with multiple solutions and the third of the three videos looks at linear congruence equations with no solutions.

The YouTube channel blackpenredpen has some videos on modular arithmetic as well:

• What does a ≡ b (mod n) mean? Basic Modular Arithmetic, Congruence? Link
• Solving congruences, 3 introductory examples. Link
• System of congruences, modular arithmetic. Link

Inconsummate Numbers

A number \(n\) is inconsummate if there is no number \(k\) which divided by its sum of digits gives \(n\). An example of a number that is not inconsummate would be 47 because there exists the number 846 such that 846 / (8+4+6) = 47. Notice that in order to confirm this, it was necessary to test numbers far ahead of 47 on the number line. The following SAGE code is designed to identify inconsummate numbers:

INPUT
# Generate a list of inconsummate numbers
# n is an inconsummate if no number k divided by its sum of digits gives n
# The larger the value of limit, the more accurate the list
limit=1000
inconsummate=[]
for number in range(1, limit):
    digits=(str(number))
    sum=0
    for n in range(0,len(digits)):
        sum+=Integer(digits[n])
    if number % sum == 0:
        inconsummate.append(number/sum)
X=Set(inconsummate)
Y=Set(range(1, limit+1))
print(Y.difference(X)) 

OUTPUT
{62, 63, 65, 66, 71, 72, 74, 75, ... }

Unfortunately after 65, the next inconsummate number is 75 and so 66, 71, 72 and 74 should be not included in the output. Notice that the limit was set to 1000. If we increase the limit to 10000, the output goes wrong here { ..., 432, 437, 438, ... } because 437 is not an inconsummate number. By setting the limit to 12000, this particular error is corrected. It would seem that by choosing a particular range (from 1 to the chosen limit), the resultant output is only reliable for about 1/30th of that range.

It's interesting to look at the frequency of inconsummate numbers:

up to \(10^1\) --> 0 --> 0%
up to \(10^2\) --> 6 --> 6%
up to \(10^3\) --> 111 --> 11.1%
up to \(10^4\) --> 1437 --> 14.37%
up to \(10^5\) --> 16430 --> 16.43%
up to \(10^6\) --> 183089 --> 18.3089%
up to \(10^7\) --> 1905285 --> 19.05285% seems to peak here
up to \(10^8\) --> 18907944 --> 18.907944%
up to \(10^9\) --> 183706706 --> 18.3706706%

Figure 1: 3 x 3 magic square 

For the ranges shown, the frequency reaches a peak of 19.05285% and then slowly decreases. The smallest 3 × 3 magic square whose entries are consecutive inconsummate numbers is shown in Figure 1. Source.

The initial inconsummate numbers are listed in OEIS A003635 and the initial members are:

62, 63, 65, 75, 84, 95, 161, 173, 195, 216, 261, 266, 272, 276, 326, 371, 372, 377, 381, 383, 386, 387, 395, 411, 416, 422, 426, 431, 432, 438, 441, 443, 461, 466, 471, 476, 482, 483, 486, 488, 491, 492, 493, 494, 497, 498, 516, 521, 522, 527, 531, 533, 536, ... 

Here is a link to a Google Document containing the full list.

on 29th August 2021, 19th February 2023 and 26th April 2026

Sums of Cubes and Squares of Sums

Today I turned 25321 days old, a prime number of days, and a prime with the property that the sum of the cubes of its digits equals the square of the sum of its digits. In other words:$$2^3+5^3+3^3+2^3+1^3 = (2+5+3+2+1)^2 = 169$$This prime is a member of OEIS A225567: Primes with nonzero digits such that sum of cubes of digits equal to square of sums. It's initial members are:$$1423, 2143, 2341, 4231, 12253, 21523, 22153, 22531, 23251, 25321, ...$$Serendipitously I then rediscovered the old and famous connection between the sum of the cubes of the first n natural numbers and the square of the sum of these some numbers, specifically:$$ \sum_1^n i^3=\big (\sum_1^n i \big)^2$$It's easy to see why the first four members of the sequence (1423, 2143, 2341 and 4231 are members) because these are simply instances of:$$1^3+2^3+3^3+4^3=(1+2+3+4)^2$$I also discovered a similar relationship involving the divisors \(d_i \) of any natural number with \(n\) divisors, namely that:$$ \sum_1^n (\sigma_0(d_i))^3=\big(\sum_1^n \sigma_0(d_i) \big) ^2$$The previous looks more difficult than it actually is and a simple example will assist. Let's consider the number 10. It has four divisors 1, 2, 5 and 10. Each of these divisors has 1, 2, 2 and 4 divisors respectively. We find that:$$1^3+2^3+2^3+4^3=(1+2+2+4)^2 = 81$$It's as simple as that and it applies to every natural number. Of course 1224 and permutations of these digits can be found in OEIS A227073: Positive numbers without the digit 0 such that sum of cubes of the digits equals the square of the sum of the digits. The initial members of this sequence are:$$ 1, 12, 21, 22, 123, 132, 213, 231, 312, 321, 333, 1224, ...$$