DouKaku?beta

challenge 立方根の計算

xは0以上1000未満の実数です。 y * y * y = xになるような実数y(立方根)を小数点以下12桁以上の正確さで 求める関数cube_rootを作って下さい。

ただし、このお題の趣旨は実数区間での探索なので、 立方根関数があっても使ってはいけません。 指数関数と対数関数も禁止します。

Pythonで表現した入出力の例： 

>>> cube_root(10.0)
2.1544346900318834
>>> _ ** 3
9.9999999999999947
>>> cube_root(100.0)
4.6415888336127793
>>> _ ** 3
100.00000000000003
post comment

Posted feedbacks - Python

print文でフォーマット表示した時は、少し違う出力になるけど、
doctestはパスしました。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
def cube_root(num, guess=1.0):
  """
  >>> cube_root(10.0)
  2.1544346900318838
  >>> _ ** 3
  10.000000000000002
  >>> cube_root(100.0)
  4.6415888336127793
  >>> _ ** 3
  100.00000000000003
  """
  improve = lambda x: ((x*2)+(num/x**2)) / 3.0
  good_enough = lambda x: abs(x**3-num) < 0.00000000001

  while True:
    if good_enough(guess):
      return guess
    guess = improve(guess)


手続き版は最初に投稿したので、

・再帰
・Yコンビネータ
・アキュームレータ

を使ったものを投稿します。
アルゴリズムは全部同じ、Newton Raphson's method です。

※ Yコンビネータは、one-linerなどの式中で再帰呼び出しをする方法。
手続きなしでの実装は#2518に投稿しました。
 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
def cube_root(x,y=1.0):
  if abs(y**3-x) < 1e-13:
    return y
  else:
    return cube_root(x, (((x/y**2)+(y*2)) / 3.0))



def cube_root(x):
  def Y(f):
    f_ = lambda v: f(f_,v)
    return f_

  good_enough = lambda y: abs(y**3-x) < 1e-13
  improve = lambda y: ((x/y**2)+(y*2)) / 3.0
  return Y(lambda f,y: y if good_enough(y) else f(improve(y)))(1.0)


def cube_root(x):
  def accgen(func, n):
    while True:
      yield n 
      n = func(n)

  good_enough = lambda y: abs(y**3-x)<1e-13
  improve = lambda y: (((x/y**2)+(y*2))/3.0)
  return (y for y in accgen(improve,1.0) if good_enough(y)).next()


なんとか法とかよく知らないので、再帰クイックソート的に
範囲を狭めてやってみました。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def cube_root(n):
  def f(a,b):
    if abs(n-(b**3)) < 0.000000000001: return b
    if max(n, ((a+b)/2)**3) == n:
      return f((a+b)/2, b)
    else:
      return f(a, (a+b)/2)
  return f(0.0, float(n))

i = cube_root(10)
print '%.13f ** 3 = %.13f' % (i, i**3)
i = cube_root(100)
print '%.13f ** 3 = %.13f' % (i, i**3)


Edu-sigメーリングリストのアーカイブより（参考リンク）

Halley法です。高階関数にしてみました。

see: [Edu-sig] Halley's Method for nth root of P 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
def gen_nth_root_func(nth, max_depth=-1, initial_guess=1.0, accuracy=1e-12):
  is_good_enough = lambda x,y: abs(y**nth-x) < accuracy
  improve = lambda x,y: ((nth+1)*x*y + (nth-1)*(y**(nth+1))) / ((nth-1)*x + (nth+1)*y**nth)

  def nth_root(num, depth=max_depth, guess=initial_guess):
    if num == 0:
        return num
    elif is_good_enough(num,guess) or depth == 0:
        return guess
    else:
        return nth_root(num, depth-1, improve(num,guess))
  return nth_root

square_root = gen_nth_root_func(nth=2) 
cube_root = gen_nth_root_func(nth=3, max_depth=-1, accuracy=1e-13)

むむむ、平方根だと何かとかぶるかと思って立方根にしたのですが、それでもかぶってましたか…orz

この問題は「線形に探索すると明らかに時間がかかりすぎるけども、二分探索なら十分な速度で答えが出る」というラインを狙って出題したので、下のようなアルゴリズムでも十分答えが出ます。 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
def cube_root(x):
	c = 8.0
	d = 4.0
	for i in range(100):
		v = c * c * c
		if v < x:
			c += d
			d /= 2
		elif v > x:
			c -= d
			d /= 2
		else:
			break
	return c


ニュートン・ラプソン法で書いてみる．
1
2
3
4
5
6
7
def cube_root(x):
  y = 0.0 
  y_next = x 
  while y != y_next:
    y = y_next
    y_next = y - (y * y * y - x) / (3.0 * y * y)
  return y

あ，ホントですね． 

y != y_next

を収束条件にしてるせいみたいです． テストの手抜きがバレますね...

もう少しまじめに収束判定した修正版を貼っておきます． テストももう少しまじめにしておいたので大丈夫なはず..． 

1
2
3
4
5
6
7
def cube_root(x):
  y = 0.0 
  y_next = x 
  while abs(y - y_next) >= 1e-12:
    y = y_next
    y_next = y - (y * y * y - x) / (3.0 * y * y)
  return y

サーバ移転で忙しくてほったらかしになっていました。

先日チャットで議論した結果、親記事がYコンビネータと呼んでいるコードは「不動点関数として有名なYコンビネータ」とは別物だという結論に達しました。というわけでYコンビネータを使ったバージョンを載せておきます。2行目と3行目がYコンビネータの定義です。

参考文献：Fixed point combinator - Wikipedia, the free encyclopedia , Pythonで階乗を求める(Yコンビネータ編)

簡単に説明すると、Yコンビネータというのは与えられた関数の不動点(与えられた関数をfとするとf(x) = xになるようなx)を求める関数の一つです。

で、Yコンビネータに与えられる関数は上のコードではどこにあるのかというと、 4行目から8行目の(1.0)の直前までです。 名前が付いていなくて説明しにくいのでこの関数を仮にstepと呼ぶことにします。 ここまでのところはY(step)という形になっているわけです。 で、Yが不動点を求める関数なのでY(step)は関数stepの不動点になります。 つまりY(step) = step(Y(step))です。 ここでstepの定義に目を移してみると、これは関数fを受け取って 「yを受け取ってそのまま返したり、fを呼んだ結果を返したりする関数」 を返す関数です。 step(Y(step))は「yを受け取ってそのまま返したり、Y(step)を呼んだ結果を返したりする関数」 を返すわけです。 Y(step)を呼ぶということはstep(Y(step))を呼ぶということと同じなので、どこまで行っても呼ばれるのはY(step)です。 こうやって再帰呼び出しを(関数に名前をつけることなく)実現しているわけです。

そして最後の(1.0)でその関数に1.0を渡しています。

とまぁ、長々と説明しましたが、ここには僕より詳しい人がたくさんいると思うので、間違っていたらつっこんでもらえるかなと半分期待しつつ自分の理解を確かめるつもりで書いてみました。 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
>>> cube_root = lambda x:(
      (lambda f:((lambda g: f(lambda x: g(g)(x)))
                  (lambda g: f(lambda x: g(g)(x)))))
      (lambda f:
           lambda y:
               y
               if abs(y ** 3 - x) < 1e-13
               else f(((x/y**2)+(y*2)) / 3.0))(1.0))

>>> cube_root(10.0)
2.1544346900318838
>>> cube_root(100.0)
4.6415888336127793

Index

Feed

Other

Link

Pathtraq

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