# RADD decomposition of an integer

## Introduction

The $$\RADD(n)$/extract_tex] operation is defined as the sum of $$\n + [\$$ the number whose decimal representation are the decimal digits of $$\n\$$ in reverse order $$$\$$, see A004086. After reversal, trailing zeros are lost. See Reverse-Then-Add Sequence and A056964.

Not all numbers can be represented as the result of a $$\RADD\$$ operation.

## Examples

$$\RADD(11) = 11 + 11 = 22\\ RADD(23) = 23 + 32 = 55\\ RADD(70) = 70 + 7 = 77\\ RADD(85) = 85 + 58 = 143\\ RADD(1100)= 1100 + 11 = 1111\$$

## Task

We are looking for a method that determines for a given number $$\m\$$, whether it can be represented with $$\RADD\$$, and if this is the case, then determines an argument $$\n\$$ such that $$\m = RADD(n)\$$, symbolically $$\n = RADD^{-1}(m)\$$

As a convention, the larger of the two summands shall be used as result, e.g., $$\RADD^{-1}(55) = 32, RADD^{-1}(143) = 85\$$

If more than one decomposition is possible, there is no preferred choice which one is used for the output.

Since the problem comes from the OEIS sequence A356648, only square numbers $$\s=n^2\$$ should be considered as a target.

## Challenge

Write a function or a program that takes an integer $$\n\gt0\$$ as input and returns the $$\RADD\$$ decompositions $$\RADD^{-1}(k_{i}^2)\$$ of as many as possible further distinct square numbers $$\k_{i}^2 \ge n^2\$$ as the result, such that there are no numbers $$\k'\$$ with $$\k_{i} whose squares $$\k'^2\$$ are also $$\RADD\$$-decomposable.

As $$\n\$$ increases, the obvious trivial method of looping through all candidates may become infeasible.

## Winning criterion

The program that delivers the most consecutive terms represented by their decompositions after $$\n=1101111\$$ $$\(1101111^2 = 220005934299 + 992439500022)\$$ in a running time of $$\300\$$ s wins. If only one term is found, the combined running time to scan the range up to this term and to find its decomposition counts.

## Tests

With $$\R2(n) := RADD^{-1}(n^2)\$$

    n = {2, 4, 11, 22, 25, 33, 101, 121, 141, 202, 222, 264, 303, 307, 451, 836, 1001}
R2(2) = 2, (2^2 = 2 + 2)
R2(22) = 341, (22^2 = 484 = 143 + 341)


for all numbers $$\n_i\le 1001\$$ not in the list no $$\RADD\$$ decomposition exists.

• Range $$\10^5\cdots 2\times 10^5\$$
n = {100001, 101101, 102201, 103801, 105270, 109901, 110011, 111111, 111221, 112211, 120021, 121121, 122221, 125092, 125129, 133431, 138259, 148489, 161619, 165269, 171959, 200002}


are the only numbers in this range for which $$\RADD\$$ decompositions exist.

• Some larger examples
n = 386221, n^2 = 149166660841 = 100166999840 + 48999661001, R2(386221)= 100166999840
n = 476311, n^2 = 226872168721 = 130002968690 + 96869200031, R2(476311) = 130002968690
n = 842336, n^2 = 709529936896 = 109600929995 + 599929006901, R2(842336) = 599929006901

• I think the standard practice here is to allow golfers define their own consistent value for the situation where no such decomposition exists. Commented Dec 31, 2022 at 19:31
• Definitely yes. Return -1 or an empty output would be options. The only condition is that it must be uniquely distinguishable from a result in which a decomposition exists. Commented Dec 31, 2022 at 20:41
• @AndersKaseorg Rephrasing the win criterion is not a problem as long as no one has submitted anything yet. What do you think: The program that produces the complete list of all 6-digit terms in the shortest time wins? Commented Dec 31, 2022 at 20:48
• Or: Returns most terms after n=100001 in a specified time, but then the problem of comparability and a running environment arises. Commented Dec 31, 2022 at 20:54
• There is no way around the problem of comparability across environments. Different computers will run a program at different speeds, even relative to other programs. You’re going to need to benchmark all the answers on the same hardware, or find someone who will do so. The “most terms” criterion is better for two reasons: it avoids the fixed output problem I described above, and it reduces the risk that you end up with an unscorable question with solutions that run so fast that they’re within the margin of error of any measurement. Commented Dec 31, 2022 at 21:40

# C++11 (gcc), O(log(n))

This algorithm can check, if a number n has a RADD-decomp and calculate the decomp in O(D)=O(log(n)), where D is the number of digits. It also includes a pre-test, with shorter loop.
In this way, a large range can be scanned without gaps:

TOTAL COUNT IN [1:8589934591]: 2219
(There are 2219 $$\\ n\in[1:2^{33}-1]\$$ such that $$\n^2\$$ has a RADD decomposition)
Took 153s on my pc.

The core routine uses byte-arrays to process large numbers.
With the new approach, squares are not explicitly calculated. However, adding odd numbers to get the next squares still creates some overhead. But it runs faster overall without conversion from int to byte-array.

Try online

#pragma GCC optimize("Ofast","unroll-loops","omit-frame-pointer","inline")
#define FRC_INLINE __attribute__((always_inline)) inline

#include <iostream>
using std::cout;

typedef uint64_t intI;
typedef __uint128_t intO;

//decimal digit
//to represent a number n, all 'digit n[i]' shall be in [0:9]
//as a temp in a calculation it may also be negative
typedef signed char digit;

//RADD pre-test
//n is an array containing L decimal digits
//returns true, if for a given n there is an R, s.t. n=R+reverse_decimal(R)
//returns false in many cases where no RADD-decomp exists
//error rate < 0.003% for n in [1:2^32]
//running time is in O(L)
template<bool check1=true>
FRC_INLINE bool preTestRADD(const digit* n, int L)
{
if (L < 1) return false;
int N = L-1;

if (n[N] == 1 && check1) {
if (n[0] == 1)
if (preTestRADD<false>(n, L)) return true;
--N;
}

for (int d = 0; N-d-d > 0; ++d) {
digit x = n[d]-n[N-d];
if (!(x==0 || x==1 || x==9 || x==-1 || x==-9)) return false;
}

return true;
}

//calc RADD-decomposition
//n is an array containing L decimal digits
//checks for a given n, if there is an R, s.t. n=R+reverse_decimal(R)
//if yes, it determines the highest such R and returns the digit count of R
//returns 0, if no RADD-decomp exists (then R is undefined)
//running time is in O(L)
template<bool check1=true>
FRC_INLINE int calcRADD(digit* R, const digit* n, int L)
{
if (L < 1) return 0;
int N = L-1;
digit cri = 0, cro = 0;

if (n[N] == 1 && check1) {
if (n[0] == 1) {
int L2;
if (L2 = calcRADD<false>(R, n, L)) return L2;
}
--N;
cro = 1;
}

int d = 0;

for (; N-d-d > 0; ++d) {
switch (n[d]-n[N-d]) {
case -9: {
if (!cri || cro) return 0;
R[d] = 0;
R[N-d] = 9;
break;
}
case -1: {
if (cri) return 0;
if (cro) {
R[d] = n[d]+1;
R[N-d] = 9;
} else {
R[d] = 0;
R[N-d] = n[d];
}
cri = cro;
cro = 1;
break;
}
case 0: {
if (cro) {
R[d] = n[d]-cri+1;
if (R[d] > 9) return 0;
R[N-d] = 9;
} else {
R[N-d] = n[d]-cri;
if (R[N-d] < 0) return 0;
R[d] = 0;
}
std::swap(cri,cro);
break;
}
case 1: {
if (!cri) return 0;
if (cro) {
R[d] = n[N-d]+1;
R[N-d] = 9;
} else {
R[d] = 0;
R[N-d] = n[N-d];
}
cri = cro;
cro = 0;
break;
}
case 9: {
if (cri || !cro) return 0;
R[d] = 0;
R[N-d] = 9;
break;
}
default: {
return 0;
}
}
}

if (!(N-d-d)) {
if ((n[d] & 1) != cri) return 0;
R[d] = n[d]/2 + 5*cro;
return N+1;
}

if (cro != cri) return 0;
return N+1;
}

constexpr int MAX_DIGITS = 40;

FRC_INLINE int getDigits(digit d[MAX_DIGITS], intO n)
{
int L = 0;
for (; L < MAX_DIGITS; ++L){
if (n <= 0xFFFFFFFFFFFFFFFF) break;
d[L] = n % 10;
n /= 10;
}
intI nt = (intI)n;
for (; L < MAX_DIGITS; ++L){
if (!nt) break;
d[L] = nt % 10;
nt /= 10;
}
return L;
}

//returns 0, if n^2 has no RADD-decomp
//otherwise it returns the highest R with n^2=R+reverse_decimal(R)
//running time is in O(log(n))
//faster for single numbers / small ranges
FRC_INLINE intO RADD2(intI n)
{
digit in[MAX_DIGITS];
digit out[MAX_DIGITS];

intO n2 = (intO)n*n;
intO nt = n2;

int L = 0;
for (; L < MAX_DIGITS; ++L){
if (nt <= 0xFFFFFFFFFFFFFFFF) break;
in[L] = nt % 10;
nt /= 10;
}
intI ntt = (intI)nt;
for (; L < MAX_DIGITS; ++L){
if (!ntt) break;
in[L] = ntt % 10;
ntt /= 10;
}

if (!preTestRADD(in, L)) return 0;

int Lout = calcRADD(out, in, L);
if (!Lout) return 0;

nt = out[Lout-1];
for (int l = Lout-2; l >= 0; --l) {
nt *= 10;
nt += out[l];
}

if (n <= 0xFFFFFFFF)
{
std::cout << n << "-> " << (intI)n2 << " R: " << (intI)nt << "\n";
} else {
std::cout << n << "-> ";
for (int l = L-1; l >= 0; --l) std::cout << (char)(in[l]+'0');
std::cout << " R: ";
for (int l = Lout-1; l >= 0; --l) std::cout << (char)(out[l]+'0');
std::cout << "\n";
}

return nt;
}

//scans range [A:B]
//no explicit calculation of squares
//faster for bigger ranges
FRC_INLINE int scan(intI A, intI B)
{
static digit n2[MAX_DIGITS];
static digit R2[MAX_DIGITS];
static digit step[MAX_DIGITS];

for (int i=0; i<MAX_DIGITS; ++i) {
step[i] = 0; n2[i] = 0;
}

intI n = A;
int n2L = getDigits(n2,(intO)n*n);
int stpL = getDigits(step,2*(intO)n-1);

int cnt = 0;
for (;;++n) {
if (!n) return cnt;

int R2L = 0;
if (preTestRADD(n2, n2L))
R2L = calcRADD(R2, n2, n2L);

if (R2L) {
++cnt;
cout << n << "-> ";
for (int l = n2L-1; l >= 0; --l) cout << (char)(n2[l]+'0');
cout << " R: ";
for (int l = R2L-1; l >= 0; --l) cout << (char)(R2[l]+'0');
cout << "\n";
}

if (n >= B) return cnt;

//n2 <- next square
int cs;
if (cs = ((*step+=2)>9 ? 1 : 0)) *step = 1;
digit cn = (*n2+=*step)>9 ? 1 : 0;
*n2 -= 10*cn;
for (int i = 1; i < MAX_DIGITS; ++i) {
if (i >= stpL && !cs && !cn) break;
if (cs = ((step[i]+=cs)>9 ? 1 : 0)) step[i] = 0;
cn = (n2[i]+=step[i]+cn)>9 ? 1 : 0;
n2[i] -= 10*cn;
}
if (step[stpL]) ++stpL;
if (n2[n2L]) ++n2L;
}
}

constexpr intI SCAN_START = 1;
//constexpr intI SCAN_START = 0xFFFFFFFF;

//constexpr intI SCAN_STOP = 1000000000;
constexpr intI SCAN_STOP = 0x1FFFFFFFF;

int main() {
int cnt = scan(SCAN_START,SCAN_STOP);

/*int cnt = 0;
for (intI n = SCAN_START; n <= SCAN_STOP; ++n) {
if (!n) break;

intO R2 = RADD2(n);

if (R2) ++cnt;
}*/

cout << "TOTAL COUNT IN [" << SCAN_START << ":" << SCAN_STOP << "]: " << cnt << std::endl;
return 0;
}

• How much do one need to modify to allow larger range, given int128 available?
– l4m2
Commented Jan 2, 2023 at 5:58
• @l4m2 thx for the suggestion. gcc __uint128_t is available on many pcs. It works, but is considerably slower. I'll try to improve it soon. Commented Jan 2, 2023 at 13:35
• If no machine-operation is available for the squaring, but uint64^2=>uint128 is a single MUL RAX
– l4m2
Commented Jan 2, 2023 at 13:59
• TOTAL COUNT IN [4294967295:8589934591]: 165， took me 10min+
– l4m2
Commented Jan 2, 2023 at 16:05
• Without int128, a potential solution is described in the answers to this question . Commented Jan 2, 2023 at 16:12

# Haskell, 4 s

sums[a]=[[a]]
sums(a:b:c)=map(a:)(sums(b:c))++map((a-1):)(filter(const(a>0))$sums((10+b):c)) pick[]=[] pick[c]=[div c 2] pick(a:b)=min 9a:pick(init b)++[max 0(a-9)] radd n=n+(read.reverse.show)n iradd n= (\s->case s of[]->0;_->(read.concat.map show.pick.head)s). filter(\x->x==reverse x).map(dropWhile(==0)).sums. map(read.(:"")).show$n

main=print\$take 1[(n,n2,i,n2-i)|n<-[1101112..],let n2=n^2,let i=iradd n2,i>0,radd i==n2]


Try it online!

26496 in 51min

Unordered but should cover [0,10^14)

#include <vector>
#include <set>
#include <bitset>
#include <iostream>
#include <math.h>

//  Sebastian

typedef uint64_t intI;
typedef __uint128_t intO;

std::set<intI> outted;

typedef signed char digit;
template<bool check1=true>
inline bool preTestRADD(const digit* n, int L)
{
if (L < 1) return false;
int N = L-1;

if (n[N] == 1 && check1) {
if (n[0] == 1)
if (preTestRADD<false>(n, L)) return true;
--N;
}

for (int d = 0; N-d-d > 0; ++d) {
digit x = n[d]-n[N-d];
if (!(x==0 || x==1 || x==9 || x==-1 || x==-9)) return false;
}

return true;
}template<bool check1=true>
inline int calcRADD(digit* R, const digit* n, int L)
{
if (L < 1) return 0;
int N = L-1;
digit cri = 0, cro = 0;

if (n[N] == 1 && check1) {
if (n[0] == 1) {
int L2;
if (L2 = calcRADD<false>(R, n, L)) return L2;
}
--N;
cro = 1;
}

int d = 0;

for (; N-d-d > 0; ++d) {
switch (n[d]-n[N-d]) {
case -9: {
if (!cri || cro) return 0;
R[d] = 0;
R[N-d] = 9;
break;
}
case -1: {
if (cri) return 0;
if (cro) {
R[d] = n[d]+1;
R[N-d] = 9;
} else {
R[d] = 0;
R[N-d] = n[d];
}
cri = cro;
cro = 1;
break;
}
case 0: {
if (cro) {
R[d] = n[d]-cri+1;
if (R[d] > 9) return 0;
R[N-d] = 9;
} else {
R[N-d] = n[d]-cri;
if (R[N-d] < 0) return 0;
R[d] = 0;
}
std::swap(cri,cro);
break;
}
case 1: {
if (!cri) return 0;
if (cro) {
R[d] = n[N-d]+1;
R[N-d] = 9;
} else {
R[d] = 0;
R[N-d] = n[N-d];
}
cri = cro;
cro = 0;
break;
}
case 9: {
if (cri || !cro) return 0;
R[d] = 0;
R[N-d] = 9;
break;
}
default: {
return 0;
}
}
}

if (!(N-d-d)) {
if ((n[d] & 1) != cri) return 0;
R[d] = n[d]/2 + 5*cro;
return N+1;
}

if (cro != cri) return 0;
return N+1;
}inline intO RADD2(intI n, intO& n2)
{
constexpr int MAX_DIGITS = 0x40;
digit in[MAX_DIGITS];
digit out[MAX_DIGITS];
n2 = (intO)n*n;

intO nt = n2;

int L;
for (L = 0; L < MAX_DIGITS; ++L){
if (nt <= 0xFFFFFFFFFFFFFFFF) break;
in[L] = nt % 10;
nt /= 10;
}
uint64_t ntt = (uint64_t)nt;
for (; L < MAX_DIGITS; ++L){
if (!ntt) break;
in[L] = ntt % 10;
ntt /= 10;
}

if (!preTestRADD(in, L)) return 0;

int Lout = calcRADD(out, in, L);
if (!Lout) return 0;

nt = out[Lout-1];
for (int l = Lout-2; l >= 0; --l) {
nt *= 10;
nt += out[l];
}

if (outted.count(n)) return 0; outted.insert(n);

if (n <= 0xFFFFFFFF)
{
std::cout << n << " -> " << (uint64_t)n2 << " R: " << (uint64_t)nt << std::endl;
} else {
std::cout << n << " -> ";
for (int l = L-1; l >= 0; --l) std::cout << (char)(in[l]+'0');
std::cout << " R: ";
for (int l = Lout-1; l >= 0; --l) std::cout << (char)(out[l]+'0');
std::cout << std::endl;
}

return nt;
}

int cnt;

void call(intI n) {
intO n2;
intO R2 = RADD2(n, n2);
if (R2) cnt++;
}

typedef unsigned int uint;
typedef unsigned long u2int;
typedef __uint128_t u4int;

std::vector<uint> ves[11111119];
//std::set<u2int> ved[10000007];

constexpr u2int pow10l(int n) {
return n?pow10l(n-1)*10:1;
}
constexpr u2int pow10g(int n) {
return n?pow10g(n-1)*10+1:0;
}

template<int P, int Q, int N>
struct ragebase { static
void rage(u2int p, uint q) {
if (ves[pow10g(Q)+q% pow10l(Q)].empty()) return;
//#pragma nounroll
for (int i=Q==0?1:0; i<19; ++i) {
ragebase<P-1, Q+1, N>::rage (p+i*pow10l(P-1), q+i*pow10l(Q));
}
}};
template<int Q, int N>
struct ragebase<0,Q,N> { static
void rage (u2int p, uint q) { //if (p==161) fprintf(stderr, "%lu/%u\n", p, q);
q = q % pow10l(N);
std::vector<uint>& r = ves[pow10g(N)+q];
//std::set<u2int>& y = ved[q];
if (r.empty()) return;
{uint pB = sqrtl(p*pow10l(N)/1000);
uint pL = pB-1, pR = pB+1;
if (pL==-1) pL=0;
for (uint i=pL; i<=pR; ++i) {
//if (y.count(i)) continue; y.insert(i);
for (uint j: r) {
call ( i*pow10l(N)+j);
}
}}
{uint pB = sqrtl(p*pow10l(N)/100);
uint pL = pB-1, pR = pB+1;
if (pL==-1) pL=0;
for (uint i=pL; i<=pR; ++i) {
//if (y.count(i)) continue; y.insert(i);
for (uint j: r) {
call ( i*pow10l(N)+j);
}
}}
{uint pB = sqrtl(p*pow10l(N)/10);
uint pL = pB-1, pR = pB+1;
if (pL==-1) pL=0;
for (uint i=pL; i<=pR; ++i) {
//if (y.count(i)) continue; y.insert(i);
for (uint j: r) {
call ( i*pow10l(N)+j);
}
}}
{uint pB = sqrtl(p*pow10l(N)/1);
uint pL = pB-1, pR = pB+3;
if (pL==-1) pL=0;
for (uint i=pL; i<=pR; ++i) {
//if (y.count(i)) continue; y.insert(i);
for (uint j: r) {
call ( i*pow10l(N)+j);
}
}}

}};
//template<int P, int Q, int N>
//using rage = ragebase<P,Q,N>::rage;

template<int N>
void run() {
std::cerr << N << ' ' << clock() << '\n';
//for (int i=0; i<pow10l(N); ++i) ves[i].clear();
//for (int i=0; i<pow10l(N); ++i) ved[i].clear();
for (u2int i=0; i<pow10l(N); ++i) ves[pow10g(N)+i*i%pow10l(N)].push_back(i);
std::cerr << N << ' ' << clock() << '\n';

ragebase<N,0,N>::rage(0,0);
}

int main() {
run<0>();
run<1>();
run<2>();
run<3>();
run<4>();
run<5>();
run<6>();
run<7>();
}


Try it online!

Old answer:

C++ (gcc), 554 in 76s

#include <stdio.h>
#include <math.h>
typedef unsigned long xint;

constexpr xint pow10(int n) {
return n?pow10(n-1)*10:1;
}

#define L(m,n,y) for(i##n=i##m; i##n<pow10(y)*19; i##n+=pow10(y)+pow10(n))
#define Lz(m,n,y) for(i##n=i##m+pow10(y)+pow10(n); i##n<pow10(y)*19; i##n+=pow10(y)+pow10(n))
#define L0(n) for(i##n=0; i##n<pow10(n)*19; i##n+=pow10(n)*2)
#define L1(n) for(i##n=0; i##n<pow10(n)*11*19; i##n+=pow10(n)*11)
#define H(s) { t=sqrtl(i0); if (t*t==i0) printf("%lu: %lu %lu %lu %lu %lu %lu %lu %lu %lu %lu\n", i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10); }

#define M1(n) Lz(1,0,n) H()
#define M2(n) L(2,1,n-1) M1(n)
#define M3(n) L(3,2,n-2) M2(n)
#define M4(n) L(4,3,n-3) M3(n)
#define M5(n) L(5,4,n-4) M4(n)
#define M6(n) L(6,5,n-5) M5(n)
#define M7(n) L(7,6,n-6) M6(n)
#define M8(n) L(8,7,n-7) M7(n)
#define M9(n) L(9,8,n-8) M8(n)

#define P(n) L0(n) M##n(2*n)  L1(n) M##n(2*n+1)

int main() {

run();
xint i0=0, i1=0, i2=0, i3=0, i4=0, i5=0, i6=0, i7=0, i8=0, i9=0, i10=0, t;

//L0(1) Lz(1,0,2) H()

//L1(1) Lz(1,0,3) H()

P(1)P(2)P(3)P(4)P(5)

}


}


