バッチファイルでProject Euler(54)

Problem 23


6が完全数で、12以上の6の倍数が過剰数であることはすぐにわかります。なので、6の剰余で分類して考えるとよいです。
例えば6の剰余が2(以下、r = 2などと表記する)なら20、r = 4なら40になります。しかし、奇数はすなわち2を因数に含めないので小さな過剰数はありません。最小の過剰数は3^3*5*7 = 945です。このように、奇数の小さな過剰数は3を含むのは明らかなので、r = 3となります。計算してみると、そのほかの奇数は28123以下にはありません。
さて、例えばr = 4なら40がその最小の過剰数なので、52以上ならば2つの過剰数の和で表されます。なので、46以下でr = 0を含まない和を考えます。ここで足してr = 4となるrのペアを考えると、(1, 3), (2, 2), (5, 5)となります。これらのうち足して46以下になるのは20 + 20 = 40のみです。このことからr = 4で2つの過剰数の和にならない自然数は、4, 10, 16, 22, 28, 34, 46となります。他のrでも同様に計算することができます。
計算するとき、当然forで回すと速くなりますが、何番目の過剰数まで見ればよいかをまず調べます。そうすればforを使いやすくなります。何番目かを調べるのには二分探索を使っています。
それでも2000までで1分以上かかりました。細かく区切って計算すれば配列が小さくなり速くなるはずですが、もう撤退します。

@echo off

setlocal enabledelayedexpansion
call :start_time
set /a N = 28123
set /a N = 2000
call :sqrt %N%
set /a L = %ERRORLEVEL%
set /a M = 500

rem // calculate prime numbers
set /a ps.size = 1
set /a ps_0 = 2
for /L %%k in (3, 2, %L%) do (
    call :is_prime %%k
    if !ERRORLEVEL! == 1 call :push_vector ps %%k
)

rem // get abundant numbers
set /a abn.size = 0
call :fill abn_0 1 12
for /L %%b in (2, %M%, %N%) do (
    set /a e = %%b + %M%
    call :sieve a b %%b !e!
    set /a e_1 = !e! - 1
    for /L %%n in (%%b, 1, !e_1!) do (
        set /a r = %%n %% 6
        if !r! NEQ 0 (
            set /a k = %%n - %%b
            set /a sd = "b_!k!" - %%n
            if !sd! GTR %%n (
                set /a size = "abn_!r!.size"
                set /a abn_!r!_!size! = %%n
                set /a abn_!r!.size += 1
            )
        )
    )
)

call :check_time

set /a s = 0
for /L %%r in (0, 1, 5) do (
    call :get_r_pairs rs1 rs2 %%r
    set /a abn_r_size = "abn_%%r.size"
    if !abn_r_size! == 0 (
        set /a mod = %N% %% 6
        if %%r LEQ !mod! (
            set /a front = %N% - !mod! + %%r - 6
        ) else (
            set /a front = %N% - !mod! + %%r - 12
        )
    ) else (
        set /a front = "abn_%%r_0"
    )
    set /a b.size = !front! / 6 + 2
    
    call :fill b !b.size! 0
    set /a max_index = "rs1.size" - 1
    for /L %%k in (0, 1, !max_index!) do (
        set /a r1 = "rs1_%%k"
        set /a r2 = "rs2_%%k"
        call :set_sum_r b !r1! !r2!
    )
    
    set /a max_b_index = !b.size! - 1
    for /L %%k in (0, 1, !max_b_index!) do (
        set /a e = "b_%%k"
        if !e! == 0 set /a s += %%k * 6 + %%r
    )
)
echo %s%
call :check_time
exit /b 0

rem // call :set_sum_r 配列 r1 r2
rem // 6の剰余がそれぞれr1とr2の過剰数の和となる数を求め
rem // それをインデックスとする配列の要素に1をセットする
:set_sum_r
    set /a size = "%1.size"
    set /a r = (%2 + %3) %% 6
    set /a size1 = "abn_%2.size"
    set /a size2 = "abn_%3.size"
    if %size1% == 0 exit /b 0
    if %size2% == 0 exit /b 0
    
    rem // 各々の剰余で最小の過剰数
    set /a front1 = "abn_%2_0"
    set /a front2 = "abn_%3_0"
    set /a front3 = "abn_%r%_0"
    
    rem // 走査しなければならない過剰数の最大のインデックス
    set /a max_num1 = %N% - %front2% - 12
    set /a max_num2 = %N% - %front1% - 12
    call :bin_search abn_%2 %max_num1%
    set /a max_index1 = %ERRORLEVEL%
    call :bin_search abn_%3 %max_num2%
    set /a max_index2 = %ERRORLEVEL%
    
    for /L %%i in (0, 1, %max_index1%) do (
        set /a num1 = "abn_%2_%%i"
        for /L %%j in (0, 1, %max_index2%) do (
            set /a num2 = "abn_%3_%%j"
            set /a num = !num1! + !num2!
            set /a index = num / 6
            if !index! LEQ %size% set /a %1_!index! = 1
        )
    )
    exit /b 0

