Talk:Rare numbers

From Rosetta Code

comments concerning   interesting observations   from an webpage

(The author's webpage, the last URL reference from this task's preamble, re-shown below:)


(a URL reference):

  •   author's  website:   rare numbers   by Shyam Sunder Gupta.     (lots of hints and some observations).


I was considering adding checks   (to the REXX program)   to assert that for:

  •   when the number of digits in a rare number is even,   the   sum   must be divisible by   11,     and
  •   when the number of digits in a rare number is   odd,   the   difference   must be divisible by   9.


n-r is divisible by 9 for all Rare numbers. n-r is also divisible by 99 when the number of digits is odd. (see Talk:Rare_numbers#30_mins_not_30_years)--Nigel Galloway (talk) 13:48, 21 September 2019 (UTC)

30 mins not 30 years

If Shyam Sunder Gupta has really spent 30 years on this he should have stayed in bed. Let me spend 30mins on it. The following took 9mins so any questions and I have 21mins to spare.

Let me consider n-r=l for a 2 digit number ng n<g. Then l=(10g+n)-(g+10n)=9(g-n) where n is 0..8 and g is 1..9.
l is one of 9 18 27 36 45 54 63 72 81. l must be a perfect square so only 9 36 and 81 are of interest.

9 -> ng=89 78 67 56 45 34 23 12 01
36-> ng=59 48 37 26 15 04
81-> ng=09

For each of these candidate ng I must determine if ng+gn is a perfect square.

09+90  99 n
59+95 154 n
48+84 132 n
37+73 110 n
26+62 114 n
15+51  66 n
04+40  44 n
89+98 187 n
78+87 165 n
67+76 143 n
56+65 121 y
45+54  99 n
34+43  77 n
23+32  55 n  
12+21  33 n
01+10  11 n

From which I see that 65 is the only Rare 2 digit number.

I love an odd number of digits. Let me call the 3 digit number nxg. l=(100g+10x+n)-(g+10x+100n).
x disappears and I am left with 99(g-n). None of 99 198 297 396 495 594 693 792 or 894 are perfect squares.
So there are no Rare 3 digit numbers.

At 4 I begin to think about using a computer. Consider nige. l=(1000e+100g+10i+n)-(e+10g+100i+1000n).
I need a table for l as above for 9(111(e-n)+10(g-i)) where n<=e and if n=e then i<g.

Before turning the computer on I'll add that for 5 digits nixge l=(10000e+1000g+100x+10i+n)-(e+10g+100x+1000i+10000n).
x disappears leaving 99(111(e-n)+10(g-i))
--Nigel Galloway (talk) 13:34, 12 September 2019 (UTC)
I have turned the computer on and produced a solution using only the above and nothing from the referenced website which completes in under a minute. The reference is rubbish, consider removing it--Nigel Galloway (talk) 10:42, 18 September 2019 (UTC)
Rubbish or not, is there anything on the referenced (Gupta's) website that is incorrect?   The properties and observations is what the REXX solution used (and others have as well) to calculate   rare   numbers,   albeit not as fast as your algorithm.   I have no idea how long Shyam Sunder Gupta's program(s) executed before it found eight rare numbers   (or how much virtual memory it needed).   Is the   F#   algorithm suitable in finding larger   rare   numbers?   I suspect (not knowing F#) that virtual memory may become a limitation.   Eight down, seventy-six more to go.     -- Gerard Schildberger (talk) 18:18, 18 September 2019 (UTC)
So is the task now to replicate [[1]]? This may not be reasonable for a RC task as I now discuss.
Why have you "no idea how long Shyam Sunder Gupta's program(s) executed before it found eight rare numbers"? On the webpage it says "the program has been made so powerful that all numbers up to 10^14 can be just checked for Rare numbers in less than a minute". This implies that it can search 10^13 in 5 secs. I believe this. The problem with the webpage is not that it is wrong, but that it is disingenuous. I estimate that the following is achievable:
10^13 -> 5 secs;
10^15 -> 60 secs;
10^17 -> 20 mins;
10^19 -> 7 hours;
10^21 -> 6 days;
10^23 -> 4 months.
I would say that 10^17 is reasonable for a RC task and is in line with the timings given on the webpage. Those who have recently obtained an 8th. generation i7 might want to observe that there is an obvious multithreading strategy and might want to prove that Goroutines are more than a bad pun on coroutines, as a suitable punishment for making me envious. (Warning, from my experience of i7s it might be wise to take it back to the shop and have water cooling installed before attempting to run it for a day and a half on full throttle). I remain to be convinced that these benchmarks are achieved by the Fortran, Ubasic programs on the webpage, or can be achieved in this task using the methods described on the webpage.
It is necessary to distinguish the algorithm I describe above from the F# implementation on the task page. The algorithm can be written to require very little memory. Obviously the F# as it stands calculates all candidates before checking them and this list grows with increasing number of digits. It is more than adequate for the current task, and I anticipate little difficulty in accommodating reasonable changes to make the task less trivial as layed out above--Nigel Galloway (talk) 14:01, 20 September 2019 (UTC)
Obviously, this task's requirement is not to replicate the list of 84   rare   numbers on  [Shyam Sunder Gupta's webpage: a list of 84 rare numbers].   The task requirements have not changed:   find and show the first   5   rare   numbers.   The last two requirements are optional.   Any hints and properties can be used as one sees fit.     -- Gerard Schildberger (talk) 17:55, 20 September 2019 (UTC)
Well, I have to admit that Nigel's (n-r) approach is considerably faster than what SSG has published as I've just added a second Go version which finds the first 25 Rare numbers (up to 15 digits) in about 42 seconds, albeit on my Core i7 which is kept cool when hitting turbo mode by a rather noisy fan.
Although I don't doubt SSG's claims (he is presumably a numerical expert), he must be using much more sophisticated methods than he has published to achieve those sort of times on antique hardware.
Incidentally, I regard it as bad form to use concurrent processing (i.e. goroutines) in RC tasks unless this is specifically asked for or is otherwise unavoidable (for processing events in some GUI package for example). It is difficult for languages which don't have this stuff built in to compete and it may disguise the usage of what are basically poor algorithms. --PureFox (talk) 23:04, 24 September 2019 (UTC)

A few more mins.

Above I considered n-r=L and investigated the nature of L. It is difficult to estimate complexity with increasing number of digits for this algorithm because: it depends on the number of L which happen to be a perfect square, which is neither random nor easy to predict; each L which is a perfect square expands into varying number of r. At first sight n+r=H does not look promising. For the number n....g n+g may take the values 1 to 18 rather than 1 to 9 for n-g and for an odd number of digits x doesn't disappear, rather it creates 10 times the number of possibilites.
For a rare number H-L=(n+r)-(n-r)=2r which implies that for a Rare number H and L must both be odd or even. From the lists H and R I can produce all the candidates for r. n is H-r therefore I must simply determine for each candidate is r H-r reversed (see C++ on the task page). The nice thing is that it is easy to determine the size of H and L and they predict timings for this algorithm. The time does not depend on the number of L which happen to be a perfect square. This simplicity means that I can predict optimizations which will make significant difference to the run time for increasing number of digits. The programming skill then is to compute the Cartesian Product of n/2 for L and (n+1)/2 for H integer ranges; then calculate the Inner Product of each with 10**x-10**y for L and 10**x+10**y for H. Efficiently!!! I especially like that this can be described as a problem in Linear Algebra subject to analysis and representation as a graph.--Nigel Galloway (talk) 12:19, 20 December 2019 (UTC)

the 1st REXX version

This is the 1st REXX version,   before all the optimizations were added:

/*REXX program to calculate and display an specified amount of   rare    numbers.       */
numeric digits 20;    w= digits() + digits() % 3 /*ensure enough decimal digs for calcs.*/
parse arg many start .                           /*obtain optional argument from the CL.*/
if  many=='' |  many==","  then  many= 3         /*Not specified?  Then use the default.*/
#= 0                                             /*the number of  rare  numbers (so far)*/
    do n=10                                      /*N=10, 'cause 1 dig #s are palindromic*/
    r= reverse(n)                                /*obtain the reverse of the number  N. */
    if r>n   then iterate                        /*Difference will be negative?  Skip it*/
    if n==r  then iterate                        /*Palindromic?   Then it can't be rare.*/
    s= n+r                                       /*obtain the    sum     of  N  and  R. */
    d= n-r                                       /*   "    "  difference  "  "   "   "  */
    if iSqrt(s)**2 \== s  then iterate           /*Not a perfect square?  Then skip it. */
    if iSqrt(d)**2 \== d  then iterate           /* "  "    "       "       "    "   "  */
    #= # + 1                                     /*bump the counter of  rare  numbers.  */
    say right( th(#), length(#) + 9)       ' rare number is:  '       right( commas(n), w)
    if #>=many  then leave                       /* [↑]  W:  the width of # with commas.*/
    end   /*n*/
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _;  do jc=length(_)-3  to 1  by -3; _=insert(',', _, jc); end;  return _
th:     parse arg th;return th||word('th st nd rd',1+(th//10)*(th//100%10\==1)*(th//10<4))
/*──────────────────────────────────────────────────────────────────────────────────────*/
iSqrt:  parse arg x;   $= 0;     q= 1;                            do while q<=x;    q= q*4
                                                                  end  /*while q<=x*/
                  do while q>1;  q= q % 4;   _= x-$-q;   $= $ % 2
                  if _>=0  then do;          x= _;       $= $ + q
                                end
                  end   /*while q>1*/;                            return $

Pretty simple,   but slow as molasses in January.

Not ready for prime time.

the 2nd REXX version

This is the 2nd REXX version,   after all of the hints   (properties of   rare   numbers)   within Shyam Sunder Gupta's
webpage   have been incorporated in this REXX program.

/*REXX program to calculate and display an specified amount of   rare    numbers.       */
numeric digits 20;    w= digits() + digits() % 3 /*ensure enough decimal digs for calcs.*/
parse arg many start .                           /*obtain optional argument from the CL.*/
if  many=='' |  many==","  then  many= 5         /*Not specified?  Then use the default.*/

@dr.=0;   @dr.2= 1; @dr.5=1 ; @dr.8= 1; @dr.9= 1 /*rare # must have these digital roots.*/
@ps.=0;   @ps.2= 1; @ps.3= 1; @ps.7= 1; @ps.8= 1 /*perfect squares    must end in these.*/
@end.=0;  @end.1=1; @end.4=1; @end.6=1; @end.9=1 /*rare # must not  end in these digits.*/
@dif.=0;  @dif.2=1; @dif.3=1; @dif.7=1; @dif.8=1; @dif.9=1 /* A─Q mustn't be these digs.*/
@noq.=0;  @noq.0=1; @noq.1=1; @noq.4=1; @noq.5=1; @noq.6=1; @noq.9=1 /*A=8, Q mustn't be*/
@149.=0;  @149.1=1; @149.4=1; @149.9=1           /*values for  Z  that need a even  Y.  */
#= 0                                             /*the number of  rare  numbers (so far)*/
@n05.=0;    do i= 1        to 9;  if i==0 | i==5  then iterate;  @n05.i= 1; end  /*¬1 ¬5*/
@eve.=0;    do i=-8  by 2  to 8;  @eve.i=1;  end /*define even   "    some are negative.*/
@odd.=0;    do i=-9  by 2  to 9;  @odd.i=1;  end /*   "   odd    "      "    "    "     */
                                                 /*N=10, 'cause 1 dig #s are palindromic*/
    do n=10;  parse var  n  a 2 b 3 '' -2 p +1 q /*get 1st\2nd\penultimate\last digits. */
    if @end.q  then iterate                      /*rare numbers can't end in: 1 4 6 or 9*/
    if q==3    then iterate

       select                                    /*weed some integers based on 1st digit*/
       when a==q  then do
                       if a==2|a==8 then nop     /*if A = Q,   then A must be  2  or 8. */
                                    else iterate /*A  not two or eight?       Then skip.*/
                       if b\==p  then iterate    /*B  not equal to  P?        Then skip.*/
                       end
       when a==2  then do; if q\==2 then iterate /*A = 2?     Then  Q  must also be  2. */
                           if b\==p then iterate /*" " "      Then  B  must equal    P. */
                  end
       when a==4  then do
                       if q\==0   then iterate   /*if Q not equal to zero, then skip it.*/
                       _= b - p                  /*calculate difference between B and P.*/
                       if @eve._  then iterate   /*Positive not even?      Then skip it.*/
                       end
       when a==6  then do
                       if @n05.q  then iterate   /*Q  not a zero or five?  Then skip it.*/
                       _= b - p                  /*calculate difference between B and P.*/
                       if @eve._  then iterate
                       end
       when a==8  then do
                       if @noq.q  then iterate   /*Q  isn't one of 2, 3, 7, 8?  Skip it.*/
                         select
                         when q==2  then            if b+p\==9                then iterate
                         when q==3  then do; if b>p         then if b-p\== 7  then iterate
                                               else if b<p  then if b-p\==-3  then iterate
                                                            else if b==p      then iterate
                                         end
                         when q==7  then do; if b>1         then if b+p\==11  then iterate
                                               else if b==0 then if b+p\== 1  then iterate
                                         end
                         when q==8  then            if b\==p                  then iterate
                         otherwise  nop
                         end   /*select*/
                       end                       /* [↓]  A  is an odd digit.            */
       otherwise  n= n + 10**(length(n) - 1) - 1 /*bump N so next N starts with even dig*/
                  iterate                        /*Now, go and use the next value of  N.*/
       end   /*select*/

    _= a - q;     if @dif._  then iterate        /*diff of A─Q must be: 0, 1, 4, 5, or 6*/
    r= reverse(n)                                /*obtain the reverse of the number  N. */
    if r>n   then iterate                        /*Difference will be negative?  Skip it*/
    if n==r  then iterate                        /*Palindromic?   Then it can't be rare.*/

    d= n-r;   parse var  d  ''  -2  y  +1  z     /*obtain the last 2 digs of difference.*/
    if @ps.z  then iterate                       /*Not 0, 1, 4, 5, 6, 9? Not perfect sq.*/
       select
       when z==0   then if y\==0    then iterate /*Does Z = 0?   Then  Y  must be zero. */
       when z==5   then if y\==2    then iterate /*Does Z = 5?   Then  Y  must be two.  */
       when z==6   then if y//2==0  then iterate /*Does Z = 6?   Then  Y  must be odd.  */
       otherwise        if @149.z   then if y//2  then iterate /*Z=1,4,9? Y must be even*/
       end   /*select*/

    s= n+r;   parse var  s  ''  -2  y  +1  z     /*obtain the last two digits of the sum*/
    if @ps.z  then iterate                       /*Not 0, 2, 5, 8, or 9? Not perfect sq.*/
       select
       when z==0   then if y\==0    then iterate /*Does Z = 0?   Then  Y  must be zero. */
       when z==5   then if y\==2    then iterate /*Does Z = 5?   Then  Y  must be two.  */
       when z==6   then if y//2==0  then iterate /*Does Z = 6?   Then  Y  must be odd.  */
       otherwise        if @149.z   then if y//2  then iterate /*Z=1,4,9? Y must be even*/
       end   /*select*/

    $= a + b                                     /*a head start on figuring digital root*/
               do k=3  for length(n) - 2         /*now, process the rest of the digits. */
               $= $ + substr(n, k, 1)            /*add the remainder of the digits in N.*/
               end   /*k*/
                                                 /*This REXX pgm uses 20 decimal digits.*/
       do while $>9                              /* [◄]  Algorithm is good for 111 digs.*/
       if $>9  then $= left($,1) + substr($,2,1)+ substr($,3,1,0)  /*>9? Reduce to a dig*/
       end   /*while*/

    if \@dr.$             then iterate           /*Doesn't have good digital root?  Skip*/
    if iSqrt(s)**2 \== s  then iterate           /*Not a perfect square?  Then skip it. */
    if iSqrt(d)**2 \== d  then iterate           /* "  "    "       "       "    "   "  */
    #= # + 1                                     /*bump the counter of  rare  numbers.  */
    say right( th(#), length(#) + 9)       ' rare number is:  '       right( commas(n), w)
    if #>=many  then leave                       /* [↑]  W:  the width of # with commas.*/
    end   /*n*/
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _;  do jc=length(_)-3  to 1  by -3; _=insert(',', _, jc); end;  return _
th:     parse arg th;return th||word('th st nd rd',1+(th//10)*(th//100%10\==1)*(th//10<4))
/*──────────────────────────────────────────────────────────────────────────────────────*/
iSqrt:  parse arg x;   $= 0;     q= 1;                            do while q<=x;    q= q*4
                                                                  end  /*while q<=x*/
                  do while q>1;  q= q % 4;   _= x-$-q;   $= $ % 2
                  if _>=0  then do;          x= _;       $= $ + q
                                end
                  end   /*while q>1*/;                            return $

Still pretty sluggish,   like molasses in March.


The above REXX program was modified to generate a group of numbers which were   AB   (two digit) numbers
concatenated with   PQ   (two digit)   numbers to yield a list of four digit numbers.

AB   are the 1st two digits of a   rare   number,   and   PQ   are the   last two digits.


This list was sorted and the duplicates removed,   and it formed a list of   (left 2 digits abutted with the right 2 digits)
numbers that every   rare   number must have   (except for the first   rare   number   (65),   which is found the  hard
(slow)   way.

Tweaks, F#

Kudos to Nigel Galloway for the F# version. I don't know the language well, but was able to cut a few corners to improve performance slightly and add some stats:

Output at Tio.run (linked):
nth        Rare Number    elapsed  completed
 1                  65      97 ms 
                           120 ms    2
                           126 ms    3
                           127 ms    4
                           140 ms    5
 2             621,770     148 ms 
                           148 ms    6
                           151 ms    7
                           253 ms    8
 3         281,089,082     261 ms 
                           283 ms    9
 4       2,022,652,202     606 ms 
 5       2,042,832,002    1162 ms 
                          2528 ms   10
                          3423 ms   11
 6     872,546,974,178   16583 ms 
 7     872,568,754,178   17427 ms 
 8     868,591,084,757   28471 ms 
                         37612 ms   12
Of course, it can't get too far in the 60 second timeout window. Sometimes it doesn't get past the 5th number, due to poor luck at Tio.run. It can do 13 digits at Tio.run, as long as you start at 13:
nth        Rare Number    elapsed  completed
 1   6,979,302,951,885   21129 ms 
                         27470 ms   13
Output on a i7 core (Visual Studio Console App):
nth        Rare Number    elapsed  completed
 1                  65      22 ms
                            24 ms    2
                            26 ms    3
                            27 ms    4
                            28 ms    5
 2             621,770      31 ms
                            31 ms    6
                            33 ms    7
                            65 ms    8
 3         281,089,082      69 ms
                            76 ms    9
 4       2,022,652,202     221 ms
 5       2,042,832,002     373 ms
                           859 ms   10
                          1220 ms   11
 6     872,546,974,178    6284 ms
 7     872,568,754,178    6486 ms
 8     868,591,084,757   10336 ms
                         13739 ms   12
 9   6,979,302,951,885   16890 ms
                         19075 ms   13
10  20,313,693,904,202  105822 ms
11  20,313,839,704,202  106337 ms
12  20,331,657,922,202  116631 ms
13  20,331,875,722,202  118324 ms
14  20,333,875,702,202  122171 ms
15  40,313,893,704,200  257026 ms
16  40,351,893,720,200  258086 ms
                        283465 ms   14
Checking up to 14 digit numbers in under 5 minutes. Given more time, it can get the correct 17th number (first 15 digit number), but has problems after that. Not sure if it's the memory requirements, or perhaps I pared down Nigel's original program too much for valid output after 14 digits.

I added this here in the discussion and did not post the revised code on the main codepage out of respect for Nigel's original contribution. All I did was tweak it, and did not add any core improvements.

I've been trying to figure out how to put some limits on the permutation of numbers generated, but don't know F# well enough to do it effectively. And when I work in languages I am familiar with, the performance is decades of magnitude worse.--Enter your username (talk) 18:37, 22 September 2019 (UTC)

As mentioned earlier in this page, I've just added a second Go version based on Nigel's approach which is infinitely faster than the first. You're welcome to try and optimize it further as you're much better at this sort of thing than I am. The perfect square checking might be capable of improvement as I'm just using a simple math.Sqrt approach here rather than having a number of preliminary filters as I did in the first version. --PureFox (talk) 23:11, 24 September 2019 (UTC)

Tweaks, F# (v2)

Thanks again Nigel Galloway, for the improved F# version. (I find it even more delightfully terse, and I understand it even less that the first version.) However, I gave it a couple of shortcuts and got this result on the core i7:

nth              Rare Number   total time   digs  (et per dig)
 1                        65         1 ms
                                     1 ms    2    (      0 ms)
                                     1 ms    3    (      0 ms)
                                     2 ms    4    (      0 ms)
                                     2 ms    5    (      0 ms)
 2                   621,770         5 ms
                                     5 ms    6    (      2 ms)
                                     6 ms    7    (      0 ms)
                                    29 ms    8    (     22 ms)
 3               281,089,082        31 ms
                                    46 ms    9    (     16 ms)
 4             2,022,652,202        83 ms
 5             2,042,832,002       128 ms
                                   461 ms   10    (    413 ms)
                                   758 ms   11    (    296 ms)
 6           872,546,974,178      1665 ms
 7           872,568,754,178      1707 ms
 8           868,591,084,757      3203 ms
                                  6932 ms   12    (   6173 ms)
 9         6,979,302,951,885      8657 ms
                                 12306 ms   13    (   5373 ms)
10        20,313,693,904,202     27181 ms
11        20,313,839,704,202     27278 ms
12        20,331,657,922,202     29120 ms
13        20,331,875,722,202     29399 ms
14        20,333,875,702,202     30104 ms
15        40,313,893,704,200     78829 ms
16        40,351,893,720,200     79112 ms
                                128604 ms   14    ( 116297 ms)
17       200,142,385,731,002    139117 ms
18       221,462,345,754,122    139142 ms
19       816,984,566,129,618    140051 ms
20       245,518,996,076,442    140431 ms
21       204,238,494,066,002    140679 ms
22       248,359,494,187,442    140687 ms
23       244,062,891,224,042    140778 ms
24       403,058,392,434,500    181008 ms
25       441,054,594,034,340    181033 ms
                                230265 ms   15    ( 101661 ms)
26     2,133,786,945,766,212    478386 ms
27     2,135,568,943,984,212    499689 ms
28     8,191,154,686,620,818    503747 ms
29     8,191,156,864,620,818    506482 ms
30     2,135,764,587,964,212    507132 ms
31     2,135,786,765,764,212    509065 ms
32     8,191,376,864,400,818    513368 ms
33     2,078,311,262,161,202    532418 ms
34     8,052,956,026,592,517    893558 ms
35     8,052,956,206,592,517    896814 ms
36     8,650,327,689,541,457    962029 ms
37     8,650,349,867,341,457    963704 ms
38     6,157,577,986,646,405    965923 ms
39     4,135,786,945,764,210   1333750 ms
40     6,889,765,708,183,410   2225840 ms
                               2251017 ms   16    (2020751 ms)
41    86,965,750,494,756,968   2446863 ms
42    22,542,040,692,914,522   2447002 ms
43    67,725,910,561,765,640   3995397 ms
                               4106524 ms   17    (1855506 ms)
44   284,684,666,566,486,482   7761433 ms
45   225,342,456,863,243,522   7887401 ms
46   225,342,458,663,243,522   7930157 ms
47   225,342,478,643,243,522   8019971 ms
48   284,684,868,364,486,482   8085227 ms
49   871,975,098,681,469,178   8439579 ms
50   865,721,270,017,296,468   9124567 ms
51   297,128,548,234,950,692   9135222 ms
52   297,128,722,852,950,692   9145493 ms
53   811,865,096,390,477,018   9249578 ms
54   297,148,324,656,930,692   9286388 ms
55   297,148,546,434,930,692   9306517 ms
56   898,907,259,301,737,498   9679708 ms
57   631,688,638,047,992,345   16317159 ms
58   619,431,353,040,136,925   16376430 ms
59   619,631,153,042,134,925   16559935 ms
60   633,288,858,025,996,145   16629165 ms
61   633,488,632,647,994,145   16685653 ms
62   653,488,856,225,994,125   18768151 ms
63   497,168,548,234,910,690   24649427 ms
                               42231937 ms   18    (38125412 ms)

Tested up to 15 digits in under 4 minutes. 16 digits in under 35 minutes, 17 digits under an hour and 10 minutes. It got to the last of the 18 digit rare numbers in under 7 hours, but it takes 11 3/4 hours to complete the block. Still not quite fast enough to go after 19, 20 or 21 digits.
Tio.run (link) version results:

nth              Rare Number   total time   digs  (et per dig)
 1                        65       156 ms 
                                   190 ms    2    (    147 ms)
                                   199 ms    3    (      0 ms)
                                   200 ms    4    (      1 ms)
                                   201 ms    5    (      0 ms)
 2                   621,770       209 ms 
                                   210 ms    6    (      9 ms)
                                   217 ms    7    (      6 ms)
                                   326 ms    8    (    109 ms)
 3               281,089,082       335 ms 
                                   386 ms    9    (     59 ms)
 4             2,022,652,202       512 ms 
 5             2,042,832,002       740 ms 
                                  1970 ms   10    (   1583 ms)
                                  3055 ms   11    (   1085 ms)
 6           872,546,974,178      6688 ms 
 7           872,568,754,178      6938 ms 
 8           868,591,084,757     13531 ms 
                                 28013 ms   12    (  24957 ms)
 9         6,979,302,951,885     34793 ms 
                                 48838 ms   13    (  20825 ms)
One serious shortcut is using [0L;1L;4L;5L;6L] instead of [0..9]. This is allowed because square numbers can't end with 2, 3, 7, or 8, and 9 doesn't count because the numbers being operated upon here have 9 as a factor. --Enter your username (talk) 21:04, 29 September 2019 (UTC)

Tweaks, Go

Thank you PureFox for the Go translation. Here are the results of some code tweaking. Feel free to re-use anything you like on the main page. Some of the corners that were cut: 1. Computed the forward number, the reverse number and the digital root from the "digits" array all at the same time, rather than going back when needed. 2. Used a Boolean array of the last 2 digits of valid squares to determine whether to do the more computationally expensive floating point square root call. 3. The dl array can be seq(-9, 8), rather than seq(-9, 9). It doesn't seem to matter on the low number of digits ( <16 ) we are covering. The digital root fail check can be either before or after the !isSquare() check, you might find it slightly faster one way or the other.

Output from a core i7:
nth               number         time   completed
 1                    65             0s
                                970.4µs  2
                                970.4µs  3
                                970.4µs  4
                                970.4µs  5
 2               621,770        970.4µs
                                970.4µs  6
                               1.9676ms  7
                               4.9833ms  8
 3           281,089,082       5.9569ms
                               6.9785ms  9
 4         2,022,652,202      19.9197ms
 5         2,042,832,002      40.8908ms
                              84.7729ms 10
                              137.632ms 11
 6       872,546,974,178     646.2443ms
 7       872,568,754,178     677.1611ms
 8       868,591,084,757     1.0701309s
                             1.2416797s 12
 9     6,979,302,951,885     1.6326327s
                             1.9797045s 13
10    20,313,693,904,202    10.5358443s
11    20,313,839,704,202    10.6046898s
12    20,331,657,922,202    12.1206037s
13    20,331,875,722,202    12.3569934s
14    20,333,875,702,202    13.0172337s
15    40,313,893,704,200    22.9626308s
16    40,351,893,720,200    23.1066218s
                            24.4769465s 14
17   200,142,385,731,002    27.4160929s
18   221,462,345,754,122    27.6274982s
19   816,984,566,129,618     30.293394s
20   245,518,996,076,442    31.7365338s
21   204,238,494,066,002    31.9389919s
22   248,359,494,187,442    32.0028151s
23   244,062,891,224,042    32.3977374s
24   403,058,392,434,500    36.6045135s
25   441,054,594,034,340    36.8109598s
                            38.1518922s 15

It can't get past 15 digits, due to the memory requirement. I purposely let the "natural" order of the results display, as I was more interested in the performance time of each result.

Output from Tio.run (linked):
nth               number         time   completed
 1                    65      104.828µs
                              131.427µs  2
                              139.331µs  3
                              178.561µs  4
                              200.489µs  5
 2               621,770      698.249µs
                              757.813µs  6
                              1.42756ms  7
                             8.927975ms  8
 3           281,089,082    10.244811ms
                            12.738426ms  9
 4         2,022,652,202     72.43969ms
 5         2,042,832,002   123.310535ms
                           218.093377ms 10
                           415.295536ms 11
 6       872,546,974,178    1.93598265s
 7       872,568,754,178   2.000165728s
 8       868,591,084,757   2.821645253s
                           3.286164167s 12
 9     6,979,302,951,885   5.924838452s
                           6.864383645s 13

It can't get past 13 digits on Tio.run, due to memory requirement. Execution time at Tio.run is often worse than this, but it always completes in the 60 second time limit. --Enter your username (talk) 22:30, 28 September 2019 (UTC)

Thanks for trying to do something here.
Curiously, it's slower than before when I run it several times on my core i7. The range to get to 15 digits is between 46 and 49 seconds whereas the previous version is steady at around 42 seconds. I've recently upgraded from Go version 1.12.9 to 1.13.1 (the latest as I post this) but I doubt whether it would affect this particular program.
As you're getting 38 seconds for your 'tweaked' version, then your core i7 is probably faster than mine though presumably that time was faster than the original version on the same machine so I'm not sure what to make of it.
As you say there are memory problems when trying to go above 15 digits, though I see that Nigel has today posted a new F# version which has managed to reach 17 digits without using Cartesian products. So I think we're going to have to take another look at it anyway. --PureFox (talk) 18:41, 29 September 2019 (UTC)
The original version of your 2nd Go program runs about 2-3 seconds slower than the version I put at the Tio.run link. So these tweaks are not much of an improvement. Thanks so much for sharing your version. Not sure why it would go slower on your core i7. I just installed Go for the first time at v 1.13.1. I am careful to keep any other programs from executing concurrently that might interfere with the timing measurements. The full name of my cpu is i7-7700 @ 3.6Ghz. Operating on Win10. It is not overclocked, I have not verified it's speed with any benchmarking software. --Enter your username (talk) 04:34, 30 September 2019 (UTC)
Mine's an Intel 8565U which has a much lower base frequency (1.8 GHz) but a higher turbo frequency (4.6 GHz) than yours (4.2 GHz I believe). I suspect yours will be faster overall but there's probably not a great deal in it.
Incidentally, I'm using Ubuntu 18.04 rather than Windows 10 but I've no reason to suppose that Go executes faster on one rather than the other nowadays.
Anyway, the good news is that I've come up with a new strategy which is significantly faster - 15 digits is processed in around 28 seconds compared to 42 seconds previously.
Basically, I'm using a combination of Nigel and Shyam's approachs which cuts down the Cartesian products quite a lot but, unfortunately, still not enough to process 16 digits before running out of memory.
To deal with the memory problem, I've added a second version which delivers the Cartesian products 100 at a time rather than en masse. Not surprisingly, the former is slower than the latter and it's back up to about 43 seconds to get to 15 digits. However, 16 digits takes less than 7 minutes and 17 digits less than 12 minutes - I couldn't be bothered to go any further than that - so it's a worthwhile trade-off.
The idiomatic way to 'yield' values in Go is to spawn a goroutine and pass it a channel (or a buffered channel in this case). However, I'd stress here that I'm not trying to parallelize the algorithm (although one certainly could) as I don't like doing this for RC tasks. --PureFox (talk) 15:14, 30 September 2019 (UTC)
I thought I'd see if I could extend the program to process 18 digit numbers but was surprised when it blew up with an OOM error after hitting the 56th rare number! Frankly, I don't understand why - there appeared to be plenty of unused memory when I profiled it. It may be due to heap fragmentation as the Go GC doesn't compact the heap after a collection and so needs to find a large enough slot for new allocations.
Anyway, I thought I'd get rid of the Cartesian product function altogether and replace it with (in effect) a nested loop and was glad I did as this has restored performance to previous levels and solved the OOM problem. 15, 16 and 17 digits are dispatched in 28 seconds, 4 minutes and 6 minutes respectively and even 18 digits completes in a tolerable 74 minutes.
To reliably go any further than this would require the use of big integers (unpleasant and relatively slow in Go) as signed 64 bit integers have a 19 digit maximum. It might be possible to use unsigned 64 bit integers (20 digit maximum) though this would require some fancy footwork to deal with negative numbers and subtraction. So I think that's my lot now :) --PureFox (talk) 20:04, 2 October 2019 (UTC)

Tweaks, Go (Turbo)

Still some performance hiding in there...

Skipping combinations 2 and 3 for the differences, and stopping at 6 instead of 9. Since no square can end in the digit 2 or 3, these can be skipped safely. 6 is the highest possible difference, so stopping at 6 is OK.

Skipping combinations 0 through 3, 7 through 9, and 12 through 14 on the sums. Since lowest possible first digit of the rare number is 2, and the sum must be greater than the rare number, and digits 2 and 3 cannot be the last digits of the sum, 4 is the lowest possible last digit of the sum. Also, no square can end in digits 7 or 8, so combinations starting with 7, 8, 12, & 13 can be eliminated. Removing combination 9 & 14 is cheating, however no solution up to 19 digits depends on the sum being a square containing the combination 9 or 14 as the first/last digit.

package main
 
import (
    "fmt"
    "math"
    "sort"
    "time"
)
 
type (
    z1 func() z2
    z2 struct {
        value    int64
        hasValue bool
    }
)
 
var pow10 [19]int64
 
func init() {
    pow10[0] = 1
    for i := 1; i < 19; i++ {
        pow10[i] = 10 * pow10[i-1]
    }
}
 
func izRev(n int, i, g uint64) bool {
    if i/uint64(pow10[n-1]) != g%10 {
        return false
    }
    if n < 2 {
        return true
    }
    return izRev(n-1, i%uint64(pow10[n-1]), g/10)
}
 
func fG(n z1, start, end, reset int, step int64, l *int64) z1 {
    i, g, e := step*int64(start), step*int64(end), step*int64(reset)
    return func() z2 {
        for i < g {
            *l += step
            i += step
            return z2{*l, true}
        }
        i = e
        *l -= (g - e)
        return n()
    }
}
 
type nLH struct{ even, odd []uint64 }
 
type zp struct {
    n z1
    g [][2]int64
}
 
func newNLH(e zp) nLH {
    var even, odd []uint64
    n, g := e.n, e.g
    for i := n(); i.hasValue; i = n() {
        for _, p := range g {
            ng, gg := p[0], p[1]
            if (ng > 0) || (i.value > 0) {
                w := uint64(ng*pow10[4] + gg + i.value)
                ws := uint64(math.Sqrt(float64(w)))
                if ws*ws == w {
                    if w%2 == 0 {
                        even = append(even, w)
                    } else {
                        odd = append(odd, w)
                    }
                }
            }
        }
    }
    return nLH{even, odd}
}
 
func makeL(n int) zp {
    g := make([]z1, n/2-3)
    g[0] = func() z2 { return z2{} }
    for i := 1; i < n/2-3; i++ {
        s := -9
        if i == n/2-4 {
            s = -10
        }
        l := pow10[n-i-4] - pow10[i+3]
        acc += l * int64(s)
        g[i] = fG(g[i-1], s, 9, -9, l, &acc)
    }
    var g0, g1, g2, g3 int64
    l0, l1, l2, l3 := pow10[n-5], pow10[n-6], pow10[n-7], pow10[n-8]
    f := func() [][2]int64 {
        var w [][2]int64
        for g0 < 7 {
            nn := g3*l3 + g2*l2 + g1*l1 + g0*l0
            gg := -1000*g3 - 100*g2 - 10*g1 - g0
            if g3 < 9 {
                g3++
            } else {
                g3 = -9
                if g2 < 9 {
                    g2++
                } else {
                    g2 = -9
                    if g1 < 9 {
                        g1++
                    } else {
                        g1 = -9
                        if g0 == 1 { g0 += 2 } 
                        g0++
                    }
                }
            }
            if bs[(pow10[10]+gg)%10000] {
                w = append(w, [2]int64{nn, gg})
            }
        }
        return w
    }
    return zp{g[n/2-4], f()}
}
 
func makeH(n int) zp {
    acc = -(pow10[n/2] + pow10[(n-1)/2])
    g := make([]z1, (n+1)/2-3)
    g[0] = func() z2 { return z2{} }
    for i := 1; i < n/2-3; i++ {
        j := 0
        if i == (n+1)/2-3 {
            j = -1
        }
        g[i] = fG(g[i-1], j, 18, 0, pow10[n-i-4]+pow10[i+3], &acc)
        if n%2 == 1 {
            g[(n+1)/2-4] = fG(g[n/2-4], -1, 9, 0, 2*pow10[n/2], &acc)
        }
    }
    g0 := int64(4)
    var g1, g2, g3 int64
    l0, l1, l2, l3 := pow10[n-5], pow10[n-6], pow10[n-7], pow10[n-8]
    f := func() [][2]int64 {
        var w [][2]int64
        for g0 < 17 {
            nn := g3*l3 + g2*l2 + g1*l1 + g0*l0
            gg := 1000*g3 + 100*g2 + 10*g1 + g0
            if g3 < 18 {
                g3++
            } else {
                g3 = 0
                if g2 < 18 {
                    g2++
                } else {
                    g2 = 0
                    if g1 < 18 {
                        g1++
                    } else {
                        g1 = 0
                        switch g0 {case 6, 9: g0 += 3 }
                        g0++
                    }
                }
            }
            if bs[gg%10000] {
                w = append(w, [2]int64{nn, gg})
            }
        }
        return w
    }
    return zp{g[(n+1)/2-4], f()}
}
 
var (
    acc  int64
    bs   = make([]bool, 10000)
    L, H nLH
)
 
func rare(n int) []uint64 {
    acc = 0
    for g := 0; g < 10000; g++ {
        bs[(g*g)%10000] = true
    }
    L = newNLH(makeL(n))
    H = newNLH(makeH(n))
    var rares []uint64
    for _, l := range L.even {
        for _, h := range H.even {
            r := (h - l) / 2
            z := h - r
            if izRev(n, r, z) {
                rares = append(rares, z)
            }
        }
    }
    for _, l := range L.odd {
        for _, h := range H.odd {
            r := (h - l) / 2
            z := h - r
            if izRev(n, r, z) {
                rares = append(rares, z)
            }
        }
    }
    if len(rares) > 0 {
        sort.Slice(rares, func(i, j int) bool {
            return rares[i] < rares[j]
        })
    }
    return rares
}
 
// Formats time in form hh:mm:ss.fff (i.e. millisecond precision).
func formatTime(d time.Duration) string {
    f := d.Milliseconds()
    s := f / 1000
    f %= 1000
    m := s / 60
    s %= 60
    h := m / 60
    m %= 60
    return fmt.Sprintf("%02d:%02d:%02d.%03d", h, m, s, f)
}
 
func commatize(n uint64) string {
    s := fmt.Sprintf("%d", n)
    le := len(s)
    for i := le - 3; i >= 1; i -= 3 {
        s = s[0:i] + "," + s[i:]
    }
    return s
}
 
func main() {
    bStart := time.Now() // block time
    tStart := bStart     // total time
    nth := 3             // i.e. count of rare numbers < 10 digits
    fmt.Println("nth             rare number    digs  block time    total time")
    for nd := 10; nd <= 19; nd++ {
        rares := rare(nd)
        if len(rares) > 0 {
            for i, r := range rares {
                nth++
                t := ""
                if i < len(rares)-1 {
                    t = "\n"
                }
                fmt.Printf("%2d  %25s%s", nth, commatize(r), t)
            }
        } else {
            fmt.Printf("%29s", "")
        }
        fbTime := formatTime(time.Since(bStart))
        ftTime := formatTime(time.Since(tStart))
        fmt.Printf("  %2d: %s  %s\n", nd, fbTime, ftTime)
        bStart = time.Now() // restart block timing
    }
}
Output:

Results on the core i7-7700 @ 3.6Ghz.

nth             rare number    digs  block time    total time
 4              2,022,652,202
 5              2,042,832,002  10: 00:00:00.001  00:00:00.001
                               11: 00:00:00.006  00:00:00.008
 6            868,591,084,757
 7            872,546,974,178
 8            872,568,754,178  12: 00:00:00.015  00:00:00.024
 9          6,979,302,951,885  13: 00:00:00.098  00:00:00.123
10         20,313,693,904,202
11         20,313,839,704,202
12         20,331,657,922,202
13         20,331,875,722,202
14         20,333,875,702,202
15         40,313,893,704,200
16         40,351,893,720,200  14: 00:00:00.269  00:00:00.392
17        200,142,385,731,002
18        204,238,494,066,002
19        221,462,345,754,122
20        244,062,891,224,042
21        245,518,996,076,442
22        248,359,494,187,442
23        403,058,392,434,500
24        441,054,594,034,340
25        816,984,566,129,618  15: 00:00:01.810  00:00:02.203
26      2,078,311,262,161,202
27      2,133,786,945,766,212
28      2,135,568,943,984,212
29      2,135,764,587,964,212
30      2,135,786,765,764,212
31      4,135,786,945,764,210
32      6,157,577,986,646,405
33      6,889,765,708,183,410
34      8,052,956,026,592,517
35      8,052,956,206,592,517
36      8,191,154,686,620,818
37      8,191,156,864,620,818
38      8,191,376,864,400,818
39      8,650,327,689,541,457
40      8,650,349,867,341,457  16: 00:00:05.122  00:00:07.325
41     22,542,040,692,914,522
42     67,725,910,561,765,640
43     86,965,750,494,756,968  17: 00:00:33.461  00:00:40.787
44    225,342,456,863,243,522
45    225,342,458,663,243,522
46    225,342,478,643,243,522
47    284,684,666,566,486,482
48    284,684,868,364,486,482
49    297,128,548,234,950,692
50    297,128,722,852,950,692
51    297,148,324,656,930,692
52    297,148,546,434,930,692
53    497,168,548,234,910,690
54    619,431,353,040,136,925
55    619,631,153,042,134,925
56    631,688,638,047,992,345
57    633,288,858,025,996,145
58    633,488,632,647,994,145
59    653,488,856,225,994,125
60    811,865,096,390,477,018
61    865,721,270,017,296,468
62    871,975,098,681,469,178
63    898,907,259,301,737,498  18: 00:01:37.823  00:02:18.611
64  2,042,401,829,204,402,402
65  2,060,303,819,041,450,202
66  2,420,424,089,100,600,242
67  2,551,755,006,254,571,552
68  2,702,373,360,882,732,072
69  2,825,378,427,312,735,282
70  6,531,727,101,458,000,045
71  6,988,066,446,726,832,640
72  8,066,308,349,502,036,608
73  8,197,906,905,009,010,818
74  8,200,756,128,308,135,597
75  8,320,411,466,598,809,138  19: 00:12:22.226  00:14:40.838
--Enter your username (talk) 01:32, 21 March 2020 (UTC)
Hey, thanks for that! I must confess I was so pleased at getting the time down to 21 minutes (and so shell-shocked by Nigel's variable naming conventions) that I hadn't even bothered to look at whether further improvement was possible.
It's a little slower on my machine (15 minutes 14 seconds), which seems to be suffering a bit of turbo lag, but a great improvement nonetheless so I'm going to update the version on the main page.
I've also updated Nigel's C++ version with the same tweaks:
#include <iostream>
#include <functional>
#include <bitset>
#include <cmath>
using Z2=std::optional<long>; using Z1=std::function<Z2()>;
constexpr std::array<const long,19> pow10{1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,10000000000,100000000000,1000000000000,10000000000000,100000000000000,1000000000000000,10000000000000000,100000000000000000,1000000000000000000};
constexpr bool izRev(int n,unsigned long i,unsigned long g){return (i/pow10[n-1]!=g%10)? false : (n<2)? true : izRev(n-1,i%pow10[n-1],g/10);}
const Z1 fG(Z1 n,int start, int end,int reset,const long step,long &l){return ([n,i{step*start},g{step*end},e{step*reset},&l,step]()mutable{
    while(i<g){l+=step; i+=step; return Z2(l);} i=e; l-=(g-e); return n();});}
struct nLH{
  std::vector<unsigned long>even{};
  std::vector<unsigned long>odd{};
  nLH(std::pair<Z1,std::vector<std::pair<long,long>>> e){auto [n,g]=e; while (auto i=n()){for(auto [ng,gg]:g){ if((ng>0)|(*i>0)){
    unsigned long w=ng*pow10[4]+gg+*i; unsigned long g=sqrt(w); if(g*g==w) (w%2==0)? even.push_back(w) : odd.push_back(w);}}}}
};
class Rare{
  long acc{0};
  const std::bitset<10000>bs;
  const std::pair<Z1,std::vector<std::pair<long,long>>> makeL(const int n){
    Z1 g[n/2-3]; g[0]=([]{return Z2{};}); 
    for(int i{1};i<n/2-3;++i){int s{(i==n/2-4)? -10:-9}; long l=pow10[n-i-4]-pow10[i+3]; acc+=l*s; g[i]=fG(g[i-1],s,9,-9,l,acc);}
    return {g[n/2-4],([g0{0},g1{0},g2{0},g3{0},l3{pow10[n-8]},l2{pow10[n-7]},l1{pow10[n-6]},l0{pow10[n-5]},this]()mutable{std::vector<std::pair<long,long>>w{}; while (g0<7){
      long n{g3*l3+g2*l2+g1*l1+g0*l0}; long g{-1000*g3-100*g2-10*g1-g0}; if(g3<9) ++g3; else{g3=-9; if(g2<9) ++g2; else{g2=-9; if(g1<9) ++g1; else{g1=-9; if(g0==1) g0=3; ++g0;}}} 
      if (bs[(pow10[10]+g)%10000]) w.push_back({n,g});} return w;})()};}
  const std::pair<Z1,std::vector<std::pair<long,long>>> makeH(const int n){ acc=-(pow10[n/2]+pow10[(n-1)/2]);
    Z1 g[(n+1)/2-3]; g[0]=([]{return Z2{};}); 
    for(int i{1};i<n/2-3;++i) g[i]=fG(g[i-1],(i==(n+1)/2-3)? -1:0,18,0,pow10[n-i-4]+pow10[i+3],acc); 
    if(n%2==1) g[(n+1)/2-4]=fG(g[n/2-4],-1,9,0,2*pow10[n/2],acc);
    return {g[(n+1)/2-4],([g0{4},g1{0},g2{0},g3{0},l3{pow10[n-8]},l2{pow10[n-7]},l1{pow10[n-6]},l0{pow10[n-5]},this]()mutable{std::vector<std::pair<long,long>>w{}; while (g0<17){
      long n{g3*l3+g2*l2+g1*l1+g0*l0}; long g{g3*1000+g2*100+g1*10+g0}; if(g3<18) ++g3; else{g3=0; if(g2<18) ++g2; else{g2=0; if(g1<18) ++g1; else{g1=0; if(g0==6||g0==9)g0+=3; ++g0;}}} 
      if (bs[g%10000]) w.push_back({n,g});} return w;})()};}
  const nLH L,H;
public: Rare(int n):L{makeL(n)},H{makeH(n)},bs{([]{std::bitset<10000>n{false}; for(int g{0};g<10000;++g) n[(g*g)%10000]=true; return n;})()}{
  std::cout<<"Rare "<<n<<std::endl;
  for(auto l:L.even) for(auto h:H.even){unsigned long r{(h-l)/2},z{(h-r)}; if(izRev(n,r,z)) std::cout<<"n="<<z<<" r="<<r<<" n-r="<<l<<" n+r="<<h<<std::endl;}
  for(auto l:L.odd)  for(auto h:H.odd) {unsigned long r{(h-l)/2},z{(h-r)}; if(izRev(n,r,z)) std::cout<<"n="<<z<<" r="<<r<<" n-r="<<l<<" n+r="<<h<<std::endl;}
  }
};
int main(){
  Rare(19);
}
Compiling this with g++ brings the overall execution time down from 30 to 21 minutes in round figures. So the figures for clang or mingw may well come in at below 10 minutes now.
Waiting to see now if Horst.h can blow us all out of the water with a super-fast Pascal version :)

--PureFox (talk) 12:55, 21 March 2020 (UTC)

I recently installed mingw and verified that the tweaked c++ 10-19 version does indeed execute in 8 2/3 minutes for the 19 digit block. The overall time for blocks 10 - 19 is 10 1/3 minutes.
I also installed gmp and tried the c++ 20+ digit version, but found issues. 10, 12, & 13 had normal output, but digits 14 and up spent the approximate calculation time without producing any solutions. Has anyone else observed this issue? (By the way, I had to change long to long long in either version. Perhaps part of my issue?) --Enter your username (talk) 02:30, 23 March 2020 (UTC)
Ah, just missed sub-10 minutes on the overall time but still impressively fast.
In C++ the width of the long type is, of course, implementation dependent and in g++ running on Ubuntu 18.04, amd64, sizeof(long) is 8 bytes. It looks like it's only 4 bytes on the version of mingw you're using (I don't know what Nigel's using), hence the need for long long.
I haven't looked at the 20+ digit C++ version at all though, at first sight, it looks similar to the 10-19 digit version but with long replaced where necessary with mpz_t. FWIW, I compiled (g++ -std=c++17 -O3 rare_ng3.cpp -lgmp -lgmpxx -o rare_ng3) and ran it on my machine and whilst it seemed to execute OK it wasn't very fast taking about 1 minute and 7.25 minutes respectively to complete up to and including the 16 and 17 digit blocks - g++ doesn't seem to be even at the races here :(
I've no idea what could be causing the issues you're having though it may be GMP related as it can be a pig to install correctly. A more fruitful approach might be to try and translate the C++ version to C# using the 128 bit integer library you found rather than System.Numerics.BigInteger though whether this will be any quicker than what you've done already I don't know. --PureFox (talk) 10:37, 23 March 2020 (UTC)

Tweaks, C++

Wringing out some more performance here, just over 5 minutes for 19 digits by themselves, and just over 6 minutes for digits 2 thru 19. See the code comments for details on what was tweaked. Curiously, I found that the g++ version executes faster than the clang++ version. Compiler arguments used: g++ rare.cpp -o rare.exe -std=c++17 -O3
Curious to see if it can go faster on a different platform - But the original executed faster on g++ than clang++ for me on my platform too. Perhaps I don't have clang++ set up the same as others?

// Rare Numbers : Nigel Galloway - December 20th., 2019

/*
speed related tweaks:
 skip certain "outer" permutations which cannot produce squares (ends with 2, 3, 7, or 8), or will not produce squares for less than 20 digits (5 for diffs, and 9 & 14 for sums).
 "outer" compute 5 digits instead of 4, that is, start at 12 digits instead of 10
 diffs pre-computation of squares does only multiples of 9
 don't compute every square up to 100_000, split diffs pre-computation of squares from sums pre-computation, and each has a separate limit for the pre-computation
 do 2 through 11 digits by only looking at "middles" (no "outers") - accomplished by adding an additional definition for nLH() that takes only the function of optional long
 when doing less that 12 digits, limit diffs to be above zero, and sums to be above a limit at which smaller squares produced are less than the forward number itself (pow10[n-1] * 4)
 makeL(), makeH() return nLH() now, instead of nLH() components.  That is, no pair<> construction
 computes "outers" to a single value (vector of long instead of a vector of pair<long, long>)
 sends one power of ten to the mutable section instead of four
 nLH.odd / nLH.even selection based on oddness of current "outer" value, instead of the sqare itself

cosmetic or unnecessary tweaks:
 start calculations at 2 digits instead of 10
 converted output to printf()
 indicate count of solutions
 added elasped time indications
 sorted output
 optional command line argument for maximum number of digits to compute
 compute the power of ten array instead of declaring it with literals
*/

#include <functional>
#include <bitset>
#include <cmath>
#include <chrono>

using namespace std;
using namespace chrono;

template <typename T> // concatenates vectors
vector<T>& operator +=(vector<T>& v, const vector<T>& w) { v.insert(v.end(), w.begin(), w.end()); return v; }

int sc = 0; // solution counter 
auto st = steady_clock::now(), st0 = st, tmp = st; double dir = 0; // for determining elasped time

using Z2 = optional<long long>;
using Z1 = function<Z2()>;
using VU = vector<unsigned long long>;
using VS = vector<string>;

constexpr array<const int,  7>  li { 1, 3, 0, 0, 2, 0, 1 };
constexpr array<const int, 17>  lu { 1, 2, 0, 0, 1, 1, 4, 0, 0, 0, 1, 4, 0, 0, 0, 1, 1 };

// powers of 10 array
constexpr auto pow10 = [] { array <long long, 19> n {1}; for (int j{0}, i{1}; i < 19; j = i++) n[i] = n[j] * 10; return n; } ();

bool izRev(int n, unsigned long long i, unsigned long long g) {
  return (i / pow10[n - 1] != g % 10) ? false : n < 2 ? true : izRev(n - 1, i % pow10[n - 1], g / 10);
}

const Z1 fG(Z1 n, int start, int end, int reset, const long long step, long long &l) {
  return [n, i{step * start}, g{step * end}, e{step * reset}, &l, step] () mutable {
    while (i < g) { i += step; return Z2(l += step); }
    l -= g - (i = e); return n(); };
}

int c = 0; // solution counter
long long acc, l, llim;

struct nLH {
  vector<unsigned long long>even{}, odd{};
  nLH(Z1 a) { unsigned long long r;
    while (auto i = a()) if (*i > llim) {
      r = sqrt(*i); if (r * r == i) *i & 1 ? even.push_back(*i) : odd.push_back(*i); } }
  nLH(Z1 a, vector<long long> b) { unsigned long long sq, r;
    while (auto i = a()) for (auto ng : b) if ((ng > 0) | (*i > 0)) {
      r = sqrt(sq = ng + *i); if (r * r == sq) ng & 1 ? even.push_back(sq) : odd.push_back(sq); } }
};

// formats elasped time
string dFmt(duration<double> et, int digs) {
  string res = ""; double dt = et.count();
  if (dt > 60.0) { int m = (int)(dt / 60.0); dt -= m * 60.0; res = to_string(m) + "m"; }
  res += to_string(dt); return res.substr(0, digs - 1) + 's';
}

// combines list of square differences with list of square sums, reports compatible results
VS dump(int nd, VU lo, VU hi) {
  VS res {};
  for (auto l : lo) for (auto h : hi) {
    auto r { (h - l) >> 1 }, z { h - r };
    if (izRev(nd, r, z)) {
      char buf[99]; sprintf(buf, "%20llu %11lu %10lu", z, (long long)sqrt(h), (long long)sqrt(l));
      res.push_back(buf); } } return res;
}

const double fac = 3.94;
const int mbs = (int)sqrt(fac * pow10[9]), mbt = (int)sqrt(fac * fac * pow10[9]) >> 3;
const bitset<100000>bs {[]{bitset<100000>n{false}; for(int g{3};g<mbs;g+=3) n[(g*g)%100000]=true; return n;}()};
const bitset<100000>bt {[]{bitset<100000>n{false}; for(int g{11};g<mbt;g++) n[(g*g)%100000]=true; return n;}()};

// reports one block of digits
void doOne(int n, nLH L, nLH H) {
  VS lines = dump(n, L.even, H.even); lines += dump(n, L.odd , H.odd); sort(lines.begin(), lines.end());
  duration<double> tet = (tmp = steady_clock::now()) - st; int ls = lines.size();
  if (ls-- > 0)
    for (int i = 0; i <= ls; i++)
      printf("%3d %s%s", ++c, lines[i].c_str(), i == ls ? "" : "\n");
  else printf("%s", string(47, ' ').c_str());
  printf("  %2d:     %s  %s\n", n, dFmt(tmp - st0, 8).c_str(), dFmt(tet, 8).c_str()); st0 = tmp;
}

class Rare {
  const nLH makeL(const int n) {
    constexpr int r = 9; acc = llim = 0; Z1 g = [] { return Z2 {}; }; int s = -r, q = (n > 11) * 5;
      for (int i = 1; i < n / 2 - q + 1; ++i) {
        l = pow10[n - i - q] - pow10[i + q - 1]; s -= i == n / 2 - q; g = fG(g, s, r, -r, l, acc += l * s); }
      return q ? nLH( g, [g0{0}, g1{0}, g2{0}, g3{0}, g4{0}, l3{pow10[n - 5]}] () mutable {
          vector<long long> w {}; long long g; while (g0 < 7) {
            if (bs[((g = -10000 * g4 -1000 * g3 - 100 * g2 - 10 * g1 - g0) + 1000000000000LL) % 100000]) w.push_back(l3 * (g4 + g3 * 10 + g2 * 100 + g1 * 1000 + g0 * 10000) + g);
            if (g4 < r) ++g4; else { g4 = -r; if (g3 < r) ++g3; else { g3 = -r; if (g2 < r) ++g2; else { g2 = -r; if (g1 < r) ++g1; else { g1 = -r; g0 += li[g0]; } } } } }
      return w; } () ) : nLH(g);
  }
  const nLH makeH(const int n) {
    constexpr int r = 18; llim = pow10[n - 1] << 2; acc = -pow10[n >> 1] - pow10[(n - 1) >> 1]; Z1 g = [] { return Z2 {}; };
    int q = (n > 11) * 5;
      for (int i = 1; i < (n >> 1) - q + 1; ++i)
        g = fG(g, 0, r, 0, pow10[n - i - q] + pow10[i + q - 1], acc); 
      if (n & 1) { l = pow10[n >> 1] << 1; g = fG(g, 0, r >> 1, 0, l, acc += l); }
      return q ? nLH(g, [g0{4}, g1{0}, g2{0}, g3{0}, g4{0}, l3{pow10[n - 5]}] () mutable {
          vector<long long> w {}; long long g; while (g0 < 17) {
            if (bt[(g = g4 * 10000 + g3 * 1000 + g2 * 100 + g1 * 10 + g0) % 100000]) w.push_back(l3 * (g4 + g3 * 10 + g2 * 100 + g1 * 1000 + g0 * 10000) + g);
            if (g4 < r) ++g4; else { g4 = 0; if (g3 < r) ++g3; else { g3 = 0; if (g2 < r) ++g2; else { g2 = 0; if (g1 < r) ++g1; else { g1 = 0; g0 += lu[g0]; } } } } }
        return w; } () ) : nLH(g);
  }
  public: Rare(int n) { doOne(n, makeL(n), makeH(n)); }
};

int main(int argc, char *argv[]) {
  int max = argc > 1 ? stoi(argv[1]) : 19; if (max < 2) max = 2; if (max > 19 ) max = 19;
  printf("%4s %19s %11s %10s %5s %11s %9s\n", "nth", "forward", "rt.sum", "rt.diff", "digs", "block.et", "total.et");
  for (int nd = 2; nd <= max; nd++) Rare(int(nd));
}
Output:

Results on the core i7-7700 @ 3.6Ghz.

 nth             forward      rt.sum    rt.diff  digs    block.et  total.et
  1                   65          11          3   2:     0.00334s  0.00334s
                                                  3:     0.00289s  0.00623s
                                                  4:     0.00244s  0.00868s
                                                  5:     0.00320s  0.01188s
  2               621770         836        738   6:     0.00308s  0.01496s
                                                  7:     0.00281s  0.01777s
                                                  8:     0.00336s  0.02113s
  3            281089082       23708        330   9:     0.00746s  0.02860s
  4           2022652202       63602        300
  5           2042832002       63602       6360  10:     0.01855s  0.04715s
                                                 11:     0.09707s  0.14422s
  6         868591084757     1275175     333333
  7         872546974178     1320616      32670
  8         872568754178     1320616      33330  12:     0.01456s  0.15879s
  9        6979302951885     3586209    1047717  13:     0.04947s  0.20826s
 10       20313693904202     6368252     269730
 11       20313839704202     6368252     270270
 12       20331657922202     6368252     329670
 13       20331875722202     6368252     330330
 14       20333875702202     6368252     336330
 15       40313893704200     6368252    6330336
 16       40351893720200     6368252    6336336  14:     0.12746s  0.33572s
 17      200142385731002    20006998      69300
 18      204238494066002    20122102    1891560
 19      221462345754122    21045662      69300
 20      244062891224042    22011022    1908060
 21      245518996076442    22140228     921030
 22      248359494187442    22206778    1891560
 23      403058392434500    20211202   19940514
 24      441054594034340    22011022   19940514
 25      816984566129618    40421606     250800  15:     0.73898s  1.07470s
 26     2078311262161202    64030648    7529850
 27     2133786945766212    65272218    2666730
 28     2135568943984212    65272218    3267330
 29     2135764587964212    65272218    3326670
 30     2135786765764212    65272218    3333330
 31     4135786945764210    65272218   63333336
 32     6157577986646405   105849161   33333333
 33     6889765708183410    83866464   82133718
 34     8052956026592517   123312255   29999997
 35     8052956206592517   123312255   30000003
 36     8191154686620818   127950856    3299670
 37     8191156864620818   127950856    3300330
 38     8191376864400818   127950856    3366330
 39     8650327689541457   127246955   33299667
 40     8650349867341457   127246955   33300333  16:     2.25107s  3.32578s
 41    22542040692914522   212329862     333300
 42    67725910561765640   269040196  251135808
 43    86965750494756968   417050956      33000  17:     13.6768s  17.0026s
 44   225342456863243522   671330638     297000
 45   225342458663243522   671330638     303000
 46   225342478643243522   671330638     363000
 47   284684666566486482   754565658      30000
 48   284684868364486482   754565658     636000
 49   297128548234950692   770186978   32697330
 50   297128722852950692   770186978   32702670
 51   297148324656930692   770186978   33296670
 52   297148546434930692   770186978   33303330
 53   497168548234910690   770186978  633363336
 54   619431353040136925  1071943279  299667003
 55   619631153042134925  1071943279  300333003
 56   631688638047992345  1083968809  297302703
 57   633288858025996145  1083968809  302637303
 58   633488632647994145  1083968809  303296697
 59   653488856225994125  1083968809  363303363
 60   811865096390477018  1273828556   33030330
 61   865721270017296468  1315452006   32071170
 62   871975098681469178  1320582934    3303300
 63   898907259301737498  1339270086   64576740  18:     42.5039s  59.5065s
 64  2042401829204402402  2021001202   18915600
 65  2060303819041450202  2020110202  199405140
 66  2420424089100600242  2200110022   19080600
 67  2551755006254571552  2259094848     693000
 68  2702373360882732072  2324811012     693000
 69  2825378427312735282  2377130742    2508000
 70  6531727101458000045  3454234451 1063822617
 71  6988066446726832640  2729551744 2554541088
 72  8066308349502036608  4016542096    2508000
 73  8197906905009010818  4046976144  133408770
 74  8200756128308135597  4019461925  495417087
 75  8320411466598809138  4079154376   36366330  19:     5m3.388s  6m2.894s

--Enter your username (talk) 23:46, 25 May 2020 (UTC)

Good to see the spirit of C is alive and well multyplying bools by integers and mysterious 3.94's. Idiotmatic? who cares? but old fashioned never, so perhaps int sc{0}; rather than int sc=0;. I have compiled the code using mingw on an Core I5 1035G1 and using g++ and clang++ on a Core I7 Q720 (now a very old machine). Interestingly the poor old Q720 takes about the same time for g++ and clang++ with this code. The time taken by the Core I5 1035G1 is the same as your i7 for 2..18 but interestingly about 10secs faster for 19. I shall look at the actual changes over the next few days, the important thing is that they mustnot be based on knowing the result. The timings I have obtained are:
 nth             forward      rt.sum    rt.diff  digs    block.et  total.et
  1                   65          11          3   2:     0.00027s  0.00027s
                                                  3:     0.00002s  0.00030s
                                                  4:     0.00001s  0.00031s
                                                  5:     0.00002s  0.00034s
  2               621770         836        738   6:     0.00005s  0.00040s
                                                  7:     0.00025s  0.00066s
                                                  8:     0.00075s  0.00142s
  3            281089082       23708        330   9:     0.00477s  0.00619s
  4           2022652202       63602        300
  5           2042832002       63602       6360  10:     0.01432s  0.02051s
                                                 11:     0.08783s  0.10835s
  6         868591084757     1275175     333333
  7         872546974178     1320616      32670
  8         872568754178     1320616      33330  12:     0.01142s  0.11977s
  9        6979302951885     3586209    1047717  13:     0.05034s  0.17011s
 10       20313693904202     6368252     269730
 11       20313839704202     6368252     270270
 12       20331657922202     6368252     329670
 13       20331875722202     6368252     330330
 14       20333875702202     6368252     336330
 15       40313893704200     6368252    6330336
 16       40351893720200     6368252    6336336  14:     0.12319s  0.29331s
 17      200142385731002    20006998      69300
 18      204238494066002    20122102    1891560
 19      221462345754122    21045662      69300
 20      244062891224042    22011022    1908060
 21      245518996076442    22140228     921030
 22      248359494187442    22206778    1891560
 23      403058392434500    20211202   19940514
 24      441054594034340    22011022   19940514
 25      816984566129618    40421606     250800  15:     0.72540s  1.01872s
 26     2078311262161202    64030648    7529850
 27     2133786945766212    65272218    2666730
 28     2135568943984212    65272218    3267330
 29     2135764587964212    65272218    3326670
 30     2135786765764212    65272218    3333330
 31     4135786945764210    65272218   63333336
 32     6157577986646405   105849161   33333333
 33     6889765708183410    83866464   82133718
 34     8052956026592517   123312255   29999997
 35     8052956206592517   123312255   30000003
 36     8191154686620818   127950856    3299670
 37     8191156864620818   127950856    3300330
 38     8191376864400818   127950856    3366330
 39     8650327689541457   127246955   33299667
 40     8650349867341457   127246955   33300333  16:     2.24472s  3.26344s
 41    22542040692914522   212329862     333300
 42    67725910561765640   269040196  251135808
 43    86965750494756968   417050956      33000  17:     13.9236s  17.1870s
 44   225342456863243522   671330638     297000
 45   225342458663243522   671330638     303000
 46   225342478643243522   671330638     363000
 47   284684666566486482   754565658      30000
 48   284684868364486482   754565658     636000
 49   297128548234950692   770186978   32697330
 50   297128722852950692   770186978   32702670
 51   297148324656930692   770186978   33296670
 52   297148546434930692   770186978   33303330
 53   497168548234910690   770186978  633363336
 54   619431353040136925  1071943279  299667003
 55   619631153042134925  1071943279  300333003
 56   631688638047992345  1083968809  297302703
 57   633288858025996145  1083968809  302637303
 58   633488632647994145  1083968809  303296697
 59   653488856225994125  1083968809  363303363
 60   811865096390477018  1273828556   33030330
 61   865721270017296468  1315452006   32071170
 62   871975098681469178  1320582934    3303300
 63   898907259301737498  1339270086   64576740  18:     42.8003s  59.9873s
 64  2042401829204402402  2021001202   18915600
 65  2060303819041450202  2020110202  199405140
 66  2420424089100600242  2200110022   19080600
 67  2551755006254571552  2259094848     693000
 68  2702373360882732072  2324811012     693000
 69  2825378427312735282  2377130742    2508000
 70  6531727101458000045  3454234451 1063822617
 71  6988066446726832640  2729551744 2554541088
 72  8066308349502036608  4016542096    2508000
 73  8197906905009010818  4046976144  133408770
 74  8200756128308135597  4019461925  495417087
 75  8320411466598809138  4079154376   36366330  19:     4m52.40s  5m52.39s
 
g++ (SUSE Linux) 9.2.1 20200109 [gcc-9-branch revision 280039]
 
 40     8650349867341457   127246955   33300333  16:     15.5694s  22.0091s
 43    86965750494756968   417050956      33000  17:     1m33.80s  1m55.81s
 63   898907259301737498  1339270086   64576740  18:     4m57.66s  6m53.47s
 75  8320411466598809138  4079154376   36366330  19:     30m31.1s  37m24.5s

clang version 9.0.1

 40     8650349867341457   127246955   33300333  16:     15.1173s  21.4643s
 43    86965750494756968   417050956      33000  17:     1m32.36s  1m53.82s
 63   898907259301737498  1339270086   64576740  18:     4m47.22s  6m41.04s
 75  8320411466598809138  4079154376   36366330  19:     29m5.71s  35m46.7s

--Nigel Galloway (talk) 14:05, 10 June 2020 (UTC)

Thanks for those comparisons, appreciated. Re 3.94, the limit is 4, I was just undercutting it a bit. --Enter your username (talk) 00:58, 17 June 2020 (UTC)
  • nLH.odd / nLH.even selection based on oddness of current "outer" value, instead of the sqare itself: Okay but the odd numbers seem to be in the even bucket and the even numbers are in the odd bucket. Probably my fault for calling them odd and even since the only requirement is that the are split, not what the bucket is called, but I couldn't sleep with them the wrong way round so I've but them in the correct named bucket.--Nigel Galloway (talk) 14:09, 4 January 2021 (UTC)
  • skip certain "outer" permutations which cannot produce squares (ends with 2, 3, 7, or 8), or will not produce squares for less than 20 digits (5 for diffs, and 9 & 14 for sums): There is no reason why 5 can not produce a rare number. 9 & 14 will not for any Rare number of any length. The situation for g0 and g1 for all rare numbers is:
Diffs (L)

g0 -> g1 
 0 ->  0
 1 -> -7..2..9
 4 -> -8..2..8
 5 -> -3;7
 6 -> -9..2..9

Sums (H)

g0 -> g1 
 4 ->  0..2..18
 6 ->  1..2..17
10 ->  9
11 ->  1..2..17
15 ->  1;11
16 ->  0..2..18

--Nigel Galloway (talk) 14:09, 4 January 2021 (UTC)

  • Do 2 through 11 digits by only looking at "middles" (no "outers") - accomplished by adding an additional definition for nLH() that takes only the function of optional long when doing less that 12 digits, limit diffs to be above zero, and sums to be above a limit at which smaller squares produced are less than the forward number itself (pow10[n-1] * 4) makeL(), makeH() return nLH() now, instead of nLH() components: Okay but no need for an extra global variable. It is really a parameter to nLH's constructor (single) as now only one constructor is required.--Nigel Galloway (talk) 14:09, 4 January 2021 (UTC)
  • The pair "construction" shouldn't incur an overhead as the compiler should optimize it away. It adds flexability to the design, the need for which I shall reserve judgement.--Nigel Galloway (talk) 14:09, 4 January 2021 (UTC)
Thank you for the code review and adopting some of the tweaks I presented. Good fix for my awkward handling of digits 2-11. I agree about the pair construction, I only removed it to speed up my C# translation, which was slightly impacting performance (perhaps due to an awkward translation on my part.)

I agree that digit '5' ought to be considered in the "diff" side. Even though it doesn't contribute to any solution in the Int64 range, if Int128 ever becomes part of C++, your program should be easily extended to the point where checking "diff 5" will contribute a solution. My decision to omit checking '5' was for performance reasons since we are constricted to Int64.--Enter your username (talk) 00:01, 25 January 2021 (UTC)

21+ digit rare numbers

Well, one anyway (so far). I tweaked the BigInteger version of the C# program to skip to start at 21 digits. Around 6 hours, I got the first one: 219,518,549,668,074,815,912, with the sum = 20,953,210,2682, and the difference = 8,877,0002. Still have no idea how long it will take to finish the block of 21 digit numbers. Since the difference found so far was a relatively low number, it probably has quite a while to go.

I am also running another instance that checks the block of 20 digit numbers (in order to verify the algorithm against the table of known rare numbers), but after 6 hours, it still hasn't come up with anything yet. A little surprising, as there are a few 20 digit rare numbers with 7 digit differences. If I don't see anything on the 20 digit run in 6 more hours, there may be some kind of issue to work out. --Enter your username (talk) 02:52, 21 October 2019 (UTC)

P.S. 5 found so far:

 Nth  Time (hours)  rare number
  85         6      219,518,549,668,074,815,912 
  86        10 1/2  837,982,875,780,054,779,738
  87        11 1/2  208,393,425,242,000,083,802 
  88        12 1/3  286,694,688,797,362,186,682
  89        13 2/3  257,661,195,832,219,326,752

--Enter your username (talk) 21:29, 22 October 2019 (UTC)

In the Go program, I took advantage of Shyam's observation that no rare number up to 10^20 ended in 3. So, as far as the 21 digit numbers are concerned, a possible fly in the ointment is that they could in theory end in 3 (with a first digit of 8) which would, of course, mean adding another key/value pair to the 'fml' dictionary with other consequent amendments. The chances are that there won't be any such numbers but you might want to bear it in mind if you have occasion to run this instance again.
As far as the 20 digit numbers are concerned, there shouldn't be any issues if you're running a BigInteger version for that too and adjusting the IsSquare method accordingly.

--PureFox (talk) 12:58, 23 October 2019 (UTC)

Regarding missing solutions at 20 digits, I found that the Math.Sqrt() function may drop a few bits of resolution at 20 digits or more. Not totally unexpected. However, the Math.Sqrt() function can be used as an initial accurate first guess in a Newton's method integer square root function.
There is an integer square root function that only uses integers and floating point multiple and division isn't used.     -- Gerard Schildberger (talk) 07:31, 12 January 2020 (UTC)
I agree the 8,3 combo should be added for 21+ digits computation. For the task as defined, (5th to 8th number), it's not needed. But for going beyond the table on Shyam Sunder Gupta's webpage, one really needs to consider that combination.
I e-mailed Shyam Sunder Gupta and he said he'd update his webpage.   He said that he figured nobody would compute rare numbers that high, so he didn't bother to enter the updated list before.   He also mentioned that he had found a rare number ending in the decimal digit 3.     -- Gerard Schildberger (talk) 07:26, 12 January 2020 (UTC)
I see this BigInteger conversion attempt as a means of reaching a few 21 digit numbers (which it did) and a means to reveal any shortcomings in the existing ulong algorithm which don't translate well to the BigInteger version (which it also did, I guess). If only it wasn't so impractically slow. If I could get it an order of magnitude faster, it would easier to persue this further. Creating a multi-tasking version would probably only go 2 to 4 times faster.--Enter your username (talk) 06:57, 24 October 2019 (UTC)

There only seem to be 5 21 digit rare numbers, so I started looking at 22 digits. I found that each odd number of digits takes less time (about 80% of the time) of the previous even number of digits, but each even/odd pair (number of digits such as 16/17 vs 14/15) takes about 20 times a long as the previous pair, so the computation time increases dramatically for results above 19 digits. Here are number of digits 21 and 22:

  S.No.			     R		 R+R1		R-R1
     85	 208393425242000083802	 204150294022	  1158663002
     86	 219518549668074815912	 209532102682	    88770002
     87	 257661195832219326752	 226998922482	  1930896002
     88	 286694688797362186682	 239452699422	  1158663002
     89	 837982875780054779738	 409384944262	   736593002
     90	2414924301133245383042	 694172869282	 33299966702
     91	2414924323311045383042	 694172869282	 33300033302
     92	2414946523311023183042	 694172869282	 33366633302
     93	2576494891793995836752	 717835697482	  3299967002
     94	2576494893971995836752	  17835697482	  3300033002
     95	2620937863931054483162	 723517958682	 26633367302
     96	2620955623931476283162	 723517958682	 26699967302
     97	2620955641393276283162	 723517958682	 26700032702
     98	2622935621573476481162	 723517958682	 33299966702
     99	2622935643751276481162	 723517958682	 33300033302
    100	2622937641933274481162	 723517958682	 33306033302
    101	2622955841933256281162	 723517958682	 33360633302
    102	2622957843751254281162	 723517958682	 33366633302
    103	2727651947516658327272	 738572306122	  6429474002
    104	2747736918335953517072	 738572306122	 63705041402
    105	2788047668617596408872	 746732524122	   266730002
    106	2788047848617776408872	 746732524122	   327330002
    107	2788047868437576408872	 746732524122	   333330002
    108	2788047888617376408872	 746732524122	   339330002
    109	2939501759705522349392	 766742069722	  2636373002
    110	2939503375709360349392	 766742069722	  2696973002
    111	2939503537707740349392	 766742069722	  2702973002
    112	2939521359525562149392	 766742069722	  3297033002
    113	2939521557527542149392	 766742069722	  3303033002
    114	2939523577527340149392	 766742069722	  3363633002
    115	2939523779525320149392	 766742069722	  3369633002
    116	2959503377707360349192	 766742069722	 63303033602
    117	6344828989519887483525	1076971535312	330300030332
    118	8045841652464561594308	1268100678462	 32999996702
    119	8045841654642561594308	1268100678462	 33000003302
    120	8655059576513659814468	1315266100062	 32969703302
    121	8655059772157639814468	1315266100062	 32970296702
    122	8655079374155679614468	1315266100062	 33029696702
    123	8655079574515659614468	1315266100062	 33030303302 

That took about 5 3/4 days of computation. Only one odd number. Lots of 2,2 combinations there. Nearly half of all rare numbers under 21 digits start and end with 2. The trend continues above 20 digits with over half starting and ending with 2. --Enter your username (talk) 06:17, 2 December 2019 (UTC)

Oops, it appears that I missed one, Shyam Sunder Gupta has published an updated list of 124 Rare numbers up to 22 digits on his webpage, back on December 15th. The one I missed is the first Rare number ending with 3.
    124	8888070771864228883913	1099179648492	754598074952
At least that verifies that there are only 5 21 digit Rare numbers. --Enter your username (talk) 22:46, 26 December 2019 (UTC)

Some 23s:

  125  (20,006,212,343,920,163,220,002)
  126  (20,404,210,361,902,143,200,402)
  127  (21,544,373,975,964,337,344,512)
  128  (22,781,275,420,027,357,218,722)
  129  (80,618,209,916,486,890,281,608)
  130  (81,313,065,142,333,312,588,218)
  131  (84,247,683,299,691,574,674,248)
  132  (89,650,295,750,128,200,205,698)

That's 8 found for 23 digits, taking ~4 1/5 days to go through the combinations. --Enter your username (talk) 17:54, 15 December 2019 (UTC)

Now that you've established that there are five 21 digit rare numbers, it might be worth mailing SSG (at guptass@rediffmail.com) to see if he'll add them to his list.
It's a pity that we don't have 128 bit integers to speed up the calculations. There's a proposal to add them to Go which has plenty of support but the Go team doesn't seem particularly keen.
C# does of course have its 'decimal' type which is 16 bit and can deal with up to 28 digit integers. Although it's much slower than 'double', it might still be faster than BigInteger.
I also found this open source library for C# which might be worth checking out. --PureFox (talk) 17:25, 2 December 2019 (UTC)
To get to 22 digits, I had to use multitasking and that C# UI128 package you spotted. The performance of the UI128 package is good, only about two and a half times worse than UI64's, whereas BigInteger is around 7 to 8 times worse. The UI128 package should go up to 38 decimal digits or so. But because it takes 20 times longer for each additional pair of digits, I don't think 38 digits a possibilty with the current algorithm. I am contemplating a new algorithm that could be more efficient.
Regarding the 5 21 digit rare numbers found, to prove that there are only 5, what is really needed is to prove that all 1021 - 5 of the possible 1021 are not rare. If the algorithm that finds 5 misses a couple, thats not a good result. I'm confident to report 5 found, but not quite absolutely sure they are the only 5 at this point in time.--Enter your username (talk) 03:04, 3 December 2019 (UTC)
Congratulations on completing the 23 digit rare numbers! You've a lot more patience (and spare computing capacity) than me :) --PureFox (talk) 18:19, 15 December 2019 (UTC)
I've entered code in C++ which I expect using clang++ on the monster beasts you have for computers will complete 21 in a day and a quarter and 22 in less than 4 days. The difference between clang++ and g++ is suprising. It would be interesting to know how MSVC does.--Nigel Galloway (talk) 12:28, 20 December 2019 (UTC)
I don't have clang but, compiling your C++ code for 10 to 19 digits on my core i7 using g++ v7.4 (-std=c++17 -O3) on Ubuntu 18.04, produces execution times of 11.7, 78, 221 and 1484 seconds for rare numbers with 16, 17, 18 and 19 digits respectively. The corresponding times for the Go entry (with the Julia entry not far behind) were 221, 355, 4532 and 6610 seconds so, even if we assume these languages are 2 or 3 times slower than C++, the algorithm you're now using is considerably more efficient than what we had before.
Thanks for that. I've timed it on a Core I5 1035G1 and obtained 17->33sec; 18->92sec; and 19->741sec. For comparison I compiled the C# and obtained 17->6m45sec (compared to 2m12sec on an i7 7700 claimed on the page).--Nigel Galloway (talk) 16:40, 13 March 2020 (UTC)
We've certainly had a wide variation in the timings for this task and it appears now that g++ is much slower than both clang and mingw for some reason. The C# time for 17 numbers looks right to me as my Go translation was only a second or two behind, again using Core I7. I'm surprised how much quicker C++ is than C#, its stable-mate VB.NET and Go for this task. Although Enter your Username has clearly tried hard to minimize the effect of GC by declaring huge swathes of variables 'static', the performance gap is still huge. --PureFox (talk) 11:35, 14 March 2020 (UTC)
I thought for good measure I'd add a Go translation of your C++ program (10 to 19 digits version) and this has cut the execution time from 54 to 21 minutes in round figures which seems more in line with expectations. For comparison, I compiled the C++ program again (using g++ 7.5.0 this time) and ran it on the same machine but the total execution time was almost identical at around 30 minutes so it's difficult to know what to make of it. Perhaps g++ is not yet quite up to speed with C++ 17 features? --PureFox (talk) 17:07, 16 March 2020 (UTC)
I thought for good measure I'd (tit for tat?) install Go and run your goTurbo code on my i5-1035G1. I obtained the following:
       40      8,650,349,867,341,457  16: 00:00:06.403  00:00:09.163
       43     86,965,750,494,756,968  17: 00:00:42.818  00:00:51.982             (mingw real    0m33.328s)
       63    898,907,259,301,737,498  18: 00:02:01.282  00:02:53.264             (mingw real    1m32.945s)
       75  8,320,411,466,598,809,138  19: 00:15:36.454  00:18:29.719             (mingw real   12m21.298s)

I think this demonstrates that the task should require at least the first 63 Rare numbers!!! --Nigel Galloway (talk) 12:22, 19 March 2020 (UTC)

LOL, the Go version is even faster on your i5 than it is on my i7, and only about 50% behind mingw! I imagine you're running it on Windows 10 whereas I'm using Ubuntu 18.04 but that shouldn't make much difference. Probably should have bought an i7 with a higher basic clock speed as performance can be a bit disappointing at times.
The stretch goal is already open-ended and Enter your username has certainly gone well beyond the call of duty there :) --PureFox (talk) 14:32, 19 March 2020 (UTC)
I find it hard to believe that clang++ could be 3 times faster than g++ (though you can get some strange results with these CPU-intensive tasks) and, although I no longer have an up to date Windows machine, on past form I'd be surprised if Visual C++ were any faster than g++ itself. Enter your username may be able to confirm the position there. --PureFox (talk) 20:53, 20 December 2019 (UTC)
I tried compiling it using MSVC. It crashed the compiler, well at least caused it to catch an internal error.--Nigel Galloway (talk) 16:40, 13 March 2020 (UTC)

Broken Link

The link to lots of facts and hints on the main page seems to be broken!