Happy 2023

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.4.0      ✔ purrr   0.3.5 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.3      ✔ forcats 0.5.2 ── Conflicts ─────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

digit <- function(n, i)
  (n %% 10^(i+1) - n %% 10^i) / 10^i

digitsum <- Vectorize(function(n) {
  digit(n, 0:floor(log10(n))) |> sum()
})

res <- tibble(n = 1:2500) |>
  mutate(div = n/digitsum(n),
         is_square = floor(sqrt(div))^2 == div) |>
  filter(is_square) |>
  mutate(i = 1:n(),
         `a(n)` = n) |>
  select(i, `a(n)`)

res |> print(n = Inf)

# A tibble: 55 × 2
       i `a(n)`
   <int>  <int>
 1     1      1
 2     2      2
 3     3      3
 4     4      4
 5     5      5
 6     6      6
 7     7      7
 8     8      8
 9     9      9
10    10     12
11    11     24
12    12     36
13    13     48
14    14     81
15    15    100
16    16    144
17    17    150
18    18    192
19    19    200
20    20    225
21    21    288
22    22    300
23    23    320
24    24    324
25    25    375
26    26    400
27    27    441
28    28    500
29    29    512
30    30    600
31    31    640
32    32    648
33    33    700
34    34    704
35    35    735
36    36    800
37    37    832
38    38    882
39    39    900
40    40    960
41    41   1014
42    42   1088
43    43   1200
44    44   1452
45    45   1458
46    46   1521
47    47   1815
48    48   2023
49    49   2025
50    50   2028
51    51   2178
52    52   2304
53    53   2312
54    54   2352
55    55   2400

res |>
  filter(`a(n)` >= 2023) |>
  pull(`a(n)`) |>
  paste(collapse = ", ")

[1] "2023, 2025, 2028, 2178, 2304, 2312, 2352, 2400"

Added bonus, candidate for A358693.

library(gmp)

is_prime <- Vectorize(function(n) {
  factors <- factorize(n)
  length(factors) == 1 && factors[1] == n
})

res <- tibble(n = 1:15000) |>
  mutate(div = n/digitsum(n),
         is_square = floor(sqrt(div))^2 == div,
         sqrt_prime = is_prime(sqrt(div)),
         prop = is_square & sqrt_prime) |>
  filter(prop) |>
  mutate(i = 1:n(),
         `a(n)` = n,
         sqrt_div = sqrt(div)) |>
  select(i, `a(n)`, sqrt_div)

res |> print(n = Inf)

# A tibble: 50 × 3
       i `a(n)` sqrt_div
   <int>  <int>    <dbl>
 1     1     12        2
 2     2     24        2
 3     3     36        2
 4     4     48        2
 5     5     81        3
 6     6    150        5
 7     7    225        5
 8     8    375        5
 9     9    441        7
10    10    735        7
11    11    882        7
12    12   1014       13
13    13   1452       11
14    14   1521       13
15    15   1815       11
16    16   2023       17
17    17   2028       13
18    18   2178       11
19    19   2312       17
20    20   2535       13
21    21   2601       17
22    22   3549       13
23    23   3610       19
24    24   4046       17
25    25   4332       19
26    26   4335       17
27    27   4624       17
28    28   4913       17
29    29   5054       19
30    30   5415       19
31    31   5491       17
32    32   5780       17
33    33   6069       17
34    34   6137       19
35    35   6358       17
36    36   6647       17
37    37   6936       17
38    38   7581       19
39    39   7942       19
40    40   8664       19
41    41   8959       17
42    42   9386       19
43    43   9522       23
44    44   9747       19
45    45  10092       29
46    46  11532       31
47    47  12321       37
48    48  12615       29
49    49  12696       23
50    50  14415       31

Someone else’s nums:

others <- c(12, 24, 36, 48, 81, 150, 225, 375, 441, 735, 882, 1014, 1452, 1521, 1815, 2023, 2028, 2178, 2312, 2535, 2601, 3549, 3610, 4046, 4332, 4335, 4624, 4913, 5054, 5415, 5491, 5780, 6069, 6137, 6358, 6647, 6936, 7581, 7942, 8664, 8959, 9386, 9522, 9747, 10092, 11532, 12321, 12615, 12696)

mine <- res |> pull("a(n)") |> head(length(others))
mine

 [1]    12    24    36    48    81   150   225   375   441   735   882  1014
[13]  1452  1521  1815  2023  2028  2178  2312  2535  2601  3549  3610  4046
[25]  4332  4335  4624  4913  5054  5415  5491  5780  6069  6137  6358  6647
[37]  6936  7581  7942  8664  8959  9386  9522  9747 10092 11532 12321 12615
[49] 12696

all(others == mine)

[1] TRUE