Rare numbers
Rare numbers are positive integers n where:
You are encouraged to solve this task according to the task description, using any language you may know.
- Definitions and restrictions
- 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#
Traditional
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
Quicker
Along the lines of the C++ version. Computing the possible squares for the sums and differences, then comparing them together and reporting the ones that have a proper forward, reverse result. To save computation time, some shortcuts have been taken to avoid generating many non-square numbers.
Update, added computation of digital root, which increased performance a few percentage points. Interestingly, the digital root is always zero for the difference list of squares, so no advantage is given by tracking it. Only the sum list of squares benefits from calculating the digital root of the partial sum and using an abbreviated sequence for the last round of permutation.
<lang csharp>using static System.Math; // for Sqrt()
using System.Collections.Generic; // for List<>, .Count
using System.Linq; // for .Last(), .ToList()
using System.Diagnostics; // for Stopwatch()
using static System.Console; // for Write(), WriteLine()
using llst = System.Collections.Generic.List<int[]>;
class Program
{
#region vars static int[] d, // permutation working array drar = new int[19], // digital root lookup array dac; // running digital root array static long[] p = new long[20], // powers of 10 ac, // accumulator array pp; // long coefficient array that combines with digits of working array static bool odd = false; // flag for odd number of digits static long sum, // calculated sum of terms (square candidate) rt; // root of sum static int cn = 0, // solution counter nd = 2, // number of digits nd1 = nd - 1, // nd helper ln, // previous value of "n" (in Recurse()) dl; // length of "d" array; static Stopwatch sw = new Stopwatch(), swt = new Stopwatch(); // for timings static List<long> sr = new List<long>(); // temporary list of squares used for building static readonly int[] tlo = new int[] { 0, 1, 4, 5, 6 }, // primary differences starting point all = Seq(-9, 9), // all possible differences odl = Seq(-9, 9, 2), // odd possible differences evl = Seq(-8, 8, 2), // even possible differences thi = new int[] { 4, 5, 6, 9, 10, 11, 14, 15, 16 }, // primary sums staring point. note: (0, 1) omitted, as any square generated will not have enough digits alh = Seq(0, 18), // all possible sums odh = Seq(1, 17, 2), // odd possible sums evh = Seq(0, 18, 2), // even possible sums ten = Seq(0, 9), // used for odd number of digits z = Seq(0, 0), // no difference, used to avoid generating a bunch of negative square candidates t7 = new int[] { -3, 7 }, // shortcut for low 5 nin = new int[] { 9 }, // shortcut for hi 10 tn = new int[] { 10 }, // shortcut for hi 0 (unused, uneeded) t12 = new int[] { 2, 12 }, // shortcut for hi 5 o11 = new int[] { 1, 11 }, // shortcut for hi 15 pos = new int[] { 0, 1, 4, 5, 6, 9 }; // shortcut for 2nd lo 0 static llst lul = new llst { z, odl, null, null, evl, t7, odl }, // shortcut lookup lo primary luh = new llst { tn, evh, null, null, evh, t12, odh, null, null, evh, nin, odh, null, null, odh, o11, evh }, // shortcut lookup hi primary l2l = new llst { pos, null, null, null, all, null, all }, // shortcut lookup lo secondary l2h = new llst { null, null, null, null, alh, null, alh, null, null, null, alh, null, null, null, alh, null, alh }, lu, l2; // shortcut lookup hi secondary static int[][] chTen = new int[][] { new int[] { 0,2,5,8,9 }, new int[] { 0,3,4,6,9 }, new int[] { 1,4,7,8 }, new int[] { 2,3,5,8 }, new int[] { 0,3,6,7,9 }, new int[] { 1,2,4,7 }, new int[] { 2,5,6,8 }, new int[] { 0,1,3,6,9 }, new int[] { 1,4,5,7 } }; static int[][] chAH = new int[][] { new int[] { 0,2,5,8,9,11,14,17,18 }, new int[] { 0,3,4,6,9,12,13,15,18 }, new int[] { 1,4,7,8,10,13,16,17 }, new int[] { 2,3,5,8,11,12,14,17 }, new int[] { 0,3,6,7,9,12,15,16,18 }, new int[] { 1,2,4,7,10,11,13,16 }, new int[] { 2,5,6,8,11,14,15,17 }, new int[] { 0,1,3,6,9,10,12,15,18 }, new int[] { 1,4,5,7,10,13,14,16 } }; #endregion vars // Returns a sequence of integers static int[] Seq(int f, int t, int s = 1) { int[] r = new int[(t - f) / s + 1]; for (int i = 0; i < r.Length; i++, f += s) r[i] = f; return r; } // Returns Integer Square Root static long ISR(long s) { return (long)Sqrt(s); }
// Recursively determines whether "r" is the reverse of "f" static bool IsRev(int nd, long f, long r) { nd--; return f / p[nd] != r % 10 ? false : (nd < 1 ? true : IsRev(nd, f % p[nd], r / 10)); } // Recursive procedure to evaluate the permutations, no shortcuts static void RecurseLE5(llst lst, int lv) { if (lv == dl) { // check if on last stage of permutation if ((sum = ac[lv - 1]) > 0) if ((rt = (long)Sqrt(sum)) * rt == sum) sr.Add(sum); } // test accumulated sum, append to result if square else foreach (int n in lst[lv]) { // set up next permutation d[lv] = n; if (lv == 0) ac[0] = pp[0] * n; else ac[lv] = ac[lv - 1] + pp[lv] * n; // update accumulated sum RecurseLE5(lst, lv + 1); } } // Recursively call next level // Recursive procedure to evaluate the hi permutations, shortcuts added to avoid generating many non-squares, digital root calc added static void Recursehi(llst lst, int lv) { int lv1 = lv - 1; if (lv == dl) { // check if on last stage of permutation if ((0x202021202030213 & (1 << (int)((sum = ac[lv1]) & 63))) != 0) // test accumulated sum, append to result if square if ((rt = (long)Sqrt(sum)) * rt == sum) sr.Add(sum); } else foreach (int n in lst[lv]) { // set up next permutation d[lv] = n; if (lv == 0) { ac[0] = pp[0] * n; dac[0] = drar[n]; } // update accumulated sum and running dr else { ac[lv] = ac[lv1] + pp[lv] * n; dac[lv] = dac[lv1] + drar[n]; if (dac[lv] > 8) dac[lv] -= 9; } switch (lv) { // shortcuts to be performed on designated levels case 0: lst[1] = lu[ln = n]; lst[2] = l2[n]; break; // primary level: set shortcuts for secondary level case 1: // secondary level: set shortcuts for tertiary level switch (ln) { // for sums case 5: case 15: lst[2] = n < 10 ? evh : odh; break; case 9: lst[2] = ((n >> 1) & 1) == 0 ? evh : odh; break; case 11: lst[2] = ((n >> 1) & 1) == 1 ? evh : odh; break; } break; } if (lv == dl - 2) lst[dl - 1] = odd ? chTen[dac[dl - 2]] : chAH[dac[dl - 2]]; // reduce last round according to dr calc Recursehi(lst, lv + 1); } } // Recursively call next level // Recursive procedure to evaluate the lo permutations, shortcuts added to avoid generating many non-squares static void Recurselo(llst lst, int lv) { int lv1 = lv - 1; if (lv == dl) { // check if on last stage of permutation if ((sum = ac[lv1]) > 0) if ((rt = (long)Sqrt(sum)) * rt == sum) sr.Add(sum); } // test accumulated sum, append to result if square else foreach (int n in lst[lv]) { // set up next permutation d[lv] = n; if (lv == 0) ac[0] = pp[0] * n; else ac[lv] = ac[lv1] + pp[lv] * n; // update accumulated sum switch (lv) { // shortcuts to be performed on designated levels case 0: lst[1] = lu[ln = n]; lst[2] = l2[n]; break; // primary level: set shortcuts for secondary level case 1: // secondary level: set shortcuts for tertiary level switch (ln) { // for difs case 1: lst[2] = (((n + 9) >> 1) & 1) == 0 ? evl : odl; break; case 5: lst[2] = n < 0 ? evl : odl; break; } break; } Recurselo(lst, lv + 1); } } // Recursively call next level // Produces a list of candidate square numbers static List<long> listEm(llst lst, llst plu, llst pl2) { d = new int[dl = lst.Count]; sr.Clear(); lu = plu; l2 = pl2; ac = new long[dl]; dac = new int[dl]; // init support vars pp = new long[dl]; for (int i = 0, j = nd1; i < dl; i++, j--) pp[i] = lst[0].Length > 6 ? p[j] + p[i] : p[j] - p[i]; // build coefficients array if (nd <= 5) RecurseLE5(lst, 0); else { if (lst[0].Length > 8) Recursehi(lst, 0); else Recurselo(lst, 0); } return sr; } // call appropriate recursive procedure // Reveals whether combining two lists of squares can produce a Rare number static void Reveal(List<long> lo, List<long> hi) { List<string> s = new List<string>(); // create temp list of results foreach (long l in lo) foreach (long h in hi) { long r = (h - l) >> 1, f = h - r; // generate all possible fwd & rev candidates from lists if (IsRev(nd, f, r)) s.Add(string.Format("{0,20} {1,11} {2,10} ", f, ISR(h), ISR(l))); } // test and append sucesses to temp list s.Sort(); if (s.Count > 0) foreach (string t in s) // if there are any, output sorted results Write("{0,2} {1}{2}", ++cn, t, t == s.Last() ? "" : "\n"); else Write("{0,48}", ""); } static void Main(string[] args) { WriteLine("{0,3}{1,20} {2,11} {3,10} {4,4}{5,16} {6, 17}", "nth", "forward", "rt.sum", "rt.dif", "digs", "block time", "total time"); p[0] = 1; for (int i = 0, j = 1; j < p.Length; i = j++) p[j] = p[i] * 10; // create powers of 10 array for (int i = 0; i < drar.Length; i++) drar[i] = (i << 1) % 9; // create digital root array llst lls = new llst { tlo }, hls = new llst { thi }; sw.Start(); swt.Start(); // initialize permutations list, timers for (; nd <= 18; nd1 = nd++, odd = !odd) { // loop through all numbers of digits if (nd > 2) if (odd) hls.Add(ten); else { lls.Add(all); hls[hls.Count - 1] = alh; } // build permutations list Reveal(listEm(lls, lul, l2l).ToList(), listEm(hls, luh, l2h)); // reveal results if (!odd && nd > 5) hls[hls.Count - 1] = alh; // restore last element of hls, so that dr shortcut doesn't mess up next nd WriteLine("{0,2}: {1} {2}", nd, sw.Elapsed, swt.Elapsed); sw.Restart(); } // 19 hls.Add(ten); Reveal(listEmU(lls, lul, l2l).ToList(), listEmU(hls, luh, l2h)); // reveal unsigned results WriteLine("{0,2}: {1} {2}", nd, sw.Elapsed, swt.Elapsed); } #region 19 static ulong usum, // unsigned calculated sum of terms (square candidate) urt; // unsigned root of sum static ulong[] acu, // unsigned accumulator array ppu; // unsigned long coefficient array that combines with digits of working array static List<ulong> sru = new List<ulong>(); // unsigned temporary list of squares used for building // Reveals whether combining two lists of unsigned squares can produce a Rare number static void Reveal(List<ulong> lo, List<ulong> hi) { List<string> s = new List<string>(); // create temp list of results foreach (ulong l in lo) foreach (ulong h in hi) { ulong r = (h - l) >> 1, f = h - r; // generate all possible fwd & rev candidates from lists if (IsRev(nd, f, r)) s.Add(string.Format("{0,20} {1,11} {2,10} ", f, ISR(h), ISR(l))); } // test and append sucesses to temp list s.Sort(); if (s.Count > 0) foreach (string t in s) // if there are any, output sorted results Write("{0,2} {1}{2}", ++cn, t, t == s.Last() ? "" : "\n"); else Write("{0,48}", ""); } // Produces a list of unsigned candidate square numbers static List<ulong> listEmU(llst lst, llst plu, llst pl2) { d = new int[dl = lst.Count]; sru.Clear(); lu = plu; l2 = pl2; acu = new ulong[dl]; dac = new int[dl]; // init support vars ppu = new ulong[dl]; for (int i = 0, j = nd1; i < dl; i++, j--) ppu[i] = (ulong)(lst[0].Length > 6 ? p[j] + p[i] : p[j] - p[i]); // build coefficients array if (lst[0].Length > 8) RecurseUhi(lst, 0); else RecurseUlo(lst, 0); return sru; } // call recursive procedure // Recursive procedure to evaluate the unsigned hi permutations, shortcuts added to avoid generating many non-squares, digital root calc added static void RecurseUhi(llst lst, int lv) { int lv1 = lv - 1; if (lv == dl) { // check if on last stage of permutation if ((0x202021202030213 & (1 << (int)((usum = acu[lv1]) & 63))) != 0) // test accumulated sum, append to result if square if ((urt = (ulong)Sqrt(usum)) * urt == usum) sru.Add(usum); } else foreach (int n in lst[lv]) { // set up next permutation d[lv] = n; if (lv == 0) { acu[0] = ppu[0] * (uint)n; dac[0] = drar[n]; } // update accumulated sum and running dr else { acu[lv] = n >= 0 ? acu[lv1] + ppu[lv] * (uint)n : acu[lv1] - ppu[lv] * (uint)-n; dac[lv] = dac[lv1] + drar[n]; if (dac[lv] > 8) dac[lv] -= 9; } switch (lv) { // shortcuts to be performed on designated levels case 0: lst[1] = lu[ln = n]; lst[2] = l2[n]; break; // primary level: set shortcuts for secondary level case 1: // secondary level: set shortcuts for tertiary level switch (ln) { // for sums case 5: case 15: lst[2] = n < 10 ? evh : odh; break; case 9: lst[2] = ((n >> 1) & 1) == 0 ? evh : odh; break; case 11: lst[2] = ((n >> 1) & 1) == 1 ? evh : odh; break; } break; } if (lv == dl - 2) lst[dl - 1] = odd ? chTen[dac[dl - 2]] : chAH[dac[dl - 2]]; // reduce last round according to dr calc RecurseUhi(lst, lv + 1); } } // Recursively call next level // Recursive procedure to evaluate the unsigned lo permutations, shortcuts added to avoid generating many non-squares static void RecurseUlo(llst lst, int lv) { int lv1 = lv - 1; if (lv == dl) { // check if on last stage of permutation if ((usum = acu[lv1]) > 0) if ((urt = (ulong)Sqrt(usum)) * urt == usum) sru.Add(usum); } // test accumulated sum, append to result if square else foreach (int n in lst[lv]) { // set up next permutation d[lv] = n; if (lv == 0) acu[0] = ppu[0] * (uint)n; else acu[lv] = n >= 0 ? acu[lv1] + ppu[lv] * (uint)n : acu[lv1] - ppu[lv] * (uint)-n; // update accumulated sum switch (lv) { // shortcuts to be performed on designated levels case 0: lst[1] = lu[ln = n]; lst[2] = l2[n]; break; // primary level: set shortcuts for secondary level case 1: // secondary level: set shortcuts for tertiary level switch (ln) { // for difs case 1: lst[2] = (((n + 9) >> 1) & 1) == 0 ? evl : odl; break; case 5: lst[2] = n < 0 ? evl : odl; break; } break; } RecurseUlo(lst, lv + 1); } } // Recursively call next level // Returns unsigned Integer Square Root static ulong ISR(ulong s) { return (ulong)Sqrt(s); } // Recursively determines whether "r" is the reverse of "f" static bool IsRev(int nd, ulong f, ulong r) { nd--; return f / (ulong)p[nd] != r % 10 ? false : (nd < 1 ? true : IsRev(nd, f % (ulong)p[nd], r / 10)); } #endregion 19
}</lang>
- Output:
Results on the core i7-7700 @ 3.6Ghz.
nth forward rt.sum rt.dif digs block time total time 1 65 11 3 2: 00:00:00.0030626 00:00:00.0030626 3: 00:00:00.0001018 00:00:00.0033254 4: 00:00:00.0000963 00:00:00.0035054 5: 00:00:00.0000928 00:00:00.0036834 2 621770 836 738 6: 00:00:00.0021741 00:00:00.0059392 7: 00:00:00.0001724 00:00:00.0061956 8: 00:00:00.0002609 00:00:00.0065384 3 281089082 23708 330 9: 00:00:00.0012672 00:00:00.0079061 4 2022652202 63602 300 5 2042832002 63602 6360 10: 00:00:00.0045628 00:00:00.0125626 11: 00:00:00.0201361 00:00:00.0328037 6 868591084757 1275175 333333 7 872546974178 1320616 32670 8 872568754178 1320616 33330 12: 00:00:00.0519065 00:00:00.0848320 9 6979302951885 3586209 1047717 13: 00:00:00.3772503 00:00:00.4622089 10 20313693904202 6368252 269730 11 20313839704202 6368252 270270 12 20331657922202 6368252 329670 13 20331875722202 6368252 330330 14 20333875702202 6368252 336330 15 40313893704200 6368252 6330336 16 40351893720200 6368252 6336336 14: 00:00:00.9416903 00:00:01.4041338 17 200142385731002 20006998 69300 18 204238494066002 20122102 1891560 19 221462345754122 21045662 69300 20 244062891224042 22011022 1908060 21 245518996076442 22140228 921030 22 248359494187442 22206778 1891560 23 403058392434500 20211202 19940514 24 441054594034340 22011022 19940514 25 816984566129618 40421606 250800 15: 00:00:07.0248881 00:00:08.4296936 26 2078311262161202 64030648 7529850 27 2133786945766212 65272218 2666730 28 2135568943984212 65272218 3267330 29 2135764587964212 65272218 3326670 30 2135786765764212 65272218 3333330 31 4135786945764210 65272218 63333336 32 6157577986646405 105849161 33333333 33 6889765708183410 83866464 82133718 34 8052956026592517 123312255 29999997 35 8052956206592517 123312255 30000003 36 8191154686620818 127950856 3299670 37 8191156864620818 127950856 3300330 38 8191376864400818 127950856 3366330 39 8650327689541457 127246955 33299667 40 8650349867341457 127246955 33300333 16: 00:00:18.1046570 00:00:26.5344137 41 22542040692914522 212329862 333300 42 67725910561765640 269040196 251135808 43 86965750494756968 417050956 33000 17: 00:02:11.8544100 00:02:38.3889020 44 225342456863243522 671330638 297000 45 225342458663243522 671330638 303000 46 225342478643243522 671330638 363000 47 284684666566486482 754565658 30000 48 284684868364486482 754565658 636000 49 297128548234950692 770186978 32697330 50 297128722852950692 770186978 32702670 51 297148324656930692 770186978 33296670 52 297148546434930692 770186978 33303330 53 497168548234910690 770186978 633363336 54 619431353040136925 1071943279 299667003 55 619631153042134925 1071943279 300333003 56 631688638047992345 1083968809 297302703 57 633288858025996145 1083968809 302637303 58 633488632647994145 1083968809 303296697 59 653488856225994125 1083968809 363303363 60 811865096390477018 1273828556 33030330 61 865721270017296468 1315452006 32071170 62 871975098681469178 1320582934 3303300 63 898907259301737498 1339270086 64576740 18: 00:05:38.5737725 00:08:16.9627994 64 2042401829204402402 2021001202 18915600 65 2060303819041450202 2020110202 199405140 66 2420424089100600242 2200110022 19080600 67 2551755006254571552 2259094848 693000 68 2702373360882732072 2324811012 693000 69 2825378427312735282 2377130742 2508000 70 6531727101458000045 3454234451 1063822617 71 6988066446726832640 2729551744 2554541088 72 8066308349502036608 4016542096 2508000 73 8197906905009010818 4046976144 133408770 74 8200756128308135597 4019461925 495417087 75 8320411466598809138 4079154376 36366330 19: 00:42:31.7490390 00:50:48.7120790
D
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
<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
Kotlin
<lang scala>import java.time.Duration import java.time.LocalDateTime import kotlin.math.sqrt
class Term(var coeff: Long, var ix1: Byte, var ix2: Byte)
const val maxDigits = 16
fun toLong(digits: List<Byte>, reverse: Boolean): Long {
var sum: Long = 0 if (reverse) { var i = digits.size - 1 while (i >= 0) { sum = sum * 10 + digits[i] i-- } } else { var i = 0 while (i < digits.size) { sum = sum * 10 + digits[i] i++ } } return sum
}
fun isSquare(n: Long): Boolean {
val root = sqrt(n.toDouble()).toLong() return root * root == n
}
fun seq(from: Byte, to: Byte, step: Byte): List<Byte> {
val res = mutableListOf<Byte>() var i = from while (i <= to) { res.add(i) i = (i + step).toByte() } return res
}
fun commatize(n: Long): String {
var s = n.toString() val le = s.length var i = le - 3 while (i >= 1) { s = s.slice(0 until i) + "," + s.substring(i) i -= 3 } return s
}
fun main() {
val startTime = LocalDateTime.now() var pow = 1L println("Aggregate timings to process all numbers up to:") // terms of (n-r) expression for number of digits from 2 to maxDigits val allTerms = mutableListOf<MutableList<Term>>() for (i in 0 until maxDigits - 1) { allTerms.add(mutableListOf()) } for (r in 2..maxDigits) { val terms = mutableListOf<Term>() pow *= 10 var pow1 = pow var pow2 = 1L var i1: Byte = 0 var i2 = (r - 1).toByte() while (i1 < i2) { terms.add(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 val fml = mapOf( 0.toByte() to listOf(listOf<Byte>(2, 2), listOf<Byte>(8, 8)), 1.toByte() to listOf(listOf<Byte>(6, 5), listOf<Byte>(8, 7)), 4.toByte() to listOf(listOf<Byte>(4, 0)), 6.toByte() to listOf(listOf<Byte>(6, 0), listOf<Byte>(8, 2)) ) // map of other digit differences for 'n' to pairs giving this value val dmd = mutableMapOf<Byte, MutableList<List<Byte>>>() for (i in 0 until 100) { val a = listOf((i / 10).toByte(), (i % 10).toByte()) val d = a[0] - a[1] dmd.getOrPut(d.toByte(), { mutableListOf() }).add(a) } val fl = listOf<Byte>(0, 1, 4, 6) val dl = seq(-9, 9, 1) // all differences val zl = listOf<Byte>(0) // zero differences only val el = seq(-8, 8, 2) // even differences only val ol = seq(-9, 9, 2) // odd differences only val il = seq(0, 9, 1) val rares = mutableListOf<Long>() val lists = mutableListOf<MutableList<List<Byte>>>() for (i in 0 until 4) { lists.add(mutableListOf()) } for (i_f in fl.withIndex()) { lists[i_f.index] = mutableListOf(listOf(i_f.value)) } var digits = mutableListOf<Byte>() var count = 0
// Recursive closure to generate (n+r) candidates from (n-r) candidates // and hence find Rare numbers with a given number of digits. fun fnpr( cand: List<Byte>, di: MutableList<Byte>, dis: List<List<Byte>>, indicies: List<List<Byte>>, nmr: Long, nd: Int, level: Int ) { if (level == dis.size) { digits[indicies[0][0].toInt()] = fml[cand[0]]?.get(di[0].toInt())?.get(0)!! digits[indicies[0][1].toInt()] = fml[cand[0]]?.get(di[0].toInt())?.get(1)!! var le = di.size if (nd % 2 == 1) { le-- digits[nd / 2] = di[le] } for (i_d in di.slice(1 until le).withIndex()) { digits[indicies[i_d.index + 1][0].toInt()] = dmd[cand[i_d.index + 1]]?.get(i_d.value.toInt())?.get(0)!! digits[indicies[i_d.index + 1][1].toInt()] = dmd[cand[i_d.index + 1]]?.get(i_d.value.toInt())?.get(1)!! } val r = toLong(digits, true) val npr = nmr + 2 * r if (!isSquare(npr)) { return } count++ print(" R/N %2d:".format(count)) val checkPoint = LocalDateTime.now() val elapsed = Duration.between(startTime, checkPoint).toMillis() print(" %9sms".format(elapsed)) val n = toLong(digits, false) println(" (${commatize(n)})") rares.add(n) } else { for (num in 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. fun fnmr(cand: MutableList<Byte>, list: List<List<Byte>>, indicies: List<List<Byte>>, nd: Int, level: Int) { if (level == list.size) { var nmr = 0L var nmr2 = 0L for (i_t in allTerms[nd - 2].withIndex()) { if (cand[i_t.index] >= 0) { nmr += i_t.value.coeff * cand[i_t.index] } else { nmr2 += i_t.value.coeff * -cand[i_t.index] if (nmr >= nmr2) { nmr -= nmr2 nmr2 = 0 } else { nmr2 -= nmr nmr = 0 } } } if (nmr2 >= nmr) { return } nmr -= nmr2 if (!isSquare(nmr)) { return } val dis = mutableListOf<List<Byte>>() dis.add(seq(0, ((fml[cand[0]] ?: error("oops")).size - 1).toByte(), 1)) for (i in 1 until cand.size) { dis.add(seq(0, (dmd[cand[i]]!!.size - 1).toByte(), 1)) } if (nd % 2 == 1) { dis.add(il) } val di = mutableListOf<Byte>() for (i in 0 until dis.size) { di.add(0) } fnpr(cand, di, dis, indicies, nmr, nd, 0) } else { for (num in list[level]) { cand[level] = num fnmr(cand, list, indicies, nd, level + 1) } } }
for (nd in 2..maxDigits) { digits = mutableListOf() for (i in 0 until nd) { digits.add(0) } if (nd == 4) { lists[0].add(zl) lists[1].add(ol) lists[2].add(el) lists[3].add(ol) } else if (allTerms[nd - 2].size > lists[0].size) { for (i in 0 until 4) { lists[i].add(dl) } } val indicies = mutableListOf<List<Byte>>() for (t in allTerms[nd - 2]) { indicies.add(listOf(t.ix1, t.ix2)) } for (list in lists) { val cand = mutableListOf<Byte>() for (i in 0 until list.size) { cand.add(0) } fnmr(cand, list, indicies, nd, 0) } val checkPoint = LocalDateTime.now() val elapsed = Duration.between(startTime, checkPoint).toMillis() println(" %2d digits: %9sms".format(nd, elapsed)) }
rares.sort() println("\nThe rare numbers with up to $maxDigits digits are:") for (i_rare in rares.withIndex()) { println(" %2d: %25s".format(i_rare.index + 1, commatize(i_rare.value))) }
}</lang>
- Output:
Aggregate timings to process all numbers up to: R/N 1: 130ms (65) 2 digits: 133ms 3 digits: 133ms 4 digits: 135ms 5 digits: 136ms R/N 2: 155ms (621,770) 6 digits: 171ms 7 digits: 176ms 8 digits: 238ms R/N 3: 243ms (281,089,082) 9 digits: 266ms R/N 4: 272ms (2,022,652,202) R/N 5: 432ms (2,042,832,002) 10 digits: 693ms 11 digits: 1037ms R/N 6: 1690ms (872,546,974,178) R/N 7: 1757ms (872,568,754,178) R/N 8: 2380ms (868,591,084,757) 12 digits: 2682ms R/N 9: 3081ms (6,979,302,951,885) 13 digits: 3612ms R/N 10: 9091ms (20,313,693,904,202) R/N 11: 9180ms (20,313,839,704,202) R/N 12: 11322ms (20,331,657,922,202) R/N 13: 11611ms (20,331,875,722,202) R/N 14: 12477ms (20,333,875,702,202) R/N 15: 26933ms (40,313,893,704,200) R/N 16: 27128ms (40,351,893,720,200) 14 digits: 29696ms R/N 17: 29759ms (200,142,385,731,002) R/N 18: 30024ms (221,462,345,754,122) R/N 19: 33577ms (816,984,566,129,618) R/N 20: 35392ms (245,518,996,076,442) R/N 21: 35662ms (204,238,494,066,002) R/N 22: 35748ms (248,359,494,187,442) R/N 23: 36108ms (244,062,891,224,042) R/N 24: 42484ms (403,058,392,434,500) R/N 25: 42760ms (441,054,594,034,340) 15 digits: 45334ms R/N 26: 106307ms (2,133,786,945,766,212) R/N 27: 130390ms (2,135,568,943,984,212) R/N 28: 134315ms (8,191,154,686,620,818) R/N 29: 137815ms (8,191,156,864,620,818) R/N 30: 139449ms (2,135,764,587,964,212) R/N 31: 141563ms (2,135,786,765,764,212) R/N 32: 146705ms (8,191,376,864,400,818) R/N 33: 163353ms (2,078,311,262,161,202) R/N 34: 204546ms (8,052,956,026,592,517) R/N 35: 209994ms (8,052,956,206,592,517) R/N 36: 251686ms (8,650,327,689,541,457) R/N 37: 254537ms (8,650,349,867,341,457) R/N 38: 256579ms (6,157,577,986,646,405) R/N 39: 307145ms (4,135,786,945,764,210) R/N 40: 347119ms (6,889,765,708,183,410) 16 digits: 350388ms 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
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
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 1st 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
via
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