立方根の計算 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;
}

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

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

先日チャットで議論した結果、親記事が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を渡しています。

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

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

単純なニュートン法。1000までならこれで十分。

my $d = 10 ** -12;

sub check($$$){
  my ($x, $y, $d) = @_;
  return abs($y - $x ** 3) < $d;
}

sub improve($$){
  my ($x, $a) = @_;
  return (2 * $x ** 3 + $a) / (3 * $x ** 2);
}

sub cube_root($){
  my $a = $_[0];
  return 0 if $a == 0;
  my $p = 1;
  my $q;
  until(check($p, $a, $d)){
    $q = $p;
    $p = improve($p, $a);
  }
  $p;
}

もっと精度が必要な場合は、Math::BigFloatを使うとよいと思います。(check, improveはまったく同じです。)

use Math::BigFloat;

my $d = Math::BigFloat->new("1E-12");

sub check($$$){
  my ($x, $y, $d) = @_;
  return abs($y - $x ** 3) < $d;
}

sub improve($$){
  my ($x, $a) = @_;
  return (2 * $x ** 3 + $a) / (3 * $x ** 2);
}

sub cube_root($){
  my $a = Math::BigFloat->new($_[0]);
  my $p = Math::BigFloat->bone();
  my $q;
  until(check($p, $a, $d)){
    $q = $p;
    $p = improve($p, $a);
  }
  $p->ffround(-12);
}

gcc version 3.4.2 (mingw-special)で動作確認しました。
終了条件のチェックはこれでいいのかな？

#include <cmath>

long double cube_root(long double x, long double min = 0.0, long double max = 1000.0)
{
    long double mid = (min + max) / 2.0;
    if((mid == min) || (mid == max))
        return mid;
    return ((std::pow(mid, 3) - x) > 0) ?  cube_root(x, min, mid) : cube_root(x, mid, max);
}

/*
#include <cstdlib>
#include <iomanip>
#include <iostream>

int main(int argc, char * argv[])
{
    long double x = std::atof(argv[1]);
    std::cout << std::setprecision(15) << cube_root(x) << std::endl;
    return EXIT_SUCCESS;
}
*/

ニュートン法で。

import std.stdio;
import std.math;

real cube_root(real x){
    real y0 = 3;
    while(true){
        real y1 = (2 * y0 + x / (y0 * y0)) / 3;
        if(abs(y1 - y0) < 1e-14) return y1;
        y0 = y1;
    }
}

void main(){
    foreach(r; [cube_root(10.0), cube_root(100.0)]){
        writefln("%.13f ^ 3 = %.13f", r, pow(r, 3));
          //=> 2.1544346900319 ^ 3 = 10.0000000000000
          //=> 4.6415888336128 ^ 3 = 100.0000000000000
    }
}

ニュートン法にしました。

public static double CubeRoot(double x)
{
  double a0 = (double)(-x);
  double[] a = { a0, 0.0, 0.0, 1.0 };
  return SolveAlgebraicEquationByNewtonMethod(a, 1e-12);
}

delegate double MathFunc(double x);

static double SolveAlgebraicEquationByNewtonMethod(double[] a, double delta)
{
  MathFunc f = CreateFunc(a);
  MathFunc fPrime = CreateFunc(Derivative(a));
  double x = 1.0;
  while (NearlyEquals(fPrime(x), 0.0, double.Epsilon))
    x += 1.0;
  double betterX = 0;
  while (true)
  {
    betterX = Newton(f, fPrime, x);
    if (NearlyEquals(x, betterX, delta))
      break;
    x = betterX;
  }
  return betterX;
}

static MathFunc CreateFunc(double[] a)
{
  return delegate(double x)
  {
    double fx = 0.0;
    for (int i = 0; i < a.Length; i++)
    {
      double d = a[i];
      if (NearlyEquals(d, 0.0, double.Epsilon))
        continue;
      for (int j = 0; j < i; j++)
        d *= x;
      fx += d;
    }
    return fx;
  };
}

static double[] Derivative(double[] a)
{
  double[] aPrime = new double[a.Length - 1];
  for (int i = 1; i < a.Length; i++)
  {
    aPrime[i - 1] = a[i] * (double)i;
  }
  return aPrime;
}

static double Newton(MathFunc f, MathFunc fPrime, double x)
{
  return x - (f(x) / fPrime(x));
}

static bool NearlyEquals(double d1, double d2, double delta)
{
  return (((d1 - d2) <= delta) && ((d2 - d1) <= delta));
}

#2868をなでしこに移植しました

●立方根(xの)
    yとは実数
    y_nextとは実数=x 
    (ABS(y-y_next)>=1e-12)の間
        y=y_next
        y_next=y-(y*y*y-x)/(3*y*y)
    yで戻る
4の立方根を表示
4^(1/3)を表示

自動微分ライブラリを使ってみました。 HackageDBにあります。

