立方根の計算 DouKaku?

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

% awk -f cuberoot.awk
8.0
2.0000000000000000
10.0
2.1544346900318931
100.0
4.6415888336127793

{
	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になってしまいますね・・・修正します

{
	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

#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はパスしました。

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ですが。

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のような小さい数でも対応できる

(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で

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

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

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

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

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

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

※ Yコンビネータは、one-linerなどの式中で再帰呼び出しをする方法。
手続きなしでの実装は#2518に投稿しました。

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()

ニュートン法で求めてみました。

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));
    }
}

範囲を半分ずつ切っていく方針で。

(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)))))

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

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 で。

| 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

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

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

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コンビネータは、one-linerなどの式中で再帰呼び出しをする方法。

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

# だがそれがいい(ぇ

あぁっ。べき乗もつかっちゃいけなかったんだ。
指数とか対数って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)

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

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

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

SICPにあった区間二分法ほとんどそのまま。

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[d2rdlar/+3/dddla-d0>sld>qlrx]sr[13k.1dddd*sddsa10/lrx]sc

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

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

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;
}

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

1	[_1]ss[q]sq[d2rdlar/+3/dddla-d0>sld>qlrx]sr[14kd0=q.1dddddsddsa10/lrx]sc

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

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

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

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

試してみたのですが、cube_root(4.0) だと
1.5874010519681996
と
1.5874010519681994
で振動して、終了しませんね・・・

ああ、すみません。振動する値というのは y です。

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

それと遅延評価ならではの数列生成と収束判定を分離した版もお願いします。

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

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

0.0だけは別判定しないとだめですね。

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

上の「こんにちはこんにちは」は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でもこんにちわこんにちわされました。

あ，ホントですね．

y != y_next

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

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

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

立方根関数、指数関数、対数関数は禁止・・・ということは
数値解析関数なら使っていいに違いないｗ
uniroot()関数は内部でBrent法を使用しています。

> cube_root(100)
[1] "4.641588833613"
> cube_root(10)
[1] "2.154434690032"

see: Brent's method - Wikipedia

1
2
3

cube_root <- function(x){
   sprintf("%.12f", uniroot(function(y)(y*y*y - x), c(0, 1000), tol=1e-12)$root)
}

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

相対誤差も絶対誤差も 1e-12 以下になるようにしてみました。浮動小数点演算の誤差は相対誤差を基本に考えた方が幸せになれるような気がしたので。

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#include<assert.h>

#define EPS 1e-12

double cubicrt(double x)
{
    double y,oldy,z=1.0;
    if(x==0.0){
        return 0.0;
    }
    while(fabs(x)>8){
        x/=8;
        z*=2;
    }
    while(fabs(x)<0.125){
        x*=8;
        z/=2;
    }
    y=x;
    do{
        oldy=y;
        y=2*oldy/3+x/(3*oldy*oldy);
    }while(fabs(y-oldy)>fabs(EPS*oldy) || fabs(y-oldy)>EPS);

    return y*z;
}

void test(double x,double y,double truey,double eps)
{
    assert(fabs(y-truey)<=fabs(eps*y) && fabs(y-truey)<=eps);
    return;
}

int main()
{
    int i=0;
    double x,y;
    srand(time(NULL));

    /* テスト */
    for(i=0;i<1000;i++){
        x=rand()*2000.0/RAND_MAX-500;
        y=cubicrt(x);
        test(x,y*y*y,x,EPS);
    }
    for(i=0;i<1000;i++){
        x=rand()*2e-310/RAND_MAX-1e-310;
        y=cubicrt(x);
        test(x,y*y*y,x,EPS);
    }
    return 0;
}

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

DouKaku?

立方根の計算

Posted feedbacks - Flatten