Rare numbers: Difference between revisions

Content added Content deleted

Inline

Revision as of 19:41, 31 December 2019

Definitions and restrictions

Rare numbers are positive integers n where:

n is expressed in base ten
r is the reverse of n (decimal digits)
n must be non-palindromic (n ≠ r)
(n+r) is the sum
(n-r) is the difference and must be positive
the sum and the difference must be perfect squares

Task

find and show the first 5 rare numbers
find and show the first 8 rare numbers (optional)
find and show more rare numbers (stretch goal)

Show all output here, on this page.

References

an OEIS entry: A035519 rare numbers.
an OEIS entry: A059755 odd rare numbers.
planetmath entry: rare numbers. (some hints)
author's website: rare numbers by Shyam Sunder Gupta. (lots of hints and some observations).

C++

I am travelling at the moment so all timings for the following are made on my laptop. Expect better than half when I run them on my desktop in the new year.

The task

The following is a simple implementation that demonstrates the principle of n-r=L and n+r=H. Interestingly it compile with both g++ and clang++, but g++ produces incorrect output. It is sufficient to meet the unambitious requirements of this task. <lang cpp> // Rare Numbers : Nigel Galloway - December 20th., 2019

include <iostream>
include <functional>
include <cmath>
include <numeric>

constexpr std::array<const long,18> pow10{1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,10000000000,100000000000,1000000000000,10000000000000,100000000000000,1000000000000000,10000000000000000,100000000000000000}; constexpr auto r1=([]{std::array<bool,10>n{}; for(auto g:{0,1,4,5,6,9}) n[g]=true; return n;})(); constexpr bool izRev( int n, long i, long g){return (i/pow10[n-1]!=g%10)? false : (n<2)? true : izRev(n-1,i%pow10[n-1],g/10);} struct nLH{

 std::vector<long>even{};
 std::vector<long>odd{};
 nLH(std::function<std::optional<long>()> g){while (auto n=g()){if (([n]{long g=sqrt(*n); return (g*g==n);})()) (*n%2==0)? even.push_back(*n) : odd.push_back(*n);}}

}; template<int z>void Rare(){

 auto L=nLH(([n=std::array<int,z/2>{},p=([]{int g{z/2}; std::array<long,z/2>n{}; for(auto& n:n){n=pow10[z-g]-pow10[g-1]; --g;} return n;})()]() mutable{
   for (auto g=n.begin();g<(n.end());++g) if (*g<9){*g+=1; while(!r1[(n[z/2-1]*9)%10]) *g+=1; return std::optional{std::inner_product(n.begin(),n.end(),p.begin(),0L)};} else *g=-9;
   return std::optional<long>{};}));
 auto H=nLH(([n=([]{std::array<int,(z+1)/2>n{}; *(n.end()-1)=1; *n.begin()-=1; return n;})(),
              p=([]{int g{z/2}; std::array<long,(z+1)/2>n{}; for(auto& i:n) if (z%2==1&n[0]==0) i=2*pow10[z/2]; else {i=pow10[z-g]+pow10[g-1]; --g;} return n;})()]() mutable{
   for (auto g=n.begin();g<(n.end());++g) if (*g<19&(z%2==0|n[0]<10)){
        *g+=1; while(!r1[n[(z+1)/2-1]]) *g+=1; return std::optional{std::inner_product(n.begin(),n.end(),p.begin(),0L)};} else *g=0;
   return std::optional<long>{};}));
 std::cout<<"Rare numbers of length "<<z<<std::endl;
 for(auto l:L.even) for(auto h:H.even){long r=(h-l)/2; if(izRev(z,r,h-r)) std::cout<<"n="<<h-r<<" r="<<r<<" n-r="<<l<<" n+r="<<h<<std::endl;}
 for(auto l:L.odd)  for(auto h:H.odd) {long r=(h-l)/2; if(izRev(z,r,h-r)) std::cout<<"n="<<h-r<<" r="<<r<<" n-r="<<l<<" n+r="<<h<<std::endl;}

} int main(){

 Rare<2>();
 Rare<3>();
 Rare<4>();
 Rare<5>();
 Rare<6>();
 Rare<7>();
 Rare<8>();
 Rare<9>();
 Rare<10>();
 Rare<11>();
 Rare<12>();
 Rare<13>();
 Rare<14>();
 Rare<15>();
 Rare<16>();

} </lang>

Output:

Rare numbers of length 2
n=65 r=56 n-r=9 n+r=121
Rare numbers of length 3
Rare numbers of length 4
Rare numbers of length 5
Rare numbers of length 6
n=621770 r=77126 n-r=544644 n+r=698896
Rare numbers of length 7
Rare numbers of length 8
Rare numbers of length 9
n=281089082 r=280980182 n-r=108900 n+r=562069264
Rare numbers of length 10
n=2022652202 r=2022562202 n-r=90000 n+r=4045214404
n=2042832002 r=2002382402 n-r=40449600 n+r=4045214404
Rare numbers of length 11
Rare numbers of length 12
n=872546974178 r=871479645278 n-r=1067328900 n+r=1744026619456
n=872568754178 r=871457865278 n-r=1110888900 n+r=1744026619456
n=868591084757 r=757480195868 n-r=111110888889 n+r=1626071280625
Rare numbers of length 13
n=6979302951885 r=5881592039796 n-r=1097710912089 n+r=12860894991681
Rare numbers of length 14
n=20313693904202 r=20240939631302 n-r=72754272900 n+r=40554633535504
n=20313839704202 r=20240793831302 n-r=73045872900 n+r=40554633535504
n=20331657922202 r=20222975613302 n-r=108682308900 n+r=40554633535504
n=20331875722202 r=20222757813302 n-r=109117908900 n+r=40554633535504
n=20333875702202 r=20220757833302 n-r=113117868900 n+r=40554633535504
n=40313893704200 r=240739831304 n-r=40073153872896 n+r=40554633535504
n=40351893720200 r=202739815304 n-r=40149153904896 n+r=40554633535504
Rare numbers of length 15
n=200142385731002 r=200137583241002 n-r=4802490000 n+r=400279968972004
n=221462345754122 r=221457543264122 n-r=4802490000 n+r=442919889018244
n=816984566129618 r=816921665489618 n-r=62900640000 n+r=1633906231619236
n=245518996076442 r=244670699815542 n-r=848296260900 n+r=490189695891984
n=204238494066002 r=200660494832402 n-r=3577999233600 n+r=404898988898404
n=248359494187442 r=244781494953842 n-r=3577999233600 n+r=493140989141284
n=244062891224042 r=240422198260442 n-r=3640692963600 n+r=484485089484484
n=403058392434500 r=5434293850304 n-r=397624098584196 n+r=408492686284804
n=441054594034340 r=43430495450144 n-r=397624098584196 n+r=484485089484484
Rare numbers of length 16
n=2133786945766212 r=2126675496873312 n-r=7111448892900 n+r=4260462442639524
n=2135568943984212 r=2124893498655312 n-r=10675445328900 n+r=4260462442639524
n=8191154686620818 r=8180266864511918 n-r=10887822108900 n+r=16371421551132736
n=8191156864620818 r=8180264686511918 n-r=10892178108900 n+r=16371421551132736
n=2135764587964212 r=2124697854675312 n-r=11066733288900 n+r=4260462442639524
n=2135786765764212 r=2124675676875312 n-r=11111088888900 n+r=4260462442639524
n=8191376864400818 r=8180044686731918 n-r=11332177668900 n+r=16371421551132736
n=2078311262161202 r=2021612621138702 n-r=56698641022500 n+r=4099923883299904
n=4135786945764210 r=124675496875314 n-r=4011111448888896 n+r=4260462442639524
n=6889765708183410 r=143818075679886 n-r=6745947632503524 n+r=7033583783863296
n=8052956026592517 r=7152956206592508 n-r=899999820000009 n+r=15205912233185025
n=8052956206592517 r=7152956026592508 n-r=900000180000009 n+r=15205912233185025
n=8650327689541457 r=7541459867230568 n-r=1108867822310889 n+r=16191787556772025
n=8650349867341457 r=7541437689430568 n-r=1108912177910889 n+r=16191787556772025
n=6157577986646405 r=5046466897757516 n-r=1111111088888889 n+r=11204044884403921

real    7m6.890s
user    7m6.614s
sys     0m0.156s

Rare numbers of length 17
n=86965750494756968 r=86965749405756968 n-r=1089000000 n+r=173931499900513936
n=22542040692914522 r=22541929604024522 n-r=111088890000 n+r=45083970296939044
n=67725910561765640 r=4656716501952776 n-r=63069194059812864 n+r=72382627063718416

real    44m23.656s
user    44m22.729s
sys     0m0.401s

10 to 19 digits

The following is a faster implementation of the algorithm. It compiles with both g++ and clang++, both produce the correct output but g++ takes 3 times as long. I expect some difference but 300%, come on!!! It will not work for lengths less than 10 or greater than 19. <lang cpp> // Rare Numbers : Nigel Galloway - December 20th., 2019

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<10){
     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; ++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{1},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; ++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);

} </lang>

Output:

Rare 10
n=2022652202 r=2022562202 n-r=90000 n+r=4045214404
n=2042832002 r=2002382402 n-r=40449600 n+r=4045214404
Rare 11
Rare 12
n=872546974178 r=871479645278 n-r=1067328900 n+r=1744026619456
n=872568754178 r=871457865278 n-r=1110888900 n+r=1744026619456
n=868591084757 r=757480195868 n-r=111110888889 n+r=1626071280625
Rare 13
n=6979302951885 r=5881592039796 n-r=1097710912089 n+r=12860894991681
Rare 14
n=20331657922202 r=20222975613302 n-r=108682308900 n+r=40554633535504
n=20331875722202 r=20222757813302 n-r=109117908900 n+r=40554633535504
n=40351893720200 r=202739815304 n-r=40149153904896 n+r=40554633535504
n=20313693904202 r=20240939631302 n-r=72754272900 n+r=40554633535504
n=20313839704202 r=20240793831302 n-r=73045872900 n+r=40554633535504
n=20333875702202 r=20220757833302 n-r=113117868900 n+r=40554633535504
n=40313893704200 r=240739831304 n-r=40073153872896 n+r=40554633535504
Rare 15
n=245518996076442 r=244670699815542 n-r=848296260900 n+r=490189695891984
n=204238494066002 r=200660494832402 n-r=3577999233600 n+r=404898988898404
n=248359494187442 r=244781494953842 n-r=3577999233600 n+r=493140989141284
n=200142385731002 r=200137583241002 n-r=4802490000 n+r=400279968972004
n=221462345754122 r=221457543264122 n-r=4802490000 n+r=442919889018244
n=441054594034340 r=43430495450144 n-r=397624098584196 n+r=484485089484484
n=403058392434500 r=5434293850304 n-r=397624098584196 n+r=408492686284804
n=244062891224042 r=240422198260442 n-r=3640692963600 n+r=484485089484484
n=816984566129618 r=816921665489618 n-r=62900640000 n+r=1633906231619236
Rare 16
n=2135568943984212 r=2124893498655312 n-r=10675445328900 n+r=4260462442639524
n=2078311262161202 r=2021612621138702 n-r=56698641022500 n+r=4099923883299904
n=8191154686620818 r=8180266864511918 n-r=10887822108900 n+r=16371421551132736
n=8191156864620818 r=8180264686511918 n-r=10892178108900 n+r=16371421551132736
n=6889765708183410 r=143818075679886 n-r=6745947632503524 n+r=7033583783863296
n=2135764587964212 r=2124697854675312 n-r=11066733288900 n+r=4260462442639524
n=2135786765764212 r=2124675676875312 n-r=11111088888900 n+r=4260462442639524
n=2133786945766212 r=2126675496873312 n-r=7111448892900 n+r=4260462442639524
n=4135786945764210 r=124675496875314 n-r=4011111448888896 n+r=4260462442639524
n=8191376864400818 r=8180044686731918 n-r=11332177668900 n+r=16371421551132736
n=8650327689541457 r=7541459867230568 n-r=1108867822310889 n+r=16191787556772025
n=8650349867341457 r=7541437689430568 n-r=1108912177910889 n+r=16191787556772025
n=8052956026592517 r=7152956206592508 n-r=899999820000009 n+r=15205912233185025
n=8052956206592517 r=7152956026592508 n-r=900000180000009 n+r=15205912233185025
n=6157577986646405 r=5046466897757516 n-r=1111111088888889 n+r=11204044884403921

real    0m42.311s
user    0m42.297s
sys     0m0.004s

Rare 17
n=67725910561765640 r=4656716501952776 n-r=63069194059812864 n+r=72382627063718416
n=86965750494756968 r=86965749405756968 n-r=1089000000 n+r=173931499900513936
n=22542040692914522 r=22541929604024522 n-r=111088890000 n+r=45083970296939044

real    3m22.854s
user    3m22.742s
sys     0m0.020s

