Harshad Numbers Revisited

On the 10th February 2017, I posted about Harshad numbers, also known as Niven numbers, that are defined to be those numbers that are divisible by the sum of their digits in a given base. The post featured 24786 which was the last of the right-truncatable Harshad numbers (in base 10). Since then I've not been reminded of such numbers until I only recently came across the Numbers Aplenty website (see recent post). The reason for the disappearance of the Harshad numbers from the OEIS is that they are relatively frequent and so do not appear in a search of its database. Only the following numbers are listed:

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 18, 20, 21, 24, 27, 30, 36, 40, 42, 45, 48, 50, 54, 60, 63, 70, 72, 80, 81, 84, 90, 100, 102, 108, 110, 111, 112, 114, 117, 120, 126, 132, 133, 135, 140, 144, 150, 152, 153, 156, 162, 171, 180, 190, 192, 195, 198, 200, 201, 204

However, the OEIS entry does contain a link to all the Harshad numbers less or equal to 100,000. There are 11,872 of them and so there density is slightly less than 12% in that range. My number for today, 25290, appears in that list in position number 3771. However, I was only aware of this because of the Numbers Aplenty entry for 25290. The OEIS is definitely limited when it comes to searching for properties of numbers in this region of the number line and beyond.

The Numbers Aplenty reference to Harshad numbers contains some interesting additional information to what I posted in my earlier blog. For example, it is possible to have at most twenty consecutive Niven numbers and this happens infinitely often. Below is a list of runs of up to 13:

Wikipedia states that more generally there are 2b but not 2b + 1 consecutive b-harshad numbers, so that for example there can be runs of up to 22 11-harshad numbers, 24 12-harshad numbers and so on. Harshad numbers fall into the category of recreational mathematics.

While it's rather easy to determine whether a number is Harshad or not, nonetheless here is the SAGE code that I developed to test "harshadness":

INPUT
#Test for Harshad Number, divisible by sum of digits
number=25290
sum=0
nu_str=str(number)
for d in nu_str:
    sum+=Integer(d)
if number % sum == 0:
    print(number, "is Harshad Number")
else:
    print(number, "is not a Harshad Number") 

OUTPUT
(25290, 'is Harshad Number')

Happy Numbers

Today I turned 25286 days old and 25286 happens to be a happy number, defined as:

Let us define a function \( s(n) \), for \( n>0 \), which gives the sum of the squares of the digits of \( n \), so, for example, \( s(37)=3^2+7^2=58 \).

