TAOCP
基本は「The Art of Computer Programming」の中のアルゴリズム。 元はGaucheの数値演算パッケージを読んでて。。。
基数変換
多倍長整数の実装では、まず計算する前にリーダが取り込んで内部にデータオブジェクトを生成する必要がある。 その場合まず必要になるのがおそらく基数変換。
基数変換の項を読む際の基礎事項として、 Vol.1の1.2.4 Integer Functions and Elementary Number Theory を読んでおくのが良い。 ここにfloorとceilingの定義というか記号が指示されている。
- floor
floor xはxを以下の最大整数。xの左右に└と┘で括られる。
- celing
ceiling xはx以上の最小整数。xの左右に┌と┐で括られる。
- xが整数の場合だけ ┌x┐ = └x┘
- xが整数で無い場合だけ ┌x┐ = └x┘ + 1
- x-1 ＜ └x┘ ≦ x ≦ ┌x┐ ＜ x+1なら └ -x┘ = -┌x┐
剰余の定義は x mod y = x - y└x/y┘ (y ≠ 0)であり、 x mod 0 = x である。
基数変換には基本になる4種の方法がある。 最初の2種は整数の変換で残りの2種は小数の変換である。
なお、一般に基数bの小数から基数Bの小数へ変換する際には厳密に同じ値には変換できない。 たとえば、基数3の0.1は1/3であり、基数10に変換したなら0.33333333...となることは明白。 この場合には適当に丸め演算を行う必要がある。
基数変換基本方法1
基数変換基本方法2
基数変換基本方法3
基数変換基本方法4
多倍長整数
多倍長整数の算術演算はマシンワードの算術演算を基礎において処理できる。
加算
アルゴリズムA
ただしx,yはリストで表現しており、リストの長さのチェックや基数bに対して
x,yが正しい入力かどうかなどのチェックはしてない。
gosh> (define (add x y b)
(let loop
((u (reverse x))
(v (reverse y))
(w '())
(k 0))
(if (and (null? u) (null? v))
(cons k w)
(loop (cdr u)
(cdr v)
(cons (modulo (+ (car u) (car v) k) b) w)
(floor (/ (+ (car u) (car v) k) b))))))
add
以下は実行例だ。 10進での18782+18782=37564。あと2進での11+10=101。
gosh> (add '(1 8 7 8 2) '(1 8 7 8 2) 10) (0.0 3.0 7.0 5.0 6.0 4) gosh> (add '(1 1) '(1 0) 2) (1.0 0.0 1)
ちなみに下は内部ではベクタで動作するバージョン。かつ、floorした値をinexact->exactで確定数に整形したもの。
(define (add x y b)
(let ((u (list->vector (reverse x)))
(v (list->vector (reverse y)))
(w (make-vector (+ (length x) 1))))
(do ((j 0 (+ j 1))
(k 0 (inexact->exact
(floor (/ (+ (vector-ref u j)
(vector-ref v j) k)
b)))))
((>= j (length x))
(vector-set! w j k)
(reverse (vector->list w)))
(vector-set! w j (modulo
(+ (vector-ref u j)
(vector-ref v j)
k)
b)))))
下は実行例。
gosh> (add '(8 7 6 5 4) '(9 1 2 4 7) 10) (1 7 8 9 0 1) gosh> (add '(16383 16383 16383 16383) '(0 0 0 1) 16384) (1 0 0 0 0)
減算
アルゴリズムS
ただしx,yはリスト表現しており、リストの長さのチェックや基数bに対して x,yが正しい入力かどうかなどのチェックはしていないし、 xの方がyより大きいかなどもチェックしていない。
gosh> (define (sub x y b)
(let loop
((u (reverse x))
(v (reverse y))
(w '())
(k 0))
(if (and (null? u) (null? v))
w
(loop (cdr u)
(cdr v)
(cons (modulo (+ (- (car u) (car v)) k) b) w)
(floor (/ (+ (- (car u) (car v)) k) b))))))
sub
以下は実行例だ。 10進での98765432-87654321=11111111。あと16進でのcde-bcd=111。
gosh> (sub '(9 8 7 6 5 4 3 2) '(8 7 6 5 4 3 2 1) 10) (1.0 1.0 1.0 1.0 1.0 1.0 1.0 1) gosh> (sub '(12 13 14) '(11 12 13) 16) (1.0 1.0 1)
なお、u＞vなら最後のkは0になるが、逆なら-1になる。 だから必ず最初に大小比較をしておくこと。
下はやはり内部ではベクタで動作するバージョン。こちらもinexact->exactで整形。 なお、除算に備えてkをconsして返すことにした。 減算の解はcadrで取得。
(define (sub x y b)
(let ((u (list->vector (reverse x)))
(v (list->vector (reverse y)))
(w (make-vector (length x))))
(do ((j 0 (+ j 1))
(k 0 (inexact->exact
(floor (/ (+ (- (vector-ref u j)
(vector-ref v j))
k)
b)))))
((>= j (length x))
(if (< k 0) #?=k) ;;; borrow check こいつでオーバーフローの検出が可能
;;; (reverse (vector->list w)))
(list k (reverse (vector->list w))))
(vector-set! w j (modulo
(+ (- (vector-ref u j)
(vector-ref v j))
k)
b)))))
下は実行例。
gosh> (sub '(1 0 0 0) '(0 9 9 9) 10) (0 (0 0 0 1)) gosh> (sub '(1 3 4 9 2 3 5 1) '(0 9 7 5 4 9 7 1) 10) (0 (0 3 7 3 7 3 8 0)) gosh> (sub '(1 2 3 4 5 6 7) '(1 2 3 4 5 6 8) 10) (-1 (9 9 9 9 9 9 9))
乗算
アルゴリズムM
内部をリスト版で組むのが大変そうだったので最初からベクタ版で作成。 ただしx,yはリスト表現しており、基数bに対して x,yが正しい入力かどうかなどのチェックはしていない。
(define (mul x y b)
(let* ((u (list->vector (reverse x)))
(v (list->vector (reverse y)))
(m (length x))
(n (length y))
(w (make-vector (+ m n)))
(t 0))
;初期化
(do ((x 0 (+ x 1)))
((>= x m))
(vector-set! w x 0))
(do ((j 0 (+ j 1)))
((>= j n)
(reverse (vector->list w)))
(if (= (vector-ref v j) 0)
(vector-set! w (+ j m) 0)
(do ((i 0 (+ i 1))
(k 0 (inexact->exact (floor (/ t b)))))
((>= i m) (vector-set! w (+ j m) k))
(set! t (+ (* (vector-ref u i)
(vector-ref v j))
(vector-ref w (+ i j)) k))
(vector-set! w (+ i j) (modulo t b)))))))
以下は実行例。
gosh> (mul '(1 2 3) '(7 8 9) 10) (0 9 7 0 4 7) gosh> (mul '(1 2 3 4 5 6 7 8 9) '(9 8 7 6 5 4 3 2 1) 10) (1 2 1 9 3 2 6 3 1 1 1 2 6 3 5 2 6 9)
除算
アルゴリズムD
一番面倒なやつ。
除算は自分で紙と鉛筆でやるときにも分かるけど、商を求めるのに推測を繰り返しつつ求めている。 つまり３でいけるかな〜？お、いける、じゃぁ４は？いけるねぇ、じゃぁ５でも割れるか？・・・あ、ダメだオーバーしちゃうよ〜、じゃぁ４だな。って感じで探索しているわけだ。
この推測値をいかに機械的に求められ、かつ高速に適正な値を割り出せる様にするかが、Normalizeと推定値q^の出し方に掛かっている模様。
「基本的には x/y をした時の商と (n*x)/(n*y) をした時の商は同じ」なのだが、n倍した時には上位の桁に対して下位の桁からの繰り上がりが反映されるため、推定値を割り出す際に最上位桁だけでなく、何割りか下位桁の情報を持ち込めることにより高速に適正な推定値から始められるということのようだ。(と思う)
細かい理屈は全く分からんが、大雑把にいってそんなところだろう。 で、nは？っていうと分母の方が桁数が増えない範囲でそれなりに大きくしており、それ以上のでかい数を掛けないのは何故？って聞かれると分からない。 効率からかそれ以上の値を掛けても意味が無いからか、はたまた論理的にそれでいいのか。。。
とにかくコードを書いてみる。 ただし、Add Back のところは省いている。 減算時に負になるってのをどう検出していいのか今のところ理解できてない。(基数bの補数表現になるんだけど)
まぁ、確率的にも低いらしいので、まずコードを書いていろいろ試して発生する具体例をみつけてからでないとデバッグできないよね〜ってのが一応の言い訳。:-P
あと、余りが・・・理屈は簡単だけどうまくコードに書けてない。後述。
(define (div x y b)
(let* ((n (length y))
(m (- (length x) n))
(d (inexact->exact (floor (/ b (+ (car y) 1))))))
;exercises 16
;vが1桁の場合を特殊扱い
;short division
(if (= n 1)
(let* ((m (length x))
(w (make-vector m))
(r 0)
(u (list->vector (reverse x)))
(v (car y)))
(do ((j (- m 1) (- j 1)))
((< j 0) (if (= r 0)
(reverse (vector->list w))
(list (reverse (vector->list w)) r)))
(vector-set! w j (inexact->exact (floor (/ (+ (* r b)
(vector-ref u j))
v))))
(set! r (modulo (+ (* r b)
(vector-ref u j))
v))))
;vが1桁より大きい場合。
(let ((u (list->vector (reverse (mul x (list d) b))))
(v (list->vector (reverse (cdr (mul y (list d) b)))))
(q (make-vector (+ m 1))))
(do ((j m (- j 1)))
((< j 0) (list (reverse (vector->list q))
(div (reverse (vector->list u)) (list d) b)))
(let loop ((q^ (inexact->exact
(floor (/ (+ (* (vector-ref u (+ j n)) b)
(vector-ref u (+ j n -1)))
(vector-ref v (- n 1))))))
(r^ (modulo (+ (* (vector-ref u (+ j n)) b)
(vector-ref u (+ j n -1)))
(vector-ref v (- n 1)))))
(if (or (= q^ b)
(> (* q^ (vector-ref v (- n 2)))
(+ (* b r^) (vector-ref u (+ j n -2)))))
(begin
(set! q^ (- q^ 1))
(set! r^ (+ r^ (vector-ref v (- n 1))))
(if (< r^ b)
(loop q^ r^))))
(let ((u-temp (reverse (vector->list u)))
(qv (append (make-list 0 (- m j))
(mul (reverse (vector->list v)) (list q^) b)
(make-list 0 j))))
(set! u (list->vector (reverse (sub u-temp qv b)))))
(vector-set! q j q^)
))))))
(define (make-list mem n)
(do ((r '() (cons mem r))
(c n (- c 1)))
((= c 0) r)))
とりあえずの実行例はこんな感じ。
gosh> (div '(1 2 3 4 5 6 7 8 9) '(1 2 3) 10) (1 0 0 3 7 1 3) gosh> (/ 123456789 123) 1003713.731707317 gosh> (div '(9 8 7 6 5 4 3 2 1 9 8 7 6 5 4 3 2 1) '(1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9) 10) (8) gosh> (/ 987654321987654321 123456789123456789) 8.0000000729
問題はアルゴリズムDでは剰余も出しているのだが、これがどうやったらいいもんだか。。。
というのも余りは上でいうところの n、コード中のd倍されており、 最終的に余った多倍長整数uをn(つまりd)で割ってやる必要がある。 これって再帰的に呼び出すことになるから商と余りをconsして返そうとすると・・・おいおい。。。 余りを出すためにdivを呼び出すと無限に・・・多分割り切れるから余りの全桁が0なら余り無しってことで、条件分けするしかないのだろうか？・・・でも元の計算でたまたま余りが無い場合はなんか返り値が中途半端だな〜。
除算デバッグ
ダメだしから。 上記のコードはバグってました。 どうにもおかしくなって来たので完全に書き直すことにする。。。あ〜あ。
;;; 作り直し除算
(define (div x y b)
(let* ((n (length y))
(m (- (length x) n))
(u (list->vector (reverse x)))
(v (list->vector (reverse y)))
(d (inexact->exact
(floor (/ b
(+ (vector-ref v (- n 1)) 1))))))
(if (= n 1)
;;; yが1桁
;;; exercises 16 short division
(let* ((m (length x))
(w (make-vector m))
(r 0)
(u (list->vector (reverse x)))
(v (car y)))
(do ((j (- m 1) (- j 1)))
((< j 0) (if (= r 0)
(reverse (vector->list w))
(list (reverse (vector->list w)) r)))
(vector-set! w j(inexact->exact (floor (/ (+ (* r b)
(vector-ref u j))
v))))
(set! r (modulo (+ (* r b)
(vector-ref u j))
v))))
;;; yが2桁以上
(let ((u (list->vector (reverse (mul x (list d) b))))
(v (list->vector (reverse (cdr (mul y (list d) b)))))
(q (make-vector (+ m 1))))
(do ((j m (- j 1))
(q^ 0 0)
(r^ 0 0))
((< j 0)
(list (reverse (vector->list q))
(div (reverse (vector->list u)) (list d) b)))
;;; q^ r^ の初期推定値
(set! q^ (inexact->exact
(floor (/ (+ (* (vector-ref u (+ j n)) b)
(vector-ref u (+ j n -1)))
(vector-ref v (- n 1))))))
(set! r^ (modulo (+ (* (vector-ref u (+ j n)) b)
(vector-ref u (+ j n -1)))
(vector-ref v (- n 1))))
;;; q^ r^の決定
(let loop ()
(if (or (= q^ b)
(> (* q^ (vector-ref v (- n 2)))
(+ (* b r^) (vector-ref u (+ j n -2)))))
(begin
(set! q^ (- q^ 1))
(set! r^ (+ r^ (vector-ref v (- n 1))))
(if (< r^ b)
(loop)))))
;;; #?=j #?=q^ #?=r^ ;;; for debug
;;; u の書き換え
(set! u (list->vector
(reverse (cadr
(sub
(reverse (vector->list u))
(append (make-list 0 (- m j))
(mul (reverse (vector->list v)) (list q^) b)
(make-list 0 j))
b)))))
(vector-set! q j q^)
)))))
(define (check x y)
(list (quotient x y) (modulo x y)))
実行結果
gosh> (load "./multi-prec.scm") #t gosh> (div '(1 2 3 4 5 6 7 8 9) '(6 2) 10) ((0 1 9 9 1 2 3 8) (0 0 0 0 0 0 0 0 3 3)) gosh> (check 123456789 62) (1991238 33) gosh> (check 12345 321) (38 147) gosh> (div '(1 2 3 4 5) '(3 2 1) 10) ((0 3 8) (0 0 0 1 4 7)) gosh> (div '(9 8 7 6 5 4 3 2 1 9 8 7 6 5 4 3 2 1) '(1 2 3 4 5 6 7 8 9) 10) ((8 0 0 0 0 0 0 0 8 0) (0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 2 0 1)) gosh> (check 987654321987654321 123456789) (8000000080 111111201) gosh> (div '(1 2 3 4 5 6 7) '(5 4 3 2 1) 10) ((0 2 2) (0 0 0 3 9 5 0 5)) gosh> (check 1234567 54321) (22 39505)
ようやく安定したかな。 add back発生がなかなか見付からない。cut-sea:2003/12/06 23:07:30 PST
add backの追加／デバッグ
gosh> (div '(1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0) '(9 8 7 6 5 4 3 2 1) 10) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 9 8 7 6 5 4 3 2 1 0 0 0 0 0 0 0 0 0 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (0 0 2 4 6 9 1 3 5 6 9 1 2 3 4 5 6 7 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 1 9 7 5 3 0 8 6 4 2 0 0 0 0 0 0 0 0 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (0 0 0 4 9 3 8 2 7 0 4 9 2 3 4 5 6 7 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 4 9 3 8 2 7 1 6 0 5 0 0 0 0 0 0 0 0) #?=k ↑ここ！結構下の方の桁でちょろっとオーバーしちゃってる #?- -1 ←これはsubから出しているデバッグプリント #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (9 9 9 9 9 9 9 9 9 8 8 8 7 3 4 5 6 7 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 0 8 8 8 8 8 8 8 8 8 9 0 0 0 0 0 0 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (9 9 9 9 1 1 1 1 0 9 9 9 8 4 4 5 6 7 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 0 0 0 9 8 7 6 5 4 3 2 1 0 0 0 0 0 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (9 9 9 9 1 0 1 2 3 3 4 5 5 2 3 5 6 7 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 0 0 0 0 9 8 7 6 5 4 3 2 1 0 0 0 0 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (9 9 9 9 1 0 0 2 4 5 8 0 0 9 1 4 6 7 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 0 0 0 0 1 9 7 5 3 0 8 6 4 2 0 0 0 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (9 9 9 9 1 0 0 0 4 8 2 7 0 0 5 0 4 7 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 0 0 0 0 0 3 9 5 0 6 1 7 2 8 4 0 0 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (9 9 9 9 1 0 0 0 0 8 7 6 3 8 7 7 6 3 8 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 0 0 0 0 0 0 7 9 0 1 2 3 4 5 6 8 0 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (9 9 9 9 1 0 0 0 0 0 8 6 2 6 4 3 0 7 0 9 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 0 0 0 0 0 0 0 7 9 0 1 2 3 4 5 6 8 0) #?="./multi-prec.scm":128:(reverse (vector->list u)) #?- (9 9 9 9 1 0 0 0 0 0 0 7 2 5 1 9 6 1 4 1 0) #?="./multi-prec.scm":129:(append (make-list 0 (- m j)) (mul (reverse (v ... #?- (0 0 0 0 0 0 0 0 0 0 0 6 9 1 3 5 8 0 2 4 7) ((0 1 2 5 9 1 1 2 4 8 8 7) (9 9 9 9 1 0 0 0 0 0 0 0 3 3 8 3 8 1 1 6 3)) gosh> (check 12345678901234567890 987654321) (12499999887 339506163)
これがテストケースかな。今度こそ大丈夫か？
でもってadd backを実装したコードがこちら。
;;; 作り直し除算
(define (div x y b)
(let* ((n (length y))
(m (- (length x) n))
(u (list->vector (reverse x)))
(v (list->vector (reverse y)))
(d (inexact->exact
(floor (/ b
(+ (vector-ref v (- n 1)) 1))))))
(if (= n 1)
;;; yが1桁
;;; exercises 16 short division
(let* ((m (length x))
(w (make-vector m))
(r 0)
(u (list->vector (reverse x)))
(v (car y)))
(do ((j (- m 1) (- j 1)))
((< j 0) (if (= r 0)
(reverse (vector->list w))
(list (reverse (vector->list w)) r)))
(vector-set! w j(inexact->exact (floor (/ (+ (* r b)
(vector-ref u j))
v))))
(set! r (modulo (+ (* r b)
(vector-ref u j))
v))))
;;; yが2桁以上
(let ((u (list->vector (reverse (mul x (list d) b))))
(v (list->vector (reverse (cdr (mul y (list d) b)))))
(q (make-vector (+ m 1))))
(do ((j m (- j 1))
(q^ 0 0)
(r^ 0 0))
((< j 0)
(list (reverse (vector->list q))
(div (reverse (vector->list u)) (list d) b)))
;;; q^ r^ の初期推定値
(set! q^ (inexact->exact
(floor (/ (+ (* (vector-ref u (+ j n)) b)
(vector-ref u (+ j n -1)))
(vector-ref v (- n 1))))))
(set! r^ (modulo (+ (* (vector-ref u (+ j n)) b)
(vector-ref u (+ j n -1)))
(vector-ref v (- n 1))))
;;; q^ r^の決定
(let loop ()
(if (or (= q^ b)
(> (* q^ (vector-ref v (- n 2)))
(+ (* b r^) (vector-ref u (+ j n -2)))))
(begin
(set! q^ (- q^ 1))
(set! r^ (+ r^ (vector-ref v (- n 1))))
(if (< r^ b)
(loop)))))
;;; #?=j #?=q^ #?=r^ ;;; for debug
;;; u の書き換え
;;; (set! u (list->vector
;;; (reverse (cadr
;;; (sub
;;; #?=(reverse (vector->list u))
;;; #?=(append (make-list 0 (- m j))
;;; (mul (reverse (vector->list v)) (list q^) b)
;;; (make-list 0 j))
;;; b)))))
(let* ((sub-return (sub
(reverse (vector->list u))
(append (make-list 0 (- m j))
(mul (reverse (vector->list v)) (list q^) b)
(make-list 0 j))
b))
(sub-flag (car sub-return))
(sub-result (cadr sub-return)))
;;; (set! u (list->vector (reverse sub-result)))
;;; add back
(if (< sub-flag 0)
(begin
(set! q^ (- q^ 1))
(set! sub-result (cadr (sub
(reverse (vector->list u))
(append (make-list 0 (- m j))
(mul (reverse (vector->list v)) (list q^) b)
(make-list 0 j))
b)))))
;;; q^ のストア
(vector-set! q j q^)
(set! u (list->vector (reverse sub-result)))
))))))
(define (check x y)
(list (quotient x y) (modulo x y)))
実行確認。
gosh> (load "./multi-prec.scm") #t gosh> (div '(1 2 3 4 5 6 7 8 9 0) '(9 8 7 6 5 4 3 2 1) 10) ((0 1) (0 0 2 4 6 9 1 3 5 6 9)) gosh> (div '(1 2 3 4 5 6 7 8 9 0) '(9 8 7) 10) ((0 1 2 5 0 8 2 8) (0 0 0 0 0 0 0 0 6 5 4)) gosh> (check 1234567890 987654321) (1 246913569) gosh> (div '(1 2 3 4 5 6 7 8 9 0) '(7 6 5 4 3 2 1) 10) ((0 1 6 1) (0 0 0 0 2 2 2 2 2 0 9)) ;;;本命のこれがさっきまでadd back不在により間違った値を返していた gosh> (div '(1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0) '(9 8 7 6 5 4 3 2 1) 10) ((0 1 2 4 9 9 9 9 9 8 8 7) (0 0 0 0 0 0 0 0 0 0 0 0 3 3 9 5 0 6 1 6 3)) gosh> (check 12345678901234567890 987654321) (12499999887 339506163)
コメント
- Shiro: Schemeのfloorはinexactな数を渡すとinexactな数値が帰って来るので、
せっかく多倍長でexactな計算をしても桁がinexactになってしまって気持ち悪いです。
(inexact->exact (floor (/ X Y)))
とか、いっそのこと(quotient X Y)
とするとよろしいのでは。
- 確定／非確定がまだ良く分かってないもので・・・。(^^;; 教科書にならうなら上ので書いた方が自分としては分かりよかったかもしれないですね。cut-sea:2003/12/05 06:56:42 PST
- とりあえずfloorに対してinexact->exactで確定数に。あとx->integerをマップして返してた部分はそのまま返すように。cut-sea:2003/12/05 20:11:29 PST
- 除算アルゴリズムでadd backを通るテストケースを見付けるのは結構面倒です。 ランダムにやってもなかなかあたりません。途中の数値をダンプしながら 色々な入力を試していると、パターンがあるので、それを元に少しづつ近づけて ゆくしか。ちなみに、Gaucheのbignum.cではhalf word (ILP32なら16bit) を 1桁で計算していますが、その場合にadd backを通るテストケースは test/number.scmにあります ("add back"で検索してください)。
- それから、アルゴリズムDでのプリスケール/ポストスケールですが、 うまくdを選んでやれば再帰にmulやdivを呼ぶ必要はなくて、シフトだけで済みます。
- あっ、そうか！dが基数bの冪乗倍にとればよいわけですね。
d = b ^ n
ex.) b = 2でn = 3を選んだ場合
d = 2 ^ 3 = 8
となるから、コード中にでてくる(mul y d b) はyを基数b=2のシステム(標準的な2進マシン)でd=8倍するってのはn=3桁左シフト！ あとはbのn乗のnを決定する時にどういった条件で決めてやるかベンチマークするなりしてコーディングしてやればよいわけだ。 ってことは別にTAOCPにあるq^で推定値を設定する必要はやはりないわけですね。cut-sea:2003/12/05 06:27:20 PST
- あれ？でもなんか変だな。最初に除数と被除数を8倍するのも単なる左シフトだから掛けるという意味がないような気がする。 最終的にシフト演算だけですませようとしているわけだからd=b^nには違いないはずだけど。
- 基数bに関係なく2の冪乗でかつvの桁が増えないで済む範囲のdを求めて やれば確かにビットシフトで済むな。 それなら基数の冪乗と違って意味があるように思える。2003/12/05 19:30:58 PST