Rare 18
n=865721270017296468 r=864692710072127568 n-r=1028559945168900 n+r=1730413980089424036
n=297128548234950692 r=296059432845821792 n-r=1069115389128900 n+r=593187981080772484
n=297128722852950692 r=296059258227821792 n-r=1069464625128900 n+r=593187981080772484
n=898907259301737498 r=894737103952709898 n-r=4170155349027600 n+r=1793644363254447396
n=811865096390477018 r=810774093690568118 n-r=1091002699908900 n+r=1622639190081045136
n=284684666566486482 r=284684665666486482 n-r=900000000 n+r=569369332232972964
n=225342456863243522 r=225342368654243522 n-r=88209000000 n+r=450684825517487044
n=225342458663243522 r=225342366854243522 n-r=91809000000 n+r=450684825517487044
n=225342478643243522 r=225342346874243522 n-r=131769000000 n+r=450684825517487044
n=284684868364486482 r=284684463868486482 n-r=404496000000 n+r=569369332232972964
n=297148324656930692 r=296039656423841792 n-r=1108668233088900 n+r=593187981080772484
n=297148546434930692 r=296039434645841792 n-r=1109111789088900 n+r=593187981080772484
n=871975098681469178 r=871964186890579178 n-r=10911790890000 n+r=1743939285572048356
n=497168548234910690 r=96019432845861794 n-r=401149115389048896 n+r=593187981080772484
n=633488632647994145 r=541499746236884336 n-r=91988886411109809 n+r=1174988378884878481
n=631688638047992345 r=543299740836886136 n-r=88388897211106209 n+r=1174988378884878481
n=653488856225994125 r=521499522658884356 n-r=131989333567109769 n+r=1174988378884878481
n=633288858025996145 r=541699520858882336 n-r=91589337167113809 n+r=1174988378884878481
n=619631153042134925 r=529431240351136916 n-r=90199912690998009 n+r=1149062393393271841
n=619431353040136925 r=529631040353134916 n-r=89800312687002009 n+r=1149062393393271841

real    9m23.649s
user    9m23.510s
sys     0m0.005s

Rare 19
n=6988066446726832640 r=462386276446608896 n-r=6525680170280223744 n+r=7450452723173441536
n=2060303819041450202 r=2020541409183030602 n-r=39762409858419600 n+r=4080845228224480804
n=2702373360882732072 r=2702372880633732072 n-r=480249000000 n+r=5404746241516464144
n=2551755006254571552 r=2551754526005571552 n-r=480249000000 n+r=5103509532260143104
n=8066308349502036608 r=8066302059438036608 n-r=6290064000000 n+r=16132610408940073216
n=2825378427312735282 r=2825372137248735282 n-r=6290064000000 n+r=5650750564561470564
n=8320411466598809138 r=8319088956641140238 n-r=1322509957668900 n+r=16639500423239949376
n=2042401829204402402 r=2042044029281042402 n-r=357799923360000 n+r=4084445858485444804
n=2420424089100600242 r=2420060019804240242 n-r=364069296360000 n+r=4840484108904840484
n=8197906905009010818 r=8180109005096097918 n-r=17797899912912900 n+r=16378015910105108736
n=8200756128308135597 r=7955318038216570028 n-r=245438090091565569 n+r=16156074166524705625
n=6531727101458000045 r=5400008541017271356 n-r=1131718560440728689 n+r=11931735642475271401

real    65m40.434s
user    65m35.421s
sys     0m2.127s

20+ digits

<lang cpp> // Rare Numbers : Nigel Galloway - December 20th., 2019

