Sub-unit squares
A sub-unit square is a square number (product of two identical non-negative integers) that remains a square after having a 1 subtracted from each digit in the square.
- E.G.
The number 1 is a sub-unit square. 1 - 1 is 0, which is also a square, though it's kind-of a degenerate case.
The number 3136 is a sub-unit square. 3136 (56²) with unit 1 subtracted from each digit is 2025 (45²).
A sub-unit square cannot contain a digit zero (0) since zero minus one is negative one. Every known sub-unit square, with the exception of 1, ends with the digits 36.
- Task
- Find and display at least the first five sub-unit squares.
- See also
ALGOL 68[edit]
Assumes LONG INT is at least 64 bit (as in ALGOL 68G version 2 or 3). Uses the optimisations in the J sample, based on Wherra's observation that the squares (other than 1) must end in 36 (see the Discussion page).
BEGIN # find some sub-unit squares: squares which are still squares when 1 #
# is subtracted from each of their digits #
# the sub-unit squares must end in 36 (see task discussion) and so #
# the square formed by subtracting 1 from each digit must end in 25 #
# so, as with the J sample, we start with the "25" square and test #
# the incremented square #
# find the sub unit squares - note that 1 is a special case #
print( ( "1" ) );
INT su count := 1;
LONG INT add ones := 1;
LONG INT power of 10 := 10;
FOR n FROM 5 BY 10 WHILE su count < 7 DO
LONG INT sq = ( LENG n ) * ( LENG n );
WHILE sq > power of 10 DO
# the square has more digits than the previous one #
power of 10 *:= 10;
add ones *:= 10 +:= 1
OD;
IF LONG INT add unit sq = sq + add ones;
add unit sq < power of 10
THEN
# squaring the number with one added to the digits has the #
# same number of digits as the original square #
IF LONG INT root = ENTIER long sqrt( add unit sq );
root * root = add unit sq
THEN
# the add unit square is actually a square #
# make sure there are no 0s in the add unit square as they #
# can't be decremented to a positive digit #
IF BOOL no zeros := TRUE;
LONG INT v := add unit sq;
WHILE v > 0 AND no zeros DO
no zeros := v MOD 10 /= 0;
v OVERAB 10
OD;
no zeros
THEN
print( ( " ", whole( add unit sq, 0 ) ) );
su count +:= 1
FI
FI
FI
OD
END
- Output:
1 36 3136 24336 5973136 71526293136 318723477136
Arturo[edit]
allSquares: map 1..1000000 'x -> x^2
subUnits: allSquares | select 'z -> not? some? digits z 'w -> w=0
| select 'x -> contains? allSquares to :integer join map digits x 'k -> k-1
print ["The first" size subUnits "sub-unit squares:"]
print subUnits
- Output:
The first 10 sub-unit squares: 36 3136 24336 118336 126736 5973136 71526293136 113531259136 137756776336 318723477136
Delphi[edit]
function IsSubUnitSquare(SQ: integer): boolean;
{Returns true if you subtract one from each digit}
{and it is still square. Assume SQ is square}
var I,J: integer;
var IA: TIntegerDynArray;
begin
Result:=False;
{Get all digits}
GetDigits(SQ,IA);
{Subtract one from each digit}
for J:=0 to High(IA) do
begin
{Zeros not allowed = they would cause negative digits}
if IA[J]=0 then exit;
IA[J]:=IA[J]-1;
end;
{Turn digits into number again}
SQ:=DigitsToInt(IA);
{Test if it is square}
if Frac(Sqrt(SQ))<>0 then exit;
Result:=True;
end;
procedure ShowSubUnitSquares(Memo: TMemo);
var I,SQ,Cnt: integer;
begin
Cnt:=0;
for I:=1 to high(Integer) do
begin
SQ:=I*I;
if IsSubUnitSquare(SQ) then
begin
Inc(Cnt);
Memo.Lines.Add(Format('%d %8.0n',[Cnt,SQ+0.0]));
if Cnt>=7 then break;
end;
end;
end;
- Output:
1 1 2 36 3 3,136 4 24,336 5 118,336 6 126,736 7 5,973,136 Elapsed Time: 6.899 ms.
J[edit]
digits=: 10&#.inv
idigs=: >:&.digits"0
(#~ 1-(0 e. digits)"0) (#~ (=<.)@%:)idigs*:5*i.2e5
1 36 3136 24336 5973136 71526293136 318723477136
Here we work in the reverse direction. Since these sub-unit squares always end in 36, their decremented digit version must always end in 25 and the square roots must always end in 5. So, we start with a sequence of multiples of 5, square them, increment their digits and discard those which are not exact squares. Finally, we discard those which contain a 0 digit. (The method used to increment digits here maps 99 to 1010. Arguably, that should instead be 110, but it would get discarded, regardless. (With a little more complexity we could discard anything with a 9 digit while generating the initial sequence, this would provide only a minor gain.))
jq[edit]
The following jq program closely follows the task description and thus includes squares greater than 1 that begin with "1", even though A061844 excludes them.
Since gojq, the Go implementation of jq, uses unbounded-precision integer arithmetic, the program can be expected to produce reliable results when used with that implementation.
By contrast, the C implementation of jq uses IEEE 754 64-bit arithmetic, and so cannot be expected to produce reliable values beyond the first few, though in fact the results using the C-implementation are accurate through at least the first 11 sub-unit squares.
# The following is sufficiently accurate for the task:
def is_square: sqrt | . == floor;
# Emit a stream of "sub-unit squares"
def subunitSquares:
def sub1:
tostring | explode
| if any(. == 48) then null
else map(.-1) | implode | tonumber
end;
foreach range(1; infinite) as $i ( null;
($i * $i) as $sq
| ($sq|sub1)
| if . == null or (tostring[-1:] | IN("3", "7", "8")) then null
elif is_square then $sq
else null
end;
select(.) ) ;
11 as $count
| "The first \($count) sub-unit squares are:",
limit($count; subunitSquares)
- Output:
"The first 11 sub-unit squares are:" 1 36 3136 24336 118336 126736 5973136 71526293136 113531259136 137756776336 318723477136
Julia[edit]
""" Rosetta Code task rosettacode.org/wiki/Sub-unit_squares """
issquare(n::Integer) = isqrt(n)^2 == n
from_subsquare(n) = (d = digits(n); all(<(9), d) && issquare(evalpoly(10, d .+ 1)))
println("Subunit squares up to 10^18: ",
pushfirst!([evalpoly(10, digits(i^2) .+ 1) for i in 25:5:10^9 if from_subsquare(i^2)], 1))
- Output:
Subunit squares up to 10^18: [1, 3136, 24336, 5973136, 71526293136, 318723477136, 264779654424693136]
Perl[edit]
Like Raku, but without the laziness.
use v5.36;
sub U ($n) { $n>1 ? (2*$n)**2 : 1 }
sub L ($n) { $n**2 }
sub R ($n) { 1 x $n }
my($l,$u,$c) = (-1);
while (++$u) {
next if U($u) =~ /0/;
my $chars = length U($u);
while ($l++) {
next if U($u) - L($l) > R($chars);
last if U($u) - L($l) < R($chars);
say U($u) and ++$c == 7 and exit
}
}
- Output:
1 36 3136 24336 5973136 71526293136 318723477136
Phix[edit]
Uses the much faster method from the OEIS page.
with javascript_semantics function sub_unit_squares(integer digits, atom d1s) -- d1s is eg 1111 for digits==4 sequence res = {} for f in factors(d1s,1) do atom g = d1s/f, a = (f+g)/2, b = (f-g)/2 if a*a-b*b=d1s then string a2 = sprintf("%d",a*a) if length(a2)=digits and not find('0',a2) then res = append(res,a2) end if end if end for res = unique(res) return res end function atom n = 1, t0 = time() sequence res = {} for digits=1 to iff(machine_bits()=32?12:18) do res &= sub_unit_squares(digits,n) n = n*10+1 end for printf(1,"Found %d sub-unit squares in %s:\n%s\n",{length(res),elapsed(time()-t0),join(res,"\n")})
As per the Julia jq entry, while A061844 explicitly excludes squares beginning with 1 apart from 1 itself, this task does not.
It would of course be trivial to exclude them using either and a2[1]!='1'
, or and length(sprintf("%d",b*b))=digits
.
Technically 32 bit will go to 16 digits and 64 bit to 19, but it won't find any more answers, so there's no point.
- Output:
(64 bit)
Found 15 sub-unit squares in 7.1s: 1 36 3136 24336 118336 126736 5973136 71526293136 113531259136 137756776336 318723477136 118277763838638336 118394723955598336 126825269186126736 264779654424693136
gmp version[edit]
I've now added an mpz_factors() routine, not yet shipped but you can run this online here.
with javascript_semantics include mpfr.e constant match_A061844 = true function sub_unit_squares(integer digits) mpz d1s= mpz_init(repeat('1',digits)) sequence res = {} mpz {g,a,b,t} = mpz_inits(4) for f in mpz_factors({"mpz",d1s},false) do mpz_fdiv_q(g,d1s,f) -- g = d1s/f mpz_add(a,f,g) assert(mpz_fdiv_q_ui(a,a,2)=0) -- a = (f+g)/2 mpz_sub(b,f,g) assert(mpz_fdiv_q_ui(b,b,2)=0) -- b = (f-g)/2 mpz_mul(a,a,a) mpz_mul(b,b,b) mpz_sub(b,a,b) if mpz_cmp(b,d1s)=0 then -- a*a-b*b = d1s string a2 = mpz_get_str(a) -- (a*a by now) if length(a2)=digits and (not match_A061844 or digits=1 or a2[1]!='1') and not find('0',a2) then res = append(res,a2) end if end if end for res = unique(res) return res end function atom t = time() sequence res = {} for digits=1 to 32 do -- (next up would be 40, 38 exceeded my patience) res &= sub_unit_squares(digits) end for printf(1,"Found %d sub-unit squares in %s:\n%s\n",{length(res),elapsed(time()-t),join(res,"\n")})
- Output:
Found 23 sub-unit squares in 0.9s: 1 36 3136 24336 5973136 71526293136 318723477136 264779654424693136 24987377153764853136 31872399155963477136 58396845218255516736 517177921565478376336 252815272791521979771662766736 518364744896318875336864648336 554692513628187865132829886736 658424734191428581711475835136 672475429414871757619952152336 694688876763154697414122245136 711197579293752874333735845136 975321699545235187287523246336 23871973274358556957126877486736 25347159162241162461433882565136 34589996454813135961785697637136
Raku[edit]
First seven take about 5 seconds with this implementation. The eighth would take several hours at least.
my @upper = 1,|(1 .. *).map((* × 2)²);
my @lower = (^∞).map: *²;
my \R = [\+] 0, 1, 10, 100 … *;
my $l = 0;
.say for (gather {
(^∞).map: -> $u {
next if @upper[$u].contains: 0;
my $chars = @upper[$u].chars;
loop {
$l++ and next if @upper[$u] - @lower[$l] > R[$chars];
take @upper[$u] if @upper[$u] - @lower[$l] == R[$chars];
last;
}
}
})[^7]
- Output:
1 36 3136 24336 5973136 71526293136 318723477136
RPL[edit]
Since brute force is not an option for 4-bit computer languages, we use the method described in the OEIS page, based on the resolution of x² - y² = 11..1
The line OVER →STR 1 1 SUB "1" ≠ AND
can be removed to obtain the OEIS sequence.
FACTS
is defined at Factors of an integer
RPL code | Comment |
---|---|
≪ IFERR FP NOT THEN DROP 0 END ≫ ‘INT?’ STO ≪ ALOG 1 - 9 / DUP FACTS 1 OVER SIZE 2 / CEIL SUB → ones factors ≪ { } 1 factors SIZE FOR j factors j GET ones OVER / + 2 / SQ IF DUP XPON ones XPON == OVER →STR "0" POS NOT AND OVER →STR 1 1 SUB "1" ≠ AND OVER ones - √ INT? AND THEN + ELSE DROP END NEXT ≫ ≫ ‘SUBUNITS’ STO |
INT? ( x → boolean ) SUBUNITS ( n → { SUS } ) ones = (10^n-1)/9 factors = half of the full FACTS list is enough for each factor x = ((factor+factor/ones)/2)^2 if x has same number of digits as ones and x does not contain any zero and x does not start with 1 and x's sub-unit square is an integer then append x to list next return list |
≪ 2 { 1 } WHILE DUP SIZE 5 < REPEAT OVER SUBUNITS + SWAP 1 + SWAP END SWAP DROP ≫ EVAL
- Output:
1: { 1 36 3136 24336 5973136 }
Returned in 41 seconds on a basic HP-48SX... but finding the following SUS (71526293136) took one hour and a half, essentially to factor 11111111111 = 21649 x 513239
Wren[edit]
Apart from the number '1' itself, it looks like the first digit can't be '1' either. In other words, the number of digits in the transformed number must remain unchanged.
import "./math" for Int
System.print("The first 7 sub-unit squares are:")
System.print(1)
var i = 2
var count = 1
while (count < 7) {
var sq = i * i
var digits = Int.digits(sq)
if (digits[0] != 1 && !digits.contains(0)) {
var sum = digits[0] - 1
for (i in 1...digits.count) sum = sum * 10 + digits[i] - 1
if (Int.isSquare(sum)) {
System.print(sq)
count = count + 1
}
}
i = i + 1
}
- Output:
The first 7 sub-unit squares are: 1 36 3136 24336 5973136 71526293136 318723477136
XPL0[edit]
func HasZero(N); \Return 'true' if N contains a zero digit
int N;
[repeat N:= N/10;
if rem(0) = 0 then return true;
until N=0;
return false;
];
func DecDigits(N); \Decrement each digit of N
int N, M, T;
[M:= 1_111_111_111;
repeat T:= N-M;
M:= M/10;
until T >= 0;
return T;
];
int N, C, N2, M, M2;
[N:= 1; C:= 0;
repeat N2:= N*N;
if not HasZero(N2) then
[M2:= DecDigits(N2);
M:= sqrt(M2);
if M*M = M2 then
[IntOut(0, N2);
ChOut(0, ^ );
C:= C+1;
];
];
N:= N+1;
until C >=7;
]
- Output:
1 36 3136 24336 118336 126736 5973136