自然数の桁数つづき

shiroさんから鋭いツッコミが入った。

shiro『log10の実装は計算を浮動小数点数で行うものが多いと思われるので、nが10^18近辺を越えたあたりから誤差が入る可能性があります。

gosh> (define (log10 x) (/ (log x) (log 10)))
log10
gosh> (log10 (expt 10 17))
17.0
gosh> (log10 (expt 10 18))
17.999999999999996

従って、大学3年生によるコード:

(define (digits n) (cond ((= n 0) 0) ((< n 10) 1) (else (+ (digits (remainder n 10)) 1))))

』

10^18あたりから誤差がはいってくるのは、倍精度浮動小数点(IEEE754)の仮数部が52bitで、10進数に直したら最小単位が大体 2.2 * 10^-16 なので、有効桁数が大体17桁だから、ってことでしょうか？んー？違う気がする。

---

とりあえず、IEEE754の2を基数とした仮数部の小数点以下ビット列を求める関数を作ってみた。

(define (significand val len)
  (cond
   ((= val 0) 0)
   ((<= len 0) (error "INVALID ARGUMENT"))
   ((< val 1) (significand (* 2 val) len))
   ((> val 2) (significand (/ val 2) len))
   (else
    (let loop((val (- val 1)) (len len) (ret ()))
      (if (= len 0)
	  (reverse ret)
	  (receive (n-val n-ret)
	      (if (>= val 0.5)
		  (values (* 2 (- val 0.5)) (cons 1 ret))
		  (values (* 2 val) (cons 0 ret)))
	    (loop n-val (- len 1) n-ret)))))))

gosh> (significand 10 10)
(0 1 0 0 0 0 0 0 0 0)       # 10 = 1.25 * 2^3 = 1 + 0.01(2) * 2^3
gosh> (significand 1.25 10)
(0 1 0 0 0 0 0 0 0 0)       # 1.25 = 1.25 * 2^0 = 1 + 0.01(2) * 2^0
gosh> (significand 0.75 10)
(1 0 0 0 0 0 0 0 0 0)       # 0.75 = 1.5 * 2^-1 = 1 + 0.1(2) * 2^-1

あれ、よく考えたら val * 2^len を丸めた整数を2進表示すればいいだけじゃねorz

感覚としては、誤差の絶対値が1より大きくなりはじめるあたりで誤差がでてきそうな感じがするんだけど、対数にその誤差が効いてくるのを数式で表現できんなー。

とりあえず、制御工学のテスト勉強でもする。