Elliptical Curves

I recently read an interesting article in Quanta Magazine that involved elliptic curves. It's a fairly high level article but an excerpt read:

Their starting point was one of the central objects in number theory: elliptic curves. Just like circles and lines, elliptic curves are both numbers and shapes. They are pairs of numbers, \(x\) and \(y\), that serve as solutions to an algebraic equation like \(y^2 = x^3 − 2x\). The graph of those solutions creates a geometric shape that looks vaguely like a vertical line extruding a bubble.

This was followed by the diagram shown in Figure 1:

Figure 1

It turns out that these elliptic curves are very important in cryptography as this following YouTube video explained:

Many elliptic curves can be written as \(y^2=x^3+ax+b\) and all the curves shown in Figure 1 are of this type. What's of interest to me is that there are integer solutions to these sorts of equations, in other words there are points \( (x,y) \) such that \(x\) and \(y\) are both integers. These could be termed elliptic Diophantine equations.

Let's consider \(y^2=x^3-x+1\) as a particular example (the graph is shown in Figure 1, bottom left). When \(x=0\), \(y \pm 1\) but what about when \(y=0\)? Here we have \(x^3-x+1=0\) and there are three solutions, two complex and one real. What is the real solution? $${\left(\frac{1}{18} \, \sqrt{23} \sqrt{3} - \frac{1}{2}\right)}^{\frac{1}{3}} + \frac{1}{3 \, {\left(\frac{1}{18} \, \sqrt{23} \sqrt{3} - \frac{1}{2}\right)}^{\frac{1}{3}}} \approx -1.3247$$So let's check what integral solutions there are between -1 and 100. Figure 2 shows the points that satisfy (permalink):

Figure 2

Extending the range to one million does not throw up any new points so these seem to be the only points that are integers. I have a book that deals with this very subject that is described as follows:

This book presents in a unified and concrete way the beautiful and deep mathematics - both theoretical and computational - on which the explicit solution of an elliptic Diophantine equation is based. It collects numerous results and methods that are scattered in the literature. Some results are hidden behind a number of routines in software packages, like Magma and Maple; professional mathematicians very often use these routines just as a black-box, having little idea about the mathematical treasure behind them. Almost 20 years have passed since the first publications on the explicit solution of elliptic Diophantine equations with the use of elliptic logarithms. The "art" of solving this type of equation has now reached its full maturity. The author is one of the main persons that contributed to the development of this art. The monograph presents a well-balanced combination of a variety of theoretical tools (from Diophantine geometry, algebraic number theory, theory of linear forms in logarithms of various forms - real/complex and p-adic elliptic - and classical complex analysis), clever computational methods and techniques (LLL algorithm and de Weger's reduction technique, AGM algorithm, Zagier's technique for computing elliptic integrals), ready-to-use computer packages. A result is the solution in practice of a large general class of Diophantine equations.

This is rather too detailed for me but the mention of Magma and Maple make me wonder what tools SageMath has to offer in terms of the solution of elliptic Diophantine equations. The book was published in August of 2013. As it turns out SageMath offers two commands. For the equation \(y^2=x^3+ax+b\), the commands are:

E=EllipticCurve([\(a,b\)])

E.integral_points()

For the case of \(y^2=x^3-x+1\), this produces the output (\(a\)=-1 and \(b\)=1):

[(-1 : 1 : 1), (0 : 1 : 1), (1 : 1 : 1), (3 : 5 : 1), (5 : 11 : 1), (56 : 419 : 1)]

The output is a little confusing but it would seem that the final 1 in each triplet signifies that the middle number can be positive or negative. It is thus equivalent to:

[(-1, ±1), ((0, ±1)), (1, ±1), (3, ±5), (5, ±11), (56, ±419)]

This is what is got earlier using my longer method. Elliptic curves is a big topic and I've only scratched the surface here. Figure 3 shows the points on the curve except for (56, ±419).

Figure 3: created using GeoGebra Classic

Sphenic Numbers Revisited