Try it online!

Search all numbers that is A+rev(A). Slower than creent record but I still post it

• Your initial #define P(n) L0(n) M##n(2*n) L1(n) M##n(2*n+1) macro magic looked more interesting to read. First version is mostly best. I also returned to my initial algorithm after quietly experimenting with asymmetry cuts and discarding those failures. Commented Jan 4, 2023 at 20:21

# Scala

I was trying to rewrite @Roman Czyborra's Haskell solution in scala.

$$\ \color{red}{\text{But the Scala code outputs None}}\$$

Any help would be appreciated.

TIO

object Main {

def sums(lst: List[Int]): List[List[Int]] = lst match {
case Nil      => Nil
case a :: Nil => List(List(a))
case a :: rest =>
sums(rest)
.flatMap(l => List((a :: l), (a - 1) :: l))
.filter(l => l.head > 0)
}

def pick(lst: List[Int]): List[Int] = lst match {
case Nil      => Nil
case c :: Nil => List(c / 2)
case a :: rest =>
val updatedRest = pick(rest.init)
(math.min(9, a) :: updatedRest) :+ math.max(0, a - 9)
}

def radd(n: Int): Int = n + (n.toString.reverse.toInt)

def iradd(n: Int): Int = {
val sumsList = sums(n.toString.map(_.asDigit).toList)
val filteredList =
sumsList.filter(x => x == x.reverse).map(_.dropWhile(_ == 0))
if (filteredList.isEmpty) 0
else pick(filteredList.head).mkString("").toInt
}

def main(args: Array[String]): Unit = {
val nums = Stream.from(1101112).map(n => (n, n * n))
val output = nums.collectFirst {
case (n, n2) => {
val i = iradd(n2)
if (i > 0 && radd(i) == n2)
Some((n, n2, i, n2 - i))
else None
}
}.flatten
println(output)
}

}