include <iostream>
include <functional>
include <bitset>
include <gmpxx.h>

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}; const bool izRev(const mpz_class n,const mpz_class i,const mpz_class g){return (i/n!=g%10)? false : (n<2)? true : izRev(n/10,i%n,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<mpz_class>even{};
 std::vector<mpz_class>odd{};
 nLH(std::pair<Z1,std::vector<std::pair<long,long>>> e){auto [n,g]=e; mpz_t w,l,y; mpz_inits(w,l,y,NULL); mpz_set_si(w,pow10[4]); 
  while (auto i=n()){for(auto [ng,gg]:g){if((ng>0)|(*i>0)){mpz_set_si(y,gg+*i); mpz_addmul_ui(y,w,ng);
  if(mpz_perfect_square_p(y)) (gg%2==0)? even.push_back(mpz_class(y)) : odd.push_back(mpz_class(y));}}} mpz_clears(w,l,y,NULL);}

}; class Rare{

 mpz_class r,z,p;
 long acc{0};
 const std::bitset<10000>bs;
 const std::pair<Z1,std::vector<std::pair<long,long>>> makeL(const int n){ //std::cout<<"Making L"<<std::endl;
   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<10){
     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; ++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]); //std::cout<<"Making H"<<std::endl;
   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{1},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; ++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;})()}{

 mpz_ui_pow_ui(p.get_mpz_t(),10,n-1);
 std::cout<<"Rare "<<n<<std::endl;
 for(auto l:L.even) for(auto h:H.even){r=(h-l)/2; z=h-r; if(izRev(p,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) {r=(h-l)/2; z=h-r; if(izRev(p,r,z)) std::cout<<"n="<<z<<" r="<<r<<" n-r="<<l<<" n+r="<<h<<std::endl;}

}}; int main(){

 Rare(20);

} </lang>

Output:

Rare 10
n=2022652202 r=2022562202 n-r=90000 n+r=4045214404
n=2042832002 r=2002382402 n-r=40449600 n+r=4045214404
Rare 11
Rare 12
n=872546974178 r=871479645278 n-r=1067328900 n+r=1744026619456
n=872568754178 r=871457865278 n-r=1110888900 n+r=1744026619456
n=868591084757 r=757480195868 n-r=111110888889 n+r=1626071280625
Rare 13
n=6979302951885 r=5881592039796 n-r=1097710912089 n+r=12860894991681
Rare 14
n=20331657922202 r=20222975613302 n-r=108682308900 n+r=40554633535504
n=20331875722202 r=20222757813302 n-r=109117908900 n+r=40554633535504
n=40351893720200 r=202739815304 n-r=40149153904896 n+r=40554633535504
n=20313693904202 r=20240939631302 n-r=72754272900 n+r=40554633535504
n=20313839704202 r=20240793831302 n-r=73045872900 n+r=40554633535504
n=20333875702202 r=20220757833302 n-r=113117868900 n+r=40554633535504
n=40313893704200 r=240739831304 n-r=40073153872896 n+r=40554633535504
Rare 15
n=245518996076442 r=244670699815542 n-r=848296260900 n+r=490189695891984
n=204238494066002 r=200660494832402 n-r=3577999233600 n+r=404898988898404
n=248359494187442 r=244781494953842 n-r=3577999233600 n+r=493140989141284
n=200142385731002 r=200137583241002 n-r=4802490000 n+r=400279968972004
n=221462345754122 r=221457543264122 n-r=4802490000 n+r=442919889018244
n=441054594034340 r=43430495450144 n-r=397624098584196 n+r=484485089484484
n=403058392434500 r=5434293850304 n-r=397624098584196 n+r=408492686284804
n=244062891224042 r=240422198260442 n-r=3640692963600 n+r=484485089484484
n=816984566129618 r=816921665489618 n-r=62900640000 n+r=1633906231619236
Rare 16
n=2135568943984212 r=2124893498655312 n-r=10675445328900 n+r=4260462442639524
n=2078311262161202 r=2021612621138702 n-r=56698641022500 n+r=4099923883299904
n=8191154686620818 r=8180266864511918 n-r=10887822108900 n+r=16371421551132736
n=8191156864620818 r=8180264686511918 n-r=10892178108900 n+r=16371421551132736
n=6889765708183410 r=143818075679886 n-r=6745947632503524 n+r=7033583783863296
n=2135764587964212 r=2124697854675312 n-r=11066733288900 n+r=4260462442639524
n=2135786765764212 r=2124675676875312 n-r=11111088888900 n+r=4260462442639524
n=2133786945766212 r=2126675496873312 n-r=7111448892900 n+r=4260462442639524
n=4135786945764210 r=124675496875314 n-r=4011111448888896 n+r=4260462442639524
n=8191376864400818 r=8180044686731918 n-r=11332177668900 n+r=16371421551132736
n=8650327689541457 r=7541459867230568 n-r=1108867822310889 n+r=16191787556772025
n=8650349867341457 r=7541437689430568 n-r=1108912177910889 n+r=16191787556772025
n=8052956026592517 r=7152956206592508 n-r=899999820000009 n+r=15205912233185025
n=8052956206592517 r=7152956026592508 n-r=900000180000009 n+r=15205912233185025
n=6157577986646405 r=5046466897757516 n-r=1111111088888889 n+r=11204044884403921

real    2m10.594s
user    2m10.563s
sys     0m0.004s

Rare 17
n=67725910561765640 r=4656716501952776 n-r=63069194059812864 n+r=72382627063718416
n=86965750494756968 r=86965749405756968 n-r=1089000000 n+r=173931499900513936
n=22542040692914522 r=22541929604024522 n-r=111088890000 n+r=45083970296939044

real    10m8.790s
user    10m8.648s
sys     0m0.009s

Rare 18
n=865721270017296468 r=864692710072127568 n-r=1028559945168900 n+r=1730413980089424036
n=297128548234950692 r=296059432845821792 n-r=1069115389128900 n+r=593187981080772484
n=297128722852950692 r=296059258227821792 n-r=1069464625128900 n+r=593187981080772484
n=898907259301737498 r=894737103952709898 n-r=4170155349027600 n+r=1793644363254447396
n=811865096390477018 r=810774093690568118 n-r=1091002699908900 n+r=1622639190081045136
n=284684666566486482 r=284684665666486482 n-r=900000000 n+r=569369332232972964
n=225342456863243522 r=225342368654243522 n-r=88209000000 n+r=450684825517487044
n=225342458663243522 r=225342366854243522 n-r=91809000000 n+r=450684825517487044
n=225342478643243522 r=225342346874243522 n-r=131769000000 n+r=450684825517487044
n=284684868364486482 r=284684463868486482 n-r=404496000000 n+r=569369332232972964
n=297148324656930692 r=296039656423841792 n-r=1108668233088900 n+r=593187981080772484
n=297148546434930692 r=296039434645841792 n-r=1109111789088900 n+r=593187981080772484
n=871975098681469178 r=871964186890579178 n-r=10911790890000 n+r=1743939285572048356
n=497168548234910690 r=96019432845861794 n-r=401149115389048896 n+r=593187981080772484
n=633488632647994145 r=541499746236884336 n-r=91988886411109809 n+r=1174988378884878481
n=631688638047992345 r=543299740836886136 n-r=88388897211106209 n+r=1174988378884878481
n=653488856225994125 r=521499522658884356 n-r=131989333567109769 n+r=1174988378884878481
n=633288858025996145 r=541699520858882336 n-r=91589337167113809 n+r=1174988378884878481
n=619631153042134925 r=529431240351136916 n-r=90199912690998009 n+r=1149062393393271841
n=619431353040136925 r=529631040353134916 n-r=89800312687002009 n+r=1149062393393271841

real    30m37.306s
user    30m36.730s
sys     0m0.000s

Rare 20
n=22134434735752443122 r=22134425753743443122 n-r=8982009000000 n+r=44268860489495886244
n=22134434753752443122 r=22134425735743443122 n-r=9018009000000 n+r=44268860489495886244
n=22134436953532443122 r=22134423535963443122 n-r=13417569000000 n+r=44268860489495886244
n=65459144877856561700 r=716565877844195456 n-r=64742579000012366244 n+r=66175710755700757156
n=22136414517954423122 r=22132445971541463122 n-r=3968546412960000 n+r=44268860489495886244
n=22136414971554423122 r=22132445517941463122 n-r=3969453612960000 n+r=44268860489495886244
n=22136456771730423122 r=22132403717765463122 n-r=4053053964960000 n+r=44268860489495886244
n=61952807156239928885 r=58882993265170825916 n-r=3069813891069102969 n+r=120835800421410754801
n=61999171315484316965 r=56961348451317199916 n-r=5037822864167117049 n+r=118960519766801516881

real    627m18.505s
user    626m57.689s
sys     0m1.796s

C#

Translation of: Go

Converted to unsigned longs in order to reach 19 digits. <lang csharp>using System; using System.Collections.Generic; using System.Linq; using static System.Console; using UI = System.UInt64; using LST = System.Collections.Generic.List<System.Collections.Generic.List<sbyte>>; using Lst = System.Collections.Generic.List<sbyte>; using DT = System.DateTime;

class Program {

   const sbyte MxD = 19;

   public struct term { public UI coeff; public sbyte a, b;
       public term(UI c, int a_, int b_) { coeff = c; a = (sbyte)a_; b = (sbyte)b_; } }

   static int[] digs;   static List<UI> res;   static sbyte count = 0;
   static DT st; static List<List<term>> tLst; static List<LST> lists;
   static Dictionary<int, LST> fml, dmd; static Lst dl, zl, el, ol, il;
   static bool odd; static int nd, nd2; static LST ixs;
   static int[] cnd, di; static LST dis; static UI Dif;

   // converts digs array to the "difference"
   static UI ToDif() { UI r = 0; for (int i = 0; i < digs.Length; i++)
           r = r * 10 + (uint)digs[i]; return r; }
   
   // converts digs array to the "sum"
   static UI ToSum() { UI r = 0; for (int i = digs.Length - 1; i >= 0; i--)
           r = r * 10 + (uint)digs[i]; return Dif + (r << 1); }

   // determines if the nmbr is square or not
   static bool IsSquare(UI nmbr) { if ((0x202021202030213 & (1 << (int)(nmbr & 63))) != 0)
       { UI r = (UI)Math.Sqrt((double)nmbr); return r * r == nmbr; } return false; }

   // returns sequence of sbytes
   static Lst Seq(sbyte from, int to, sbyte stp = 1) { Lst res = new Lst();
       for (sbyte item = from; item <= to; item += stp) res.Add(item); return res; }

   // Recursive closure to generate (n+r) candidates from (n-r) candidates
   static void Fnpr(int lev) { if (lev == dis.Count) { digs[ixs[0][0]] = fml[cnd[0]][di[0]][0];
           digs[ixs[0][1]] = fml[cnd[0]][di[0]][1]; int le = di.Length, i = 1;
           if (odd) digs[nd >> 1] = di[--le]; foreach (sbyte d in di.Skip(1).Take(le - 1)) {
               digs[ixs[i][0]] = dmd[cnd[i]][d][0]; digs[ixs[i][1]] = dmd[cnd[i++]][d][1]; }
           if (!IsSquare(ToSum())) return; res.Add(ToDif()); WriteLine("{0,16:n0}{1,4}   ({2:n0})",
               (DT.Now - st).TotalMilliseconds, ++count, res.Last()); }
       else foreach (var n in dis[lev]) { di[lev] = n; Fnpr(lev + 1); } }

   // Recursive closure to generate (n-r) candidates with a given number of digits.
   static void Fnmr (LST list, int lev) { if (lev == list.Count) { Dif = 0; sbyte i = 0;
           foreach (var t in tLst[nd2]) { if (cnd[i] < 0) Dif -= t.coeff * (UI)(-cnd[i++]);
               else Dif += t.coeff * (UI)cnd[i++]; } if (Dif <= 0 || !IsSquare(Dif)) return;
           dis = new LST { Seq(0, fml[cnd[0]].Count - 1) };
           foreach (int ii in cnd.Skip(1)) dis.Add(Seq(0, dmd[ii].Count - 1));
           if (odd) dis.Add(il); di = new int[dis.Count]; Fnpr(0);
       } else foreach(sbyte n in list[lev]) { cnd[lev] = n; Fnmr(list, lev + 1); } }

   static void init() { UI pow = 1;
       // terms of (n-r) expression for number of digits from 2 to maxDigits
       tLst = new List<List<term>>(); foreach (int r in Seq(2, MxD)) {
           List<term> terms = new List<term>(); pow *= 10; UI p1 = pow, p2 = 1;
           for (int i1 = 0, i2 = r - 1; i1 < i2; i1++, i2--) {
               terms.Add(new term(p1 - p2, i1, i2)); p1 /= 10; p2 *= 10; }
           tLst.Add(terms); }
       //  map of first minus last digits for 'n' to pairs giving this value
       fml = new Dictionary<int, LST> {
           [0] = new LST { new Lst { 2, 2 }, new Lst { 8, 8 } },
           [1] = new LST { new Lst { 6, 5 }, new Lst { 8, 7 } },
           [4] = new LST { new Lst { 4, 0 } },
           [6] = new LST { new Lst { 6, 0 }, new Lst { 8, 2 } } };
       // map of other digit differences for 'n' to pairs giving this value
       dmd = new Dictionary<int, LST>();
       for (sbyte i = 0; i < 10; i++) for (sbyte j = 0, d = i; j < 10; j++, d--) {
               if (dmd.ContainsKey(d)) dmd[d].Add(new Lst { i, j });
               else dmd[d] = new LST { new Lst { i, j } }; }
       dl = Seq(-9, 9);    // all differences
       zl = Seq( 0, 0);    // zero differences only
       el = Seq(-8, 8, 2); // even differences only
       ol = Seq(-9, 9, 2); // odd differences only
       il = Seq( 0, 9); lists = new List<LST>();
       foreach (sbyte f in fml.Keys) lists.Add(new LST { new Lst { f } }); }

   static void Main(string[] args) { init(); res = new List<UI>(); st = DT.Now; count = 0;
       WriteLine("{0,5}{1,12}{2,4}{3,14}", "digs", "elapsed(ms)", "R/N", "Unordered Rare Numbers");
       for (nd = 2, nd2 = 0, odd = false; nd <= MxD; nd++, nd2++, odd = !odd) { digs = new int[nd];
           if (nd == 4) { lists[0].Add(zl); lists[1].Add(ol); lists[2].Add(el); lists[3].Add(ol); }
           else if (tLst[nd2].Count > lists[0].Count) foreach (LST list in lists) list.Add(dl);
           ixs = new LST(); 
           foreach (term t in tLst[nd2]) ixs.Add(new Lst { t.a, t.b });
           foreach (LST list in lists) { cnd = new int[list.Count]; Fnmr(list, 0); }
           WriteLine("  {0,2}  {1,10:n0}", nd, (DT.Now - st).TotalMilliseconds); }
       res.Sort();
       WriteLine("\nThe {0} rare numbers with up to {1} digits are:", res.Count, MxD);
       count = 0; foreach (var rare in res) WriteLine("{0,2}:{1,27:n0}", ++count, rare);
       if (System.Diagnostics.Debugger.IsAttached) ReadKey(); }

}</lang>

Output:

Results from a core i7-7700 @ 3.6Ghz. This C# version isn't as fast as the Go version using the same hardware. C# computes up to 17, 18 and 19 digits in under 9 minutes, 1 2/3 hours and over 2 1/2 hours respectively. (Go is about 6 minutes, 1 1/4 hours, and under 2 hours).

The long-to-ulong conversion isn't causing the reduced performance, C# has more overhead as compared to Go. This C# version can easily be converted to use BigIntegers to go beyond 19 digits, but becomes around eight times slower. (ugh!)

 digs elapsed(ms) R/N  Rare Numbers
              27   1   (65)
   2          28
   3          28
   4          29
   5          29
              29   2   (621,770)
   6          29
   7          30
   8          34
              34   3   (281,089,082)
   9          36
              36   4   (2,022,652,202)
              61   5   (2,042,832,002)
  10         121
  11         176
             448   6   (872,546,974,178)
             481   7   (872,568,754,178)
             935   8   (868,591,084,757)
  12       1,232
           1,577   9   (6,979,302,951,885)
  13       2,087
           6,274  10   (20,313,693,904,202)
           6,351  11   (20,313,839,704,202)
           8,039  12   (20,331,657,922,202)
           8,292  13   (20,331,875,722,202)
           9,000  14   (20,333,875,702,202)
          21,212  15   (40,313,893,704,200)
          21,365  16   (40,351,893,720,200)
  14      23,898
          23,964  17   (200,142,385,731,002)
          24,198  18   (221,462,345,754,122)
          27,241  19   (816,984,566,129,618)
          28,834  20   (245,518,996,076,442)
          29,074  21   (204,238,494,066,002)
          29,147  22   (248,359,494,187,442)
          29,476  23   (244,062,891,224,042)
          35,481  24   (403,058,392,434,500)
          35,721  25   (441,054,594,034,340)
  15      38,231
          92,116  26   (2,133,786,945,766,212)
         113,469  27   (2,135,568,943,984,212)
         116,787  28   (8,191,154,686,620,818)
         119,647  29   (8,191,156,864,620,818)
         120,912  30   (2,135,764,587,964,212)
         122,735  31   (2,135,786,765,764,212)
         127,126  32   (8,191,376,864,400,818)
         141,793  33   (2,078,311,262,161,202)
         179,832  34   (8,052,956,026,592,517)
         184,647  35   (8,052,956,206,592,517)
         221,279  36   (8,650,327,689,541,457)
         223,721  37   (8,650,349,867,341,457)
         225,520  38   (6,157,577,986,646,405)
         273,238  39   (4,135,786,945,764,210)
         312,969  40   (6,889,765,708,183,410)
  16     316,349
         322,961  41   (86,965,750,494,756,968)
         323,958  42   (22,542,040,692,914,522)
         502,805  43   (67,725,910,561,765,640)
  17     519,583
         576,058  44   (284,684,666,566,486,482)
         707,530  45   (225,342,456,863,243,522)
         756,188  46   (225,342,458,663,243,522)
         856,346  47   (225,342,478,643,243,522)
         928,546  48   (284,684,868,364,486,482)
       1,311,170  49   (871,975,098,681,469,178)
       2,031,664  50   (865,721,270,017,296,468)
       2,048,209  51   (297,128,548,234,950,692)
       2,057,281  52   (297,128,722,852,950,692)
       2,164,878  53   (811,865,096,390,477,018)
       2,217,508  54   (297,148,324,656,930,692)
       2,242,999  55   (297,148,546,434,930,692)
       2,576,805  56   (898,907,259,301,737,498)
       3,169,675  57   (631,688,638,047,992,345)
       3,200,223  58   (619,431,353,040,136,925)
       3,482,517  59   (619,631,153,042,134,925)
       3,550,566  60   (633,288,858,025,996,145)
       3,623,653  61   (633,488,632,647,994,145)
       4,605,503  62   (653,488,856,225,994,125)
       5,198,241  63   (497,168,548,234,910,690)
  18   6,028,721
       6,130,826  64   (2,551,755,006,254,571,552)
       6,152,283  65   (2,702,373,360,882,732,072)
       6,424,945  66   (2,825,378,427,312,735,282)
       6,447,566  67   (8,066,308,349,502,036,608)
       6,677,925  68   (2,042,401,829,204,402,402)
       6,725,119  69   (2,420,424,089,100,600,242)
       6,843,016  70   (8,320,411,466,598,809,138)
       7,161,527  71   (8,197,906,905,009,010,818)
       7,198,112  72   (2,060,303,819,041,450,202)
       7,450,028  73   (8,200,756,128,308,135,597)
       7,881,502  74   (6,531,727,101,458,000,045)
       9,234,318  75   (6,988,066,446,726,832,640)
  19   9,394,513

The 75 rare numbers with up to 19 digits are:
 1:                         65
 2:                    621,770
 3:                281,089,082
 4:              2,022,652,202
 5:              2,042,832,002
 6:            868,591,084,757
 7:            872,546,974,178
 8:            872,568,754,178
 9:          6,979,302,951,885
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
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
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
41:     22,542,040,692,914,522
42:     67,725,910,561,765,640
43:     86,965,750,494,756,968
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
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

D

Translation of: Go

Scaled down from the full duration showed in the go example because I got impatient and have not spent time figuring out where the inefficeny is. <lang d>import std.algorithm; import std.array; import std.conv; import std.datetime.stopwatch; import std.math; import std.stdio;

struct Term {

   ulong coeff;
   byte ix1, ix2;

}

enum maxDigits = 16;

ulong toUlong(byte[] digits, bool reverse) {

   ulong sum = 0;
   if (reverse) {
       for (int i = digits.length - 1; i >= 0; --i) {
           sum = sum * 10 + digits[i];
       }
   } else {
       for (size_t i = 0; i < digits.length; ++i) {
           sum = sum * 10 + digits[i];
       }
   }
   return sum;

}

bool isSquare(ulong n) {

   if ((0x202021202030213 & (1 << (n & 63))) != 0) {
       auto root = cast(ulong)sqrt(cast(double)n);
       return root * root == n;
   }
   return false;

}

byte[] seq(byte from, byte to, byte step) {

   byte[] res;
   for (auto i = from; i <= to; i += step) {
       res ~= i;
   }
   return res;

}

string commatize(ulong n) {

   auto s = n.to!string;
   auto le = s.length;
   for (int i = le - 3; i >= 1; i -= 3) {
       s = s[0..i] ~ "," ~ s[i..$];
   }
   return s;

}

void main() {

   auto sw = StopWatch(AutoStart.yes);
   ulong pow = 1;
   writeln("Aggregate timings to process all numbers up to:");
   // terms of (n-r) expression for number of digits from 2 to maxDigits
   Term[][] allTerms = uninitializedArray!(Term[][])(maxDigits - 1);
   for (auto r = 2; r <= maxDigits; r++) {
       Term[] terms;
       pow *= 10;
       ulong pow1 = pow;
       ulong pow2 = 1;
       byte i1 = 0;
       byte i2 = cast(byte)(r - 1);
       while (i1 < i2) {
           terms ~= Term(pow1 - pow2, i1, i2);

           pow1 /= 10;
           pow2 *= 10;

           i1++;
           i2--;
       }
       allTerms[r - 2] = terms;
   }
   //  map of first minus last digits for 'n' to pairs giving this value
   byte[][][byte] fml = [
       0: [[2, 2], [8, 8]],
       1: [[6, 5], [8, 7]],
       4: 4, 0,
       6: [[6, 0], [8, 2]]
   ];
   // map of other digit differences for 'n' to pairs giving this value
   byte[][][byte] dmd;
   for (byte i = 0; i < 100; i++) {
       byte[] a = [i / 10, i % 10];
       auto d = a[0] - a[1];
       dmd[cast(byte)d] ~= a;
   }
   byte[] fl = [0, 1, 4, 6];
   auto dl = seq(-9, 9, 1);    // all differences
   byte[] zl = [0];            // zero diferences only
   auto el = seq(-8, 8, 2);    // even differences only
   auto ol = seq(-9, 9, 2);    // odd differences only
   auto il = seq(0, 9, 1);
   ulong[] rares;
   byte[][][] lists = uninitializedArray!(byte[][][])(4);
   foreach (i, f; fl) {
       lists[i] = f;
   }
   byte[] digits;
   int count = 0;

   // Recursive closure to generate (n+r) candidates from (n-r) candidates
   // and hence find Rare numbers with a given number of digits.
   void fnpr(byte[] cand, byte[] di, byte[][] dis, byte[][] indicies, ulong nmr, int nd, int level) {
       if (level == dis.length) {
           digits[indicies[0][0]] = fml[cand[0]][di[0]][0];
           digits[indicies[0][1]] = fml[cand[0]][di[0]][1];
           auto le = di.length;
           if (nd % 2 == 1) {
               le--;
               digits[nd / 2] = di[le];
           }
           foreach (i, d; di[1..le]) {
               digits[indicies[i + 1][0]] = dmd[cand[i + 1]][d][0];
               digits[indicies[i + 1][1]] = dmd[cand[i + 1]][d][1];
           }
           auto r = toUlong(digits, true);
           auto npr = nmr + 2 * r;
           if (!isSquare(npr)) {
               return;
           }
           count++;
           writef("     R/N %2d:", count);
           auto ms = sw.peek();
           writef("  %9s", ms);
           auto n = toUlong(digits, false);
           writef("  (%s)\n", commatize(n));
           rares ~= n;
       } else {
           foreach (num; dis[level]) {
               di[level] = num;
               fnpr(cand, di, dis, indicies, nmr, nd, level + 1);
           }
       }
   }

   // Recursive closure to generate (n-r) candidates with a given number of digits.
   void fnmr(byte[] cand, byte[][] list, byte[][] indicies, int nd, int level) {
       if (level == list.length) {
           ulong nmr, nmr2;
           foreach (i, t; allTerms[nd - 2]) {
               if (cand[i] >= 0) {
                   nmr += t.coeff * cand[i];
               } else {
                   nmr2 += t.coeff * -cast(int)(cand[i]);
                   if (nmr >= nmr2) {
                       nmr -= nmr2;
                       nmr2 = 0;
                   } else {
                       nmr2 -= nmr;
                       nmr = 0;
                   }
               }
           }
           if (nmr2 >= nmr) {
               return;
           }
           nmr -= nmr2;
           if (!isSquare(nmr)) {
               return;
           }
           byte[][] dis;
           dis ~= seq(0, cast(byte)(fml[cand[0]].length - 1), 1);
           for (auto i = 1; i < cand.length; i++) {
               dis ~= seq(0, cast(byte)(dmd[cand[i]].length - 1), 1);
           }
           if (nd % 2 == 1) {
               dis ~= il;
           }
           byte[] di = uninitializedArray!(byte[])(dis.length);
           fnpr(cand, di, dis, indicies, nmr, nd, 0);
       } else {
           foreach (num; list[level]) {
               cand[level] = num;
               fnmr(cand, list, indicies, nd, level + 1);
           }
       }
   }

   for (int nd = 2; nd <= maxDigits; nd++) {
       digits = uninitializedArray!(byte[])(nd);
       if (nd == 4) {
           lists[0] ~= zl;
           lists[1] ~= ol;
           lists[2] ~= el;
           lists[3] ~= ol;
       } else if (allTerms[nd - 2].length > lists[0].length) {
           for (int i = 0; i < 4; i++) {
               lists[i] ~= dl;
           }
       }
       byte[][] indicies;
       foreach (t; allTerms[nd - 2]) {
           indicies ~= [t.ix1, t.ix2];
       }
       foreach (list; lists) {
           byte[] cand = uninitializedArray!(byte[])(list.length);
           fnmr(cand, list, indicies, nd, 0);
       }
       auto ms = sw.peek();
       writefln("  %2d digits:  %9s", nd, ms);
   }

   rares.sort;
   writefln("\nThe rare numbers with up to %d digits are:", maxDigits);
   foreach (i, rare; rares) {
       writefln("  %2d:  %25s", i + 1, commatize(rare));
   }

}</lang>

Output:

Aggregate timings to process all numbers up to:
     R/N  1:  183 ╬╝s and 2 hnsecs  (65)
   2 digits:  193 ╬╝s and 8 hnsecs
   3 digits:  301 ╬╝s and 3 hnsecs
   4 digits:  380 ╬╝s and 1 hnsec
   5 digits:  447 ╬╝s
     R/N  2:  732 ╬╝s and 9 hnsecs  (621,770)
   6 digits:  767 ╬╝s and 1 hnsec
   7 digits:  1 ms, 291 ╬╝s, and 5 hnsecs
   8 digits:  5 ms, 602 ╬╝s, and 2 hnsecs
     R/N  3:  5 ms, 900 ╬╝s, and 2 hnsecs  (281,089,082)
   9 digits:  7 ms, 537 ╬╝s, and 1 hnsec
     R/N  4:  7 ms, 869 ╬╝s, and 5 hnsecs  (2,022,652,202)
     R/N  5:  32 ms, 826 ╬╝s, and 4 hnsecs  (2,042,832,002)
  10 digits:  96 ms, 422 ╬╝s, and 3 hnsecs
  11 digits:  161 ms, 218 ╬╝s, and 4 hnsecs
     R/N  6:  468 ms, 23 ╬╝s, and 9 hnsecs  (872,546,974,178)
     R/N  7:  506 ms, 702 ╬╝s, and 3 hnsecs  (872,568,754,178)
     R/N  8:  1 sec, 39 ms, 845 ╬╝s, and 6 hnsecs  (868,591,084,757)
  12 digits:  1 sec, 437 ms, 602 ╬╝s, and 8 hnsecs
     R/N  9:  1 sec, 835 ms, 95 ╬╝s, and 6 hnsecs  (6,979,302,951,885)
  13 digits:  2 secs, 487 ms, 165 ╬╝s, and 9 hnsecs
     R/N 10:  7 secs, 241 ms, 437 ╬╝s, and 1 hnsec  (20,313,693,904,202)
     R/N 11:  7 secs, 330 ms, 171 ╬╝s, and 2 hnsecs  (20,313,839,704,202)
     R/N 12:  9 secs, 290 ms, 907 ╬╝s, and 3 hnsecs  (20,331,657,922,202)
     R/N 13:  9 secs, 582 ms, 920 ╬╝s, and 5 hnsecs  (20,331,875,722,202)
     R/N 14:  10 secs, 383 ms, 769 ╬╝s, and 1 hnsec  (20,333,875,702,202)
     R/N 15:  25 secs, 835 ms, and 933 ╬╝s  (40,313,893,704,200)
     R/N 16:  26 secs, 14 ms, 774 ╬╝s, and 4 hnsecs  (40,351,893,720,200)
  14 digits:  30 secs, 110 ms, 971 ╬╝s, and 7 hnsecs
     R/N 17:  30 secs, 216 ms, 437 ╬╝s, and 3 hnsecs  (200,142,385,731,002)
     R/N 18:  30 secs, 489 ms, 719 ╬╝s, and 2 hnsecs  (221,462,345,754,122)
     R/N 19:  34 secs, 83 ms, 642 ╬╝s, and 9 hnsecs  (816,984,566,129,618)
     R/N 20:  35 secs, 971 ms, 413 ╬╝s, and 3 hnsecs  (245,518,996,076,442)
     R/N 21:  36 secs, 250 ms, 787 ╬╝s, and 8 hnsecs  (204,238,494,066,002)
     R/N 22:  36 secs, 332 ms, 714 ╬╝s, and 2 hnsecs  (248,359,494,187,442)
     R/N 23:  36 secs, 696 ms, 902 ╬╝s, and 2 hnsecs  (244,062,891,224,042)
     R/N 24:  44 secs, 896 ms, and 665 ╬╝s  (403,058,392,434,500)
     R/N 25:  45 secs, 181 ms, 141 ╬╝s, and 5 hnsecs  (441,054,594,034,340)
  15 digits:  49 secs, 315 ms, 407 ╬╝s, and 4 hnsecs
     R/N 26:  1 minute, 55 secs, 748 ms, 43 ╬╝s, and 4 hnsecs  (2,133,786,945,766,212)
     R/N 27:  2 minutes, 21 secs, 484 ms, 683 ╬╝s, and 7 hnsecs  (2,135,568,943,984,212)
     R/N 28:  2 minutes, 25 secs, 438 ms, 771 ╬╝s, and 7 hnsecs  (8,191,154,686,620,818)
     R/N 29:  2 minutes, 28 secs, 883 ms, 999 ╬╝s, and 6 hnsecs  (8,191,156,864,620,818)
     R/N 30:  2 minutes, 30 secs, 410 ms, and 831 ╬╝s  (2,135,764,587,964,212)
     R/N 31:  2 minutes, 32 secs, 594 ms, 842 ╬╝s, and 1 hnsec  (2,135,786,765,764,212)
     R/N 32:  2 minutes, 37 secs, 880 ms, 100 ╬╝s, and 5 hnsecs  (8,191,376,864,400,818)
     R/N 33:  2 minutes, 55 secs, 943 ms, 190 ╬╝s, and 5 hnsecs  (2,078,311,262,161,202)
     R/N 34:  3 minutes, 49 secs, 750 ms, 39 ╬╝s, and 5 hnsecs  (8,052,956,026,592,517)
     R/N 35:  3 minutes, 55 secs, 554 ms, 720 ╬╝s, and 1 hnsec  (8,052,956,206,592,517)
     R/N 36:  4 minutes, 41 secs, 59 ms, 309 ╬╝s, and 4 hnsecs  (8,650,327,689,541,457)
     R/N 37:  4 minutes, 43 secs, 951 ms, and 206 ╬╝s  (8,650,349,867,341,457)
     R/N 38:  4 minutes, 46 secs, 85 ms, 249 ╬╝s, and 7 hnsecs  (6,157,577,986,646,405)
     R/N 39:  5 minutes, 59 secs, 80 ms, 228 ╬╝s, and 5 hnsecs  (4,135,786,945,764,210)
     R/N 40:  7 minutes, 10 secs, 573 ms, 592 ╬╝s, and 2 hnsecs  (6,889,765,708,183,410)
  16 digits:  7 minutes, 16 secs, 827 ms, 76 ╬╝s, and 4 hnsecs

The rare numbers with up to 16 digits are:
   1:                         65
   2:                    621,770
   3:                281,089,082
   4:              2,022,652,202
   5:              2,042,832,002
   6:            868,591,084,757
   7:            872,546,974,178
   8:            872,568,754,178
   9:          6,979,302,951,885
  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
  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
  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

F#

The Function

This solution demonstrates the concept described in Talk:Rare_numbers#30_mins_not_30_years. It doesn't use Cartesian_product_of_two_or_more_lists#Extra_Credit <lang fsharp> // Find all Rare numbers with a digits. Nigel Galloway: September 18th., 2019. let rareNums a=

 let tN=set[1L;4L;5L;6L;9L]
 let izPS g=let n=(float>>sqrt>>int64)g in n*n=g
 let n=[for n in [0..a/2-1] do yield ((pown 10L (a-n-1))-(pown 10L n))]|>List.rev
 let rec fN i g e=seq{match e with 0->yield g |e->for n in i do yield! fN [-9L..9L] (n::g) (e-1)}|>Seq.filter(fun g->let g=Seq.map2(*) n g|>Seq.sum in g>0L && izPS g)
 let rec fG n i g e l=seq{
   match l with
    h::t->for l in max 0L (0L-h)..min 9L (9L-h) do if e>1L||l=0L||tN.Contains((2L*l+h)%10L) then yield! fG (n+l*e+(l+h)*g) (i+l*g+(l+h)*e) (g/10L) (e*10L) t
   |_->if n>(pown 10L (a-1)) then for l in (if a%2=0 then [0L] else [0L..9L]) do let g=l*(pown 10L (a/2)) in if izPS (n+i+2L*g) then yield (i+g,n+g)} 
 fN [0L..9L] [] (a/2) |> Seq.collect(List.rev >> fG 0L 0L (pown 10L (a-1)) 1L)

</lang>

43 down

<lang fsharp> let test n=

 let t = System.Diagnostics.Stopwatch.StartNew()
 for n in (rareNums n) do printfn "%A" n
 t.Stop()
 printfn "Elapsed Time: %d ms for length %d" t.ElapsedMilliseconds n

[2..17] |> Seq.iter test </lang>

Output:

(56L, 65L)
Elapsed Time: 31 ms for length 2
Elapsed Time: 0 ms for length 3
Elapsed Time: 0 ms for length 4
Elapsed Time: 0 ms for length 5
(77126L, 621770L)
Elapsed Time: 6 ms for length 6
Elapsed Time: 6 ms for length 7
Elapsed Time: 113 ms for length 8
(280980182L, 281089082L)
Elapsed Time: 72 ms for length 9
(2022562202L, 2022652202L)
(2002382402L, 2042832002L)
Elapsed Time: 1525 ms for length 10
Elapsed Time: 1351 ms for length 11
(871479645278L, 872546974178L)
(871457865278L, 872568754178L)
(757480195868L, 868591084757L)
Elapsed Time: 27990 ms for length 12
(5881592039796L, 6979302951885L)
Elapsed Time: 26051 ms for length 13
(20240939631302L, 20313693904202L)
(20240793831302L, 20313839704202L)
(20222975613302L, 20331657922202L)
(20222757813302L, 20331875722202L)
(20220757833302L, 20333875702202L)
(240739831304L, 40313893704200L)
(202739815304L, 40351893720200L)
Elapsed Time: 552922 ms for length 14
(200137583241002L, 200142385731002L)
(221457543264122L, 221462345754122L)
(816921665489618L, 816984566129618L)
(244670699815542L, 245518996076442L)
(200660494832402L, 204238494066002L)
(244781494953842L, 248359494187442L)
(240422198260442L, 244062891224042L)
(5434293850304L, 403058392434500L)
(43430495450144L, 441054594034340L)
Elapsed Time: 512282 ms for length 15
(2126675496873312L, 2133786945766212L)
(2124893498655312L, 2135568943984212L)
(8180266864511918L, 8191154686620818L)
(8180264686511918L, 8191156864620818L)
(2124697854675312L, 2135764587964212L)
(2124675676875312L, 2135786765764212L)
(8180044686731918L, 8191376864400818L)
(2021612621138702L, 2078311262161202L)
(7152956206592508L, 8052956026592517L)
(7152956026592508L, 8052956206592517L)
(7541459867230568L, 8650327689541457L)
(7541437689430568L, 8650349867341457L)
(5046466897757516L, 6157577986646405L)
(124675496875314L, 4135786945764210L)
(143818075679886L, 6889765708183410L)
Elapsed Time: 11568713 ms for length 16
(86965749405756968L, 86965750494756968L)
(22541929604024522L, 22542040692914522L)
(4656716501952776L, 67725910561765640L)
Elapsed Time: 11275839 ms for length 17

Go

This uses many of the hints within Shyam Sunder Gupta's webpage combined with Nigel Galloway's general approach (see Talk page) of working from (n-r) and deducing the Rare numbers with various numbers of digits from there.

As the algorithm used does not generate the Rare numbers in order, a sorted list is also printed. <lang go>package main

import (

   "fmt"
   "math"
   "sort"
   "time"

)

type term struct {

   coeff    uint64
   ix1, ix2 int8

}

const maxDigits = 19

func toUint64(digits []int8, reverse bool) uint64 {

   sum := uint64(0)
   if !reverse {
       for i := 0; i < len(digits); i++ {
           sum = sum*10 + uint64(digits[i])
       }
   } else {
       for i := len(digits) - 1; i >= 0; i-- {
           sum = sum*10 + uint64(digits[i])
       }
   }
   return sum

}

func isSquare(n uint64) bool {

   if 0x202021202030213&(1<<(n&63)) != 0 {
       root := uint64(math.Sqrt(float64(n)))
       return root*root == n
   }
   return false

}

func seq(from, to, step int8) []int8 {

   var res []int8
   for i := from; i <= to; i += step {
       res = append(res, i)
   }
   return res

}

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() {

   start := time.Now()
   pow := uint64(1)
   fmt.Println("Aggregate timings to process all numbers up to:")
   // terms of (n-r) expression for number of digits from 2 to maxDigits
   allTerms := make([][]term, maxDigits-1)
   for r := 2; r <= maxDigits; r++ {
       var terms []term
       pow *= 10
       pow1, pow2 := pow, uint64(1)
       for i1, i2 := int8(0), int8(r-1); i1 < i2; i1, i2 = i1+1, i2-1 {
           terms = append(terms, term{pow1 - pow2, i1, i2})
           pow1 /= 10
           pow2 *= 10
       }
       allTerms[r-2] = terms
   }
   //  map of first minus last digits for 'n' to pairs giving this value
   fml := map[int8][][]int8{
       0: {{2, 2}, {8, 8}},
       1: {{6, 5}, {8, 7}},
       4: Template:4, 0,
       6: {{6, 0}, {8, 2}},
   }
   // map of other digit differences for 'n' to pairs giving this value
   dmd := make(map[int8][][]int8)
   for i := int8(0); i < 100; i++ {
       a := []int8{i / 10, i % 10}
       d := a[0] - a[1]
       dmd[d] = append(dmd[d], a)
   }
   fl := []int8{0, 1, 4, 6}
   dl := seq(-9, 9, 1) // all differences
   zl := []int8{0}     // zero differences only
   el := seq(-8, 8, 2) // even differences only
   ol := seq(-9, 9, 2) // odd differences only
   il := seq(0, 9, 1)
   var rares []uint64
   lists := make([][][]int8, 4)
   for i, f := range fl {
       lists[i] = [][]int8Template:F
   }
   var digits []int8
   count := 0

   // Recursive closure to generate (n+r) candidates from (n-r) candidates
   // and hence find Rare numbers with a given number of digits.
   var fnpr func(cand, di []int8, dis [][]int8, indices [][2]int8, nmr uint64, nd, level int)
   fnpr = func(cand, di []int8, dis [][]int8, indices [][2]int8, nmr uint64, nd, level int) {
       if level == len(dis) {
           digits[indices[0][0]] = fml[cand[0]][di[0]][0]
           digits[indices[0][1]] = fml[cand[0]][di[0]][1]
           le := len(di)
           if nd%2 == 1 {
               le--
               digits[nd/2] = di[le]
           }
           for i, d := range di[1:le] {
               digits[indices[i+1][0]] = dmd[cand[i+1]][d][0]
               digits[indices[i+1][1]] = dmd[cand[i+1]][d][1]
           }
           r := toUint64(digits, true)
           npr := nmr + 2*r
           if !isSquare(npr) {
               return
           }
           count++
           fmt.Printf("     R/N %2d:", count)
           ms := uint64(time.Since(start).Milliseconds())
           fmt.Printf("  %9s ms", commatize(ms))
           n := toUint64(digits, false)
           fmt.Printf("  (%s)\n", commatize(n))
           rares = append(rares, n)
       } else {
           for _, num := range dis[level] {
               di[level] = num
               fnpr(cand, di, dis, indices, nmr, nd, level+1)
           }
       }
   }

   // Recursive closure to generate (n-r) candidates with a given number of digits.
   var fnmr func(cand []int8, list [][]int8, indices [][2]int8, nd, level int)
   fnmr = func(cand []int8, list [][]int8, indices [][2]int8, nd, level int) {
       if level == len(list) {
           var nmr, nmr2 uint64
           for i, t := range allTerms[nd-2] {
               if cand[i] >= 0 {
                   nmr += t.coeff * uint64(cand[i])
               } else {
                   nmr2 += t.coeff * uint64(-cand[i])
                   if nmr >= nmr2 {
                       nmr -= nmr2
                       nmr2 = 0
                   } else {
                       nmr2 -= nmr
                       nmr = 0
                   }
               }
           }
           if nmr2 >= nmr {
               return
           }
           nmr -= nmr2
           if !isSquare(nmr) {
               return
           }
           var dis [][]int8
           dis = append(dis, seq(0, int8(len(fml[cand[0]]))-1, 1))
           for i := 1; i < len(cand); i++ {
               dis = append(dis, seq(0, int8(len(dmd[cand[i]]))-1, 1))
           }
           if nd%2 == 1 {
               dis = append(dis, il)
           }
           di := make([]int8, len(dis))
           fnpr(cand, di, dis, indices, nmr, nd, 0)
       } else {
           for _, num := range list[level] {
               cand[level] = num
               fnmr(cand, list, indices, nd, level+1)
           }
       }
   }

   for nd := 2; nd <= maxDigits; nd++ {
       digits = make([]int8, nd)
       if nd == 4 {
           lists[0] = append(lists[0], zl)
           lists[1] = append(lists[1], ol)
           lists[2] = append(lists[2], el)
           lists[3] = append(lists[3], ol)
       } else if len(allTerms[nd-2]) > len(lists[0]) {
           for i := 0; i < 4; i++ {
               lists[i] = append(lists[i], dl)
           }
       }
       var indices [][2]int8
       for _, t := range allTerms[nd-2] {
           indices = append(indices, [2]int8{t.ix1, t.ix2})
       }
       for _, list := range lists {
           cand := make([]int8, len(list))
           fnmr(cand, list, indices, nd, 0)
       }
       ms := uint64(time.Since(start).Milliseconds())
       fmt.Printf("  %2d digits:  %9s ms\n", nd, commatize(ms))
   }

   sort.Slice(rares, func(i, j int) bool { return rares[i] < rares[j] })
   fmt.Printf("\nThe rare numbers with up to %d digits are:\n", maxDigits)
   for i, rare := range rares {
       fmt.Printf("  %2d:  %25s\n", i+1, commatize(rare))
   }

}</lang>

Output:

Timings are for an Intel Core i7-8565U machine with 32GB RAM running Go 1.13.3 on Ubuntu 18.04.

Aggregate timings to process all numbers up to:
     R/N  1:          0 ms  (65)
   2 digits:          0 ms
   3 digits:          0 ms
   4 digits:          0 ms
   5 digits:          0 ms
     R/N  2:          0 ms  (621,770)
   6 digits:          0 ms
   7 digits:          0 ms
   8 digits:          3 ms
     R/N  3:          3 ms  (281,089,082)
   9 digits:          4 ms
     R/N  4:          5 ms  (2,022,652,202)
     R/N  5:         31 ms  (2,042,832,002)
  10 digits:         71 ms
  11 digits:        109 ms
     R/N  6:        328 ms  (872,546,974,178)
     R/N  7:        355 ms  (872,568,754,178)
     R/N  8:        697 ms  (868,591,084,757)
  12 digits:        848 ms
     R/N  9:      1,094 ms  (6,979,302,951,885)
  13 digits:      1,406 ms
     R/N 10:      5,121 ms  (20,313,693,904,202)
     R/N 11:      5,187 ms  (20,313,839,704,202)
     R/N 12:      6,673 ms  (20,331,657,922,202)
     R/N 13:      6,887 ms  (20,331,875,722,202)
     R/N 14:      7,479 ms  (20,333,875,702,202)
     R/N 15:     17,112 ms  (40,313,893,704,200)
     R/N 16:     17,248 ms  (40,351,893,720,200)
  14 digits:     18,338 ms
     R/N 17:     18,356 ms  (200,142,385,731,002)
     R/N 18:     18,560 ms  (221,462,345,754,122)
     R/N 19:     21,181 ms  (816,984,566,129,618)
     R/N 20:     22,491 ms  (245,518,996,076,442)
     R/N 21:     22,674 ms  (204,238,494,066,002)
     R/N 22:     22,734 ms  (248,359,494,187,442)
     R/N 23:     22,994 ms  (244,062,891,224,042)
     R/N 24:     26,868 ms  (403,058,392,434,500)
     R/N 25:     27,063 ms  (441,054,594,034,340)
  15 digits:     28,087 ms
     R/N 26:     74,120 ms  (2,133,786,945,766,212)
     R/N 27:     92,245 ms  (2,135,568,943,984,212)
     R/N 28:     94,972 ms  (8,191,154,686,620,818)
     R/N 29:     97,313 ms  (8,191,156,864,620,818)
     R/N 30:     98,361 ms  (2,135,764,587,964,212)
     R/N 31:     99,971 ms  (2,135,786,765,764,212)
     R/N 32:    103,603 ms  (8,191,376,864,400,818)
     R/N 33:    115,711 ms  (2,078,311,262,161,202)
     R/N 34:    140,972 ms  (8,052,956,026,592,517)
     R/N 35:    145,099 ms  (8,052,956,206,592,517)
     R/N 36:    175,023 ms  (8,650,327,689,541,457)
     R/N 37:    177,145 ms  (8,650,349,867,341,457)
     R/N 38:    178,693 ms  (6,157,577,986,646,405)
     R/N 39:    205,564 ms  (4,135,786,945,764,210)
     R/N 40:    220,480 ms  (6,889,765,708,183,410)
  16 digits:    221,485 ms
     R/N 41:    226,520 ms  (86,965,750,494,756,968)
     R/N 42:    227,431 ms  (22,542,040,692,914,522)
     R/N 43:    345,314 ms  (67,725,910,561,765,640)
  17 digits:    354,815 ms
     R/N 44:    387,879 ms  (284,684,666,566,486,482)
     R/N 45:    510,788 ms  (225,342,456,863,243,522)
     R/N 46:    556,239 ms  (225,342,458,663,243,522)
     R/N 47:    652,051 ms  (225,342,478,643,243,522)
     R/N 48:    718,148 ms  (284,684,868,364,486,482)
     R/N 49:  1,095,093 ms  (871,975,098,681,469,178)
     R/N 50:  1,785,243 ms  (865,721,270,017,296,468)
     R/N 51:  1,800,799 ms  (297,128,548,234,950,692)
     R/N 52:  1,809,196 ms  (297,128,722,852,950,692)
     R/N 53:  1,913,085 ms  (811,865,096,390,477,018)
     R/N 54:  1,965,104 ms  (297,148,324,656,930,692)
     R/N 55:  1,990,018 ms  (297,148,546,434,930,692)
     R/N 56:  2,296,972 ms  (898,907,259,301,737,498)
     R/N 57:  2,691,110 ms  (631,688,638,047,992,345)
     R/N 58:  2,716,487 ms  (619,431,353,040,136,925)
     R/N 59:  2,984,451 ms  (619,631,153,042,134,925)
     R/N 60:  3,047,183 ms  (633,288,858,025,996,145)
     R/N 61:  3,115,724 ms  (633,488,632,647,994,145)
     R/N 62:  3,978,143 ms  (653,488,856,225,994,125)
     R/N 63:  4,255,985 ms  (497,168,548,234,910,690)
  18 digits:  4,531,846 ms
     R/N 64:  4,606,094 ms  (2,551,755,006,254,571,552)
     R/N 65:  4,624,539 ms  (2,702,373,360,882,732,072)
     R/N 66:  4,873,160 ms  (2,825,378,427,312,735,282)
     R/N 67:  4,893,810 ms  (8,066,308,349,502,036,608)
     R/N 68:  5,109,513 ms  (2,042,401,829,204,402,402)
     R/N 69:  5,152,863 ms  (2,420,424,089,100,600,242)
     R/N 70:  5,263,434 ms  (8,320,411,466,598,809,138)
     R/N 71:  5,558,356 ms  (8,197,906,905,009,010,818)
     R/N 72:  5,586,801 ms  (2,060,303,819,041,450,202)
     R/N 73:  5,763,382 ms  (8,200,756,128,308,135,597)
     R/N 74:  6,008,475 ms  (6,531,727,101,458,000,045)
     R/N 75:  6,543,047 ms  (6,988,066,446,726,832,640)
  19 digits:  6,609,905 ms

The rare numbers with up to 19 digits are:
   1:                         65
   2:                    621,770
   3:                281,089,082
   4:              2,022,652,202
   5:              2,042,832,002
   6:            868,591,084,757
   7:            872,546,974,178
   8:            872,568,754,178
   9:          6,979,302,951,885
  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
  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
  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
  41:     22,542,040,692,914,522
  42:     67,725,910,561,765,640
  43:     86,965,750,494,756,968
  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
  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

Julia

Translation of: Go

<lang julia>using Formatting, Printf

struct Term

   coeff::UInt64
   ix1::Int8
   ix2::Int8

end

function toUInt64(dgits, reverse)

   return reverse ? foldr((i, j) -> i + 10j, UInt64.(dgits)) :
                    foldl((i, j) -> 10i + j, UInt64.(dgits))

end

function issquare(n)

   if 0x202021202030213 & (1 << (UInt64(n) & 63)) != 0
       root = UInt64(floor(sqrt(n)))
       return root * root == n
   end
   return false

end

seq(from, to, step) = Int8.(collect(from:step:to))

commatize(n::Integer) = format(n, commas=true)

const verbose = true const count = [0]

""" Recursive closure to generate (n+r) candidates from (n-r) candidates and hence find Rare numbers with a given number of digits. """ function fnpr(cand, di, dis, indices, nmr, nd, level, dgits, fml, dmd, start, rares, il)

   if level == length(dis)
       dgits[indices[1][1] + 1] = fml[cand[1]][di[1] + 1][1]
       dgits[indices[1][2] + 1] = fml[cand[1]][di[1] + 1][2]
       le = length(di)
       if nd % 2 == 1
           le -= 1
           dgits[nd ÷ 2 + 1] = di[le + 1]
       end
       for (i, d) in enumerate(di[2:le])
           dgits[indices[i+1][1] + 1] = dmd[cand[i+1]][d + 1][1]
           dgits[indices[i+1][2] + 1] = dmd[cand[i+1]][d + 1][2]
       end
       r = toUInt64(dgits, true)
       npr = nmr + 2 * r
       !issquare(npr) && return
       count[1] += 1
       verbose && @printf("     R/N %2d:", count[1])
       !verbose && print("$count rares\b\b\b\b\b\b\b\b\b")
       ms = UInt64(time() * 1000 - start)
       verbose && @printf("  %9s ms", commatize(Int(ms)))
       n = toUInt64(dgits, false)
       verbose && @printf("  (%s)\n", commatize(BigInt(n)))
       push!(rares, n)
   else
       for num in dis[level + 1]
           di[level + 1] = num
           fnpr(cand, di, dis, indices, nmr, nd, level + 1, dgits, fml, dmd, start, rares, il)
       end
   end

end # function fnpr

Recursive closure to generate (n-r) candidates with a given number of digits.
var fnmr func(cand []int8, list [][]int8, indices [][2]int8, nd, level int)

function fnmr(cand, list, indices, nd, level, allterms, fml, dmd, dgits, start, rares, il)

   if level == length(list)
       nmr, nmr2 = zero(UInt64), zero(UInt64)
       for (i, t) in enumerate(allterms[nd - 1])
           if cand[i] >= 0
               nmr += t.coeff * UInt64(cand[i])
           else
               nmr2 += t.coeff * UInt64(-cand[i])
               if nmr >= nmr2
                   nmr -= nmr2
                   nmr2 = zero(nmr2)
               else
                   nmr2 -= nmr
                   nmr = zero(nmr)
               end
           end
       end
       nmr2 >= nmr && return
       nmr -= nmr2
       !issquare(nmr) && return
       dis = [[seq(0, Int8(length(fml[cand[1]]) - 1), 1)] ;
           [seq(0, Int8(length(dmd[c]) - 1), 1) for c in cand[2:end]]]
       isodd(nd) && push!(dis, il)
       di = zeros(Int8, length(dis))
       fnpr(cand, di, dis, indices, nmr, nd, 0, dgits, fml, dmd, start, rares, il)
   else
       for num in list[level + 1]
           cand[level + 1] = num
           fnmr(cand, list, indices, nd, level + 1, allterms, fml, dmd, dgits, start, rares, il)
       end
   end

end # function fnmr

function findrare(maxdigits = 19)

   start = time() * 1000.0
   pow = one(UInt64)
   verbose && println("Aggregate timings to process all numbers up to:")
   # terms of (n-r) expression for number of digits from 2 to maxdigits
   allterms = Vector{Vector{Term}}()
   for r in 2:maxdigits
       terms = Term[]
       pow *= 10
       pow1, pow2, i1, i2 = pow, one(UInt64), zero(Int8), Int8(r - 1)
       while i1 < i2
           push!(terms, Term(pow1 - pow2, i1, i2))
           pow1, pow2, i1, i2 = pow1 ÷ 10, pow2 * 10, i1 + 1, i2 - 1
       end
       push!(allterms, terms)
   end
   #  map of first minus last digits for 'n' to pairs giving this value
   fml = Dict(
       0 => [2 => 2, 8 => 8],
       1 => [6 => 5, 8 => 7],
       4 => [4 => 0],
       6 => [6 => 0, 8 => 2],
   )
   # map of other digit differences for 'n' to pairs giving this value
   dmd = Dict{Int8, Vector{Vector{Int8}}}()
   for i in 0:99
       a = [Int8(i ÷ 10), Int8(i % 10)]
       d = a[1] - a[2]
       v = get!(dmd, d, [])
       push!(v, a)
   end
   fl = Int8[0, 1, 4, 6]
   dl = seq(-9, 9, 1)  # all differences
   zl = Int8[0]        # zero differences only
   el = seq(-8, 8, 2)  # even differences only
   ol = seq(-9, 9, 2)  # odd differences only
   il = seq(0, 9, 1)
   rares = UInt64[]
   lists = [[[f]] for f in fl]
   dgits = Int8[]
   count[1] = 0

   for nd = 2:maxdigits
       dgits = zeros(Int8, nd)
       if nd == 4
           push!(lists[1], zl)
           push!(lists[2], ol)
           push!(lists[3], el)
           push!(lists[4], ol)
       elseif length(allterms[nd - 1]) > length(lists[1])
           for i in 1:4
               push!(lists[i], dl)
           end
       end
       indices = Vector{Vector{Int8}}()
       for t in allterms[nd - 1]
           push!(indices, Int8[t.ix1, t.ix2])
       end
       for list in lists
           cand = zeros(Int8, length(list))
           fnmr(cand, list, indices, nd, 0, allterms, fml, dmd, dgits, start, rares, il)
       end
       ms = UInt64(time() * 1000 - start)
       verbose && @printf("  %2d digits:  %9s ms\n", nd, commatize(Int(ms)))
   end

   sort!(rares)
   @printf("\nThe rare numbers with up to %d digits are:\n", maxdigits)
   for (i, rare) in enumerate(rares)
       @printf("  %2d:  %25s\n", i, commatize(BigInt(rare)))
   end

end # findrare function

findrare()

</lang>

Output:

Timings are on a 2.9 GHz i5 processor, 16 Gb RAM, under Windows 10.

Aggregate timings to process all numbers up to:
     R/N  1:          5 ms  (65)
   2 digits:        132 ms
   3 digits:        133 ms
   4 digits:        134 ms
   5 digits:        134 ms
     R/N  2:        135 ms  (621,770)
   6 digits:        135 ms
   7 digits:        136 ms
   8 digits:        140 ms
     R/N  3:        141 ms  (281,089,082)
   9 digits:        143 ms
     R/N  4:        144 ms  (2,022,652,202)
     R/N  5:        168 ms  (2,042,832,002)
  10 digits:        209 ms
  11 digits:        251 ms
     R/N  6:        443 ms  (872,546,974,178)
     R/N  7:        467 ms  (872,568,754,178)
     R/N  8:        773 ms  (868,591,084,757)
  12 digits:        919 ms
     R/N  9:      1,178 ms  (6,979,302,951,885)
  13 digits:      1,510 ms
     R/N 10:      4,662 ms  (20,313,693,904,202)
     R/N 11:      4,722 ms  (20,313,839,704,202)
     R/N 12:      6,028 ms  (20,331,657,922,202)
     R/N 13:      6,223 ms  (20,331,875,722,202)
     R/N 14:      6,753 ms  (20,333,875,702,202)
     R/N 15:     15,475 ms  (40,313,893,704,200)
     R/N 16:     15,594 ms  (40,351,893,720,200)
  14 digits:     16,749 ms
     R/N 17:     16,772 ms  (200,142,385,731,002)
     R/N 18:     17,006 ms  (221,462,345,754,122)
     R/N 19:     20,027 ms  (816,984,566,129,618)
     R/N 20:     21,669 ms  (245,518,996,076,442)
     R/N 21:     21,895 ms  (204,238,494,066,002)
     R/N 22:     21,974 ms  (248,359,494,187,442)
     R/N 23:     22,302 ms  (244,062,891,224,042)
     R/N 24:     27,158 ms  (403,058,392,434,500)
     R/N 25:     27,405 ms  (441,054,594,034,340)
  15 digits:     28,744 ms
     R/N 26:     79,350 ms  (2,133,786,945,766,212)
     R/N 27:     99,360 ms  (2,135,568,943,984,212)
     R/N 28:    102,426 ms  (8,191,154,686,620,818)
     R/N 29:    105,135 ms  (8,191,156,864,620,818)
     R/N 30:    106,334 ms  (2,135,764,587,964,212)
     R/N 31:    108,038 ms  (2,135,786,765,764,212)
     R/N 32:    112,142 ms  (8,191,376,864,400,818)
     R/N 33:    125,607 ms  (2,078,311,262,161,202)
     R/N 34:    154,417 ms  (8,052,956,026,592,517)
     R/N 35:    159,075 ms  (8,052,956,206,592,517)
     R/N 36:    192,323 ms  (8,650,327,689,541,457)
     R/N 37:    194,651 ms  (8,650,349,867,341,457)
     R/N 38:    196,344 ms  (6,157,577,986,646,405)
     R/N 39:    227,492 ms  (4,135,786,945,764,210)
     R/N 40:    244,502 ms  (6,889,765,708,183,410)
  16 digits:    245,658 ms
     R/N 41:    251,178 ms  (86,965,750,494,756,968)
     R/N 42:    252,157 ms  (22,542,040,692,914,522)
     R/N 43:    382,883 ms  (67,725,910,561,765,640)
  17 digits:    393,371 ms
     R/N 44:    427,555 ms  (284,684,666,566,486,482)
     R/N 45:    549,740 ms  (225,342,456,863,243,522)
     R/N 46:    594,392 ms  (225,342,458,663,243,522)
     R/N 47:    688,221 ms  (225,342,478,643,243,522)
     R/N 48:    753,385 ms  (284,684,868,364,486,482)
     R/N 49:  1,108,538 ms  (871,975,098,681,469,178)
     R/N 50:  1,770,255 ms  (865,721,270,017,296,468)
     R/N 51:  1,785,243 ms  (297,128,548,234,950,692)
     R/N 52:  1,793,571 ms  (297,128,722,852,950,692)
     R/N 53:  1,892,872 ms  (811,865,096,390,477,018)
     R/N 54:  1,941,208 ms  (297,148,324,656,930,692)
     R/N 55:  1,964,502 ms  (297,148,546,434,930,692)
     R/N 56:  2,267,616 ms  (898,907,259,301,737,498)
     R/N 57:  2,677,207 ms  (631,688,638,047,992,345)
     R/N 58:  2,702,836 ms  (619,431,353,040,136,925)
     R/N 59:  2,960,274 ms  (619,631,153,042,134,925)
     R/N 60:  3,019,846 ms  (633,288,858,025,996,145)
     R/N 61:  3,084,695 ms  (633,488,632,647,994,145)
     R/N 62:  3,924,801 ms  (653,488,856,225,994,125)
     R/N 63:  4,229,162 ms  (497,168,548,234,910,690)
  18 digits:  4,563,375 ms
     R/N 64:  4,643,118 ms  (2,551,755,006,254,571,552)
     R/N 65:  4,662,645 ms  (2,702,373,360,882,732,072)
     R/N 66:  4,925,324 ms  (2,825,378,427,312,735,282)
     R/N 67:  4,947,368 ms  (8,066,308,349,502,036,608)
     R/N 68:  5,170,716 ms  (2,042,401,829,204,402,402)
     R/N 69:  5,216,832 ms  (2,420,424,089,100,600,242)
     R/N 70:  5,329,680 ms  (8,320,411,466,598,809,138)
     R/N 71:  5,634,991 ms  (8,197,906,905,009,010,818)
     R/N 72:  5,665,799 ms  (2,060,303,819,041,450,202)
     R/N 73:  5,861,019 ms  (8,200,756,128,308,135,597)
     R/N 74:  6,136,091 ms  (6,531,727,101,458,000,045)
     R/N 75:  6,770,242 ms  (6,988,066,446,726,832,640)
  19 digits:  6,846,705 ms

The rare numbers with up to 19 digits are:
   1:                         65
   2:                    621,770
   3:                281,089,082
   4:              2,022,652,202
   5:              2,042,832,002
   6:            868,591,084,757
   7:            872,546,974,178
   8:            872,568,754,178
   9:          6,979,302,951,885
  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
  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
  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
  41:     22,542,040,692,914,522
  42:     67,725,910,561,765,640
  43:     86,965,750,494,756,968
  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
  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

Phix

naieve

Ridiculously slow, 90s just for the first 3. <lang Phix>function revn(atom n, integer nd)

   atom r = 0
   for i=1 to nd do
       r = r*10+remainder(n,10)
       n = floor(n/10)
   end for
   return r

end function

integer nd = 2, count = 0 atom lim = 99, n = 9, t0 = time() while true do

   n += 1
   atom r = revn(n,nd)
   if r<n then
       atom s = n+r,
            d = n-r
       if s=power(floor(sqrt(s)),2)
       and d=power(floor(sqrt(d)),2) then
           count += 1
           ?{count,n,elapsed(time()-t0)}
           if count=3 then exit end if
       end if
   end if
   if n=lim then

-- ?{"lim",lim,elapsed(time()-t0)}

       lim = lim*10+9
       nd += 1
   end if

end while</lang>

Output:

{1,65,"0s"}
{2,621770,"0.2s"}
{3,281089082,"1 minute and 29s"}

advanced

Translation of: VB.NET

Translation of: Go

Quite slow: over 10 mins for first 25, vs 53 secs of Go... easily the worst such (ie Phix vs. Go) comparison I have yet seen.
This task exposes some of the weaknesses of Phix, specifically subscripting speed (suprisingly not usually that much of an issue), and the fact that functions are never inlined. The relatively innoculous-looking and dirt simple to_atom() routine is responsible for over 30% of the run time. Improvements welcome: run p -d test, examine the list.asm produced from this source, and discuss or make suggestions on my talk page. <lang Phix>constant maxDigits = 15

enum COEFF, TDXA, TDXB -- struct term = {atom coeff, integer idxa, idxb}

                       -- (see allTerms below)

integer nd, -- number of digits

       count   -- of solutions found earlier, for lower nd

sequence rares -- (cleared after sorting/printing for each nd)

function to_atom(sequence digits) -- convert digits array to an atom value

   atom r = 0
   for i=1 to length(digits) do
       r = r * 10 + digits[i]
   end for
   return r

end function

-- psq eliminates 52 out of 64 of numbers fairly cheaply, which translates -- to approximately 66% of numbers, or around 10% off the overall time. -- NB: only tested to 9,007,199,254,740,991, then again I found no more new -- bit patterns after just 15^2.

constant psq = int_to_bits(#02030213,32)& -- #0202021202030213 --> bits,

              int_to_bits(#02020212,32)    -- in 32/64-bit compatible way.

function isSquare(atom n) -- determine if n is a perfect square or not

   if psq[and_bits(n,63)+1] then
       atom r = floor(sqrt(n))
       return r * r = n
   end if
   return false

end function

procedure fnpr(integer level, atom nmr, sequence di, dis, candidates, indices, fml, dmd) -- generate (n+r) candidates from (n-r) candidates

   if level>length(dis) then
       sequence digits = repeat(0,nd)
       -- (the precise why of how this populates digits has eluded me...)
       integer {a,b} = indices[1],
               c = candidates[1]+1,
               d = di[1]+1
       digits[a] = fml[c][d][1]
       digits[b] = fml[c][d][2]
       integer le = length(di)
       if remainder(nd,2) then
           d = floor(nd/2)+1
           digits[d] = di[le]
           le -= 1
       end if
       for dx=2 to le do
           {a,b} = indices[dx]
           c = candidates[dx]+10
           d = di[dx]+1
           digits[a] = dmd[c][d][1]
           digits[b] = dmd[c][d][2]
       end for
       atom npr = nmr + to_atom(reverse(digits))*2 -- (npr == 'n + r')
       if isSquare(npr) then
           rares &= to_atom(digits)
           -- (note this gets overwritten by sorted set:)
           printf(1,"working...  %2d: %,d\r", {count+length(rares),rares[$]})
       end if
   else
       for n=0 to dis[level] do
           di[level] = n
           fnpr(level+1, nmr, di, dis, candidates, indices, fml, dmd)
       end for
   end if

end procedure

procedure fnmr(sequence terms, list, candidates, indices, fml, dmd, integer level) -- generate (n-r) candidates with a given number of digits.

   if level>length(list) then
       atom nmr = 0 -- (nmr == 'n - r')
       for i=1 to length(terms) do
           nmr += terms[i][COEFF] * candidates[i]
       end for
       if nmr>0 and isSquare(nmr) then
           integer c = candidates[1]+1,
                   l = length(fml[c])-1
           sequence dis = {l}
           for i=2 to length(candidates) do
               c = candidates[i]+10
               l = length(dmd[c])-1
               dis = append(dis,l)
           end for
           if remainder(nd,2) then dis = append(dis,9) end if
           -- (above generates dis of eg {1,4,7,9} for nd=7, which as far
           --  as I (lightly) understand it scans for far fewer candidate 
           --  pairs than a {9,9,9,9} would, or something like that.)
           sequence di = repeat(0,length(dis))
           -- (di is the current "dis-scan", eg {0,0,0,0} to {1,4,7,9})
           fnpr(1, nmr, di, dis, candidates, indices, fml, dmd)
       end if
   else
       for n=1 to length(list[level]) do
           candidates[level] = list[level][n]
           fnmr(terms, list, candidates, indices, fml, dmd, level+1)
       end for
   end if

end procedure

constant dl = tagset(9,-9), -- all differences (-9..+9 by 1)

        zl = tagset(0, 0),    -- zero difference  (0 only)
        el = tagset(8,-8, 2), -- even differences (-8 to +8 by 2)
        ol = tagset(9,-9, 2), -- odd  differences (-9..+9 by 2)
        il = tagset(9, 0)     -- all integers (0..9 by 1)

procedure main()

   atom start = time()

   -- terms of (n-r) expression for number of digits from 2 to maxdigits
   sequence allTerms = {}
   atom pow = 1
   for r=2 to maxDigits do
       sequence terms = {}
       pow *= 10
       atom p1 = pow, p2 = 1
       integer tdxa = 0, tdxb = r-1
       while tdxa < tdxb do
           terms = append(terms,{p1-p2, tdxa, tdxb}) -- {COEFF,TDXA,TDXB}
           p1 /= 10
           p2 *= 10
           tdxa += 1
           tdxb -= 1
       end while
       allTerms = append(allTerms,terms)
   end for

--/* --(This is what the above loop creates:) --pp(allTerms,{pp_Nest,1,pp_StrFmt,3,pp_IntCh,false,pp_IntFmt,"%d",pp_FltFmt,"%d",pp_Maxlen,148}) {Template:9,0,1,

Template:99,0,2,
{{999,0,3}, {90,1,2}},
{{9999,0,4}, {990,1,3}},
{{99999,0,5}, {9990,1,4}, {900,2,3}},
{{999999,0,6}, {99990,1,5}, {9900,2,4}},
{{9999999,0,7}, {999990,1,6}, {99900,2,5}, {9000,3,4}},
{{99999999,0,8}, {9999990,1,7}, {999900,2,6}, {99000,3,5}},
{{999999999,0,9}, {99999990,1,8}, {9999900,2,7}, {999000,3,6}, {90000,4,5}},
{{9999999999,0,10}, {999999990,1,9}, {99999900,2,8}, {9999000,3,7}, {990000,4,6}},
{{99999999999,0,11}, {9999999990,1,10}, {999999900,2,9}, {99999000,3,8}, {9990000,4,7}, {900000,5,6}},
{{999999999999,0,12}, {99999999990,1,11}, {9999999900,2,10}, {999999000,3,9}, {99990000,4,8}, {9900000,5,7}},
{{9999999999999,0,13}, {999999999990,1,12}, {99999999900,2,11}, {9999999000,3,10}, {999990000,4,9}, {99900000,5,8}, {9000000,6,7}},
{{99999999999999,0,14}, {9999999999990,1,13}, {999999999900,2,12}, {99999999000,3,11}, {9999990000,4,10}, {999900000,5,9}, {99000000,6,8}}}

--*/

   -- map of first minus last digits for 'n' to pairs giving this value
   sequence fml = repeat({},10) -- (aka 0..9)
   --      (fml == 'first minus last')
   fml[1] = {{2, 2}, {8, 8}}
   fml[2] = {{6, 5}, {8, 7}}
   fml[5] = Template:4, 0

-- fml[6] = Template:8, 3 -- (um? - needs longer lists, & that append(lists[4],dl) below)

   fml[7] = {{6, 0}, {8, 2}}

-- sequence lists = {Template:0,Template:1,Template:4,Template:5,Template:6}

   sequence lists = {Template:0,Template:1,Template:4,Template:6}
   -- map of other digit differences for 'n' to pairs giving this value
   sequence dmd = repeat({},19) -- (aka -9..+9, so add 10 when indexing dmd)
   --      (dmd == 'digit minus digit')
   for tens=0 to 9 do
       integer d = tens+10
       for ones=0 to 9 do
           dmd[d] = append(dmd[d], {tens,ones})
           d -= 1
       end for
   end for

--/* --(This is what the above loop creates:) --pp(dmd,{pp_Nest,1,pp_StrFmt,3,pp_IntCh,false}) {Template:0,9,

{{0,8}, {1,9}},
{{0,7}, {1,8}, {2,9}},
{{0,6}, {1,7}, {2,8}, {3,9}},
{{0,5}, {1,6}, {2,7}, {3,8}, {4,9}},
{{0,4}, {1,5}, {2,6}, {3,7}, {4,8}, {5,9}},
{{0,3}, {1,4}, {2,5}, {3,6}, {4,7}, {5,8}, {6,9}},
{{0,2}, {1,3}, {2,4}, {3,5}, {4,6}, {5,7}, {6,8}, {7,9}},
{{0,1}, {1,2}, {2,3}, {3,4}, {4,5}, {5,6}, {6,7}, {7,8}, {8,9}},
{{0,0}, {1,1}, {2,2}, {3,3}, {4,4}, {5,5}, {6,6}, {7,7}, {8,8}, {9,9}},
{{1,0}, {2,1}, {3,2}, {4,3}, {5,4}, {6,5}, {7,6}, {8,7}, {9,8}},
{{2,0}, {3,1}, {4,2}, {5,3}, {6,4}, {7,5}, {8,6}, {9,7}},
{{3,0}, {4,1}, {5,2}, {6,3}, {7,4}, {8,5}, {9,6}},
{{4,0}, {5,1}, {6,2}, {7,3}, {8,4}, {9,5}},
{{5,0}, {6,1}, {7,2}, {8,3}, {9,4}},
{{6,0}, {7,1}, {8,2}, {9,3}},
{{7,0}, {8,1}, {9,2}},
{{8,0}, {9,1}},
Template:9,0}

--*/

   count = 0
   printf(1,"digits time  nth rare numbers:\n")
   nd = 2
   while nd <= maxDigits do
       rares = {}
       sequence terms = allTerms[nd-1]
       if nd=4 then
           lists[1] = append(lists[1],zl)
           lists[2] = append(lists[2],ol)
           lists[3] = append(lists[3],el)

-- lists[4] = append(lists[4],dl) -- if fml[6] = Template:8, 3 -- lists[5] = append(lists[5],ol) -- ""

           lists[4] = append(lists[4],ol)  -- else
       elsif length(terms)>length(lists[1]) then
           for i=1 to length(lists) do
               lists[i] = append(lists[i],dl)
           end for
       end if
       sequence indices = {}
       for t=1 to length(terms) do
           sequence term = terms[t]
           -- (we may as well make this 1-based while here)
           indices = append(indices,{term[TDXA]+1,term[TDXB]+1})
       end for
       for i=1 to length(lists) do
           sequence list = lists[i],
                    candidates = repeat(0,length(list))
           fnmr(terms, list, candidates, indices, fml, dmd, 1)
       end for
       -- (re-)output partial results for this nd-set in sorted order:
       rares = sort(rares)
       for i=1 to length(rares) do
           count += 1
           printf(1,"%12s %2d: %,19d \n", {"",count,rares[i]})
       end for
       printf(1,"  %2d  %5s\n", {nd, elapsed_short(time()-start)})
       nd += 1
   end while

end procedure main()|</lang>

Output:

digits time  nth rare numbers:
              1:                  65
   2     0s
   3     0s
   4     0s
   5     0s
              2:             621,770
   6     0s
   7     0s
   8     0s
              3:         281,089,082
   9     0s
              4:       2,022,652,202
              5:       2,042,832,002
  10     0s
  11     1s
              6:     868,591,084,757
              7:     872,546,974,178
              8:     872,568,754,178
  12    15s
              9:   6,979,302,951,885
  13    29s
             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   6:07
             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  10:42

REXX

(See the discussion page for a simplistic 1^st version that computes rare numbers only using the task's basic rules).

Most of the hints (properties of rare numbers) within Shyam Sunder Gupta's webpage have been incorporated in this
REXX program and the logic is now expressed within the list of AB...PQ (abutted numbers within the @g list).

These improvements made this REXX version around 25% faster than the previous version (see the discussion page). <lang rexx>/*REXX program calculates and displays a specified amount of rare numbers. */ numeric digits 20; w= digits() + digits() % 3 /*use enough dec. digs for calculations*/ parse arg many . /*obtain optional argument from the CL.*/ if many== | many=="," then many= 5 /*Not specified? Then use the default.*/ @g= 2002 2112 2222 2332 2442 2552 2662 2772 2882 2992 4000 4010 4030 4050 4070 4090 4100 ,

   4110 4120 4140 4160 4180 4210 4230 4250 4270 4290 4300 4320 4340 4360 4380 4410 4430 ,
   4440 4450 4470 4490 4500 4520 4540 4560 4580 4610 4630 4650 4670 4690 4700 4720 4740 ,
   4760 4780 4810 4830 4850 4870 4890 4900 4920 4940 4960 4980 4990 6010 6015 6030 6035 ,
   6050 6055 6070 6075 6090 6095 6100 6105 6120 6125 6140 6145 6160 6165 6180 6185 6210 ,
   6215 6230 6235 6250 6255 6270 6275 6290 6295 6300 6305 6320 6325 6340 6345 6360 6365 ,
   6380 6385 6410 6415 6430 6435 6450 6455 6470 6475 6490 6495 6500 6505 6520 6525 6540 ,
   6545 6560 6565 6580 6585 6610 6615 6630 6635 6650 6655 6670 6675 6690 6695 6700 6705 ,
   6720 6725 6740 6745 6760 6765 6780 6785 6810 6815 6830 6835 6850 6855 6870 6875 6890 ,
   6895 6900 6905 6920 6925 6940 6945 6960 6965 6980 6985 8007 8008 8017 8027 8037 8047 ,
   8057 8067 8077 8087 8092 8097 8107 8117 8118 8127 8137 8147 8157 8167 8177 8182 8187 ,
   8197 8228 8272 8297 8338 8362 8387 8448 8452 8477 8542 8558 8567 8632 8657 8668 8722 ,
   8747 8778 8812 8837 8888 8902 8927 8998      /*4 digit abutted numbers for AB and PQ*/

@g#= words(@g)

        /* [↓]─────────────────boolean arrays are used for checking for digit presence.*/

@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.*/ @149.=0; @149.1=1; @149.4=1; @149.9=1 /*values for Z that need an even Y. */ @odd.=0; do i=-9 by 2 to 9; @odd.i=1 /* " " N " " " " A. */

         end   /*i*/

@gen.=0; do i=1 for words(@g); parse value word(@g,i) with a 2 b 3 p 4 q; @gen.a.b.p.q=1

              /*# AB···PQ  could be a good rare value*/
         end   /*i*/

div9= 9 /*dif must be ÷ 9 when N has even #digs*/ evenN= \ (10 // 2) /*initial value for evenness of N. */

= 0 /*the number of rare numbers (so far)*/

   do n=10                                      /*Why 10?  All 1 dig #s are palindromic*/
   parse var   n   a  2  b  3    -2  p  +1  q /*get 1st\2nd\penultimate\last digits. */
   if @odd.a  then do;  n=n+10**(length(n)-1)-1 /*bump N so next N starts with even dig*/
                        evenN=\(length(n+1)//2) /*flag when N has an even # of digits. */
                        if evenN  then div9=  9 /*when dif isn't divisible by   9  ... */
                                  else div9= 99 /*  "   "    "        "     "  99   "  */
                        iterate                 /*let REXX do its thing with  DO  loop.*/
                   end                          /* {it's allowed to modify a DO index} */
   if \@gen.a.b.p.q  then iterate               /*can  N  not be a rare AB···PQ number?*/
   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.*/
   dif= n-r;   parse var  dif    -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*/                          /* [↑]  the OTHERWISE handles Z=8 case.*/
   if dif//div9\==0  then iterate               /*Difference isn't ÷ by div9? Then skip*/
   sum= n+r;   parse var  sum    -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*/                          /* [↑]  the OTHERWISE handles Z=8 case.*/
   if evenN  then if sum//11 \==0  then iterate /*N has even #digs? Sum must be ÷ by 11*/
   $= 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*/
      do while $>9                              /* [◄]  Algorithm is good for 111 digs.*/
      if $>9  then $= left($,1) + substr($,2,1) + substr($,3,1,0)     /*>9?  Reduce it.*/
      end   /*while*/
   if \@dr.$                 then iterate       /*Doesn't have good digital root?  Skip*/
   if iSqrt(sum)**2 \== sum  then iterate       /*Not a perfect square?  Then skip it. */
   if iSqrt(dif)**2 \== dif  then iterate       /* "  "    "       "       "    "   "  */
   #= # + 1;                 call tell          /*bump rare number counter;  display #.*/
   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 _ tell: say right(th(#),length(#)+9) ' rare number is:' right(commas(n),w); 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

         do while q>1; q=q%4; _= x-$-q;  $= $%2;  if _>=0  then do;      x=_;  $=$+q; end
         end   /*while q>1*/;                     return $</lang>

output when using the input of: 8

       1st  rare number is:                           65
       2nd  rare number is:                      621,770
       3rd  rare number is:                  281,089,082
       4th  rare number is:                2,022,652,202
       5th  rare number is:                2,042,832,002
       6th  rare number is:              868,591,084,757
       7th  rare number is:              872,546,974,178
       8th  rare number is:              872,568,754,178

Visual Basic .NET

Translation of: C#

via

Translation of: Go

Surprisingly slow, I expected performance to be a little slower than C#, but this is quite a bit slower. This vb.net version takes 1 2/3 minutes to do what the C# version can do in 2/3 of a minute.

<lang vbnet>Imports System.Console Imports DT = System.DateTime Imports Lsb = System.Collections.Generic.List(Of SByte) Imports Lst = System.Collections.Generic.List(Of System.Collections.Generic.List(Of SByte)) Imports UI = System.UInt64

Module Module1

   Const MxD As SByte = 15

   Public Structure term
       Public coeff As UI : Public a, b As SByte
       Public Sub New(ByVal c As UI, ByVal a_ As Integer, ByVal b_ As Integer)
           coeff = c : a = CSByte(a_) : b = CSByte(b_)
       End Sub
   End Structure

   Dim nd, nd2, count As Integer, digs, cnd, di As Integer()
   Dim res As List(Of UI), st As DT, tLst As List(Of List(Of term))
   Dim lists As List(Of Lst), fml, dmd As Dictionary(Of Integer, Lst)
   Dim dl, zl, el, ol, il As Lsb, odd As Boolean, ixs, dis As Lst, Dif As UI

   ' converts digs array to the "difference"
   Function ToDif() As UI
       Dim r As UI = 0 : For i As Integer = 0 To digs.Length - 1 : r = r * 10 + digs(i)
       Next : Return r
   End Function

   ' converts digs array to the "sum"
   Function ToSum() As UI
       Dim r As UI = 0 : For i As Integer = digs.Length - 1 To 0 Step -1 : r = r * 10 + digs(i)
       Next : Return Dif + (r << 1)
   End Function

   '  determines if the nmbr is square or not
   Function IsSquare(nmbr As UI) As Boolean
       If (&H202021202030213 And (1UL << (nmbr And 63))) <> 0 Then _
           Dim r As UI = Math.Sqrt(nmbr) : Return r * r = nmbr Else Return False
   End Function

   '// returns sequence of SBbytes
   Function Seq(from As SByte, upto As Integer, Optional stp As SByte = 1) As Lsb
       Dim res As Lsb = New Lsb()
       For item As SByte = from To upto Step stp : res.Add(item) : Next : Return res
   End Function

   ' Recursive closure to generate (n+r) candidates from (n-r) candidates
   Sub Fnpr(ByVal lev As Integer)
       If lev = dis.Count Then
           digs(ixs(0)(0)) = fml(cnd(0))(di(0))(0) : digs(ixs(0)(1)) = fml(cnd(0))(di(0))(1)
           Dim le As Integer = di.Length, i As Integer = 1
           If odd Then le -= 1 : digs(nd >> 1) = di(le)
           For Each d As SByte In di.Skip(1).Take(le - 1)
               digs(ixs(i)(0)) = dmd(cnd(i))(d)(0)
               digs(ixs(i)(1)) = dmd(cnd(i))(d)(1) : i += 1 : Next
           If Not IsSquare(ToSum()) Then Return
           res.Add(ToDif()) : count += 1
           WriteLine("{0,16:n0}{1,4}   ({2:n0})", (DT.Now - st).TotalMilliseconds, count, res.Last())
       Else
           For Each n In dis(lev) : di(lev) = n : Fnpr(lev + 1) : Next
       End If
   End Sub

   ' Recursive closure to generate (n-r) candidates with a given number of digits.
   Sub Fnmr(ByVal list As Lst, ByVal lev As Integer)
       If lev = list.Count Then
           Dif = 0 : Dim i As SByte = 0 : For Each t In tLst(nd2)
               If cnd(i) < 0 Then Dif -= t.coeff * CULng(-cnd(i)) _
                             Else Dif += t.coeff * CULng(cnd(i))
               i += 1 : Next
           If Dif <= 0 OrElse Not IsSquare(Dif) Then Return
           dis = New Lst From {Seq(0, fml(cnd(0)).Count - 1)}
           For Each i In cnd.Skip(1) : dis.Add(Seq(0, dmd(i).Count - 1)) : Next
           If odd Then dis.Add(il)
           di = New Integer(dis.Count - 1) {} : Fnpr(0)
       Else
           For Each n As SByte In list(lev) : cnd(lev) = n : Fnmr(list, lev + 1) : Next
       End If
   End Sub

   Sub init()
       Dim pow As UI = 1
       ' terms of (n-r) expression for number of digits from 2 to maxDigits
       tLst = New List(Of List(Of term))() : For Each r As Integer In Seq(2, MxD)
           Dim terms As List(Of term) = New List(Of term)()
           pow *= 10 : Dim p1 As UI = pow, p2 As UI = 1
           Dim i1 As Integer = 0, i2 As Integer = r - 1
           While i1 < i2 : terms.Add(New term(p1 - p2, i1, i2))
               p1 = p1 / 10 : p2 = p2 * 10 : i1 += 1 : i2 -= 1 : End While
           tLst.Add(terms) : Next
       ' map of first minus last digits for 'n' to pairs giving this value
       fml = New Dictionary(Of Integer, Lst)() From {
           {0, New Lst() From {New Lsb() From {2, 2}, New Lsb() From {8, 8}}},
           {1, New Lst() From {New Lsb() From {6, 5}, New Lsb() From {8, 7}}},
           {4, New Lst() From {New Lsb() From {4, 0}}},
           {6, New Lst() From {New Lsb() From {6, 0}, New Lsb() From {8, 2}}}}
       ' map of other digit differences for 'n' to pairs giving this value
       dmd = New Dictionary(Of Integer, Lst)()
       For i As SByte = 0 To 10 - 1 : Dim j As SByte = 0, d As SByte = i
           While j < 10 : If dmd.ContainsKey(d) Then dmd(d).Add(New Lsb From {i, j}) _
               Else dmd(d) = New Lst From {New Lsb From {i, j}}
               j += 1 : d -= 1 : End While : Next
       dl = Seq(-9, 9)    ' all  differences
       zl = Seq(0, 0)     ' zero difference
       el = Seq(-8, 8, 2) ' even differences
       ol = Seq(-9, 9, 2) ' odd  differences
       il = Seq(0, 9)
       lists = New List(Of Lst)()
       For Each f As SByte In fml.Keys : lists.Add(New Lst From {New Lsb From {f}}) : Next
   End Sub

   Sub Main(ByVal args As String())
       init() : res = New List(Of UI)() : st = DT.Now : count = 0
       WriteLine("{0,5}{1,12}{2,4}{3,14}", "digs", "elapsed(ms)", "R/N", "Rare Numbers")
       nd = 2 : nd2 = 0 : odd = False : While nd <= MxD
           digs = New Integer(nd - 1) {} : If nd = 4 Then
               lists(0).Add(zl) : lists(1).Add(ol) : lists(2).Add(el) : lists(3).Add(ol)
           ElseIf tLst(nd2).Count > lists(0).Count Then
               For Each list As Lst In lists : list.Add(dl) : Next : End If
           ixs = New Lst() : For Each t As term In tLst(nd2) : ixs.Add(New Lsb From {t.a, t.b}) : Next
           For Each list As Lst In lists : cnd = New Integer(list.Count - 1) {} : Fnmr(list, 0) : Next
           WriteLine("  {0,2}  {1,10:n0}", nd, (DT.Now - st).TotalMilliseconds)
           nd += 1 : nd2 += 1 : odd = Not odd : End While
       res.Sort() : WriteLine(vbLf & "The {0} rare numbers with up to {1} digits are:", res.Count, MxD)
       count = 0 : For Each rare In res : count += 1 : WriteLine("{0,2}:{1,27:n0}", count, rare) : Next
       If System.Diagnostics.Debugger.IsAttached Then ReadKey()
   End Sub

End Module</lang>

Output:

 digs elapsed(ms) R/N  Rare Numbers
              25   1   (65)
   2          26
   3          26
   4          27
   5          27
              28   2   (621,770)
   6          29
   7          30
   8          41
              42   3   (281,089,082)
   9          46
              47   4   (2,022,652,202)
             116   5   (2,042,832,002)
  10         273
  11         422
           1,363   6   (872,546,974,178)
           1,476   7   (872,568,754,178)
           2,937   8   (868,591,084,757)
  12       3,584
           4,560   9   (6,979,302,951,885)
  13       5,817
          18,234  10   (20,313,693,904,202)
          18,471  11   (20,313,839,704,202)
          23,626  12   (20,331,657,922,202)
          24,454  13   (20,331,875,722,202)
          26,599  14   (20,333,875,702,202)
          60,784  15   (40,313,893,704,200)
          61,246  16   (40,351,893,720,200)
  14      65,387
          65,465  17   (200,142,385,731,002)
          66,225  18   (221,462,345,754,122)
          76,417  19   (816,984,566,129,618)
          81,727  20   (245,518,996,076,442)
          82,461  21   (204,238,494,066,002)
          82,694  22   (248,359,494,187,442)
          83,729  23   (244,062,891,224,042)
          99,241  24   (403,058,392,434,500)
         100,009  25   (441,054,594,034,340)
  15     104,207

The 25 rare numbers with up to 15 digits are:
 1:                         65
 2:                    621,770
 3:                281,089,082
 4:              2,022,652,202
 5:              2,042,832,002
 6:            868,591,084,757
 7:            872,546,974,178
 8:            872,568,754,178
 9:          6,979,302,951,885
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
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