DouKaku?beta

challenge exp(pi * sqrt(n))が整数に近くなるnを探す

1以上200未満の整数nのうち、 exp(pi * sqrt(n))がほとんど整数であるようなnを求めるコードを書いてください。 なお、expは底がeである指数関数 - Wikipedia、 piは円周率、sqrtは平方根です。また「ほとんど整数である」とは 整数からプラスマイナス0.0001の範囲にあることとします。

Pythonで34行のスクリプトを書いて得られた出力の例が下のようになります。 

37 199148647.999978
58 24591257752.000000
67 147197952743.999999
163 262537412640768744.000000
この問題は光成さんに教えて頂いた e^{π*sqrt{163}}≒26253741640768744 が元になっています。ご協力ありがとうございました。
post comment

Posted feedbacks - Python

Decimalを使って多倍長計算。sqrtは組み込みのものを使い、expは安直にマクローリン展開で収束するまで回しています。psycoを使ってますが、それでもかなり遅いです。

see: Pythonで10進数計算したり数値を丸めたりする 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
from decimal import *
import itertools

def decimal_exp(x):
    curr, bunsi, bunbo = Decimal("1.0"), Decimal("1.0"), 1
    for k in itertools.count(1):
        bunsi *= x
        bunbo *= k
        next = curr + bunsi / bunbo
        if curr == next:
            return curr
        curr = next

def main():
    pi = Decimal("3.141592653589793238462643383279")
    for n in xrange(1, 200):
        x = decimal_exp(pi * Decimal(n).sqrt()) # result
        i = x.quantize(Decimal("1.0"), rounding=ROUND_HALF_UP) # nearest integer
        if (x - i).quantize(Decimal("0.00001"), rounding=ROUND_DOWN) == 0:
            print n, x

if __name__ == '__main__':
    try:
        import psyco
    except ImportError:
        pass
    else:
        psyco.full()
    main()


上記のコードをctypesを使ってpythonで実装してみました。

Cのヘッダファイル内で#ifdefやら#defineが多用されているので、もしかしたら
環境依存になっているかもしれません。

呼び出し毎にインスタンスを生成するexp関数と、静的なインスタンスを持つ
exp2関数の速度を比較すると、exp2の方が20%ほどはやいようです。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
from ctypes import *

gmp = cdll.LoadLibrary('libgmp.so')

for s in ['set_default_prec', 'init', 'clear', 'set', 'set_ui', 'set_str', 'mul'
, 'mul_ui', 'div', 'add', 'cmp', 'sqrt', 'printf', 'snprintf']:
  try:
    gmp.__setattr__(s, gmp.__getattr__('__gmpf_'+s))
  except:
    gmp.__setattr__(s, gmp.__getattr__('__gmp_'+s))

def mpf():
  class mpf_t(Structure):
    _fields_ = [("_mp_prec", c_int), ("_mp_size", c_int), ("_mp_exp", c_int), ("
_mp_d", c_void_p)]

    def __init__(self):
      Structure.__init__(self)
      gmp.init(byref(self))

    def __del__(self):
      try:
        gmp.clear(byref(self))
      except:
        pass

  return pointer(mpf_t())

def exp(rop, op):
  curr, bunsi, bunbo, next = mpf(), mpf(), mpf(), mpf()

  gmp.set_ui(curr, c_ulong(1))
  gmp.set_ui(bunsi, c_ulong(1))
  gmp.set_ui(bunbo, c_ulong(1))

  k = 1
  while True:
    gmp.mul(bunsi, bunsi, op)
    gmp.mul_ui(bunbo, bunbo, c_ulong(k))
    gmp.div(next, bunsi, bunbo)
    gmp.add(next, next, curr)
    if gmp.cmp(curr, next) == 0:
      gmp.set(rop, curr)
      return
    gmp.set(curr, next)
    k += 1

def exp2(rop, op):
  try:
    curr, bunsi, bunbo, next = exp2.curr, exp2.bunsi, exp2.bunbo, exp2.next
  except:
    exp2.curr, exp2.bunsi, exp2.bunbo, exp2.next = mpf(), mpf(), mpf(), mpf()
    curr, bunsi, bunbo, next = exp2.curr, exp2.bunsi, exp2.bunbo, exp2.next

  gmp.set_ui(curr, c_ulong(1))
  gmp.set_ui(bunsi, c_ulong(1))
  gmp.set_ui(bunbo, c_ulong(1))

  k = 1
  while True:
    gmp.mul(bunsi, bunsi, op)
    gmp.mul_ui(bunbo, bunbo, c_ulong(k))
    gmp.div(next, bunsi, bunbo)
    gmp.add(next, next, curr)
    if gmp.cmp(curr, next) == 0:
      gmp.set(rop, curr)
      return
    gmp.set(curr, next)
    k += 1

def frac_f(x):
  buf = c_buffer(256)
  gmp.snprintf(buf, len(buf)-1, "%.Ff", x)
  return float(buf.value[buf.value.index('.'):])

def main():
  from time import time
  t = time()

  gmp.set_default_prec(128)

  x, pi = mpf(), mpf()

  gmp.set_str(pi, "3.141592653589793238462643383279502884197", 10)

  for i in range(1, 200):
    gmp.set_ui(x, c_ulong(i))
    gmp.sqrt(x, x)
    gmp.mul(x, x, pi)
#    exp(x, x)
    exp2(x, x)
    f = frac_f(x)
    if f < 0.0001 or (1-f) < 0.0001:
      gmp.printf("%d %.6Ff\n", i, x)

  print time() - t

if __name__ == '__main__':
  main()

#3130の匿名さんのアイデアを拝借して組みなおしてみたところ、2.6倍早くなりました。

see: これを組んでいて、Python2.5のバグを踏んでしまった（汗） 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
from decimal import *
import itertools

def divmod(x, y): # workaround for python2.5 decimal mod bug.
    div = x // y
    mod = x - div * y
    return div, mod

def decimal_exp(x):
    div, mod = divmod(x, Decimal("2.3025850929940456840179914546844"))
    curr, bunsi, bunbo = Decimal("1.0"), Decimal("1.0"), 1
    for k in itertools.count(1):
        bunsi *= mod
        bunbo *= k
        next = curr + bunsi / bunbo
        if curr == next:
            return (10 ** div) * curr
        curr = next

def main():
    pi = Decimal("3.141592653589793238462643383279")
    for n in xrange(1, 200):
        x = decimal_exp(pi * Decimal(n).sqrt()) # result
        i = x.quantize(Decimal("1.0"), rounding=ROUND_HALF_UP) # nearest integer
        if (x - i).quantize(Decimal("0.0001"), rounding=ROUND_DOWN) == 0:
            print n, x

if __name__ == '__main__':
    try:
        import psyco
    except ImportError:
        pass
    else:
        psyco.full()
    main()

Index

Feed

Other

Link

Pathtraq

loading...
Creative Commons License このサイトの内容は、 クリエイティブ・コモンズ・ライセンスの下でライセンスされています。 [詳細]