Jump to content

Sub-unit squares

From Rosetta Code
Revision as of 23:06, 3 November 2023 by Xing216 (talk | contribs) (Added Python implementation for Sub-unit squares task using gmpy2 library)
Task
Sub-unit squares
You are encouraged to solve this task according to the task description, using any language you may know.

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

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

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

Works with: Delphi version 6.0


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.

EasyLang

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
.

FreeBASIC

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.

Haskell

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]

J

   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

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

""" Rosetta Code task rosettacode.org/wiki/Sub-unit_squares """

using Polynomials
evalpoly( base, d ) = polyval( Poly( d ), base )
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.

Translation of: Lua
# 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]

Lua

Translation of: ALGOL 68

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

MiniScript

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

Nim

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

Perl

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

Library: Phix/online

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 7.1s:
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)
    mpz d1s= mpz_init(repeat('1',digits))
    sequence res = {}
    mpz {g,a,b,t} = mpz_inits(4)
    for f in mpz_factors({"mpz",d1s},1,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

Should you set match_A061844 to false, ie include numbers >1 that begin with 1, it'll find fifty of them.

Cheating slightly, but at least matching/verifying A061844 all bar the last:

--for digits in tagset(32)&{40,42,54,60,64,66,72,84} do -- Found 73 sub-unit squares in 4 minutes and 17s
--for digits in {90} do -- Found 3 sub-unit squares in 4 minutes and 18s
--for digits in {96} do -- Found 1 sub-unit squares in 8 minutes and 30s
--for digits in {180} do -- no dice, eventually blows up mpz_factors()

Python

Library: gmpy2
# 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]

Raku

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

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

Rust

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]

Wren

Library: Wren-math

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

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 
Cookies help us deliver our services. By using our services, you agree to our use of cookies.