rem // call :get_r_pairs rs1 rs2 r
rem // 足して6の剰余がrになるペアを列挙する
rem // ただし0は含まない
:get_r_pairs
    set /a k = 0
    set /a end = %3 / 2
    for /L %%r in (1, 1, %end%) do (
        set /a %1_!k! = %%r
        set /a %2_!k! = %3 - %%r
        set /a k += 1
    )
    
    set /a begin = %3 + 1
    set /a end = (%3 + 6) / 2
    for /L %%r in (%begin%, 1, %end%) do (
        set /a %1_!k! = %%r
        set /a %2_!k! = %3 + 6 - %%r
        set /a k += 1
    )
    
    set /a %1.size = %k%
    set /a %2.size = %k%
    exit /b 0

rem // bin_search 配列 値
rem // 与えられた値を超えない最大の要素のインデックスを返す
:bin_search
    setlocal
    set /a b = 0
    set /a e = "%1.size"
    :loop_bin_search
        set /a m = (%b% + %e%) / 2
        if %b% == %m% exit /b %m%
        set /a vm = "%1_%m%"
        if %vm% LEQ %2 (
            set /a b = %m%
        ) else (
            set /a e = %m%
        )
        goto :loop_bin_search

:sqrt
    setlocal
    set /a r = %1
    :loop_sqrt
        set /a r1 = (%r% + %1 / %r%) / 2
        if %r1% GEQ %r% exit /b %r%
        set /a r = %r1%
        goto :loop_sqrt

:is_prime
    set /a max_index = "ps.size" - 1
    for /L %%k in (0, 1, %max_index%) do (
        set /a r = %1 %% "ps_%%k"
        if !r! == 0 exit /b 0
    )
    exit /b 1

:sieve
    set /a _e = %4 - 1
    call :range %1 %3 %4
    call :fill %2 %M% 1
    set /a k = 0
    set /a p = 2
    :loop_sieve
        set /a _b = (%3 + %p% - 1) / %p% * %p%
        if %_b% == %p% set /a _b += %p%
        for /L %%k in (%_b%, %p%, %_e%) do (
            set /a i = %%k - %3
            set /a a_!i! /= %p%
            set /a _s = 1 + %p%
            set /a r = "a_!i!" %% %p%
            if !r! == 0 call :div_exp
            set /a b_!i! *= !_s!
        )
        set /a k += 1
        set /a p = "ps_%k%"
        if %p% GTR 0 (
            set /a p_sq = %p% * %p%
            if !p_sq! LSS %4 goto :loop_sieve
        )
    
    set /a M_1 = %M% - 1
    for /L %%k in (0, 1, %M_1%) do (
        set /a ak = "%1_%%k"
        if !ak! GTR 1 set /a %2_%%k *= !ak! + 1
    )
    exit /b 0

:div_exp
    :loop_div_exp
        set /a a_!i! /= %p%
        set /a _s = !_s! * %p% + 1
        set /a r = a_!i! %% %p%
        if %r% NEQ 0 exit /b 0
        goto :loop_div_exp

:sum_divs
    setlocal
    set /a s = 1
    call :sqrt %1
    set /a max = %ERRORLEVEL%
    for /L %%d in (2, 1, %max%) do (
        set /a r = %1 %% %%d
        if !r! == 0 (
            set /a s += %%d
            set /a q = %1 / %%d
            if %%d NEQ !q! set /a s += !q!
        )
    )
    exit /b %s%

:divs
    :loop_divs
        set /a q /= !d!
        set /a _s = !_s! * !d! + 1
        set /a r = !q! %% !d!
        if !r! NEQ 0 exit /b 0
        goto :loop_divs

:range
    if "%3" == "" call :range %1 0 %2 1 & exit /b 0
    if "%4" == "" call :range %1 %2 %3 1 & exit /b 0
    
    set /a max_index = (%3 - %2) / %4 - 1
    set /a %1.size = %max_index% + 1
    for /L %%i in (0, 1, %max_index%) do (
        set /a %1_%%i = %2 + %%i * %4
    )
    exit /b 0

:fill
    set /a %1.size = %2
    set /a max_index = %2 - 1
    for /L %%k in (0, 1, %max_index%) do set /a %1_%%k = %3
    exit /b 0

:push_vector
    set /a size = "%1.size"
    set /a %1_%size% = %2
    set /a %1.size += 1
    exit /b 0

:print_vector
    setlocal
    set /a size = "%1.size"
    if %size% == 0 echo %2 & exit /b 0
    set /a e = "%1_0"
    set s=%e%
    set /a k = 1
    :loop_print_vector
        if %k% == %size% echo %s% & exit /b 0
        set /a e = "%1_%k%"
        set s=%s% %e%
        set /a k += 1
        goto :loop_print_vector

:start_time
    call :get_time
    set /a __t0 = %ERRORLEVEL%
    exit /b 0

:check_time
    setlocal
    call :get_time
    set /a t = %ERRORLEVEL% - %__t0%
    if %t% LSS 10 (
        echo 0.0%t%s
    ) else (
        if %t% LSS 100 (
            echo 0.%t%s
        ) else (
            echo %t:~0,-2%.%t:~-2%s
        )
    )
    exit /b 0

:get_time
    setlocal
    set t=%TIME%
    set /a h = %t:~0,2%
    set /a m = 1%t:~3,2% %% 100
    set /a s = 1%t:~6,2% %% 100
    set /a ss = 1%t:~-2% %% 100
    set /a ret = ((%h% * 60 + %m%) * 60 + %s%) * 100 + %ss%
    exit /b %ret%