CodeIq:「デジタル・ルート」問題の解答

Kawazoe(@riverplus)さんのCodeIQの問題です。
問題はこちらの御世話になりました。解答例もあります。
~~->https://blog.goo.ne.jp/r-de-r/e/bad4f9b4a5c4fecd67db04fe31c51c78~~
(リンク切れしていました2025.5.2)
URLが変わっていましたが見つかりました。(2025.5.3追加)
->https://sangaku0418.hatenablog.com/entry/2016/11/03/100100

ーーー2025..5.2追加
ある自然数の各桁の数をすべて足す操作を行います。(199->19)
更に結果が一桁になるまで、計算結果の数字に同じことを繰り返します。
kという数に上記の操作を行うのにかかった回数をF(k)とします。(F(199)=3)
1<=k<=nのすべてのkに対するF(ｋ)の和をG(n)とします。
1<=n<=10^8なるnのG(n)を求める問題です。
ーーー

こちらにも情報があります。
->https://togetter.com/li/1040869

同じ数字を含む数は同じF値となるので、
組合せを網羅する解き方はすぐ浮かびます。
場合の数が大きくなりそうなので、他の方法をいろいろ考えましたが、
その時はうまくいかず、結局組み合わせを利用しました。
数字の組み合わせのみ直接求めて、その組み合わせの並び方の数は計算しました。
そしてすべての組み合わせに対して、
その組み合わせのF()に並び方の数を乗じました。

00..0~99..9(k桁)のすべての数と、0..9から重複を許してk個並べる場合の数が
一対一に対応しますので、k個の組み合わせごとに並び方の数を計算します。
ある組み合わせの数の和を取りそれに上位桁の数字の和を足すと
各桁の和が求まります。

例によって、abcd..tuというn桁の数なら、
nけた目の0~a-1とn-1けた目の000..0~999..9から
Σ(0<=x<=a\*10^n-1)F(x)を計算し、
aとn-1桁目の0~b-1とn-2桁の00..0~99..9から
Σ(a\*10^n<=x<=a\*10^n+b\*10^(n-1)-1)F(x)を計算し
a+bとn-2桁目の0~c-1とn-3桁の0..0~9..9からΣ(a\*10^n+b\*10^(n-1)<=x<a\*10^n+b\*10^(n-1)+c\*10^(n-2)-1)F(x)
を計算し...というように小さい桁に移動していきます。

collect/5で0~10^i-1の組み合わせを求め、上位桁の情報とともにf/3に送り、
f/3ではその組み合わせの並べ方の数を求めるとともに、
その組み合わせの最初の合算を行い、sums/1に送ってF()を計算し、
F()に並べ方の数を乗じます。
メインですべての結果を合計します。
最後はa+b+c+..+tと0..uからF()を計算して足します。
nけたの00..00~00..09の和は1と計算されますが実際は0なので補正します。

