Source code in the picat programming language

import sat. 
to_num(List, Base, Num) => 
        Len = length(List),
        Num #= sum([List[I] * Base**(Len-I) : I in 1..Len]).

palindrom(S) =>
    N = len(S), 
    Start :: 1..N, % start at the first non-zero position:
    foreach(I in 1..N)
        I1 #= max(1, min(N, N-(I-Start))), % I1 is the symmetry index partner of I (if relevant)
        element(I1, S, S1),                % S1 is the respective digit
        I #<  Start #=> S[I] #= 0,         % skip leading 0´s
        I #=  Start #=> S[I] #> 0,         % Start points to the first non-zero digit
        I #>= Start #=> S[I] #= S1         % palindromic symmetry
    end.

constrain(Max, B, X) =>
    Len = floor(log(Max) / log(B)) + 1, % length of Max in Base B representation
    Digits = new_list(Len), Digits :: 0..B-1, 
    to_num(Digits, B, X), % Digits show the Base B representation of X
    palindrom(Digits).

main =>
    N = 11, % maximum number of decimal digits for search, can be set freely
    Max = 10**N - 1, % maximum number
    X :: 2..Max,
    constrain(Max, 2, X), 
    constrain(Max, 3, X),
    Pnumbers = solve_all([X]),
    foreach([Y] in [[0], [1]] ++ Pnumbers.sort()) % start with 0 and 1, then show solutions > 1
        printf("%w %s %s%n", Y, to_radix_string(Y,2), to_radix_string(Y,3))
    end.