If we start from a number \(n\)  and we repeatedly apply \( s(\cdot) \), we obtain a sequence \(S_n\) of numbers \(n\), \(s(n) \), \(s(s(n)\), \(\dots \), and so on.

A number \(n\) is called happy if \(S_n \) contains the number 1.

Note that \(s(1)=1 \), so in that case the sequence \( S_n \) has an infinite tail of \(1\)'s.

If a number is not happy then it is easy to see that at a certain point \(S_n\) will enter the infinite loop$$ \dots,4, 16, 37, 58, 89, 145, 42, 20, 4,\dots $$So, for example, starting from  \(94\)  we obtain \(94\rightarrow97\rightarrow130\rightarrow10\rightarrow1\), so \(94\) is happy. 

On the contrary, starting from 61 we obtain \(61\rightarrow37\rightarrow58\rightarrow89 \) and thus 61 is not happy, since 89 belongs to the unhappy loop.

Source: http://www.numbersaplenty.com/set/happy_number/

According to Wikipedia: 

by inspection of the first million or so happy numbers, it appears they have a natural density of around 0.15. Perhaps surprisingly, then, the happy numbers do not have an asymptotic density. The upper density of the happy numbers is greater than 0.18577, and the lower density is less than 0.1138. After 25286, the next happy number is 25294 (followed by 25295).

As usual I tried to write some SAGE code to determine whether a number was happy or not. Here is what I came up with:

INPUT (using a happy number 25286):
entered_number=25286
number=str(entered_number)
while sum!= 1 and sum!=4:
    sum=0
    for x in range(len(number)):
        sum+=Integer(number[x])^2
    number=str(sum)
    print(sum) 

OUTPUT
133
19
82
68
100
1

INPUT (using an unhappy number 89)
entered_number=89
number=str(entered_number)
while sum!= 1 and sum!=4:
    sum=0
    for x in range(len(number)):
        sum+=Integer(number[x])^2
    number=str(sum)
    print(sum) 

OUTPUT
145
42
20
4

The 4 is used to prevent the program going into an infinite loop and, at the same time, to mark a number as happy. Unhappy numbers will always enter the following loop:

... 4, 16, 37, 58, 89, 145, 42, 20, 4, ...

Any of the above numbers could be used to identify an unhappy number and break the loop.

The Wikipedia article on unhappy numbers includes the following Python code for determining whether a number is happy or unhappy (rather than displaying the trajectory as I did in my SAGE program):

INPUT (using happy number 25286)
def square(x):
    return int(x) * int(x)
def happy(number):
    return sum(map(square, list(str(number))))
def is_happy(number):
    seen_numbers = set()
    while number > 1 and (number not in seen_numbers):
        seen_numbers.add(number)
        number = happy(number)
    return number == 1
is_happy(25286) 

OUTPUT
True

INPUT (using unhappy number 89)
def square(x):
    return int(x) * int(x)
def happy(number):
    return sum(map(square, list(str(number))))
def is_happy(number):
    seen_numbers = set()
    while number > 1 and (number not in seen_numbers):
        seen_numbers.add(number)
        number = happy(number)
    return number == 1
is_happy(89) 

OUTPUT
False 

As can be seen, the Python code involves a quite different approach. It defines three functions, the second building on the first and the third building on the second : square(x), happy(number) and is_happy(number).

The command sum(map(square, list(str(number)))) caught my eye. This is a very useful command that takes two inputs:

• a function (in this case, square)
• a list (in this case, str(number))

It outputs a list (in this case, the sum of the number's digits squared).

Sphenic Numbers

Today I turned 25285 days old and, as I discovered in Numbers Aplenty, 25285 is a sphenic number. The definition given on that site is:

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 20,710,806 sphenic numbers. 

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  \((P(2)^3-3\,P(2)\,P(4)+2\,P(6))/6\), 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

The sum of the reciprocals of the squares of the sphenic numbers does indeed approach 0.003696244 as can be seen by taking the numbers from 30 to 310 and applying the following SAGE code:

INPUT: 

sphenic=[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]
sum=0
for n in sphenic:
    sum+=1/n^2
print(sum.n())

OUTPUT: 

0.00320320889263633 

The following is a modification of the Wikipedia entry for sphenic numbers:

In number theory, a sphenic number is a natural number|positive integer that is the product of three distinct prime numbers. Thus a sphenic number is a product ''pqr'' where ''p'', ''q'', and ''r'' are three distinct prime numbers. This definition is more stringent than simply requiring the integer to have exactly three prime factors. For instance, \(60 = 2^2 × 3 × 5 \) has exactly 3 prime factors, but is not sphenic. 

The smallest sphenic number is 30 = 2 × 3 × 5, the product of the smallest three primes. The largest known sphenic number is
$$(2^{77232917}− 1) \cdot (2^{74,207,281} − 1) \cdot (2^{57,885,161} − 1)$$It is the product of the three largest known primes. 

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 Möbius function of any sphenic number is -1. 

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 square-free. 

The numbers 2013 (3×11×61), 2014 (2×19×53), and 2015 (5×13×31) are all sphenic. The next three consecutive sphenic years will be 2665 (5×13×41), 2666 (2×31×43) and 2667 (3×7×127) (see OEIS A165936).

The OEIS sequence A007304 (sphenic numbers, products of 3 distinct primes) mentions a that a sphenic brick is a rectangular parallelopiped whose sides are components of a sphenic number, namely whose sides are three distinct primes. Example: The distinct prime triple (3,5,7) produces a 3x5x7 unit brick which has volume 105 cubic units. 

From my early days of investigating numbers, I had considered triprimes (number that factor into three, not necessarily distinct, primes) as having unique representations as rectangular prisms. I was interested in the relationship between a prisms volume and its surface area, coining the term "cubicity". Here is a snapshot of a February 2014 Instagram post:

Numbers Aplenty

Today I just stumbled upon a new site that lists information about natural numbers. It's called Numbers Aplenty.

The categories under which numbers are classified looks like this:

So for instance, today I am 25284 days old and entering that number into the search box yields:

• 25284 has 36 divisors (see below), whose sum is σ = 70224. Its totient is φ = 7056. 

• The previous prime is 25261. The next prime is 25301. The reversal of 25284 is 48252 

• Adding to 25284 its reverse (48252), we get a triangular number (73536 = T383). 

• It is a Harshad number since it is a multiple of its sum of digits (21). 

• 25284 is a Rhonda number in base 10. 

• Its product of digits (640) is a multiple of the sum of its prime factors (64). 

• It is a nialpdrome in base 14. 

• It is a self number, because there is not a number n which added to its sum of digits gives 25284. 

• It is an unprimeable number. 

• It is a polite number, since it can be written in 11 ways as a sum of consecutive naturals, for example, 567 + ... + 609. 

• 225284 is an apocalyptic number. 

• It is an amenable number. 

• It is a practical number, because each smaller number is the sum of distinct divisors of 25284, and also a Zumkeller number, because its divisors can be partitioned in two sets with the same sum (35112). 

• 25284 is an abundant number, since it is smaller than the sum of its proper divisors (44940). 

• It is a pseudoperfect number, because it is the sum of a subset of its proper divisors. 

• 25284 is a wasteful number, since it uses less digits than its factorization. 

• 25284 is an evil number, because the sum of its binary digits is even. 

• The sum of its prime factors is 64 (or 55 counting only the distinct ones). 

• The product of its digits is 640, while the sum is 21. 

• The square root of 25284 is about 159.0094336824. The cubic root of 25284 is about 29.3504835430. 

• The spelling of 25284 in words is "twenty-five thousand, two hundred eighty-four".

Many of the terms mentioned above I'd never heard before. The site will be a useful adjunct to the Online Encyclopedia of Integer Sequences or OEIS.

Generating Lucky Numbers in Python

I got to thinking how I could generate lucky numbers in Python (for an earlier post explaining what lucky numbers are: click here). I experimented a little but wasn't making any real progress which is not surprising given that I'm really just a beginner. A quick search revealed the following exercise: Python Math: Print the first n Lucky Numbers along with the required code. What surprised me was the brevity of the code:

I've worked through the code now and understand it, although I struggled initially, especially. Let's suppose n is given a value of 10. This means:

• List = range(-1,109,2) = (-1, 1, 3, 5, ..., 107, 108)
• while List[i:] = List[2:] = (3, 5, ..., 107, 108) because initially i has a value of 2
• List[List[i]::List[i]] = List[List[2]::List[2]] = List[3, 5, ..., 107, 108]
• i+=1 means that value of i, will increment by 1, as the while loop is traversed

The while loop will continue as long as List[i:] is True but eventually the i will be greater than the length of List and List[:i] will be FALSE. In SageMathCell, things look like this:

The above code is purely Python and doesn't used any SAGE code at all. As far as I know, there's no command in the latter that will generate lucky numbers (as can be done for prime numbers). There are 82 exercises with solutions at the W3resource site. It would be instructive to attempt some of these. 

Solving Frobenius Equations and Computing Frobenius Numbers

A Frobenius equation is an equation of the form:$$a_1x_n + \ldots + a_nx_n=m$$where \( a_1, \dots, a_n \) are positive integers, \( m \) is an integer and the coordinates \(x_1, \dots, x_n \) of solutions are required to be non-negative integers.

The Frobenius number of \( a_1, \ldots, a_n \) is the largest \( m \) for which the Frobenius equation \(a_1x_n + \ldots + a_nx_n=m \) has no solutions.

These definitions were taken from the Wolfram Language & System Documentation Center Tutorial. A Frobenius equation and number are both seemingly abstract concepts but there are practical, concrete applications. In the tutorial, the following example is given: in how many ways can a total of 42 cents be created using 1, 5, 10, and 25 cent coins? The problem can be expressed as a Frobenius equation:$$x_1+5x_2+10x_3+25x_4=42$$WolframAlpha solves this using the command: FrobeniusSolve[{1, 5, 10, 25}, 42] and the following output is produced:

{{2, 0, 4, 0}, {2, 1, 1, 1}, {2, 2, 3, 0}, {2, 3, 0, 1}, {2, 4, 2, 0}, {2, 6, 1, 0}, {2, 8, 0, 0}, {7, 0, 1, 1}, {7, 1, 3, 0}, {7, 2, 0, 1}, {7, 3, 2, 0}, {7, 5, 1, 0}, {7, 7, 0, 0}, {12, 0, 3, 0}, {12, 1, 0, 1}, {12, 2, 2, 0}, {12, 4, 1, 0}, {12, 6, 0, 0}, {17, 0, 0, 1}, {17, 1, 2, 0}, {17, 3, 1, 0}, {17, 5, 0, 0}, {22, 0, 2, 0}, {22, 2, 1, 0}, {22, 4, 0, 0}, {27, 1, 1, 0}, {27, 3, 0, 0}, {32, 0, 1, 0}, {32, 2, 0, 0}, {37, 1, 0, 0}, {42, 0, 0, 0}}

The Frobenius number is found using FrobeniusNumber[{1, 5, 10, 25}] which yields in this case -1. This is presumably WolframAlpha's way of saying that for this situation there is no such number because one of the building blocks is 1 and that can be used to build any positive integer. Let's look at another Frobenius equation:$$12x_1+16x_2+20x_3+27x_4=123$$Here FrobeniusSolve[{12, 16, 20, 27}, 123] yields:

{{0, 1, 4, 1}, {0, 6, 0, 1}, {1, 4, 1, 1}, {2, 2, 2, 1}, {3, 0, 3, 1}, {4, 3, 0, 1}, {5, 1, 1, 1}, {8, 0, 0, 1}}

The Frobenius number in this case (found using FrobeniusNumber[{12, 16, 20, 27}]) is 89. The actual output is {}. It can be noted that FrobeniusSolve[{12, 16, 20, 27},73]) also produces {} as output but 73<89 and so 73 is not the Frobenius number for the equation. What this means is that the Frobenius equation \(12x_1+16x_2+20x_3+27x_4=m \) always has a solution provided \(m>89\).

Here is a link to the MacTutor archive page about Ferdinand Georg Frobenius. I haven't yet found how to solve Frobenius equations or calculate a Frobenius number in SageMath. The only reference I've found is the following (link) which is clearly not the same:

We can compute the Frobenius coordinates and go back and forth: 

sage: Partition([7,3,1]).frobenius_coordinates()
([6, 1], [2, 0]) 

sage: Partition(frobenius_coordinates=([6,1],[2,0]))
[7, 3, 1] 

sage: all(mu == Partition(frobenius_coordinates=mu.frobenius_coordinates())
....:     for n in range(30) for mu in Partitions(n))
True