```Julia
#julia ver 1.41,1.53,1.10.3

function solve(n0,a)  
    function ntor(k)
        u = Int(floor(log10(k)))
        ar = []
        x = k
        for i in 1 : u+1
            (x,y) = divrem(x,10)
            push!(ar,y)
        end
        return ar
    end
    function ncr(k,h)            
        mul = 1
        for i in 1 : h
            mul = div(mul*k,i)
            k -= 1
        end
        return mul
    end
    function sums(n)　　　　#ある組み合わせのF(),最初の計算は済んでいる
        if n < 10
            return 1　　　#f/3での計算分　
        else
            sum = 0
            while n > 0
                (n,r) = divrem(n,10)
                sum += r
            end
            return 1 + sums(sum)　　#1はf/3での計算分
        end
    end
    function f(a,b,ar) #ある組み合わせのF()の総計を求める
        ssum=0
        for i in 1:10
            ssum += ar[i] * (i-1)  #ある組み合わせの最初の和
        end
        com = 1　#ある組み合わせの並び替えの数
        k = length(ar0)
        for i in 1 : length(ar)
            com *= ncr(k,ar[i])
            k -= ar[i]
        end
        ret = 0
        for i in 0 : b
            ret += sums(a+i+ssum) * com　
        end
        return ret
    end
  function collect(n,m,a,b,ar)　#組み合わせを作ってf/3に送る
    if m > 9 return 0   
    elseif n == 0
        return f(a,b,ar)  #ar->0~9それぞれの数->0~99..99に対応
    else
        sum = 0
        ar1 = copy(ar)
        ar1[m+1] += 1
        sum += collect(n-1,m,a,b,ar1) + collect(n,m+1,a,b,ar)
        return sum
    end
  end
  msum = 0
  if n0 > 9
    ar0 = ntor(n0)
    x,y = 0,0
    for i in length(ar0)-1 : -1 : 1
      x = pop!(ar0)
      if x != 0
    　  msum += collect(i,0,y,x-1,zeros(Int,10)) #上位桁と0~x*10^(i-1)-1
        y += x　　　#iより上位の桁の数の和
      end
    end
    for i in 0 : ar0[1]
      msum += sums(y+i)   #1のけた
    end
    msum -= 10　　　　#最初の000..00~00..09を補正
  end
  s = a == msum ? "ok " : "no "
  println(s,n0,"->",msum)
 end

solve(9, 0)
solve(10, 1)
solve(20, 12)
solve(48, 48)
solve(100, 136)
solve(149, 200)
solve(5432, 10539)
solve(9876, 19951) 
solve(54321, 113023)
solve(90817263, 207841726)
solve(99003157, 227032799)
```
Prologで書くにあたり、もっと簡便で早そうな方法がないか考えていた時、
和が高々二桁(9と入力値の桁数の積)なのに気づき、和で分類することにしました。
一桁位大きい数を扱えます。
ar[k,i+1]にk桁以下の各桁の和がiになる場合の数を入れます。（配列は1から）
k-1桁目までの和がxでa(0<=a<=9)を加えてy=x+aになったとすると
ar[k,y+1]+=ar[k-1,x+1]　
一桁目の数字はar[1,_]に入れます。
Prologではput/5でリストに入れて行います。

fun/2でput/5を呼び出しています。
NをM+1桁の数としますとdivmod(N,10^M,D,Y)なら
0~D\*10^M-1とD\*10^M~D\*10^M+Yに分け
再帰的に0~D-1とfun(10^M-1,R1)及びDとfun(Y,R2)実行し、
その中でput/5を呼びリストを更新します。
fn/2でF(N)を、gn/4でG(N)を計算します。
最初に合算するときが回数に入っていませんので
2桁以上の時は最後にn0-9を足します。

