平方根の計算、再び

30分プログラム、その652。平方根の求め方の方法で、平方根の計算方法を試してみました。
「なんで奇数を引くんだよ。わけ解んねーよ」と騒いでたら、id:Roccoさん、別名@hi_saitoさんに、平方根の求め方というサイトを紹介してもらいました。

というわけで、さっそくコードに落してみよう。
でも、なかなかややこしくて、まだ整数の範囲しか扱えません。でも、前回の平方根の計算 - みずぴー日記よりは高速だと思います。

使い方

gosh> (sqrt 529)
(2 3)

ソースコード

#! /opt/local/bin/gosh
;; -*- mode:scheme; coding:utf-8 -*-
;;
;; sqrt.scm -
;;
;; Copyright(C) 2009 by mzp
;; Author: MIZUNO Hiroki / mzpppp at gmail dot com
;; http://howdyworld.org
;;
;; Timestamp: 2009/09/02 21:21:30
;;
;; This program is free software; you can redistribute it and/or
;; modify it under MIT Lincence.
;;
(use gauche.collection)

(define (prev odd)
  (- odd 2))

(define (next odd)
  (+ odd 2))

(define (sqrt-digit n i)
  (let loop [(n n)
	     (i i)
	     (count 0)
	     ]
    (if (< n i)
	(values count n (prev i))
	(loop (- n i) (next i)  (+ 1 count)))))

(define (split-digit n unit)
  (let loop [(n n)
	     (xs '())]
    (if (< n unit)
	(cons n xs)
	(receive (q r)
		 (quotient&remainder n unit)
		 (loop q (cons r xs))))))

(define (sqrt n)
  (receive (n _)
	   (map-accum (lambda (n accum)
			(receive (ans r minus)
				 (sqrt-digit (+ n (* 100 (car accum)))
					     (cdr accum))
				 (values ans
					 (cons r (+ (* (+ minus 1) 10) 1)))))
		      (cons 0 1)
		      (split-digit n 100))
	   n))

参考

• 過去の30分プログラム
• 平方根の求め方