Last active
September 28, 2022 22:38
-
-
Save demotomohiro/caef89a2a94500e235a6570c77af1d18 to your computer and use it in GitHub Desktop.
Calculate an automorphic number
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# This program calculate an automorphic number. | |
# 何乗しても下N桁が同じになる数(自己同形数 Automorphic number)を計算して出力します。 | |
# | |
# How to compile and run | |
# 1. Install Nim (https://nim-lang.org/install.html) | |
# 2. $ nimble refresh | |
# 3. $ nimble install cligen | |
# 4. $ nimble install https://github.com/demotomohiro/bigints@#add-per-bit-op | |
# 5. $ nim c -d:danger --passC:-flto --passL:"-flto -s" amn.nim | |
# 6. $ ./amn run --numDigits=10000 > result.txt | |
import bigints | |
import std/[times, monotimes, math] | |
func verify(x: BigInt; numDigits: Natural): bool = | |
when false: | |
((x * x) mod (pow(10'bi, numDigits))) == x | |
else: | |
let | |
base = pow(10'bi, ceilDiv(numDigits, 2)) | |
(xhigh, xlow) = divmod(x, base) | |
((xhigh * xlow * 2.initBigInt * base + xlow * xlow) mod pow(10'bi, numDigits)) == x | |
func amn(numDigits: Natural): BigInt = | |
# 5^n*l + 1 ≡ 1 (mod 5^n) | |
# l ∈ [1, 2^n]; 5^n*l + 1 ≡ 0 (mod 2^n) | |
# (5^n*l + 1)*2^n ≡ 2^n (mod 5^n*2^n) | |
# (5^n*l + 1)*2^n*(5^n*l + 1)/(2^n) ≡ 2^n*(5^n*l + 1)/(2^n) (mod 10^n) | |
# (5^n*l + 1)*(5^n*l + 1) ≡ (5^n*l + 1) (mod 10^n) | |
let | |
n5 = pow(5'bi, numDigits) | |
n2 = pow(2'bi, numDigits) | |
mask = n2 - 1.initBigInt | |
l = block: | |
var ret: BigInt | |
when false: | |
for l in countUp(1'bi, n2, 2): | |
if ((n5 * l + 1'bi) and mask) == 0'bi: | |
ret = l | |
break | |
elif false: | |
# l ∈ [1, 2^n]; 5^n*l + 1 ≡ 0 (mod 2^n) | |
# 5^n*l ≡ -1 (mod 2^n) | |
# -5^n*l ≡ 1 (mod 2^n) | |
ret = invmod(-n5, n2) | |
elif false: | |
# l ∈ [1, 2^n]; 5^n*l + 1 ≡ 0 (mod 2^n) | |
# 5^n*l ≡ -1 (mod 2^n) | |
assert n5.testBit(0) | |
ret = 1'bi | |
var | |
tmp = n5 and mask | |
shifted = tmp | |
for i in 1 ..< numDigits: | |
shifted.clearBit(numDigits - 1) | |
shifted = shifted shl 1 | |
if not tmp.testBit(i): | |
tmp += shifted | |
#tmp.clearBit(numDigits) | |
ret.setBit(i) | |
assert (tmp and (n2 - 1'bi)) == n2 - 1'bi | |
else: | |
assert n5.testBit(0) | |
ret = 1'bi | |
var | |
tmp = n5 and mask | |
a = tmp | |
for i in 1 ..< numDigits: | |
tmp = tmp shr 1 | |
a.clearBit(numDigits - i) | |
if not tmp.testBit(0): | |
tmp += a | |
ret.setBit(i) | |
assert ret != 0'bi | |
ret | |
#debugEcho "l: ", l | |
result = n5 * l + 1'bi | |
#assert verify(result, numDigits) | |
proc test = | |
doAssert amn(1) == 6'bi | |
doAssert amn(2) == 76'bi | |
doAssert amn(3) == 376'bi | |
doAssert amn(4) == 9376'bi | |
doAssert amn(5) == 9376'bi | |
doAssert amn(6) == 109376'bi | |
doAssert amn(7) == 7109376'bi | |
doAssert amn(8) == 87109376'bi | |
doAssert amn(9) == 787109376'bi | |
doAssert amn(10) == 1787109376'bi | |
var prev = 1'bi | |
for i in 11 .. 256: | |
let x = amn(i) | |
doAssert verify(x, i) | |
doAssert x >= prev | |
prev = x | |
echo "All tests completed" | |
proc bench(numDigits = 1024 * 256) = | |
let | |
s = cpuTime() | |
sm = getMonoTime() | |
se = epochTime() | |
ret = amn(numDigits) | |
e = cpuTime() | |
em = getMonoTime() | |
ee = epochTime() | |
echo "cpuTime: ", e - s | |
echo "monoTime: ", em - sm | |
echo "epochTime: ", ee - se | |
echo ret and (1 shl 62 - 1).initBigInt | |
proc run(numDigits = 1024 * 256; verify: bool = false) = | |
let | |
s = cpuTime() | |
res = amn(numDigits) | |
e = cpuTime() | |
stderr.write "Time: ", e - s, '\n' | |
if verify: | |
stderr.write "Verifying result.." | |
doAssert verify(res, numDigits) | |
stderr.write "done\n" | |
echo res | |
when isMainModule: | |
import cligen | |
dispatchMulti([run], [test], [bench]) |
This file has been truncated, but you can view the full file.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
$ ./amn bench
cpuTime: 8.135876248999999
monoTime: 8 seconds, 168 milliseconds, 770 microseconds, and 49 nanoseconds
epochTime: 8.16877007484436
0
$ time ./amn run --numDigits=4000000 --verify > result.txt
Time: 1483.781442576
Verifying result..done
real 1828m47.330s
user 1824m34.068s
sys 0m2.456s