Skip to content

Instantly share code, notes, and snippets.

@demotomohiro

demotomohiro/amn.nim

Last active September 28, 2022 22:38
Show Gist options
  • Save demotomohiro/caef89a2a94500e235a6570c77af1d18 to your computer and use it in GitHub Desktop.
Save demotomohiro/caef89a2a94500e235a6570c77af1d18 to your computer and use it in GitHub Desktop.
Download ZIP
Calculate an automorphic number
Raw
amn.nim
# 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])
Raw
amn4000000digits.txt
This file has been truncated, but you can view the full file.
@demotomohiro
Copy link
Author

$ ./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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment