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.
You are encouraged to solve this task according to the task description, using any language you may know.
- 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
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
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 %
integer suCount, addOnes, powerOf10, n;
% find the sub unit squares - note that 1 is a special case %
writeon( "1" );
suCount := 1;
addOnes := 1;
powerOf10 := 10;
n := -5;
while suCount < 5 do begin
integer sq, addUnitSq;
n := n + 10;
sq := n * n;
while sq > powerOf10 do begin
% the square has more digits than the previous one %
powerOf10 := powerOf10 * 10;
addOnes := ( addOnes * 10 ) + 1
end;
addUnitSq := sq + addOnes;
if addUnitSq < powerOf10 then begin
integer root;
% squaring the number with one added to the digits has the %
% same number of digits as the original square %
root := entier( sqrt( addUnitSq ) );
if root * root = addUnitSq then begin
logical noZeros;
integer v;
% 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 %
noZeros := true;
v := addUnitSq;
while noZeros and v > 0 do begin
noZeros := v rem 10 not = 0;
v := v div 10
end;
if noZeros then begin
writeon( i_w := 1, s_w := 0, " ", addUnitSq );
suCount := suCount + 1
end if_noZeros
end if_addUnitSq_is_a_square
end if_addUnitSq_lt_powerOf10
end while_suCount_lt_5
end.- Output:
1 36 3136 24336 5973136
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
# syntax: GAWK -f SUB-UNIT_SQUARES.AWK
# converted from FreeBASIC
BEGIN {
n = 1
while (count <= 7) {
n2 = n * n
if (!has_zero(n2)) {
m2 = dec_digits(n2)
m = int(sqrt(m2))
if (m * m == m2) {
printf("%d ",n2)
count++
}
}
n++
}
printf("\n")
exit(0)
}
function dec_digits(n, m,t) {
m = 11111111111
while (1) {
t = n - m
m = int(m/10)
if (t >= 0) {
break
}
}
return(t)
}
function has_zero(n) {
return(n ~ "0" ? 1 : 0)
}
- Output:
1 36 3136 24336 118336 126736 5973136 71526293136
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>
std::vector<uint32_t> find_digits(uint64_t number) {
std::vector<uint32_t> result;
while ( number != 0 ) {
result.emplace(result.begin(), number % 10);
number /= 10;
}
return result;
}
uint64_t find_sub_unit(const std::vector<uint32_t>& digits) {
uint64_t result = 0;
for ( const uint32_t& digit : digits ) {
result = 10 * result + ( digit - 1 );
}
return result;
}
bool is_square(const uint64_t& number) {
const uint64_t squareRoot = std::sqrt(number);
return number == squareRoot * squareRoot;
}
bool contains(const std::vector<uint32_t>& digits, uint32_t digit) {
return std::find(digits.begin(), digits.end(), digit) != digits.end();
}
int main() {
const uint32_t sub_units_required = 7;
std::cout << "The first " << sub_units_required << " sub-unit squares are:" << std::endl;
std::cout << 1 << std::endl;
uint64_t number = 2;
uint32_t count = 1;
while ( count < sub_units_required ) {
const uint64_t square = number * number;
std::vector<uint32_t> digits = find_digits(square);
if ( ! contains(digits, 0) && digits[0] != 1 ) {
const uint64_t sub_unit = find_sub_unit(digits);
if ( is_square(sub_unit) ) {
std::cout << square << std::endl;
count++;
}
}
number++;
}
}
- Output:
The first 7 sub-unit squares are: 1 36 3136 24336 5973136 71526293136 318723477136
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.
func has9 v .
while v > 0
if v mod 10 = 9
return 1
.
v = v div 10
.
return 0
.
ones = 1
pow10 = 10
while count < 5
while n * n > pow10
pow10 *= 10
ones = ones * 10 + 1
.
if has9 (n * n) = 0
sq = n * n + ones
if sqrt sq mod 1 = 0
write " " & sq
count += 1
.
.
n += 1
.- Output:
1 36 3136 24336 5973136
see Talk:Sub-unit_squares#Our_friend_the_digital_root
//Sub-unit squares. Nigel Galloway: March 4th., 2026
let dr=[|[1UL;4UL;7UL;9UL];[1UL];[9UL];[1UL;4UL;7UL];[4UL];[9UL];[1UL;4UL;7UL];[7UL];[9UL]|]
let start numD dr=((pown 10UL numD)-1UL)/9UL+dr-uint64(numD%9)
let rec zFree n=match System.Math.DivRem(n,10UL)
with _,0UL->false |0UL,_->true |n,_->zFree n
let izSq (n:uint64)=let g=sqrt(float n)|>uint64 in n=g*g
let rec fG n g=seq{let r=n*100UL+36UL in if zFree n && izSq r then yield r
let n=n+9UL in if n<g then yield! fG n g}
let fN numD=
let s=((pown 10UL numD)-1UL)/9UL
dr[numD%9]|>Seq.collect(fun n->fG(start(numD-2) n)(pown 10UL(numD-2))|>Seq.filter(fun n->izSq(n-s)))
seq{yield 1UL; yield 36UL; yield! Seq.initInfinite((+)3)|>Seq.collect(fun g->fN g)}|>Seq.iter(printfn "%d")
- Output:
1 36 3136 24336 118336 126736 5973136 71526293136 113531259136 137756776336 318723477136
Function HasZero(n As Uinteger) As Boolean
Return Iif(Instr(Str(n), "0"), True, False)
End Function
Function DecDigits(n As Uinteger) As Uinteger
Dim As Integer t, m = 1111111111
Do
t = n - m
m /= 10
Loop Until t >= 0
Return t
End Function
Dim As Uinteger n = 1, c = 0
Dim As Uinteger n2, m, m2
Do
n2 = n*n
If Not HasZero(n2) Then
m2 = DecDigits(n2)
m = Sqr(m2)
If m*m = m2 Then
Print n2 & " ";
c += 1
End If
End If
n += 1
Loop Until c >= 7
Sleep
- Output:
Same as XPL0 entry.
import Data.Char ( digitToInt , intToDigit )
isSquareNumber :: Int -> Bool
isSquareNumber n = root * root == n
where
root :: Int
root = floor $ sqrt $ fromIntegral n
isSubunitSquare :: Int -> Bool
isSubunitSquare n
|elem '0' numstr = False
|otherwise = isSquareNumber n && isSquareNumber subunitnum
where
numstr :: String
numstr = show n
subunitnum :: Int
subunitnum = read $ map (intToDigit . pred . digitToInt ) numstr
solution :: [Int]
solution = take 7 $ filter isSubunitSquare [1,2..]
- Output:
[1,36,3136,24336,118336,126736,5973136]
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.))
import java.util.ArrayList;
import java.util.List;
public class Sub_unitSquares {
public static void main(String[] args) {
final int sub_unitsRequired = 7;
System.out.println("The first " + sub_unitsRequired + " sub-unit squares are:");
System.out.println(1);
long number = 2;
int count = 1;
while ( count < sub_unitsRequired ) {
final long square = number * number;
List<Integer> digits = digits(square);
if ( ! digits.contains(0) && digits.get(0) != 1 ) {
final long sub_unit = sub_unit(digits);
if ( isSquare(sub_unit) ) {
System.out.println(square);
count += 1;
}
}
number += 1;
}
}
private static boolean isSquare(long number) {
final long squareRoot =(long) Math.sqrt(number);
return number == squareRoot * squareRoot;
}
private static long sub_unit(List<Integer> digits) {
long result = 0;
for ( int digit : digits ) {
result = 10 * result + ( digit - 1 );
}
return result;
}
private static List<Integer> digits(long number) {
List<Integer> result = new ArrayList<Integer>();
while ( number != 0 ) {
result.add(0, (int) ( number % 10 ) );
number /= 10;
}
return result;
}
}
- Output:
The first 7 sub-unit squares are: 1 36 3136 24336 5973136 71526293136 318723477136
This doesn't include 1, as if numbers with leading ones are disallowed then 1 being an exception doesn't really make sense to me. Uses the fact that all real sub-unit squares (i.e. not 1) end in 36.
function is_square(n) {
return Number.isInteger(Math.sqrt(n));
}
function is_subunit(n) s (i.e.{
const s = [...String(n)].map(Number);
if (s.includes(0) || s[0] === 1) {
return false;
}
const nsub = Number(s.map(x => x - 1).join(""));
if (is_square(n) && is_square(nsub)) {
console.log(`${n} = ${Math.sqrt(n)}^2 and ${nsub} = ${Math.sqrt(nsub)}^2`);
return true;
}
return false;
}
let count = 0, n = 36;
while (count < 4) {
if (is_subunit(n)) {
count++;
}
n += 100;
}
- Output:
36 = 6^2 and 25 = 5^2 3136 = 56^2 and 2025 = 45^2 24336 = 156^2 and 13225 = 115^2 5973136 = 2444^2 and 4862025 = 2205^2
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
""" Rosetta Code task rosettacode.org/wiki/Sub-unit_squares """
issquare(n::Integer) = isqrt(n)^2 == n
from_subsquare(n) = (d = digits(n); all( x -> x < 9, d) && issquare(evalpoly(10, d .+ 1)))
println("Subunit squares up to 10^12: ",
pushfirst!([evalpoly(10, digits(i^2) .+ 1) for i in 5:10:10^6 if from_subsquare(i^2)], 1))
- Output:
Subunit squares up to 10^12: [1, 36, 3136, 24336, 5973136, 71526293136, 318723477136]
Alternative procedural version.
# 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
begin # find the sub unit squares - note that 1 is a special case
local maxCount = 7
local suSquares = collect( 1 : maxCount ) # array initially 1, 2, 3, ...
local suCount = 1
local addOnes = 1
local powerOf10 = 10
local n = -5
while suCount < maxCount
n += 10
sq = n^2
while sq > powerOf10
# the square has more digits than the previous one
powerOf10 *= 10
addOnes *= 10
addOnes += 1
end
addUnitSq = sq + addOnes
if addUnitSq < powerOf10
# squaring the number with one added to the digits has the
# same number of digits as the original square
local root = isqrt( addUnitSq )
if root^2 == addUnitSq
# the addUnitSquare is actually a square
# make sure there are no 0s in the addUnitSquare as they
# can't be decremented to a positive digit
if all( x -> x != 0, digits( addUnitSq ) )
suCount += 1
suSquares[ suCount ] = addUnitSq
end
end
end
end
println( suSquares )
end
- Output:
[1, 36, 3136, 24336, 5973136, 71526293136, 318723477136]
Tested with Lua 5.2.4
do --[[ 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
io.write( "1" )
local suCount = 1
local addOnes = 1
local powerOf10 = 10
local n = -5
while suCount < 7 do
n = n + 10
local sq = n * n
while sq > powerOf10 do
-- the square has more digits than the previous one
powerOf10 = powerOf10 * 10
addOnes = ( addOnes * 10 ) + 1
end
local addUnitSq = sq + addOnes
if addUnitSq < powerOf10
then
-- squaring the number with one added to the digits has the
-- same number of digits as the original square
local root = math.floor( math.sqrt( addUnitSq ) )
if root * root == addUnitSq
then
-- the addUnitSquare is actually a square
-- make sure there are no 0s in the addUnitSquare as they
-- can't be decremented to a positive digit
local noZeros = true
local v = addUnitSq
while v > 0 and noZeros do
noZeros = v % 10 ~= 0
v = math.floor( v / 10 )
end
if noZeros
then
io.write( " ", addUnitSq )
suCount = suCount + 1
end
end
end
end
end
- Output:
1 36 3136 24336 5973136 71526293136 318723477136
subUnitSquare[1] = True;
subUnitSquare[n_] := With[{dig = IntegerDigits[n]},
IntegerQ[Sqrt[n]] &&
Min[dig] > 0 &&
First[dig] != 1 &&
IntegerQ[Sqrt[FromDigits[dig - 1]]]];
res = {};
n = 1;
While[Length[res] < 5,
If[subUnitSquare[n^2], AppendTo[res, n^2]];
n++;
];
res // TableForm
- Output:
1 36 3136 24336 5973136
ctr = 1
tt = time
while ctr < 1000000
a = ctr ^ 2
aa = str(a)
if a == 1 or (a % 100 == 36 and aa.indexOf("0") == null) then
b = a - ("1" * aa.len).val
c = b ^ 0.5
if floor(c) == c then print a
end if
ctr +=1
while not (ctr % 10 == 4 or ctr % 10 == 6)
ctr += 1
end while
end while
- Output:
1 36 3136 24336 118336 126736 5973136 71526293136 113531259136 137756776336 318723477136
import std/[algorithm, math, strutils]
func digits(n: Positive): seq[int] =
## Return the sequence of digits of "n".
var n = n.Natural
while n != 0:
result.add n mod 10
n = n div 10
result.reverse()
func toInt(digits: seq[int]): int =
## Convert a sequence of digits to an integer.
for d in digits:
result = 10 * result + d
func isSquare(n: int): bool =
## Return true if "n" is square.
let r = sqrt(n.toFloat).int
result = r * r == n
echo "First eight sub-unit squares:"
echo 1
var n = 0
var count = 1
while count < 8:
inc n, 5
block Check:
var digs = digits(n * n)
for d in digs.mitems:
if d == 9: break Check
inc d
let s = digs.toInt
if s.isSquare:
inc count
echo s
- Output:
First eight sub-unit squares: 1 36 3136 24336 5973136 71526293136 318723477136 264779654424693136
...but only finding 5 sub-unit squares as integers are 32-bit in oberonc and OBNC.
MODULE SubUnitSquares;
(* 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 *)
IMPORT Out, Math;
VAR suCount, addOnes, powerOf10, n, sq, addUnitSq, root, v : INTEGER;
noZeros : BOOLEAN;
BEGIN
(* find the sub unit squares - note that 1 is a special case *)
Out.String( "1" );
suCount := 1;
addOnes := 1;
powerOf10 := 10;
n := -5;
WHILE suCount < 5 DO
INC( n, 10 );
sq := n * n;
WHILE sq > powerOf10 DO
(* the square has more digits than the previous one *)
powerOf10 := powerOf10 * 10;
addOnes := ( addOnes * 10 ) + 1
END;
addUnitSq := sq + addOnes;
IF addUnitSq < powerOf10 THEN
(* squaring the number with one added to the digits has the *)
(* same number of digits as the original square *)
root := FLOOR( Math.sqrt( FLT( addUnitSq ) ) );
IF root * root = addUnitSq 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 *)
noZeros := TRUE;
v := addUnitSq;
WHILE noZeros & ( v > 0 ) DO
noZeros := v MOD 10 # 0;
v := v DIV 10
END;
IF noZeros THEN
Out.String( " " );Out.Int( addUnitSq, 0 );
INC( suCount )
END
END
END
END
END SubUnitSquares.
- Output:
1 36 3136 24336 5973136
simple brute force.
program SubUnitSquares;
{$OPTIMIZATION ON,ALL}
procedure CheckForZero(n,i,j:NativeUint);
var
q,dgt: NativeUint;
begin
repeat
q := n DIV 10;
dgt := n-10*q;
if dgt = 0 then
EXIT;
n := q
until q = 0;
j := (j-1) SHR 1;i := (i-1) SHR 1;
writeln(j:10,sqr(j):20,' -> ',sqr(i):20,i:10);
end;
var
sqI,sqJ,dSqrI,dSqrJ,Pot10Limit,RepUnit : NativeUInt;
BEGIN
Pot10Limit := 10;
RepUnit := 1;
sqj := 1; // 1*1
sqI := RepUnit;// 0+RepUnit
dSqrI := 1;
dSqrJ := 2*1+1;
repeat
if sqJ = sqI then
CheckForZero(sqj,dSqrI,dSqrJ);
repeat
sqJ += dSqrJ;
inc(dSqrJ,2);
until sqJ>=sqI;
//one more digit
if sqJ >= Pot10Limit then
begin
// choose dSqrI and sqJ so, that next sqJ starts with digit 2
RepUnit += Pot10Limit;
dSqrI := trunc(sqrt(2*Pot10Limit-RepUnit));
dSqrJ := trunc(sqrt(2*Pot10Limit));
Pot10Limit *= 10;
sqI := dSqrI*dSqrI+RepUnit;
dSqrI := 2*dSqrI+1;
sqJ := dSqrJ*dSqrJ;
dSqrJ := 2*dSqrJ+1;
end;
repeat
sqI += dSqrI;
inc(dSqrI,2);
until sqI>=sqJ;
until dSqrJ >2*514567445; //2*2*1000*1000*1000;//
END.
- @TIO.RUN:
1 1 -> 0 0
6 36 -> 25 5
56 3136 -> 2025 45
156 24336 -> 13225 115
2444 5973136 -> 4862025 2205
267444 71526293136 -> 60415182025 245795
564556 318723477136 -> 207612366025 455645
514567444 264779654424693136 -> 153668543313582025 392005795
@home 0.22user before 0.34
@TIO.RUN Real time: 0.979 s before Real time: 0.871 s
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
A more efficient approach:
use 5.036;
use bigint try => 'GMP';
use ntheory qw(:all);
sub difference_of_two_squares_solutions ($n) {
my @solutions;
my $limit = sqrtint($n);
foreach my $divisor (divisors($n)) {
last if $divisor > $limit;
my $p = $divisor;
my $q = $n / $divisor;
($p + $q) % 2 == 0 or next;
my $x = ($q + $p) >> 1;
my $y = ($q - $p) >> 1;
unshift @solutions, [$x, $y];
}
return @solutions;
}
my $N = 20; # how many terms to compute
my %seen = (1 => 1);
my $index = 1;
say($index, ': ', 1);
OUTER: for (my $n = 1 ; ; ++$n) {
my $r = (10**$n - 1) / 9;
foreach my $xy (difference_of_two_squares_solutions($r)) {
my $xsqr = $xy->[0]**2;
my @d = todigits($xsqr);
next if $d[0] == 1;
next if !vecall { $_ } @d;
next if !is_square(fromdigits([map { $_ - 1 } @d]));
if (!$seen{$xsqr}++) {
say(++$index, ': ', $xsqr);
last OUTER if ($index >= $N);
}
}
}
- Output:
1: 1 2: 36 3: 3136 4: 24336 5: 5973136 6: 71526293136 7: 318723477136 8: 264779654424693136 9: 24987377153764853136 10: 31872399155963477136 11: 58396845218255516736 12: 517177921565478376336 13: 252815272791521979771662766736 14: 518364744896318875336864648336 15: 554692513628187865132829886736 16: 658424734191428581711475835136 17: 672475429414871757619952152336 18: 694688876763154697414122245136 19: 711197579293752874333735845136 20: 975321699545235187287523246336
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 (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 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 - 32bit only shows the first 11, in 0s)
Found 15 sub-unit squares in 2.2s: 1 36 3136 24336 118336 126736 5973136 71526293136 113531259136 137756776336 318723477136 118277763838638336 118394723955598336 126825269186126736 264779654424693136
Should you uncomment the additional test so the output matches A061844, it'll only show 8 (or 7 under 32-bit).
gmp version
You can run this online here.
with javascript_semantics
include mpfr.e
constant match_A061844 = true
function sub_unit_squares(integer digits)
progress("sub_unit_squares(%d)\r",{digits})
mpz d1s = mpz_init(repeat('1',digits))
sequence res = {}
mpz {g,a,b,t} = mpz_inits(4)
atom t1 = time()+1
sequence fs = mpz_factors({"mpz",d1s},1,false)
for i,f in fs do
mpz_fdiv_q(g,d1s,f) -- g = d1s/f
mpz_add(a,f,g)
mpz_sub(b,f,g)
assert(mpz_fdiv_q_ui(a,a,2)=0) -- a = (f+g)/2
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
if time()>t1 then
progress("sub_unit_squares(%d): factor %d/%d...\r",{digits,i,length(fs)})
t1 = time()+1
end if
end for
progress("")
res = unique(res)
return res
end function
atom t = time()
sequence res = {}, tests = tagset(32)
--tests &= {40,42,54,60,64,66,72,84,90,96} -- see below
for digits in tests do
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.5s: 1 36 3136 24336 5973136 71526293136 318723477136 264779654424693136 24987377153764853136 31872399155963477136 58396845218255516736 517177921565478376336 252815272791521979771662766736 518364744896318875336864648336 554692513628187865132829886736 658424734191428581711475835136 672475429414871757619952152336 694688876763154697414122245136 711197579293752874333735845136 975321699545235187287523246336 23871973274358556957126877486736 25347159162241162461433882565136 34589996454813135961785697637136
Should you set match_A061844 to false, ie include numbers >1 that begin with 1, it'll find another 27 bringing the total to 50.
While sub_unit_squares(38) exceeded my patience, let alone sub_unit_squares(180), at a push and cheating slightly but at least matching/verifying A061844, extending tests with {40,42,54,60,64,66,72,84,90,96} gets (or 156 when match_A061844 is false):
Found 77 sub-unit squares in 11 minutes and 26s: 1 36 3136 24336 5973136 71526293136 318723477136 264779654424693136 24987377153764853136 31872399155963477136 58396845218255516736 517177921565478376336 252815272791521979771662766736 518364744896318875336864648336 554692513628187865132829886736 658424734191428581711475835136 672475429414871757619952152336 694688876763154697414122245136 711197579293752874333735845136 975321699545235187287523246336 23871973274358556957126877486736 25347159162241162461433882565136 34589996454813135961785697637136 2858541763747552538199941619545257144336 214785886789716796533667464535274377236736 233292528132679183463629157143235636286736 244671849793441155421899813243325528686736 271571567929448516411695557685613529966736 322388381596588665613523969581347191316736 385414415625146742626881165526237149942336 494827714874767379344736911473964125592336 729191918879671448289782722539515523333136 739265858539339252384919139328667324488336 451616391374794616993675837721511769881724292768597136 225155131198959884498111695266319246921397575282262981916736 235533149167838594417145126422746861139269181233838789763136 247165591495471317848341556972927924777146487462821199752336 253185584977741571193838716134142968698828866536879387589136 255118517389891411837338761327472341665993165261943933973136 258794681328961385589665886866946129315831982462313647221136 264638279877394985436962113298797369525794994152638329716736 274993879679611537524346229892753979256782689558714439941136 317227344555757949695298635797266762327152673658761948316736 322778753177461749385975279423679631448585894831889496126736 359789363125538627844318913392816245152758765579763311486736 393441389391373434369197686452869535953373734343532842366736 414994447929493177987216182314238968483589993433234743387136 425694769821125433165871546269244331365369797668913817432336 432644594437123781553749497955531461492791463379194125966736 437562512527813432113133967677338524499138828662399714286736 449852946468284983783199424761436133724851432965597753717136 482134751872773928421842374421611988974497596294794558126736 521975939684451673253113662113934173792645168416634454516736 527522135343484614696922441512961219362733539657782899395136 562641631377938237248318531478861249592277376289475543144336 571185124453583434863858314331338695842867554344647955413136 572684988569547784165849427621699664494716729842793127896336 614898423328697552174926774454594881158921253578236236792336 633696482563387284528138567949516822355885534948846913486736 643289923123766413527889349998824883237881684711852761436736 667191497411655246112479559487494956367164427966476793774336 687676275633214498821343918455779943946116536644394167797136 738581581372137137454777164721328295342861142971169656536336 762742925841272953426291121719451613531524131158969765243136 775543441492587483465982446818139435456241975693619925238336 848683965199869214645864196533185358672353998728624638126736 947138642349237674483499271295327433171899196661387174341136 964891172337149147689484881922848848795484622547192848916736 5569531161785335471834187656758668584186145889911322362616238336 392218699139394917429267256144352916588419663779199636243639414336 848357125461439331518387328465823399571373158999455173816177975229891136 416423361382974375988119796395668527245964353555414577638838888367828457388415992336 479474784622574743882185337518915188918922331343493318921536469363248612788886766736 241835624857156131887688175431621247512542273949817527242228179562863761282638896774325136 398574359223975871554875716496874244751857937743666637999272258371878627497383572366715136 917299545663587142783558195917968989124276866721326391286476151415422729817578418459766336 689912968699291934123339418659749956627118749536638166817937181829315824413919788889194395285136
I will note that in sub_unit_squares(96), mpz_factors() takes 7 mins, vs. mpz_pollard_rho() 3 mins, and builds a list of 4,194,304 factors, of which we're only really interested in a middle block of just 171,710, so there may be a potential two-fold or more speedup in there, somewhere. Alas the needed mpz_pollard_rho() [as occurs within mpz_factors] of sub_unit_squares(38|180) also exceeded my patience, so it wouldn't help with that.
local int = require "int"
print("The first 7 sub-unit squares are:")
print(1)
local i = 2
local count = 1
while count < 7 do
local sq = i * i
local digits = int.digits(sq)
if digits[1] != 1 and !digits:contains(0) then
local sum = digits[1] - 1
for j = 2, #digits do sum = sum * 10 + digits[j] - 1 end
if int.issquare(sum) then
print(sq)
count += 1
end
end
i +=1
end
- Output:
The first 7 sub-unit squares are: 1 36 3136 24336 5973136 71526293136 318723477136
# sub-unit_squares.py by Xing216
from gmpy2 import is_square
def digits(n: int) -> list[int]:
return [int(d) for d in str(n)]
def get_sub_unit(digits_n: list[int], n: int) -> int:
return int(''.join([str(d-1) for d in digits_n]))
def is_sub_unit_square(n: int) -> bool:
if not is_square(n): return False
digits_n = digits(n)
if 0 in digits_n: return False
sub_unit = get_sub_unit(digits_n, n)
if len(str(sub_unit)) != len(str(n)): return False
return is_square(sub_unit)
res = []
i = 1
while len(res) < 5:
if is_sub_unit_square(i): res.append(i)
i += 1
print(res)
- Output:
[1, 36, 3136, 24336, 5973136]
digits <- function(n, nd) rev((n %/% (10^(1:nd - 1))) %% 10)
recombine <- function(v) Reduce(function(x, y) 10*x + y, v)
is_square <- function(n) sqrt(n) %% 1 == 0
is_subunit <- function(n, nd) {
digs <- digits(n, nd)
if (0 %in% digs || digs[1] == 1) return(list(result = FALSE))
nsub <- recombine(digs - 1)
res <- is_square(n) && is_square(nsub)
list(result = res, n = n, rn = sqrt(n), nsub = nsub, rnsub = sqrt(nsub))
}
nd <- 2
count <- 0
n <- 36
while (count < 4) {
test <- is_subunit(n, nd)
if (test$result) {
do.call(sprintf, c("%i = %i^2 and %i = %i^2", tail(test, -1))) |> cat("\n")
count <- count + 1
}
n <- n + 100
if (n > 10^(nd)) nd <- nd + 1
}
- Output:
36 = 6^2 and 25 = 5^2 3136 = 56^2 and 2025 = 45^2 24336 = 156^2 and 13225 = 115^2 5973136 = 2444^2 and 4862025 = 2205^2
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
Rebol [
title: "Rosetta code: Sub-unit squares"
file: %Sub-unit_squarese.r3
url: https://rosettacode.org/wiki/Sub-unit_squares
]
;; A "sub-unit square" is a perfect square where:
;; 1. No digit is 0 (since subtracting 1 would produce a -1 digit)
;; 2. Its last two digits are "36" (necessary condition — see note below)
;; 3. Subtracting 1 from every digit yields another perfect square
sub-unit-square?: function [
"Return true if the number is a sub-unit square."
num [integer!]
][
;; Convert the number to a mutable string of digit characters
str: append clear "" num
all [
;; Reject any number containing a 0 digit — subtracting 1 would be invalid
not find str #"0"
;; Fast pre-filter: sub-unit squares must end in "36"
;; (only squares ending in 36 can produce a square when each digit is decremented)
"36" == skip tail str -2
;; Decrement every digit character in-place by 1
;; (e.g. "149" becomes "038", which as an integer is 38)
forall str [str/1: str/1 - 1]
;; Parse the modified string back into an integer
num: to integer! str
;; Compute the integer square root of the decremented number
sqr: to integer! square-root num
;; Confirm it is a perfect square (guards against floating-point rounding)
sqr * sqr == num
]
]
;; Collect the first 10 sub-unit squares
sub-unit-squares: copy []
n: 1
while [10 > length? sub-unit-squares][
num: n * n
;; Test the current square and collect it if it qualifies
if sub-unit-square? num [ append sub-unit-squares num ]
;; Advance to the next integer so we test the next perfect square
++ n
]
print "The first 10 sub-unit squares:"
print sub-unit-squares
- Output:
The first 10 sub-unit squares: 36 3136 24336 118336 126736 5973136 71526293136 113531259136 137756776336 318723477136
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
fn is_square_number( number : u64 ) -> bool {
let root : f64 = (number as f64).sqrt( ).floor( ) ;
(root as u64) * (root as u64) == number
}
fn is_subunit_square( number : u64 ) -> bool {
let numberstring = number.to_string( ) ;
let numstring : &str = numberstring.as_str( ) ;
if numstring.contains( '0' ) {
return false ;
}
else {
if number < 10 {
is_square_number( number ) && is_square_number( number - 1 )
}
else {
let mut digits : Vec<u64> = Vec::new( ) ;
let mut num : u64 = number ;
while num != 0 {
digits.push( num % 10 ) ;
num /= 10 ;
}
for n in digits.iter_mut( ) {
*n -= 1 ;
}
let mut sum : u64 = 0 ;
let mut factor : u64 = 1 ;
for d in digits {
sum += d * factor ;
factor *= 10 ;
}
is_square_number( sum ) && is_square_number( number )
}
}
}
fn main() {
let mut solution : Vec<u64> = Vec::new( ) ;
let mut current : u64 = 1 ;
while solution.len( ) != 7 {
if is_subunit_square( current ) {
solution.push( current ) ;
}
current += 1 ;
}
println!("{:?}" , solution ) ;
}
- Output:
[1, 36, 3136, 24336, 118336, 126736, 5973136]
Using the method described in the OEIS entry, representing (10^n - 1)/9 as a difference of squares x^2 - y^2 and checking x^2 to see if it satisfies the conditions.
var N = 20 # how many terms to compute
var arr = Set(1)
for n in (1..Inf) {
var r = (10**n - 1)/9
arr << r.diff_of_squares.map{.head}.map{.sqr.digits}.grep {|d|
(d[-1] != 1) && d.none{.is_zero} && d.map{.dec}.digits2num.is_square
}.map{.digits2num}...
break if (arr.len >= N)
}
arr.sort.first(N).each_kv {|k,n|
say "#{'%2d' % k+1}: #{n}"
}
- Output:
1: 1 2: 36 3: 3136 4: 24336 5: 5973136 6: 71526293136 7: 318723477136 8: 264779654424693136 9: 24987377153764853136 10: 31872399155963477136 11: 58396845218255516736 12: 517177921565478376336 13: 252815272791521979771662766736 14: 518364744896318875336864648336 15: 554692513628187865132829886736 16: 658424734191428581711475835136 17: 672475429414871757619952152336 18: 694688876763154697414122245136 19: 711197579293752874333735845136 20: 975321699545235187287523246336
Square candidates, check for zeros, check for sub-unit squares.
▽=⊸⁅√⊸⍜°⋕-₁▽⊸≡(≤˜⨂@0⟜⧻°⋕)ⁿ2↘1⇡1000000- Output:
[1 36 3136 24336 118336 126736 5973136 71526293136 113531259136 137756776336 318723477136]
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
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
import "std/vec.zc"
import "std/math.zc"
fn digits<T>(n: T) -> Vec<int> {
let digs = Vec<int>::new();
if n == 0 {
digs << 0;
return digs;
}
while n > 0 {
digs << (n % 10);
n /= 10;
}
digs.reverse();
return digs;
}
fn is_square<T>(n: T) -> bool {
let sqrt = (T)Math::round(Math::sqrt((f64)n));
return sqrt * sqrt == n;
}
fn main() {
println "The first 7 sub-unit squares are:";
println "1";
let count = 1;
for let i: u64 = 2; count < 7; i++ {
let sq = i * i;
let digs = digits(sq);
if digs[0] != 1 && !digs.contains(0) {
let sum: u64 = digs[0] - 1;
for j in 1..digs.length() { sum = sum * 10 + digs[j] - 1; }
if is_square(sum) {
println "{sq}";
count++;
}
}
}
}
- Output:
The first 7 sub-unit squares are: 1 36 3136 24336 5973136 71526293136 318723477136