```Prolog
%SWI-Prolog ver 7.42,7.64,9.25
%:-initialization(start).       %Online Prolog Compiler,ideone.    
%start.

for(-1,_,R,R):-!. %0~Y-1までloopを繰り返し
for(Y,AR1,AR0,R):-put(1,Y,AR1,AR0,AR2),Y1 is Y-1,for(Y1,AR1,AR2,R).

put(K,_,AR1,R,R):-
  K1 is K-1,length(AR1,K1),!.
put(K,J,AR1,R1,R):-  
  K1 is J+K,nth1(K,AR1,X),nth1(K1,R1,Y,AR2),Z is X+Y,
  nth1(K1,R2,Z,AR2),K2 is K+1,put(K2,J,AR1,R2,R).

cal(0,R,R):-!. 　　%各桁の和
cal(N,R1,R):-divmod(N,10,N1,M),R2 is R1+M,cal(N1,R2,R).
　
fn(N,0):-N<10,!.　　%F(n)
fn(N,R):-cal(N,0,R1),fn(R1,R2),R is 1+R2.

fun(N,AR):-
  N<10,!,N1 is N+1,length(L1,N1),maplist(=(1),L1),N2 is 10-N1,
  length(L2,N2),maplist(=(0),L2),append(L1,L2,AR).
fun(N,AR):-
  M is floor(log10(N)),K0 is (M+1)*9+1,length(AR0,K0),maplist(=(0),AR0),    %AR0=zeros(K0)
  X is 10^M,X1 is X-1,divmod(N,X,Y,RE),fun(X1,AR1),fun(RE,AR2),
  Y1 is Y-1,for(Y1,AR1,AR0,AR3),put(1,Y,AR2,AR3,AR).

gn(N,_,R,R):-
  N<11,!.     % G(n)
gn(N,AR,R1,R):-
  N1 is N-1,fn(N1,X),nth1(N,AR,Y),R2 is R1+X*Y,N2 is N-1,gn(N2,AR,R2,R).

solve(N,R):-
  fun(N,AR),length(AR,M),gn(M,AR,0,R1),(N>9->R is R1+N-9;R = R1).

start:-
  f(N,A),(solve(N,R)->(R==A->Str=" ok  ";Str=" no  ");
  write(fail)),write(" "),write(Str),writeln(N->R),fail.
start.


f(7,0).
f(9, 0).
f(10, 1).
f(20, 12).
f(48, 48).
f(100, 136).
f(149, 200).
f(5432, 10539).
f(9876, 19951).
f(54321, 113023).
f(90817263, 207841726).
f(99003157, 227032799).
```
最後に上の方法で再帰でなく下の位から実行するようにしたコードです。
断然早いです。

```Julia
function solve703(n0,a)
    function ntor(k)
      ar = []
      x = k
      while x > 0
          (x,y) = divrem(x,10)
          push!(ar,y)
      end
      return ar
  end
  function fn(x)
    if x < 10
      return 0
    else
      sum = 0
      while x > 0
        (x,r) = divrem(x,10)
        sum += r
      end
      return 1 + fn(sum)
    end
  end
  if n0 < 10
    sum = 0
  else
    ar0 = ntor(n0)　　　　#入力値を数列化
    m = length(ar0)
    k0 = (m+1)*9+1
    ar1 = zeros(Int,m,k0,2)　　
    ar2 = zeros(Int,m,k0,2)
    for i in 1:10
      ar1[1,i,1] = 1
    end
    for i in 1 : ar0[1] + 1
      ar2[1,i,1] = 1
    end
    for i in 2 : m-1
      for j0 in 0 : k0-10
          j = j0 + 1
        for k in 0 : 9            
           x = j0+k
             ar1[i,x+1,1] += ar1[i-1,j,1]  #1~iけた目の0~9すべてを含む
        end
        for k in 0 : ar0[i] - 1　#i桁目をaとすると0~a-1と0~i-1桁の0~9すべて
            x = j0+k
            ar2[i,x+1,1] += ar1[i-1,j,1]
        end
        x = j0 + ar0[i]
          ar2[i,x+1,1] += ar2[i-1,j,1]  #i桁目のaとi-1桁目のar2
       end
    end
    sum=0
    for j0 in 0 : k0 - 10
      j = j0 + 1
      for k in 0 : ar0[m] - 1
           x = j0 + k
           ar1[m,x+1,1] += ar1[m-1,j,1]
      end
      x = ar0[m] + j0
      ar2[m,x+1,1] += ar2[m-1,j,1]
    end
    for i in 9 : k0-1
      sum += fn(i) * (ar1[m,i+1,1] + ar2[m,i+1,1])
    end
    sum += (n0 - 9)
  end
  s = a == sum ? "ok " : "no "
  println(s,n0,"->",sum)
 end
```
追記
　最後のコードが断然速いと書いて投稿したのですが、
  togetter.comを見に行ったら、まさに桁違いの速さのコードが載っていました。
　私は他人のコードを読めませんし、コメントも理解できていないのですが、
　よろしかったら見に行かれてください。