After this post, I discovered that I'd already made an earlier post about sphenic numbers. No matter but it alerted me to the fact that I've made so many posts to this mathematics blog that I'm losing track of what I've posted.

Today I turned 25474 days old. This number factors to 2 * 47 * 271. Yesterday's number, 25473, factors to 3 * 7 * 1213. Both are sphenic numbers, described by Numbers Aplenty as follows:

A number \(n\) is called sphenic if it is the product of 3 distinct primes. For example, 370 is a sphenic number because it is the product of the 3 primes 2, 5 and 37. Sphenic numbers are quite common: up to \(10^8\) there are 20710806 sphenic numbers (that's about 20%). 

The sum of the reciprocals of the sphenic numbers diverges, while the sum of the reciprocal of their squares converges to \(0.003696244...\), which can be expressed as: $$ \frac{(P(2)^3-3 \, P(2) \, P(4)+2 \,P(6))}{6}\\ \text { where }P(s)=\sum_{p\mathrm{\ prime}}\frac{1}{p^s}$$ is the so-called prime Zeta function. 

The first sphenic numbers are 30, 42, 66, 70, 78, 102, 105, 110, 114, 130, 138, 154, 165, 170, 174, 182, 186, 190, 195, 222, 230, 231, 238, 246, 255, 258, 266, 273, 282, 285, 286, 290, 310

Wikipedia adds that:

All sphenic numbers have exactly eight divisors. If we express the sphenic number as \( n = p \cdot q \cdot r\) where \(p\), \(q\), and \(r\) are distinct primes, then the set of divisors of \(n\) will be \({1, p, q, r, pq, pr, qr, n}\). The converse does not hold. For example, 24 is not a sphenic number, but it has exactly eight divisors. 

All sphenic numbers are by definition squarefree, because the prime factors must be distinct. 

The first case of two consecutive sphenic integers is 230 = 2×5×23 and 231 = 3×7×11. The first case of three is 1309 = 7×11×17, 1310 = 2×5×131, and 1311 = 3×19×23. There is no case of more than three, because every fourth consecutive positive integer is divisible by 4 = 2×2 and therefore not squarefree. 

The numbers 2013 (3×11×61), 2014 (2×19×53), and 2015 (5×13×31) are all sphenic. It's interesting that these very recent calendar years formed a sphenic triplet, although I didn't know it at the time I was living through them. The next three consecutive sphenic years will be 2665 (5×13×41), 2666 (2×31×43) and 2667 (3×7×127) (see OEIS A248202 for a list of the central number of such triples). 

Sphenic Brick

In terms of geometry, each sphenic number can be considered to represent the volume of a unique and "primitive" rectangular prism (sometimes called a sphenic brick) whose dimensions are given by its three prime factors. I'm using primitive here in the same sense as "primitive Pythagorean triad" such as 3, 4 and 5 (as opposed to 6, 8 and 10). By its definition however, a sphenic number can never represent the volume of a cube or a rectangular prism with a square cross-section.

Each sphenic number \( n = p \cdot q \cdot r\) can be associated with another number, namely the surface area of the rectangular prism  \( 2 \, (p \cdot q + p \cdot r+q \cdot r) \). For example, the sphenic number  \( 7429 = 17 \cdot 19 \cdot 23\) can be viewed as a rectangular prism with an associated surface area of \(2302\) square units. The ratio between area and volume can then be explored. The table below shows the values of such ratios for sphenic numbers between 25400 and 25500:

It is possible for the volume and surface area to be equal. In a range of numbers between 1 and 1000, the only such dimensions that produce this are:

• 3 x 7 x 42   --> 882 
• 3 x 8 x 24   --> 576
• 3 x 9 x 18   --> 486
• 3 x 10 x 15 --> 450
• 4 x 5 x 20   --> 400
• 4 x 6 x 12   --> 288

Whether these are the only values with this property I don't know but none of the above numbers (288, 400, 450, 486, 576 and 882) are sphenic so it's likely that there are no sphenic numbers with this property.

To determine the sphenic numbers within a given range, this SageMath code (link to SageMathCell server) can be used or the box below (sometimes temperamental) can be experimented with:

Prime Producing Linear Polynomials

An example of a famous prime producing polynomial is \(n^2-n+41 \) which produces primes for values of \(n\) from 1 to 40. These primes are shown below:

41, 43, 47, 53, 61, 71, 83, 97, 113, 131, 151, 173, 197, 223, 251, 281, 313, 347, 383, 421, 461, 503, 547, 593, 641, 691, 743, 797, 853, 911, 971, 1033, 1097, 1163, 1231, 1301, 1373, 1447, 1523, 1601

Figure 1 is an excerpt relating to this polynomial from an interesting article on prime-generating polynomials:

Figure 1

Figure 2 shows a table of record prime producing polynomials taken from MathWorld:

Figure 2

However, prime producing linear polynomials get less attention than the above quadratic, cubic and higher order polynomials. Today I turned 25410 days old and discovered that the linear polynomial \(25410 \, n+1 \) has some interesting prime producing qualities. Specifically, it produces primes for values of \(n\) from 1 to 6: 25411, 50821, 76231, 101641, 127051, 152461. As it turns out, 25410 is the first coefficient of \(\ a*n+1\) to produce six primes in a row. The next such number is 26040.

OEIS A237190 displays a list of coefficients of of \(\ a*n+1\) that produce five primes in a row:

10830, 25410, 26040, 88740, 165900, 196560, 211050, 224400, 230280, 247710, 268500, 268920, 375480, 377490, 420330, 451410, 494340, 512820, 592620, 604170, 735750, 751290, 765780, 799170, 808080, 952680, 975660, 1053690, 1064190, 1132860, 1156170, 1532370, 1559580

As can be seen, 10830 is the first coefficient with this property. The OEIS entry also provides a list of the first 1000 such coefficients and this is a good starting point when trying to find coefficients that produce six, seven, eight or more primes in a row. With a little manipulation in a spreadsheet, the OEIS data can be pasted directly into a list set up on SageMathCell.

Of the first 1000 coefficients that produce five primes, only the following go on to produce seven primes in a row (in the range 1 to 247289070): 

512820, 8224860, 22240680, 24462900, 26486460, 62871480, 93784530, 99597960, 139819680, 196474950

Of these only the following go on to produce eight primes in a row: 

512820, 22240680, 26486460, 99597960

All four of the above fail to produce a ninth prime. Thus 512820 is the smallest coefficient that will produce eight primes in a row. In terms of primes produced over a given range and not sequentially, it would seem that 15213870 is the most prolific. It produces 40 primes between \(n+1 \) and \(100 \,n+1\) whereas 512820 produces only 33 in the same range.

What do these numbers have in common. To begin with they are all divisible by 30 and are highly factorable. The following table shows the factorisation and number of divisors for the numbers in OEIS A237190:

Figure 3

It would seem that by adding 1 to multiples of such highly factorisable numbers increases the odds of generating a prime number. Interestingly, for the 1000 coefficients mentioned earlier, the polynomial \(25410 \, n+41\) produces a record 41 primes in the range from 1 to 100, more than for any other coefficient. This naturally leads to an investigation of polynomials of the form \(25410 \, n+p\) where \(p\) is prime. Here is a table of what showed up for primes from 2 to 97:

Figure 4

Clearly \(25410 \, n+97\) wins that race and in fact holds the record for all primes below 1000. For primes up to 100000, \(25410 \, n+5471\) takes the lead with 48 primes and only surrenders it to \(25410 \, n+982231\) with 49 primes when considering primes up to one million. Below are the 49 primes generated by \(25410 \, n+982231\) as n ranges from 1 to 100:

Figure  5

Figure 6

I must say that SageMath makes it incredibly easy to investigate matters such as this. Interestingly, of all the possible coefficients from 1 to 999,999, it is \(25410 \, n+982231\) that produces the highest score of 49 primes out of 100. I was beginning to start to think that this might be something of a record until I went back to \(\ a*n+41\) and thought I'd test all the numbers up to 999,999. The table in FIGURE 6 shows the surprising result. 

The factorisations again show the same pattern as with the earlier numbers:

• 6 = 2 * 3
• 12 = 2 * 2 * 3
• 30 = 2 * 3 * 5
• 60 = 2 * 2 * 3 *5
• 210 = 2 * 3 * 5 * 7

Testing with other primes between 2 and 9999 and coefficients from 1 to 9999 (trying upper bounds of 999999 seemed to overwhelm the SAGE server) showed the best result was with the unpretentious \( 6 \,n+5 \) that produces a total of 56 primes between \(n=1\) and \(n=100\).

However, this 56% hit rate for primes when n is between 1 and 100 is not sustained. For example, setting lower and upper bounds of 4000 and 5000 for \( 6 \,n+5 \) generates 289 primes between 24005 and 30005. This is an average of 289/1000 or about 28.9%.

on August 30th 2021

See related post titled Linear Prime Chains uploaded on December 21st 2022.

Self Numbers and Junction Numbers

Today I turned 25407 days old and discovered from Numbers Aplenty that it is a self number or Columbian number. Such a number is defined as one that cannot be written as \(x+\mathrm{sod}(x) \) where \(\mathrm{sod}(\cdot) \) denotes the sum of digits. One is reminded immediately of an inconsummate number that cannot be written as \(x ÷ \mathrm{sod}(x) \) where \(\mathrm{sod}(\cdot) \) denotes the sum of digits. View my blog post about these latter numbers here.

Self numbers are base dependent. From Wikipedia:

In general, for even bases, all odd numbers below the base number are self numbers, since any number below such an odd number would have to also be a 1-digit number which when added to its digit would result in an even number. For odd bases, all odd numbers are self numbers. The set of self numbers in a given base q is infinite and has a positive asymptotic density: when q is odd, this density is 1/2.

D.R.Kaprekar, who introduced these numbers in 1949, gave the following test to determine quickly if a number is self.

Define the functions
\[
\begin{array}{lcl}
d(n)& = & \mathrm{the\ number\ of\ digits\ of\ }n\,,\\[1mm]
r(n)\vphantom{M^{M^M}a_{M_M}}&=& 1+({{(n-1)}\pmod9})\,, \\[1mm]
h(n)&=&\left\{%
\begin{array}{ll}
r(n)/2&\mathrm{\ if\ }r(n)\mathrm{\ is\ even}\,,\\
(r(n)+9)/2&\mathrm{\ if\ }r(n)\mathrm{\ is\ odd}\,.\\
\end{array}\right.\\
\end{array}
\]Then a number \(n\)  is self if, for every \(k=0,1,\dots,d(n)\), it holds
\[
\mathrm{sod}(|n-h(n)-9\cdot k|) \neq h(n)+9\cdot k\,
\]where \(|x|\)  denotes the absolute value of \(x\).

Kaprekar's algorithm of course makes it very easy to write the necessary SAGE code to identify self numbers within a given range of numbers. Here is the code that I developed:

If the SAGE code runs in the SAGE CELL window (sometimes there are glitches), you'll see the output 25407 25418 25420 25431 25442 and so 25418 marks the next occurrence of a self number. The first self numbers (OEIS A003052) are:

1, 3, 5, 7, 9, 20, 31, 42, 53, 64, 75, 86, 97, 108, 110, 121, 132, 143, 154, 165, 176, 187, 198, 209, 211, 222, 233, 244, 255, 266, 277, 288, 299, 310, 312, 323, 334, 345, 356, 367, 378, 389, 400, 411, 413, 424, 435, 446, 457, 468, 479, 490, 501, 512, 514, 525

A self prime is a self number that is prime (OEIS A006378). The first few self primes are:

3, 5, 7, 31, 53, 97, 211, 233, 277, 367, 389, 457, 479, 547, 569, 613, 659, 727, 839, 883, 929, 1021, 1087, 1109, 1223, 1289, 1447, 1559, 1627, 1693, 1783, 1873, ... 

Again from Wikipedia: "In October 2006 Luke Pebody demonstrated that the largest known Mersenne prime that is at the same time a self number is \(2^{24036583}−1 \). This is then the largest known self prime as of 2006."

By contrast, a junction number is one that can be written in more than one way as  \(x+\mathrm{sod}(x) \). For example, \(818 \)  is a junction number because it has two generators, \(796 \) and \(805 \). Indeed \[796+7+9+6 = 805+8+0+5 = 818.\]The first junction numbers are:

101, 103, 105, 107, 109, 111, 113, 115, 117, 202, 204, 206, 208, 210, 212, 214, 216, 218, 303, 305, 307, 309, 311, 313, 315, 317, 319, 404, 406, 408, 410, 412, 414, 416, 418, 420, 505, 507

It's easy enough to modify the previous SAGE code so that a list of junction numbers within a given range is generated:

The output generated from this code should be 101 103 105 107 109 111 113 115 117 together with the two smaller numbers that produce them. The junction numbers occur in runs of nine odd numbers and this is clearly seen in Figure 1 that shows the graph of such numbers between 1 and 1000.

Figure 1

Kaprekar's original paper can be viewed here.

Lotto Simulations

I've written about Lotto before in this post (Oz Lotto) and another (Losing at Lotto). The current post is not about probabilities but simulation. How can we use SageMath to simulate a lottery draw? Well as it turns out it's not too difficult.

The key is to use the ZZ.random_element(lowerbound, upperbound+1) function as in the following code. The trick is not to allow duplicate numbers and that is the reason for the inclusion of Set. This function takes our lotto list which might look like this the list [20, 37, 6, 37, 1, 19, 16] and turns it, after the sort, into the ordered set {1, 6, 16, 19, 20, 37}.

Here is the code I developed. I'm really just a tyro at both Python and SageMath so there may well be much better ways to achieve my goals but this is what I came up with:

# Generate a set of numbers for lotto games
# let N represent how many numbers are required in the game
# let G represent the number of games
# let U represent the upper bound of numbers
G=8
N=7
U=40
for i in range(G):
    lotto=[]
    while len(Set(lotto))<N:
        lotto.append(ZZ.random_element(1,U+1))
    print sorted(Set(lotto))

The following is an example of the output:

[11, 23, 27, 30, 37, 38, 40]
[4, 7, 8, 15, 20, 29, 38]
[8, 9, 13, 14, 21, 28, 36]
[7, 11, 15, 21, 28, 33, 40]
[1, 15, 17, 18, 28, 33, 39]
[4, 7, 11, 14, 18, 29, 32]
[1, 10, 12, 16, 21, 31, 33]
[2, 7, 8, 12, 13, 19, 27]

How really random this output is I don't know. I'm sure SageMath has other ways to generate more truly random numbers but this is sufficient for the modest purposes of the exercise.

The code has the added advantage of being able to accommodate "fixed" numbers. For example, let's say we always want the numbers 23 and 32 included in the output. All that needs to be done is to start with lotto = [23, 32] instead of it being an empty list.

For example:

# Generate a set of numbers for lotto games
# let N represent how many numbers are required in the game
# let G represent the number of games
# let U represent the upper bound of numbers
G=8
N=7
U=40
for i in range(G):
    lotto=[23, 32]
    while len(Set(lotto))<N:
        lotto.append(ZZ.random_element(1,U+1))
    print sorted(Set(lotto))

Here an example of the output:

[10, 11, 21, 23, 32, 35, 37]
[5, 6, 17, 23, 31, 32, 39]
[6, 11, 17, 23, 31, 32, 39]
[3, 15, 18, 23, 26, 30, 32]
[2, 12, 13, 23, 27, 28, 32]
[7, 18, 23, 26, 32, 35, 38]
[2, 11, 17, 19, 23, 32, 33]
[3, 16, 23, 26, 27, 31, 32] 

Not that I expect it will make any difference to my chances of winning but I'll test this out in some coming Lotto games. Oddly enough what got me thinking about generating some random lotto numbers was this article titled JERRY AND MARGE GO LARGE in the Huffington Post. It's about an American couple who exploited a feature in their state's lottery to make millions of dollars. It's quite an interesting story.

As I've been doing I'll include a SageMathCell with the code so that anyone curious can experiment:

Semiprime Chains

Today I turned 25399 days old. This number is a semiprime, meaning that it has two prime factors. In the case of 25399, the prime factors are 11 and 2309. However, the number is a member of OEIS A177216: numbers n that are the product of two distinct primes such that \(2n-1, 4n-3, 8n-7, 16n-15, 32n-31, 64n-63 \text{ and } 128n-127\) are also products of two distinct primes. Below is a screenshot of the mathematical tweet to celebrate the occasion.

Evaluating this leads to the following table:

The chain is thus of the form \(2^a \times n- (2^a-1) \) with the integer \(a \geq 0 \).

The initial members of OEIS A177216 are:

11293, 12139, 25399, 31261, 36199, 44869, 49471, 62521, 72397, 83086, 89737, 91705, 98941, 124846, 125041, 134023, 138994, 144793, 164041, 166171, 170431, 173311, 182527, 199543, 224962, 244294, 258169, 259891, 263086, 275281, 277987 

I was interested to see how many of these terms would survive an extension of the process to \(256 \times n - 255 \), \(512 \times n - 511 \), \(1024 \times n - 1023 \) etc.

Extending the process in the same range to \(256 \times n - 255 \) leads to a considerable thinning of the numbers:

31261 36199 44869 49471 62521 72397 83086 138994 173311 182527 224962 259891 277987

Extending the range to \(512 \times n - 511 \) gives:

31261 36199 138994 173311 259891

Extending the range to \(1024 \times n - 1023 \) gives: 173311

So 173311 is the last man standing. Will it stand any further? As it turns out it will. It holds up to \(2048 \times n - 2047 \) but, at \(4096 \times n - 4095 \), there are three prime factors (419, 601 and 2819) and not two. For \(4096 \times n - 4095 \) to hold, we need to extend the range and this leads to 2212801). Below I've listed the semiprime chain associated with the number 173311.

173311 also has the property that is formed by a concatenation of two primes 173 and 311.

The following SAGE code shows that 2212801 is the only number up to 10 million with a prime chain up to \(4096 \times n - 4095 \):

INPUT

for n in range(1,10000000):

    if len(list(divisors(n)))==4:

        if len(list(divisors(2*n-1)))==4:

            if len(list(divisors(4*n-3)))==4:

                if len(list(divisors(8*n-7)))==4:

                    if len(list(divisors(16*n-15)))==4:

                        if len(list(divisors(32*n-31)))==4:

                            if len(list(divisors(64*n-63)))==4:

                                if len(list(divisors(128*n-127)))==4:

                                    if len(list(divisors(256*n-255)))==4:

                                        if len(list(divisors(512*n-511)))==4:

                                            if len(list(divisors(1024*n-1023)))==4:

                                                if len(list(divisors(2048*n-2047)))==4:

                                                    if len(list(divisors(4096*n-4095)))==4:

                                                        print(n),

OUTPUT

2212801

You can try running the above SAGE code using the evaluate button below the SAGE cell. Sometimes it's a little slow but generally it should work. You can modify the code by adding a # in front of a line e.g. # if len(list(divisors(4096*n-4095)))==4: will remove the 4096 condition and apply only the 2048 and lower conditions. Try it out. You may need to refresh the page after each evaluation.

Rhonda Numbers

Today I turned 25398 days old and I made Rhonda numbers the focus of my mathematical tweet of the day:

In base 8, 25398 is represented as 61466 and the product of these base 8 digits is 864. The sum of the prime factors (including multiplicity) of 25398 is 2+3*2+17+83 = 108 and 108 * 8 = 864. It's only the number itself that needs to be changed into a base other than 10. The prime factors and their sum all remain in base 10. So 25398 is a Rhonda number to base 8.

OEIS A100970 lists the initial (base 10) numbers that are Rhonda numbers in base 8 and they are:

1836, 6318, 6622, 10530, 14500, 14739, 17655, 18550, 25398, ...

Numbers Aplenty has this to say about the Rhonda numbers in general:

Rhonda numbers exist only in composite bases. Indeed, the product of the digits of a number in a prime base \(b\) cannot be divisible by \(b\), since every digit is smaller than \(b\). 

Kevin Brown has proved that there are infinite Rhonda numbers. 

The first base-10 Rhonda numbers are 1568, 2835, 4752, 5265, 5439, 5664, 5824, 5832, 8526, 12985, 15625, 15698, 19435, 25284. 

1000 is the smallest Rhonda number in two bases, namely 16 and 36, since: 

\(1000=2^3\cdot 5^3\), \(1000=(3,14,8)_{16}=(27,28)_{36}\) 

\(3\cdot14\cdot8 = 16\cdot(2\cdot3+3\cdot5)\) for base 16   

\(27\cdot28=36\cdot(2\cdot3+3\cdot5)\) for base 36. 

The first smallest Rhonda numbers with respect to one, two, ... , ten bases are 560, 1000, 10200, 5670, 63945, 158400, 322920, 140800, 1200420, 889200.

From MathWorld, the first few Rhonda numbers for small composite bases b are summarised in the following table:

	OEIS	Rhonda numbers
4	A100968	10206, 11935, 12150, 16031, 45030, 94185, ...
6	A100969	855, 1029, 3813, 5577, 7040, 7304, 15104, 19136, ...
8	A100970	1836, 6318, 6622, 10530, 14500, 14739, 17655, 18550, 25398, ...
9	A100973	15540, 21054, 25331, 44360, 44660, 44733, 47652, ...
10	A099542	1568, 2835, 4752, 5265, 5439, 5664, 5824, 5832, 8526, 12985, ...
12	A100971	560, 800, 3993, 4425, 4602, 4888, 7315, 8296, 9315, 11849, 12028, ...
14	A100972	11475, 18655, 20565, 29631, 31725, 45387, 58404, 58667, 59950, ...
15	A100974	2392, 2472, 11468, 15873, 17424, 18126, 19152, 20079, 24388, ...
16	A100975	1000, 1134, 6776, 15912, 19624, 20043, 20355, 23946, 26296, ...

Here is the SageMath code that I developed to generate the sequence in bases from 2 to 36 (with base 16 being used for the example):

b=16 #enter required base

L=[]

for x in range(2,27000):

    d=Integer(x).str(base=b)

    product=1

    for i in range(len(d)):

        product=product*int(d[i], base=b)

    P=list(factor(x))

    sum=0

    for p in P:

        sum+=p[0]*p[1]

    if product/sum ==b and product !=0:

        L.append(x)

print(L)

[1000, 1134, 6776, 15912, 19624, 20043, 20355, 23946, 26296] 

This algorithm will work for all bases from 2 up to 36 but there are no Rhonda numbers in prime bases, so there will be no output for these bases. 

Let's return to what was said in the Numbers Aplenty quote above, viz.:

The first smallest Rhonda numbers that can be written one base only, two bases, three bases etc. are: 

• one --> 560
• two --> 1000
• three --> 10200
• four --> 5670
• five --> 63945
• six --> 158400
• seven --> 322920
• eight --> 140800
• nine --> 1200420
• ten --> 889200

This site gives 140800 are an example. This is the smallest number that is a Rhonda number in eight different bases. See Figure 1.

Figure 1

Let's use base 198 to illustrate how this works. In base 198, 140800 can be represented as \( (3, 117, 22)_{198} \) and the product of these digits (3, 117 and 22) is 7722. Now 140800 factorises to 2^9 * 5^2 * 11 and the sum of these prime factors is 39 (remember we are working in base 10 here). If we multiply 39 by 198, we get 7722. Thus 140800 is a Rhonda number in base 198 and so it can be shown for the other bases.

on Thursday, April 1st 2021

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.

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

Practical Numbers

Today I turned 25320 days old but I accidentally entered 25230 into the OEIS and discovered that it was a practical number, specifically one that formed the central member of a triple of practical numbers. It is a number \(n\) such that \(n-2\), \(n\), \(n+2\) are all practical numbers (OEIS A287682). So in this case, the other members of the triple are 25228 and 25232. This is not a common occurrence as can be seen by the initial members of the sequence:

4, 6, 18, 30, 198, 306, 462, 702, 1482, 2550, 3330, 4422, 5778, 6102, 6498, 9042, 11178, 11778, 14418, 15498, 17298, 17442, 19458, 20862, 21582, 22878, 23322, 23550, 25230, ...

This led me to investigate what characterised a practical number and during that process I realised that I'd wrongly entered my number of the day. However, as it turns out 24320 is also a practical number but not a member of a triplet (the next entry in OEIS A287682 is 26622). I was encouraged to continue my investigations. According to Wikipedia:

In number theory, a practical number or panarithmic number is a positive integer \(n\) such that all smaller positive integers can be represented as sums of distinct divisors of \(n\). For example, 12 is a practical number because all the numbers from 1 to 11 can be expressed as sums of its divisors 1, 2, 3, 4, and 6: as well as these divisors themselves, we have 5 = 3 + 2, 7 = 6 + 1, 8 = 6 + 2, 9 = 6 + 3, 10 = 6 + 3 + 1, and 11 = 6 + 3 + 2.

The practical numbers themselves are rather frequent. OEIS A005153 lists these initial practical numbers:

1, 2, 4, 6, 8, 12, 16, 18, 20, 24, 28, 30, 32, 36, 40, 42, 48, 54, 56, 60, 64, 66, 72, 78, 80, 84, 88, 90, 96, 100, 104, 108, 112, 120, 126, 128, 132, 140, 144, 150, 156, 160, 162, 168, 176, 180, 192, 196, 198, 200, 204, 208, 210, 216, 220, 224, 228, 234, 240, 252, ...

As can be seen, all except 1 are multiples of 2 and this is one condition for a number being practical. The fully rigorous statement of what determines a practical number would be:

A positive integer greater than one with prime factorisation

\(n=p_1^{\alpha_1} \dots p_k^{\alpha_k}\) (with the primes in sorted order) is practical if and only if each of its prime factors \(p_i \) is small enough for \(p_i-1 \) to have a representation as a sum of smaller divisors. For this to be true, the first prime \(p_1 \) must equal 2 and, for every \(i \) from 2 to \( k \), each successive prime \(p_i \) must obey the inequality:$$ p_i\leq1+\sigma(p_1^{\alpha_1}p_2^{\alpha_2}\dots p_{i-1}^{\alpha_{i-1}})=1+\prod_{j=1}^{i-1}\frac{p_j^{\alpha_j+1}-1}{p_j-1}$$where \( \sigma(x) \) denotes the sum of the divisors of \(x\). For example, 2 × 3^2 × 29 × 823 = 429606 is practical, because the inequality above holds for each of its prime factors: 3 ≤ \(\ \sigma \)(2) + 1 = 4, 29 ≤ \( \sigma \)(2 × 3^2) + 1 = 40, and 823 ≤ \(\sigma \)(2 × 3^2 × 29) + 1 = 1171.

on August 1st 2021

With some difficulty I developed some SAGE code to calculate the practical numbers contained in a given range. The example at the end of this post is for the range from 26400 to 26500. It includes multiples of 2 because all powers of 2 are practical numbers along with all perfect numbers and primorials. Such numbers are not merely of interest to recreational mathematicians. They are of interest to professional mathematicians because many of their properties are similar to the properties of the prime numbers.