このライブラリのderiv関数を使うと導関数を「数値的に」ではなく「記号的に」求める事ができます。つまりf'(x) = {f(x + h) - f(x)} / hという近似を利用するのではなく、 (x^3)' = 3x^2 の様に直接導関数を生成します。

円周率と自然対数も求めてみました。

import Data.Number.Dif

newton f init eps =
  let f' = deriv f
      next = init - (unDif f init) / (f' init)
  in if abs (next - init) < eps
    then init
    else newton f next eps


-- aの立方根は x^3 - a = 0の解
cbrt a = newton (\x -> x**3 - a) 1 1e-13

main = do
  putStrLn $ "cbrt 10   = " ++ show (cbrt 10)
  putStrLn $ "cbrt 100  = " ++ show (cbrt 100)
  putStrLn $ "cbrt 1000 = " ++ show (cbrt 1000)
  putStrLn $ "pi        = " ++ show (newton sin 3 1e-13)
  putStrLn $ "e         = " ++ show (newton (\x -> log x - 1) 2 1e-13)

{-
*Main> :main
cbrt 10   = 2.154434690031884
cbrt 100  = 4.641588833612779
cbrt 1000 = 10.000000000000004
pi        = 3.141592653589793
e         = 2.7182818284590455
-}

1> c(cube_root).
{ok,cube_root}
2> io:format("~15.13f~n", [cube_root:cube_root(10.0)]).
2.1544346900319
ok
3> io:format("~15.13f~n", [cube_root:cube_root(100.0)]).
4.6415888336128
ok

-module(cube_root).
-export([cube_root/1]).

cube_root(Num) -> cube_root(Num, 0.0, 1000.0).
cube_root(Num, Min, Max) ->
    Err = 0.0000000000001,
    Mid = (Min + Max) / 2.0,
    if
        ((Mid - Err) < Min) and (Max < (Mid + Err)) -> Mid;
        ((Mid * Mid * Mid) - Num) > 0 -> cube_root(Num, Min, Mid);
        true -> cube_root(Num, Mid, Max)
    end.

電卓でルートキーを使って立方根を求める方法を試そうと思って。

!立方根の計算
FUNCTION cube_root(x)
   LET  y = x
   LET  E = 1e-13
   DO
      LET  y0 = y
      LET  y = SQR(SQR(y * x))
   LOOP  WHILE y0 - y > E
   LET  cube_root = y
END FUNCTION

FOR i = 1 TO 1000
   PRINT  i;cube_root(i)
NEXT i
END

整数部分は配列(リスト)のインデックスで、小数部分は二分探索的なアプローチで求めています。
三乗根がピッタリ整数の場合は、整数で答えを出すようにしたので、ちょっと長くなりました。

実行例:
arc> (search-cube-root 10.0)
2.1544346900318487
arc> (expt (search-cube-root 10.0) 3)
9.999999999999513
arc> (search-cube-root 100.0)
4.641588833612786
arc> (expt (search-cube-root 100.0) 3)
100.00000000000048
arc> (search-cube-root 125.0)
5

(= cube-list (map [expt _ 3] (range 1 10)))
(= acc 0.000000000001)

(def search-real (cr of rn)
  (let gosa (- rn (expt cr 3))
    (if (< (abs gosa) acc)
        cr
        (if (< 0 gosa)
            (search-real (+ cr of) (/ of 2) rn)
            (search-real (- cr of) (/ of 2) rn)))))

(def search-cube-root (rn)
  (if (is (type rn) 'int)
      (let i (trunc rn)
        (if (mem i cube-list)
            (+ (pos i cube-list) 1)
            (search-real (+ (pos i (sort < (cons i cube-list))) 0.5) 0.25 rn)))
      (search-real (+ (pos rn (sort < (cons rn cube-list))) 0.5) 0.25 rn)))

指数関数使用禁止なのを忘れとった。
1 行目は
(= cube-list (map [* _ _ _] (range 1 10)))

5 行目は
  (let gosa (- rn (* cr cr cr))

ですね。

そのまま。関係ないデスが、left、rightの値を変えれば精度も変更可能です。

def cubeRoot(num, left=2, right=-16){
    def ans = 0
    left.downto(right){ rank ->
        def findnum = (9..1).find{ (ans + it*(10**rank))**3 <= num }?:0
        ans += findnum*(10**rank)
    }
    ans
}

println cubeRoot(10)
println cubeRoot(100)

普通に

<?php

function cube_root($num)
{
    $a = $num;
    
    while(true)
    {
        $b = ($num/($a*$a) + $a*2) / 3;
        ++$i;
        $diff = abs($b - $a);
        if($diff < 1e-12)break;
        $a = $b;
    }
    return $b;
}

//test code
$ret = cube_root(10);
printf("%.16f\n%.16f\n", $ret, $ret*$ret*$ret);
$ret = cube_root(100);
printf("%.16f\n%.16f\n", $ret, $ret*$ret*$ret);

DouKaku?

立方根の計算

Posted feedbacks - Flatten

Index

Feed

Other

Link

Pathtraq