Rare numbers
You are encouraged to solve this task according to the task description, using any language you may know.
- Definitions and restrictions
Rare numbers are positive integers n where:
- n is expressed in base ten
- r is the reverse of n (decimal digits)
- n must be non-palindromic (n ≠ r)
- (n+r) is the sum
- (n-r) is the difference and must be positive
- the sum and the difference must be perfect squares
- Task
-
- find and show the first 5 rare numbers
- find and show the first 8 rare numbers (optional)
- find and show more rare numbers (stretch goal)
Show all output here, on this page.
- References
-
- an OEIS entry: A035519 rare numbers.
- an OEIS entry: A059755 odd rare numbers.
- planetmath entry: rare numbers. (some hints)
- author's website: rare numbers by Shyam Sunder Gupta. (lots of hints and some observations).
ALGOL 68
which is
(naive version)
PROC revn = ( LONG INT na, nda )LONG INT:
BEGIN
LONG INT n := na, nd := nda, r := 0, i := 0;
WHILE i +:= 1;
i <= nd
DO
r *:= 10 +:= ( n MOD 10 );
n OVERAB 10
OD;
r
END # revn # ;
LONG INT nd := 2, count := 0, lim := 90, n := 20;
DO
n +:= 1;
LONG INT r = revn( n, nd );
IF r < n THEN
LONG INT s = n + r, d = n - r;
IF IF ODD nd
THEN d MOD 1089 = 0
ELSE s MOD 121 = 0
FI
THEN
IF LONG REAL root s = long sqrt( s );
root s = ENTIER root s
THEN
IF LONG REAL root d = long sqrt( d );
root d = ENTIER root d
THEN
count +:= 1;
print( ( whole( count, 0 ), ": ", whole( n, 0 ), newline ) );
IF count >= 5 THEN stop FI
FI
FI
FI;
IF n = lim
THEN
lim *:= 10;
nd +:= 1;
n := ( lim OVER 9 ) * 2
FI
FI
OD- Output:
1: 65 2: 621770 3: 281089082 4: 2022652202 5: 2042832002
C#
Traditional
Converted to unsigned longs in order to reach 19 digits.
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(); }
}
- 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.
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
}
- 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
C++
Calculate L and H independently
// Rare Numbers : Nigel Galloway - December 20th., 2019;
// Nigel Galloway/Enter your username - January 4th., 2021 (see discussion page.
#include <functional>
#include <bitset>
#include <cmath>
using namespace std;
using Z2 = optional<long long>; using Z1 = function<Z2()>;
// powers of 10 array
constexpr auto pow10 = [] { array <long long, 19> n {1}; for (int j{0}, i{1}; i < 19; j = i++) n[i] = n[j] * 10; return n; } ();
long long acc, l;
bool izRev(int n, unsigned long long i, unsigned long long g) {
return (i / pow10[n - 1] != g % 10) ? false : n < 2 ? true : izRev(n - 1, i % pow10[n - 1], g / 10);
}
const Z1 fG(Z1 n, int start, int end, int reset, const long long step, long long &l) {
return [n, i{step * start}, g{step * end}, e{step * reset}, &l, step] () mutable {
while (i<g){i+=step; return Z2(l+=step);}
l-=g-(i=e); return n();};
}
struct nLH {
vector<unsigned long long>even{}, odd{};
nLH(const Z1 a, const vector<long long> b, long long llim){while (auto i = a()) for (auto ng : b)
if(ng>0 | *i>llim){unsigned long long sq{ng+ *i}, r{sqrt(sq)}; if (r*r == sq) ng&1 ? odd.push_back(sq) : even.push_back(sq);}}
};
const double fac = 3.94;
const int mbs = (int)sqrt(fac * pow10[9]), mbt = (int)sqrt(fac * fac * pow10[9]) >> 3;
const bitset<100000>bs {[]{bitset<100000>n{false}; for(int g{3};g<mbs;++g) n[(g*g)%100000]=true; return n;}()};
constexpr array<const int, 7>li{1,3,0,0,1,1,1},lin{0,-7,0,0,-8,-3,-9},lig{0,9,0,0,8,7,9},lil{0,2,0,0,2,10,2};
const nLH makeL(const int n){
constexpr int r{9}; acc=0; Z1 g{[]{return Z2{};}}; int s{-r}, q{(n>11)*5}; vector<long long> w{};
for (int i{1};i<n/2-q+1;++i){l=pow10[n-i-q]-pow10[i+q-1]; s-=i==n/2-q; g=fG(g,s,r,-r,l,acc+=l*s);}
if(q){long long g0{0}, g1{0}, g2{0}, g3{0}, g4{0}, l3{pow10[n-5]}; while (g0<7){const long long g{-10000*g4-1000*g3-100*g2-10*g1-g0};
if (bs[(g+1000000000000LL)%100000]) w.push_back(l3*(g4+g3*10+g2*100+g1*1000+g0*10000)+g);
if(g4<r) ++g4; else{g4= -r; if(g3<r) ++g3; else{g3= -r; if(g2<r) ++g2; else{g2= -r; if(g1<lig[g0]) g1+=lil[g0]; else {g0+=li[g0];g1=lin[g0];}}}}}}
return q ? nLH(g,w,0) : nLH(g,{0},0);
}
const bitset<100000>bt {[]{bitset<100000>n{false}; for(int g{11};g<mbt;++g) n[(g*g)%100000]=true; return n;}()};
constexpr array<const int, 17>lu{0,0,0,0,2,0,4,0,0,0,1,4,0,0,0,1,1},lun{0,0,0,0,0,0,1,0,0,0,9,1,0,0,0,1,0},lug{0,0,0,0,18,0,17,0,0,0,9,17,0,0,0,11,18},lul{0,0,0,0,2,0,2,0,0,0,0,2,0,0,0,10,2};
const nLH makeH(const int n){
acc= -pow10[n>>1]-pow10[(n-1)>>1]; Z1 g{[]{ return Z2{};}}; int q{(n>11)*5}; vector<long long> w {};
for (int i{1}; i<(n>>1)-q+1; ++i) g = fG(g,0,18,0,pow10[n-i-q]+pow10[i+q-1], acc);
if (n & 1){l=pow10[n>>1]<<1; g=fG(g,0,9,0,l,acc+=l);}
if(q){long long g0{4}, g1{0}, g2{0}, g3{0}, g4{0},l3{pow10[n-5]}; while (g0<17){const long long g{g4*10000+g3*1000+g2*100+g1*10+g0};
if (bt[g%100000]) w.push_back(l3*(g4+g3*10+g2*100+g1*1000+g0*10000)+g);
if (g4<18) ++g4; else{g4=0; if(g3<18) ++g3; else{g3=0; if(g2<18) ++g2; else{g2=0; if(g1<lug[g0]) g1+=lul[g0]; else{g0+=lu[g0];g1=lun[g0];}}}}}}
return q ? nLH(g,w,0) : nLH(g,{0},pow10[n-1]<<2);
}
#include <chrono>
using namespace chrono; using VU = vector<unsigned long long>; using VS = vector<string>;
template <typename T> // concatenates vectors
vector<T>& operator +=(vector<T>& v, const vector<T>& w) { v.insert(v.end(), w.begin(), w.end()); return v; }
int c{0}; // solution counter
auto st{steady_clock::now()}, st0{st}, tmp{st}; // for determining elasped time
// formats elasped time
string dFmt(duration<double> et, int digs) {
string res{""}; double dt{et.count()};
if (dt > 60.0) { int m = (int)(dt / 60.0); dt -= m * 60.0; res = to_string(m) + "m"; }
res += to_string(dt); return res.substr(0, digs - 1) + 's';
}
// combines list of square differences with list of square sums, reports compatible results
VS dump(int nd, VU lo, VU hi) {
VS res {};
for (auto l : lo) for (auto h : hi) {
auto r { (h - l) >> 1 }, z { h - r };
if (izRev(nd, r, z)) {
char buf[99]; sprintf(buf, "%20llu %11lu %10lu", z, (long long)sqrt(h), (long long)sqrt(l));
res.push_back(buf); } } return res;
}
// reports one block of digits
void doOne(int n, nLH L, nLH H) {
VS lines = dump(n, L.even, H.even); lines += dump(n, L.odd , H.odd); sort(lines.begin(), lines.end());
duration<double> tet = (tmp = steady_clock::now()) - st; int ls = lines.size();
if (ls-- > 0)
for (int i{0}; i <= ls; ++i)
printf("%3d %s%s", ++c, lines[i].c_str(), i == ls ? "" : "\n");
else printf("%s", string(47, ' ').c_str());
printf(" %2d: %s %s\n", n, dFmt(tmp - st0, 8).c_str(), dFmt(tet, 8).c_str()); st0 = tmp;
}
void Rare(int n) { doOne(n, makeL(n), makeH(n)); }
int main(int argc, char *argv[]) {
int max{argc > 1 ? stoi(argv[1]) : 19}; if (max < 2) max = 2; if (max > 19 ) max = 19;
printf("%4s %19s %11s %10s %5s %11s %9s\n", "nth", "forward", "rt.sum", "rt.diff", "digs", "block.et", "total.et");
for (int nd{2}; nd <= max; ++nd) Rare(nd);
}
- Output:
Processor Intel(R) Core(TM) i5-1035G1 CPU @ 1.00GHz
nth forward rt.sum rt.diff digs block.et total.et
1 65 11 3 2: 0.00003s 0.00003s
3: 0.00002s 0.00006s
4: 0.00001s 0.00008s
5: 0.00003s 0.00011s
2 621770 836 738 6: 0.00007s 0.00018s
7: 0.00035s 0.00054s
8: 0.00110s 0.00164s
3 281089082 23708 330 9: 0.00657s 0.00821s
4 2022652202 63602 300
5 2042832002 63602 6360 10: 0.02246s 0.03068s
11: 0.11082s 0.14151s
6 868591084757 1275175 333333
7 872546974178 1320616 32670
8 872568754178 1320616 33330 12: 0.00868s 0.15019s
9 6979302951885 3586209 1047717 13: 0.03915s 0.18935s
10 20313693904202 6368252 269730
11 20313839704202 6368252 270270
12 20331657922202 6368252 329670
13 20331875722202 6368252 330330
14 20333875702202 6368252 336330
15 40313893704200 6368252 6330336
16 40351893720200 6368252 6336336 14: 0.11688s 0.30624s
17 200142385731002 20006998 69300
18 204238494066002 20122102 1891560
19 221462345754122 21045662 69300
20 244062891224042 22011022 1908060
21 245518996076442 22140228 921030
22 248359494187442 22206778 1891560
23 403058392434500 20211202 19940514
24 441054594034340 22011022 19940514
25 816984566129618 40421606 250800 15: 0.69490s 1.00114s
26 2078311262161202 64030648 7529850
27 2133786945766212 65272218 2666730
28 2135568943984212 65272218 3267330
29 2135764587964212 65272218 3326670
30 2135786765764212 65272218 3333330
31 4135786945764210 65272218 63333336
32 6157577986646405 105849161 33333333
33 6889765708183410 83866464 82133718
34 8052956026592517 123312255 29999997
35 8052956206592517 123312255 30000003
36 8191154686620818 127950856 3299670
37 8191156864620818 127950856 3300330
38 8191376864400818 127950856 3366330
39 8650327689541457 127246955 33299667
40 8650349867341457 127246955 33300333 16: 2.18232s 3.18347s
41 22542040692914522 212329862 333300
42 67725910561765640 269040196 251135808
43 86965750494756968 417050956 33000 17: 13.1765s 16.3599s
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: 41.6983s 58.0583s
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: 5m1.342s 5m59.40s
20+ digits
// 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);
}
- Output:
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 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 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
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.
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));
}
}
- 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
EasyLang
fastfunc next n .
while 1 = 1
n += 1
h = n
nrev = 0
while h > 0
nrev = nrev * 10 + h mod 10
h = h div 10
.
if sqrt (n + nrev) mod 1 = 0
if n - nrev >= 1 and sqrt (n - nrev) mod 1 = 0
return n
.
.
.
.
for cnt to 5
n = next n
print n
.
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
// 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)
43 down
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
- 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
FreeBASIC
Made some changes and added a simple test to speed things up. Results in about 1 minute.
Function revn(n As ULongInt, nd As ULongInt) As ULongInt
Dim As ULongInt r
For i As UInteger = 1 To nd
r = r * 10 + n Mod 10
n = n \ 10
Next i
Return r
End Function
Dim As UInteger nd = 2, count, lim = 90, n = 20
Do
n += 1
Dim As ULongInt r = revn(n,nd)
If r < n Then
Dim As ULongInt s = n + r, d = n - r
If nd And 1 Then
If d Mod 1089 <> 0 Then GoTo jump
Else
If s Mod 121 <> 0 Then GoTo jump
End If
If Frac(Sqr(s)) = 0 And Frac(Sqr(d)) = 0 Then
count += 1
Print count; ": "; n
If count = 5 Then Exit Do : End If
End If
End If
jump:
If n = lim Then
lim = lim * 10
nd += 1
n = (lim \ 9) * 2
End If
Loop
Print
Print "Done"
Sleep
- Output:
1: 65 2: 621770 3: 281089082 4: 2022652202 5: 2042832002
Go
Traditional
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.
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: {{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] = [][]int8{{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))
}
}
- 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
Quicker
Twice as quick as the first Go version though, despite being a faithful translation, a little slower than the C# and VB.NET versions.
package main
import (
"fmt"
"math"
"sort"
"time"
)
type llst = [][]int
var (
d []int // permutation working slice
drar [19]int // digital root lookup array
dac []int // running digital root slice
p [20]int64 // powers of 10
ac []int64 // accumulator slice
pp []int64 // coefficient slice that combines with digits of working slice
sr []int64 // temporary list of squares used for building
)
var (
odd = false // flag for odd number of digits
sum int64 // calculated sum of terms (square candidate)
rt int64 // root of sum
cn = 0 // solution counter
nd = 2 // number of digits
nd1 = nd - 1 // 'nd' helper
ln int // previous value of 'n' (in recurse())
dl int // length of 'd' slice
)
var (
tlo = []int{0, 1, 4, 5, 6} // primary differences starting point
all = seq(-9, 9, 1) // all possible differences
odl = seq(-9, 9, 2) // odd possible differences
evl = seq(-8, 8, 2) // even possible differences
thi = []int{4, 5, 6, 9, 10, 11, 14, 15, 16} // primary sums starting point
alh = seq(0, 18, 1) // all possible sums
odh = seq(1, 17, 2) // odd possible sums
evh = seq(0, 18, 2) // even possible sums
ten = seq(0, 9, 1) // used for odd number of digits
z = seq(0, 0, 1) // no difference, avoids generating a bunch of negative square candidates
t7 = []int{-3, 7} // shortcut for low 5
nin = []int{9} // shortcut for hi 10
tn = []int{10} // shortcut for hi 0 (unused, unneeded)
t12 = []int{2, 12} // shortcut for hi 5
o11 = []int{1, 11} // shortcut for hi 15
pos = []int{0, 1, 4, 5, 6, 9} // shortcut for 2nd lo 0
)
var (
lul = llst{z, odl, nil, nil, evl, t7, odl} // shortcut lookup lo primary
luh = llst{tn, evh, nil, nil, evh, t12, odh, nil, nil, evh, nin, odh, nil, nil,
odh, o11, evh} // shortcut lookup hi primary
l2l = llst{pos, nil, nil, nil, all, nil, all} // shortcut lookup lo secondary
l2h = llst{nil, nil, nil, nil, alh, nil, alh, nil, nil, nil, alh, nil, nil, nil,
alh, nil, alh} // shortcut lookup hi secondary
lu, l2 llst // ditto
chTen = llst{{0, 2, 5, 8, 9}, {0, 3, 4, 6, 9}, {1, 4, 7, 8}, {2, 3, 5, 8},
{0, 3, 6, 7, 9}, {1, 2, 4, 7}, {2, 5, 6, 8}, {0, 1, 3, 6, 9}, {1, 4, 5, 7}}
chAH = llst{{0, 2, 5, 8, 9, 11, 14, 17, 18}, {0, 3, 4, 6, 9, 12, 13, 15, 18}, {1, 4, 7, 8, 10, 13, 16, 17},
{2, 3, 5, 8, 11, 12, 14, 17}, {0, 3, 6, 7, 9, 12, 15, 16, 18}, {1, 2, 4, 7, 10, 11, 13, 16},
{2, 5, 6, 8, 11, 14, 15, 17}, {0, 1, 3, 6, 9, 10, 12, 15, 18}, {1, 4, 5, 7, 10, 13, 14, 16}}
)
// Returns a sequence of integers.
func seq(f, t, s int) []int {
r := make([]int, (t-f)/s+1)
for i := 0; i < len(r); i, f = i+1, f+s {
r[i] = f
}
return r
}
// Returns Integer Square Root.
func isr(s int64) int64 {
return int64(math.Sqrt(float64(s)))
}
// Recursively determines whether 'r' is the reverse of 'f'.
func isRev(nd int, f, r int64) bool {
nd--
if f/p[nd] != r%10 {
return false
}
if nd < 1 {
return true
}
return isRev(nd, f%p[nd], r/10)
}
// Recursive function to evaluate the permutations, no shortcuts.
func recurseLE5(lst llst, lv int) {
if lv == dl { // check if on last stage of permutation
sum = ac[lv-1]
if sum > 0 {
rt = int64(math.Sqrt(float64(sum)))
if rt*rt == sum { // test accumulated sum, append to result if square
sr = append(sr, sum)
}
}
} else {
for _, n := range lst[lv] { // set up next permutation
d[lv] = n
if lv == 0 {
ac[0] = pp[0] * int64(n)
} else {
ac[lv] = ac[lv-1] + pp[lv]*int64(n) // update accumulated sum
}
recurseLE5(lst, lv+1) // recursively call next level
}
}
}
// Recursive function to evaluate the hi permutations, shortcuts added to avoid generating many non-squares, digital root calc added.
func recursehi(lst llst, lv int) {
lv1 := lv - 1
if lv == dl { // check if on last stage of permutation
sum = ac[lv1]
if (0x202021202030213 & (1 << (int(sum) & 63))) != 0 { // test accumulated sum, append to result if square
rt = int64(math.Sqrt(float64(sum)))
if rt*rt == sum {
sr = append(sr, sum)
}
}
} else {
for _, n := range lst[lv] { // set up next permutation
d[lv] = n
if lv == 0 {
ac[0] = pp[0] * int64(n)
dac[0] = drar[n] // update accumulated sum and running dr
} else {
ac[lv] = ac[lv1] + pp[lv]*int64(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: // primary level: set shortcuts for secondary level
ln = n
lst[1] = lu[ln]
lst[2] = l2[n]
case 1: // secondary level: set shortcuts for tertiary level
switch ln { // for sums
case 5, 15:
if n < 10 {
lst[2] = evh
} else {
lst[2] = odh
}
case 9:
if ((n >> 1) & 1) == 0 {
lst[2] = evh
} else {
lst[2] = odh
}
case 11:
if ((n >> 1) & 1) == 1 {
lst[2] = evh
} else {
lst[2] = odh
}
}
}
if lv == dl-2 {
// reduce last round according to dr calc
if odd {
lst[dl-1] = chTen[dac[dl-2]]
} else {
lst[dl-1] = chAH[dac[dl-2]]
}
}
recursehi(lst, lv+1) // recursively call next level
}
}
}
// Recursive function to evaluate the lo permutations, shortcuts added to avoid
// generating many non-squares.
func recurselo(lst llst, lv int) {
lv1 := lv - 1
if lv == dl { // check if on last stage of permutation
sum = ac[lv1]
if sum > 0 {
rt = int64(math.Sqrt(float64(sum)))
if rt*rt == sum { // test accumulated sum, append to result if square
sr = append(sr, sum)
}
}
} else {
for _, n := range lst[lv] { // set up next permutation
d[lv] = n
if lv == 0 {
ac[0] = pp[0] * int64(n)
} else {
ac[lv] = ac[lv1] + pp[lv]*int64(n) // update accumulated sum
}
switch lv { // shortcuts to be performed on designated levels
case 0: // primary level: set shortcuts for secondary level
ln = n
lst[1] = lu[ln]
lst[2] = l2[n]
case 1: // secondary level: set shortcuts for tertiary level
switch ln { // for difs
case 1:
if (((n + 9) >> 1) & 1) == 0 {
lst[2] = evl
} else {
lst[2] = odl
}
case 5:
if n < 0 {
lst[2] = evl
} else {
lst[2] = odl
}
}
}
recurselo(lst, lv+1) // Recursively call next level
}
}
}
// Produces a list of candidate square numbers.
func listEm(lst, plu, pl2 llst) []int64 {
dl = len(lst)
d = make([]int, dl)
sr = sr[:0]
lu = plu
l2 = pl2
ac = make([]int64, dl)
dac = make([]int, dl) // init support vars
pp = make([]int64, dl)
for i, j := 0, nd1; i < dl; i, j = i+1, j-1 {
// build coefficients array
if len(lst[0]) > 6 {
pp[i] = p[j] + p[i]
} else {
pp[i] = p[j] - p[i]
}
}
// call appropriate recursive function
if nd <= 5 {
recurseLE5(lst, 0)
} else if len(lst[0]) > 8 {
recursehi(lst, 0)
} else {
recurselo(lst, 0)
}
return sr
}
// Reveals whether combining two lists of squares can produce a Rare number.
func reveal(lo, hi []int64) {
var s []string // create temp list of results
for _, l := range lo {
for _, h := range hi {
r := (h - l) >> 1
f := h - r // generate all possible fwd & rev candidates from lists
if isRev(nd, f, r) { // test and append sucesses to temp list
s = append(s, fmt.Sprintf("%20d %11d %10d ", f, isr(h), isr(l)))
}
}
}
sort.Strings(s)
if len(s) > 0 {
for _, t := range s { // if there are any, output sorted results
cn++
tt := ""
if t != s[len(s)-1] {
tt = "\n"
}
fmt.Printf("%2d %s%s", cn, t, tt)
}
} else {
fmt.Printf("%48s", "")
}
}
/* Unsigned variables and functions for nd == 19 */
var (
usum uint64 // unsigned calculated sum of terms (square candidate)
urt uint64 // unsigned root of sum
acu []uint64 // unsigned accumulator slice
ppu []uint64 // unsigned long coefficient slice that combines with digits of working slice
sru []uint64 // unsigned temporary list of squares used for building
)
// Returns Unsigned Integer Square Root.
func isrU(s uint64) uint64 {
return uint64(math.Sqrt(float64(s)))
}
// Recursively determines whether 'r' is the reverse of 'f'.
func isRevU(nd int, f, r uint64) bool {
nd--
if f/uint64(p[nd]) != r%10 {
return false
}
if nd < 1 {
return true
}
return isRevU(nd, f%uint64(p[nd]), r/10)
}
// Recursive function to evaluate the unsigned hi permutations, shortcuts added to avoid
// generating many non-squares, digital root calc added.
func recurseUhi(lst llst, lv int) {
lv1 := lv - 1
if lv == dl { // check if on last stage of permutation
usum = acu[lv1]
if (0x202021202030213 & (1 << (int(usum) & 63))) != 0 { // test accumulated sum, append to result if square
urt = uint64(math.Sqrt(float64(usum)))
if urt*urt == usum {
sru = append(sru, usum)
}
}
} else {
for _, n := range lst[lv] { // set up next permutation
d[lv] = n
if lv == 0 {
acu[0] = ppu[0] * uint64(n)
dac[0] = drar[n] // update accumulated sum and running dr
} else {
if n >= 0 {
acu[lv] = acu[lv1] + ppu[lv]*uint64(n)
} else {
acu[lv] = acu[lv1] - ppu[lv]*uint64(-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: // primary level: set shortcuts for secondary level
ln = n
lst[1] = lu[ln]
lst[2] = l2[n]
case 1: // secondary level: set shortcuts for tertiary level
switch ln { // for sums
case 5, 15:
if n < 10 {
lst[2] = evh
} else {
lst[2] = odh
}
case 9:
if ((n >> 1) & 1) == 0 {
lst[2] = evh
} else {
lst[2] = odh
}
case 11:
if ((n >> 1) & 1) == 1 {
lst[2] = evh
} else {
lst[2] = odh
}
}
}
if lv == dl-2 {
// reduce last round according to dr calc
if odd {
lst[dl-1] = chTen[dac[dl-2]]
} else {
lst[dl-1] = chAH[dac[dl-2]]
}
}
recurseUhi(lst, lv+1) // recursively call next level
}
}
}
// Recursive function to evaluate the unsigned lo permutations, shortcuts added to avoid
// generating many non-squares.
func recurseUlo(lst llst, lv int) {
lv1 := lv - 1
if lv == dl { // check if on last stage of permutation
usum = acu[lv1]
if usum > 0 {
urt = uint64(math.Sqrt(float64(usum)))
if urt*urt == usum { // test accumulated sum, append to result if square
sru = append(sru, usum)
}
}
} else {
for _, n := range lst[lv] { // set up next permutation
d[lv] = n
if lv == 0 {
acu[0] = ppu[0] * uint64(n)
} else {
if n >= 0 {
acu[lv] = acu[lv1] + ppu[lv]*uint64(n) // update accumulated sum
} else {
acu[lv] = acu[lv1] - ppu[lv]*uint64(-n)
}
}
switch lv { // shortcuts to be performed on designated levels
case 0: // primary level: set shortcuts for secondary level
ln = n
lst[1] = lu[ln]
lst[2] = l2[n]
case 1: // secondary level: set shortcuts for tertiary level
switch ln { // for difs
case 1:
if (((n + 9) >> 1) & 1) == 0 {
lst[2] = evl
} else {
lst[2] = odl
}
case 5:
if n < 0 {
lst[2] = evl
} else {
lst[2] = odl
}
}
}
recurseUlo(lst, lv+1) // Recursively call next level
}
}
}
// Produces a list of candidate square numbers.
func listEmU(lst, plu, pl2 llst) []uint64 {
dl = len(lst)
d = make([]int, dl)
sru = sru[:0]
lu = plu
l2 = pl2
acu = make([]uint64, dl)
dac = make([]int, dl) // init support vars
ppu = make([]uint64, dl)
for i, j := 0, nd1; i < dl; i, j = i+1, j-1 {
// build coefficients array
if len(lst[0]) > 6 {
ppu[i] = uint64(p[j] + p[i])
} else {
ppu[i] = uint64(p[j] - p[i])
}
}
// call appropriate recursive functin on
if len(lst[0]) > 8 {
recurseUhi(lst, 0)
} else {
recurseUlo(lst, 0)
}
return sru
}
// Reveals whether combining two lists of unsigned squares can produce a Rare number.
func revealU(lo, hi []uint64) {
var s []string // create temp list of results
for _, l := range lo {
for _, h := range hi {
r := (h - l) >> 1
f := h - r // generate all possible fwd & rev candidates from lists
if isRevU(nd, f, r) { // test and append sucesses to temp list
s = append(s, fmt.Sprintf("%20d %11d %10d ", f, isrU(h), isrU(l)))
}
}
}
sort.Strings(s)
if len(s) > 0 {
for _, t := range s { // if there are any, output sorted results
cn++
tt := ""
if t != s[len(s)-1] {
tt = "\n"
}
fmt.Printf("%2d %s%s", cn, t, tt)
}
} else {
fmt.Printf("%48s", "")
}
}
var (
bStart time.Time // block start time
tStart time.Time // total start time
)
// Formats time in form hh:mm:ss.fff (i.e. millisecond precision).
func formatTime(d time.Duration) string {
f := d.Milliseconds()
s := f / 1000
f %= 1000
m := s / 60
s %= 60
h := m / 60
m %= 60
return fmt.Sprintf("%02d:%02d:%02d.%03d", h, m, s, f)
}
func main() {
start := time.Now()
fmt.Printf("%3s%20s %11s %10s %3s %11s %11s\n", "nth", "forward", "rt.sum", "rt.dif", "digs", "block time", "total time")
p[0] = 1
for i, j := 0, 1; j < len(p); j++ {
p[j] = p[i] * 10 // create powers of 10 array
i = j
}
for i := 0; i < len(drar); i++ {
drar[i] = (i << 1) % 9 // create digital root array
}
bStart = time.Now()
tStart = bStart
lls := llst{tlo}
hls := llst{thi}
for nd <= 18 { // loop through all numbers of digits
if nd > 2 {
if odd {
hls = append(hls, ten)
} else {
lls = append(lls, all)
hls[len(hls)-1] = alh
}
} // build permutations list
tmp1 := listEm(lls, lul, l2l)
tmp2 := make([]int64, len(tmp1))
copy(tmp2, tmp1)
reveal(tmp2, listEm(hls, luh, l2h)) // reveal results
if !odd && nd > 5 {
hls[len(hls)-1] = alh // restore last element of hls, so that dr shortcut doesn't mess up next nd
}
bTime := formatTime(time.Since(bStart))
tTime := formatTime(time.Since(tStart))
fmt.Printf("%2d: %s %s\n", nd, bTime, tTime)
bStart = time.Now() // restart block timing
nd1 = nd
nd++
odd = !odd
}
// nd == 19
hls = append(hls, ten)
tmp3 := listEmU(lls, lul, l2l)
tmp4 := make([]uint64, len(tmp3))
copy(tmp4, tmp3)
revealU(tmp4, listEmU(hls, luh, l2h)) // reveal unsigned results
fbTime := formatTime(time.Since(bStart))
ftTime := formatTime(time.Since(tStart))
fmt.Printf("%2d: %s %s\n", nd, fbTime, ftTime)
}
- Output:
nth forward rt.sum rt.dif digs block time total time
1 65 11 3 2: 00:00:00.000 00:00:00.000
3: 00:00:00.000 00:00:00.000
4: 00:00:00.000 00:00:00.000
5: 00:00:00.000 00:00:00.000
2 621770 836 738 6: 00:00:00.000 00:00:00.000
7: 00:00:00.000 00:00:00.000
8: 00:00:00.001 00:00:00.001
3 281089082 23708 330 9: 00:00:00.006 00:00:00.008
4 2022652202 63602 300
5 2042832002 63602 6360 10: 00:00:00.015 00:00:00.023
11: 00:00:00.036 00:00:00.060
6 868591084757 1275175 333333
7 872546974178 1320616 32670
8 872568754178 1320616 33330 12: 00:00:00.057 00:00:00.117
9 6979302951885 3586209 1047717 13: 00:00:00.376 00:00:00.494
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.981 00:00:01.475
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.042 00:00:08.517
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.521 00:00:27.039
41 22542040692914522 212329862 333300
42 67725910561765640 269040196 251135808
43 86965750494756968 417050956 33000 17: 00:02:13.481 00:02:40.521
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:06:17.288 00:08:57.810
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:45:42.006 00:54:39.816
Turbo
A smallish turbo anyway :)
The following, which is a translation of the C++ 10 to 19 digits program and includes improvements suggested by Enter your username (see Talk page), completes in around 15.25 minutes which is more than 3.5 times faster than the 'quicker' version above.
Curiously, the C++ program itself when compiled with g++ 7.5.0 (-std=c++17 -O3) on the same machine (and incorporating the same improvements) completes in just over 21 minutes though I understand that when using other C++ compilers (clang, mingw) it's much faster than this.
package main
import (
"fmt"
"math"
"sort"
"time"
)
type (
z1 func() z2
z2 struct {
value int64
hasValue bool
}
)
var pow10 [19]int64
func init() {
pow10[0] = 1
for i := 1; i < 19; i++ {
pow10[i] = 10 * pow10[i-1]
}
}
func izRev(n int, i, g uint64) bool {
if i/uint64(pow10[n-1]) != g%10 {
return false
}
if n < 2 {
return true
}
return izRev(n-1, i%uint64(pow10[n-1]), g/10)
}
func fG(n z1, start, end, reset int, step int64, l *int64) z1 {
i, g, e := step*int64(start), step*int64(end), step*int64(reset)
return func() z2 {
for i < g {
*l += step
i += step
return z2{*l, true}
}
i = e
*l -= (g - e)
return n()
}
}
type nLH struct{ even, odd []uint64 }
type zp struct {
n z1
g [][2]int64
}
func newNLH(e zp) nLH {
var even, odd []uint64
n, g := e.n, e.g
for i := n(); i.hasValue; i = n() {
for _, p := range g {
ng, gg := p[0], p[1]
if (ng > 0) || (i.value > 0) {
w := uint64(ng*pow10[4] + gg + i.value)
ws := uint64(math.Sqrt(float64(w)))
if ws*ws == w {
if w%2 == 0 {
even = append(even, w)
} else {
odd = append(odd, w)
}
}
}
}
}
return nLH{even, odd}
}
func makeL(n int) zp {
g := make([]z1, n/2-3)
g[0] = func() z2 { return z2{} }
for i := 1; i < n/2-3; i++ {
s := -9
if i == n/2-4 {
s = -10
}
l := pow10[n-i-4] - pow10[i+3]
acc += l * int64(s)
g[i] = fG(g[i-1], s, 9, -9, l, &acc)
}
var g0, g1, g2, g3 int64
l0, l1, l2, l3 := pow10[n-5], pow10[n-6], pow10[n-7], pow10[n-8]
f := func() [][2]int64 {
var w [][2]int64
for g0 < 7 {
nn := g3*l3 + g2*l2 + g1*l1 + g0*l0
gg := -1000*g3 - 100*g2 - 10*g1 - g0
if g3 < 9 {
g3++
} else {
g3 = -9
if g2 < 9 {
g2++
} else {
g2 = -9
if g1 < 9 {
g1++
} else {
g1 = -9
if g0 == 1 {
g0 = 3
}
g0++
}
}
}
if bs[(pow10[10]+gg)%10000] {
w = append(w, [2]int64{nn, gg})
}
}
return w
}
return zp{g[n/2-4], f()}
}
func makeH(n int) zp {
acc = -(pow10[n/2] + pow10[(n-1)/2])
g := make([]z1, (n+1)/2-3)
g[0] = func() z2 { return z2{} }
for i := 1; i < n/2-3; i++ {
j := 0
if i == (n+1)/2-3 {
j = -1
}
g[i] = fG(g[i-1], j, 18, 0, pow10[n-i-4]+pow10[i+3], &acc)
if n%2 == 1 {
g[(n+1)/2-4] = fG(g[n/2-4], -1, 9, 0, 2*pow10[n/2], &acc)
}
}
g0 := int64(4)
var g1, g2, g3 int64
l0, l1, l2, l3 := pow10[n-5], pow10[n-6], pow10[n-7], pow10[n-8]
f := func() [][2]int64 {
var w [][2]int64
for g0 < 17 {
nn := g3*l3 + g2*l2 + g1*l1 + g0*l0
gg := 1000*g3 + 100*g2 + 10*g1 + g0
if g3 < 18 {
g3++
} else {
g3 = 0
if g2 < 18 {
g2++
} else {
g2 = 0
if g1 < 18 {
g1++
} else {
g1 = 0
if g0 == 6 || g0 == 9 {
g0 += 3
}
g0++
}
}
}
if bs[gg%10000] {
w = append(w, [2]int64{nn, gg})
}
}
return w
}
return zp{g[(n+1)/2-4], f()}
}
var (
acc int64
bs = make([]bool, 10000)
L, H nLH
)
func rare(n int) []uint64 {
acc = 0
for g := 0; g < 10000; g++ {
bs[(g*g)%10000] = true
}
L = newNLH(makeL(n))
H = newNLH(makeH(n))
var rares []uint64
for _, l := range L.even {
for _, h := range H.even {
r := (h - l) / 2
z := h - r
if izRev(n, r, z) {
rares = append(rares, z)
}
}
}
for _, l := range L.odd {
for _, h := range H.odd {
r := (h - l) / 2
z := h - r
if izRev(n, r, z) {
rares = append(rares, z)
}
}
}
if len(rares) > 0 {
sort.Slice(rares, func(i, j int) bool {
return rares[i] < rares[j]
})
}
return rares
}
// Formats time in form hh:mm:ss.fff (i.e. millisecond precision).
func formatTime(d time.Duration) string {
f := d.Milliseconds()
s := f / 1000
f %= 1000
m := s / 60
s %= 60
h := m / 60
m %= 60
return fmt.Sprintf("%02d:%02d:%02d.%03d", h, m, s, f)
}
func commatize(n uint64) string {
s := fmt.Sprintf("%d", n)
le := len(s)
for i := le - 3; i >= 1; i -= 3 {
s = s[0:i] + "," + s[i:]
}
return s
}
func main() {
bStart := time.Now() // block time
tStart := bStart // total time
nth := 3 // i.e. count of rare numbers < 10 digits
fmt.Println("nth rare number digs block time total time")
for nd := 10; nd <= 19; nd++ {
rares := rare(nd)
if len(rares) > 0 {
for i, r := range rares {
nth++
t := ""
if i < len(rares)-1 {
t = "\n"
}
fmt.Printf("%2d %25s%s", nth, commatize(r), t)
}
} else {
fmt.Printf("%29s", "")
}
fbTime := formatTime(time.Since(bStart))
ftTime := formatTime(time.Since(tStart))
fmt.Printf(" %2d: %s %s\n", nd, fbTime, ftTime)
bStart = time.Now() // restart block timing
}
}
- Output:
nth rare number digs block time total time
4 2,022,652,202
5 2,042,832,002 10: 00:00:00.002 00:00:00.002
11: 00:00:00.008 00:00:00.011
6 868,591,084,757
7 872,546,974,178
8 872,568,754,178 12: 00:00:00.022 00:00:00.033
9 6,979,302,951,885 13: 00:00:00.097 00:00:00.131
10 20,313,693,904,202
11 20,313,839,704,202
12 20,331,657,922,202
13 20,331,875,722,202
14 20,333,875,702,202
15 40,313,893,704,200
16 40,351,893,720,200 14: 00:00:00.265 00:00:00.396
17 200,142,385,731,002
18 204,238,494,066,002
19 221,462,345,754,122
20 244,062,891,224,042
21 245,518,996,076,442
22 248,359,494,187,442
23 403,058,392,434,500
24 441,054,594,034,340
25 816,984,566,129,618 15: 00:00:01.774 00:00:02.170
26 2,078,311,262,161,202
27 2,133,786,945,766,212
28 2,135,568,943,984,212
29 2,135,764,587,964,212
30 2,135,786,765,764,212
31 4,135,786,945,764,210
32 6,157,577,986,646,405
33 6,889,765,708,183,410
34 8,052,956,026,592,517
35 8,052,956,206,592,517
36 8,191,154,686,620,818
37 8,191,156,864,620,818
38 8,191,376,864,400,818
39 8,650,327,689,541,457
40 8,650,349,867,341,457 16: 00:00:05.407 00:00:07.578
41 22,542,040,692,914,522
42 67,725,910,561,765,640
43 86,965,750,494,756,968 17: 00:00:35.401 00:00:42.979
44 225,342,456,863,243,522
45 225,342,458,663,243,522
46 225,342,478,643,243,522
47 284,684,666,566,486,482
48 284,684,868,364,486,482
49 297,128,548,234,950,692
50 297,128,722,852,950,692
51 297,148,324,656,930,692
52 297,148,546,434,930,692
53 497,168,548,234,910,690
54 619,431,353,040,136,925
55 619,631,153,042,134,925
56 631,688,638,047,992,345
57 633,288,858,025,996,145
58 633,488,632,647,994,145
59 653,488,856,225,994,125
60 811,865,096,390,477,018
61 865,721,270,017,296,468
62 871,975,098,681,469,178
63 898,907,259,301,737,498 18: 00:01:42.745 00:02:25.725
64 2,042,401,829,204,402,402
65 2,060,303,819,041,450,202
66 2,420,424,089,100,600,242
67 2,551,755,006,254,571,552
68 2,702,373,360,882,732,072
69 2,825,378,427,312,735,282
70 6,531,727,101,458,000,045
71 6,988,066,446,726,832,640
72 8,066,308,349,502,036,608
73 8,197,906,905,009,010,818
74 8,200,756,128,308,135,597
75 8,320,411,466,598,809,138 19: 00:12:48.590 00:15:14.316
J
The following function determines whether a given integer is "rare", as can be seen in the sample output.
To use it to find rare numbers, one would simply apply it to each integer seriatim, and keep the numbers where its output is true (I.@:rare i. AS_HIGH_AS_YOU_WANT_TO_CHECK); but since these numbers are, well, rare, actually doing so would take a very long time.
rare =: ( np@:] *. (nbrPs rr) ) b10
np =: -.@:(-: |.) NB. Not palindromic
nbrPs =: > *. sdPs NB. n is Bigger than R and the perfect square constraint is satisfied
sdPs =: + *.&:ps - NB. n > rr and both their sum and difference are perfect squares
ps =: 0 = 1 | %: NB. Perfect square (integral sqrt)
rr =: 10&#.@:|. NB. Do note we do reverse the digits twice (once here, once in np)
b10 =: 10&#.^:_1 NB. Base 10 digits
- Output:
NB. From OEIS
R =: 65 621770 281089082 2022652202 2042832002 868591084757 872546974178 872568754178 6979302951885 20313693904202 20313839704202 20331657922202 20331875722202 20333875702202 40313893704200
rare"0 R
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Java
import java.time.Duration;
import java.time.LocalDateTime;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.concurrent.atomic.AtomicInteger;
import java.util.concurrent.atomic.AtomicReference;
public class RareNumbers {
public interface Consumer5<A, B, C, D, E> {
void apply(A a, B b, C c, D d, E e);
}
public interface Consumer7<A, B, C, D, E, F, G> {
void apply(A a, B b, C c, D d, E e, F f, G g);
}
public interface Recursable5<A, B, C, D, E> {
void apply(A a, B b, C c, D d, E e, Recursable5<A, B, C, D, E> r);
}
public interface Recursable7<A, B, C, D, E, F, G> {
void apply(A a, B b, C c, D d, E e, F f, G g, Recursable7<A, B, C, D, E, F, G> r);
}
public static <A, B, C, D, E> Consumer5<A, B, C, D, E> recurse(Recursable5<A, B, C, D, E> r) {
return (a, b, c, d, e) -> r.apply(a, b, c, d, e, r);
}
public static <A, B, C, D, E, F, G> Consumer7<A, B, C, D, E, F, G> recurse(Recursable7<A, B, C, D, E, F, G> r) {
return (a, b, c, d, e, f, g) -> r.apply(a, b, c, d, e, f, g, r);
}
private static class Term {
long coeff;
byte ix1, ix2;
public Term(long coeff, byte ix1, byte ix2) {
this.coeff = coeff;
this.ix1 = ix1;
this.ix2 = ix2;
}
}
private static final int MAX_DIGITS = 16;
private static long toLong(List<Byte> digits, boolean reverse) {
long sum = 0;
if (reverse) {
for (int i = digits.size() - 1; i >= 0; --i) {
sum = sum * 10 + digits.get(i);
}
} else {
for (Byte digit : digits) {
sum = sum * 10 + digit;
}
}
return sum;
}
private static boolean isNotSquare(long n) {
long root = (long) Math.sqrt(n);
return root * root != n;
}
private static List<Byte> seq(byte from, byte to, byte step) {
List<Byte> res = new ArrayList<>();
for (byte i = from; i <= to; i += step) {
res.add(i);
}
return res;
}
private static String commatize(long n) {
String s = String.valueOf(n);
int le = s.length();
int i = le - 3;
while (i >= 1) {
s = s.substring(0, i) + "," + s.substring(i);
i -= 3;
}
return s;
}
public static void main(String[] args) {
final LocalDateTime startTime = LocalDateTime.now();
long pow = 1L;
System.out.println("Aggregate timings to process all numbers up to:");
// terms of (n-r) expression for number of digits from 2 to maxDigits
List<List<Term>> allTerms = new ArrayList<>();
for (int i = 0; i < MAX_DIGITS - 1; ++i) {
allTerms.add(new ArrayList<>());
}
for (int r = 2; r <= MAX_DIGITS; ++r) {
List<Term> terms = new ArrayList<>();
pow *= 10;
long pow1 = pow;
long pow2 = 1;
byte i1 = 0;
byte i2 = (byte) (r - 1);
while (i1 < i2) {
terms.add(new Term(pow1 - pow2, i1, i2));
pow1 /= 10;
pow2 *= 10;
i1++;
i2--;
}
allTerms.set(r - 2, terms);
}
// map of first minus last digits for 'n' to pairs giving this value
Map<Byte, List<List<Byte>>> fml = Map.of(
(byte) 0, List.of(List.of((byte) 2, (byte) 2), List.of((byte) 8, (byte) 8)),
(byte) 1, List.of(List.of((byte) 6, (byte) 5), List.of((byte) 8, (byte) 7)),
(byte) 4, List.of(List.of((byte) 4, (byte) 0)),
(byte) 6, List.of(List.of((byte) 6, (byte) 0), List.of((byte) 8, (byte) 2))
);
// map of other digit differences for 'n' to pairs giving this value
Map<Byte, List<List<Byte>>> dmd = new HashMap<>();
for (int i = 0; i < 100; ++i) {
List<Byte> a = List.of((byte) (i / 10), (byte) (i % 10));
int d = a.get(0) - a.get(1);
dmd.computeIfAbsent((byte) d, k -> new ArrayList<>()).add(a);
}
List<Byte> fl = List.of((byte) 0, (byte) 1, (byte) 4, (byte) 6);
List<Byte> dl = seq((byte) -9, (byte) 9, (byte) 1); // all differences
List<Byte> zl = List.of((byte) 0); // zero differences only
List<Byte> el = seq((byte) -8, (byte) 8, (byte) 2); // even differences only
List<Byte> ol = seq((byte) -9, (byte) 9, (byte) 2); // odd differences only
List<Byte> il = seq((byte) 0, (byte) 9, (byte) 1);
List<Long> rares = new ArrayList<>();
List<List<List<Byte>>> lists = new ArrayList<>();
for (int i = 0; i < 4; ++i) {
lists.add(new ArrayList<>());
}
for (int i = 0; i < fl.size(); ++i) {
List<List<Byte>> temp1 = new ArrayList<>();
List<Byte> temp2 = new ArrayList<>();
temp2.add(fl.get(i));
temp1.add(temp2);
lists.set(i, temp1);
}
final AtomicReference<List<Byte>> digits = new AtomicReference<>(new ArrayList<>());
AtomicInteger count = new AtomicInteger();
// Recursive closure to generate (n+r) candidates from (n-r) candidates
// and hence find Rare numbers with a given number of digits.
Consumer7<List<Byte>, List<Byte>, List<List<Byte>>, List<List<Byte>>, Long, Integer, Integer> fnpr = recurse((cand, di, dis, indicies, nmr, nd, level, func) -> {
if (level == dis.size()) {
digits.get().set(indicies.get(0).get(0), fml.get(cand.get(0)).get(di.get(0)).get(0));
digits.get().set(indicies.get(0).get(1), fml.get(cand.get(0)).get(di.get(0)).get(1));
int le = di.size();
if (nd % 2 == 1) {
le--;
digits.get().set(nd / 2, di.get(le));
}
for (int i = 1; i < le; ++i) {
digits.get().set(indicies.get(i).get(0), dmd.get(cand.get(i)).get(di.get(i)).get(0));
digits.get().set(indicies.get(i).get(1), dmd.get(cand.get(i)).get(di.get(i)).get(1));
}
long r = toLong(digits.get(), true);
long npr = nmr + 2 * r;
if (isNotSquare(npr)) {
return;
}
count.getAndIncrement();
System.out.printf(" R/N %2d:", count.get());
LocalDateTime checkPoint = LocalDateTime.now();
long elapsed = Duration.between(startTime, checkPoint).toMillis();
System.out.printf(" %9sms", elapsed);
long n = toLong(digits.get(), false);
System.out.printf(" (%s)\n", commatize(n));
rares.add(n);
} else {
for (Byte num : dis.get(level)) {
di.set(level, num);
func.apply(cand, di, dis, indicies, nmr, nd, level + 1, func);
}
}
});
// Recursive closure to generate (n-r) candidates with a given number of digits.
Consumer5<List<Byte>, List<List<Byte>>, List<List<Byte>>, Integer, Integer> fnmr = recurse((cand, list, indicies, nd, level, func) -> {
if (level == list.size()) {
long nmr = 0;
long nmr2 = 0;
List<Term> terms = allTerms.get(nd - 2);
for (int i = 0; i < terms.size(); ++i) {
Term t = terms.get(i);
if (cand.get(i) >= 0) {
nmr += t.coeff * cand.get(i);
} else {
nmr2 += t.coeff * -cand.get(i);
if (nmr >= nmr2) {
nmr -= nmr2;
nmr2 = 0;
} else {
nmr2 -= nmr;
nmr = 0;
}
}
}
if (nmr2 >= nmr) {
return;
}
nmr -= nmr2;
if (isNotSquare(nmr)) {
return;
}
List<List<Byte>> dis = new ArrayList<>();
dis.add(seq((byte) 0, (byte) (fml.get(cand.get(0)).size() - 1), (byte) 1));
for (int i = 1; i < cand.size(); ++i) {
dis.add(seq((byte) 0, (byte) (dmd.get(cand.get(i)).size() - 1), (byte) 1));
}
if (nd % 2 == 1) {
dis.add(il);
}
List<Byte> di = new ArrayList<>();
for (int i = 0; i < dis.size(); ++i) {
di.add((byte) 0);
}
fnpr.apply(cand, di, dis, indicies, nmr, nd, 0);
} else {
for (Byte num : list.get(level)) {
cand.set(level, num);
func.apply(cand, list, indicies, nd, level + 1, func);
}
}
});
for (int nd = 2; nd <= MAX_DIGITS; ++nd) {
digits.set(new ArrayList<>());
for (int i = 0; i < nd; ++i) {
digits.get().add((byte) 0);
}
if (nd == 4) {
lists.get(0).add(zl);
lists.get(1).add(ol);
lists.get(2).add(el);
lists.get(3).add(ol);
} else if (allTerms.get(nd - 2).size() > lists.get(0).size()) {
for (int i = 0; i < 4; ++i) {
lists.get(i).add(dl);
}
}
List<List<Byte>> indicies = new ArrayList<>();
for (Term t : allTerms.get(nd - 2)) {
indicies.add(List.of(t.ix1, t.ix2));
}
for (List<List<Byte>> list : lists) {
List<Byte> cand = new ArrayList<>();
for (int i = 0; i < list.size(); ++i) {
cand.add((byte) 0);
}
fnmr.apply(cand, list, indicies, nd, 0);
}
LocalDateTime checkPoint = LocalDateTime.now();
long elapsed = Duration.between(startTime, checkPoint).toMillis();
System.out.printf(" %2d digits: %9sms\n", nd, elapsed);
}
Collections.sort(rares);
System.out.printf("\nThe rare numbers with up to %d digits are:\n", MAX_DIGITS);
for (int i = 0; i < rares.size(); ++i) {
System.out.printf(" %2d: %25s\n", i + 1, commatize(rares.get(i)));
}
}
}
- Output:
Aggregate timings to process all numbers up to:
R/N 1: 17ms (65)
2 digits: 19ms
3 digits: 20ms
4 digits: 21ms
5 digits: 21ms
R/N 2: 29ms (621,770)
6 digits: 55ms
7 digits: 57ms
8 digits: 72ms
R/N 3: 76ms (281,089,082)
9 digits: 85ms
R/N 4: 89ms (2,022,652,202)
R/N 5: 237ms (2,042,832,002)
10 digits: 451ms
11 digits: 503ms
R/N 6: 833ms (872,546,974,178)
R/N 7: 869ms (872,568,754,178)
R/N 8: 1324ms (868,591,084,757)
12 digits: 1582ms
R/N 9: 1888ms (6,979,302,951,885)
13 digits: 2299ms
R/N 10: 6199ms (20,313,693,904,202)
R/N 11: 6272ms (20,313,839,704,202)
R/N 12: 7831ms (20,331,657,922,202)
R/N 13: 8070ms (20,331,875,722,202)
R/N 14: 8708ms (20,333,875,702,202)
R/N 15: 19681ms (40,313,893,704,200)
R/N 16: 19823ms (40,351,893,720,200)
14 digits: 21712ms
R/N 17: 21758ms (200,142,385,731,002)
R/N 18: 21991ms (221,462,345,754,122)
R/N 19: 24990ms (816,984,566,129,618)
R/N 20: 26552ms (245,518,996,076,442)
R/N 21: 26774ms (204,238,494,066,002)
R/N 22: 26846ms (248,359,494,187,442)
R/N 23: 27158ms (244,062,891,224,042)
R/N 24: 32349ms (403,058,392,434,500)
R/N 25: 32576ms (441,054,594,034,340)
15 digits: 34465ms
R/N 26: 85843ms (2,133,786,945,766,212)
R/N 27: 105944ms (2,135,568,943,984,212)
R/N 28: 109027ms (8,191,154,686,620,818)
R/N 29: 111677ms (8,191,156,864,620,818)
R/N 30: 112849ms (2,135,764,587,964,212)
R/N 31: 114572ms (2,135,786,765,764,212)
R/N 32: 118622ms (8,191,376,864,400,818)
R/N 33: 132237ms (2,078,311,262,161,202)
R/N 34: 164708ms (8,052,956,026,592,517)
R/N 35: 169421ms (8,052,956,206,592,517)
R/N 36: 203280ms (8,650,327,689,541,457)
R/N 37: 205555ms (8,650,349,867,341,457)
R/N 38: 207237ms (6,157,577,986,646,405)
R/N 39: 246082ms (4,135,786,945,764,210)
R/N 40: 275691ms (6,889,765,708,183,410)
16 digits: 278088ms
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
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()
- 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
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)))
}
}
- 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
Lambdatalk
Just coding the definition of a rare number without any optimization.
{def lt_israre
{lambda {:n}
{let { {:n :n}
{:r {W.reverse :n}}
} {if {and {> :n :r}
{isInt {sqrt {+ :n :r}}}
{isInt {sqrt {- :n :r}}}}
then :n
else}}}}
-> lt_israre
{S.map lt_israre {S.serie 1 700000}}
-> 65 621770 // computed in 7650ms
Testing:
{S.map lt_israre {S.serie 1 280000000}}
-> ... crushes Firefox working in my small iPad Pro.
And so I ask javascript some help:
LAMBDATALK.DICT["js_israres"] = function() {
var args = arguments[0].trim().split(" "),
i0 = Number( args[0] ),
i1 = Number( args[1] ),
a = [];
var israre = function(n) {
var r = Number( n.toString().split("").reverse().join("") );
return (n > r) && (Number.isInteger(Math.sqrt(n+r)))
&& (Number.isInteger(Math.sqrt(n-r)))
};
for (var i=i0; i < i1; i++)
if (israre(i)) a.push(i);
return a
};
Testing:
{js_israres 1 2050000000}
-> [65,621770,281089082,2022652202,2042832002]]
// computed in 784307ms ~ 13 minutes
Too slow to try to go further.
langur
not optimized
It could look something like the following (ignoring whatever optimizations the other examples are using), if it was fast enough. I did not have the time/processor to test finding the first 5. The israre() function appears to return the right answer, tested with individual numbers.
val perfectsquare = fn n: (n ^/ 2) div 1
val israre = fn(n) {
val r = reverse(n)
if n == r: return false
val sum = n + r
val diff = n - r
diff > 0 and perfectsquare(sum) and perfectsquare(diff)
}
val findfirst = fn(mx) {
for[=[]] i = 0; len(_for) < mx; i += 1 {
if israre(i) {
_for ~= [i]
}
}
}
writeln findfirst(5)Lua
do -- find the first five rare numbers
local function revn ( na )
local n, r = na, 0
while n > 0 do
r = r * 10 r = r + ( n % 10 )
n = math.floor( n / 10 )
end
return r
end -- revn
local nd, count, lim, n = 2, 0, 90, 20
local oddNd = nd % 2 == 1
while true do
n = n + 1
local r = revn( n )
if r < n then
local s, d = n + r, n - r
local tryThis = false
if oddNd
then tryThis = d % 1089 == 0
else tryThis = s % 121 == 0
end
if tryThis then
local rootS = math.sqrt( s )
if rootS == math.floor(rootS )
then
local rootD = math.sqrt( d )
if rootD == math.floor( rootD )
then
count = count + 1
io.write( count, ": ", n, "\n" )
if count >= 5 then os.exit() end
end
end
end
if n == lim
then
lim = lim * 10
nd = nd + 1
oddNd = not oddNd
n = math.floor( lim / 9 ) * 2
end
end
end
end
- Output:
1: 65 2: 621770 3: 281089082 4: 2022652202 5: 2042832002
Mathematica /Wolfram Language
c = Compile[{{k, _Integer}},
Module[{out = {0}, start = 0, stop = 0, rlist = {0}, r = 0,
sum = 0.0, diff = 0.0, imax = 8, step = 10},
Do[
If[j == k, imax = 2, imax = 8];
Do[
If[i == 2,
start = i 10^j + 2;
stop = (i + 1) 10^j - 1;
step = 10;
,
start = i 10^j;
stop = (i + 1) 10^j - 1;
step = 1;
];
Do[
rlist = IntegerDigits[n];
r = 0;
Do[
r += rlist[[ri]] 10^(ri - 1)
,
{ri, 1, Length[rlist]}
];
If[r != n,
sum = n + r;
sum = Sqrt[sum];
If[Floor[sum] == sum,
diff = n - r;
If[diff > 0,
diff = Sqrt[diff];
If[Floor[diff] == diff,
AppendTo[out, n]
]
]
]
]
,
{n, start, stop, step}
]
,
{i, 2, imax, 2}
]
,
{j, 0, k}
];
out
]
];
Rest[c[9]] (*takes about 310 sec*)
- Output:
{65, 621770, 281089082, 2022652202, 2042832002}
Nim
Translation of the second Go version, limited to 18 digits.
import algorithm, math, strformat, times
type Llst = seq[seq[int]]
const
# Powers of 10.
P = block:
var p: array[19, int64]
p[0] = 1i64
for i in 1..18: p[i] = 10 * p[i - 1]
p
# Digital root lookup array.
Drar = block:
var drar: array[19, int]
for i in 0..18: drar[i] = i shl 1 mod 9
drar
var
d: seq[int] # permutation working slice
dac: seq[int] # running digital root slice
ac: seq[int64] # accumulator slice
pp: seq[int64] # coefficient slice that combines with digits of working slice
sr: seq[int64] # temporary list of squares used for building
var
odd = false # flag for odd number of digits
sum: int64 # calculated sum of terms (square candidate)
cn = 0 # solution counter
nd = 2 # number of digits
nd1 = nd - 1 # 'nd' helper
ln: int # previous value of 'n' (in recurse())
dl: int # length of 'd' slice
func newIntSeq(f, t, s: int): seq[int] =
## Return a sequence of integers.
result = newSeq[int]((t - f) div s + 1)
var f = f
for i in 0..result.high:
result[i] = f
inc f, s
const
Tlo = @[0, 1, 4, 5, 6] # primary differences starting point
All = newIntSeq(-9, 9, 1) # all possible differences
Odl = newIntSeq(-9, 9, 2) # odd possible differences
Evl = newIntSeq(-8, 8, 2) # even possible differences
Thi = @[4, 5, 6, 9, 10, 11, 14, 15, 16] # primary sums starting point
Alh = newIntSeq(0, 18, 1) # all possible sums
Odh = newIntSeq(1, 17, 2) # odd possible sums
Evh = newIntSeq(0, 18, 2) # even possible sums
Ten = newIntSeq(0, 9, 1) # used for odd number of digits
Z = newIntSeq(0, 0, 1) # no difference, avoids generating a bunch of negative square candidates
T7 = @[-3, 7] # shortcut for low 5
Nin = @[9] # shortcut for hi 10
Tn = @[10] # shortcut for hi 0 (unused, unneeded)
T12 = @[2, 12] # shortcut for hi 5
O11 = @[1, 11] # shortcut for hi 15
Pos = @[0, 1, 4, 5, 6, 9] # shortcut for 2nd lo 0
var
lul: Llst = @[Z, Odl, @[], @[], Evl, T7, Odl] # shortcut lookup lo primary
luh: Llst = @[Tn, Evh, @[], @[], Evh, T12, Odh, @[], @[],
Evh, Nin, Odh, @[], @[], Odh, O11, Evh] # shortcut lookup hi primary
l2l: Llst = @[Pos, @[], @[], @[], All, @[], All] # shortcut lookup lo secondary
l2h: Llst = @[@[], @[], @[], @[], Alh, @[], Alh, @[], @[],
@[], Alh, @[], @[], @[], Alh, @[], Alh] # shortcut lookup hi secondary
chTen: Llst = @[@[0, 2, 5, 8, 9], @[0, 3, 4, 6, 9], @[1, 4, 7, 8],
@[2, 3, 5, 8], @[0, 3, 6, 7, 9], @[1, 2, 4, 7],
@[2, 5, 6, 8], @[0, 1, 3, 6, 9], @[1, 4, 5, 7]]
chAH: Llst = @[@[0, 2, 5, 8, 9, 11, 14, 17, 18], @[0, 3, 4, 6, 9, 12, 13, 15, 18],
@[1, 4, 7, 8, 10, 13, 16, 17], @[2, 3, 5, 8, 11, 12, 14, 17],
@[0, 3, 6, 7, 9, 12, 15, 16, 18], @[1, 2, 4, 7, 10, 11, 13, 16],
@[2, 5, 6, 8, 11, 14, 15, 17], @[0, 1, 3, 6, 9, 10, 12, 15, 18],
@[1, 4, 5, 7, 10, 13, 14, 16]]
var lu, l2: Llst
func isr(s: int64): int64 {.inline.} =
## Return integer square root.
int64(sqrt(float(s)))
proc isRev(nd: int; f, r: int64): bool =
## Recursively determines whether 'r' is the reverse of 'f'.
let nd = nd - 1
if f div P[nd] != r mod 10: return false
if nd < 1: return true
result = isRev(nd, f mod P[nd], r div 10)
proc recurseLE5(lst: Llst; lv: int) =
## Recursive function to evaluate the permutations, no shortcuts.
if lv == dl: # Check if on last stage of permutation.
sum = ac[lv - 1]
if sum > 0:
let rt = int64(sqrt(float(sum)))
if rt * rt == sum: sr.add sum
else:
for 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.
proc recursehi(lst: var Llst; lv: int) =
## Recursive function to evaluate the hi permutations.
## Shortcuts added to avoid generating many non-squares, digital root calc added.
let lv1 = lv - 1
if lv == dl: # Check if on last stage of permutation.
sum = ac[lv1]
if (0x202021202030213 and (1 shl (sum and 63))) != 0:
# Test accumulated sum, append to result if square.
let rt = int64(sqrt(float64(sum)))
if rt * rt == sum: sr.add sum
else:
for 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: dec dac[lv], 9
case lv # Shortcuts to be performed on designated levels.
of 0: # Primary level: set shortcuts for secondary level.
ln = n
lst[1] = lu[ln]
lst[2] = l2[n]
of 1: # Secondary level: set shortcuts for tertiary level.
case ln # For sums.
of 5, 15: lst[2] = if n < 10: Evh else: Odh
of 9: lst[2] = if (n shr 1 and 1) == 0: Evh else: Odh
of 11: lst[2] = if (n shr 1 and 1) == 1: Evh else: Odh
else: discard
else: discard
if lv == dl - 2:
# Reduce last round according to dr calc.
lst[dl - 1] = if odd: chTen[dac[dl - 2]] else: chAH[dac[dl - 2]]
recursehi(lst, lv + 1) # Recursively call next level.
proc recurselo(lst: var Llst; lv: int) =
## Recursive function to evaluate the lo permutations.
## Shortcuts added to avoid generating many non-squares.
let lv1 = lv - 1
if lv == dl: # Check if on last stage of permutation.
sum = ac[lv1]
if sum > 0:
let rt = int64(sqrt(float64(sum)))
if rt * rt == sum: sr.add sum
else:
for 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.
case lv # Shortcuts to be performed on designated levels.
of 0: # Primary level: set shortcuts for secondary level.
ln = n
lst[1] = lu[ln]