Hi,
Are there scaled (fixed point) versions of division and sqrt in picolisp? I
may have missed them...
In any event I wrote a couple of functions to provide that functionality
that may be useful to others as well. Performance is reasonable.
It is quite nice to be able to be able to do arbitrary precision arithmetic
in this way =)
/Lindsay
Code:
# Scaled fixed point long division. (if necessary, round last digit)
# TODO: Brute force implementation...
# There is probably a better way to do it.
: (pp '/*)
(de /* (N D S)
(let (Q NIL R NIL Acc NIL Cnt 0)
(default S *Scl)
(setq S (+ 2 S))
(do S
(T (>= Cnt S))
(setq
Q (/ N D)
R (% N D)
Acc (cons Q Acc)
Cnt (inc Cnt)
N R )
(when (and (gt0 N) (<= N D))
(while (and (gt0 N) (<= N D))
(setq
N (* N 10)
Acc (cons 0 Acc)
Cnt (inc Cnt) ) )
(pop 'Acc)
(dec 'Cnt) ) )
(setq R (pop 'Acc))
(if (<= 5 R)
(setq Acc (cons (+ 1 (pop 'Acc)) Acc)) )
(setq Acc (flip Acc))
(format
(pack (cons (car Acc) "." (cdr Acc)))
(- S 2) ) ) )
-> /*
# Scaled fixed point sqrt(); no rounding
# TODO: Converges in 4-7 iterations in tests, 13 is probably overkill
: (pp 'sqrt*)
(de sqrt* (N S)
(let (P (/ N 2) M 13)
(default S *Scl)
(setq N (* N (** 10 S)))
(while (ge0 (dec 'M))
(setq P (/ (+ P (/ N P)) 2)) )
P ) )
-> sqrt*
Some Tests:
# Arbitrary Fixed Precision Square Root
# Compare: bc.
# also: https://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil
: (scl 64) (format (sqrt* 2.0) *Scl)
-> "1.4142135623730950488016887242096980785696718753769480731766797379"
: (scl 64) (format (sqrt* 0.5) *Scl)
-> "0.7071067811865475244008443621048490392848359376884740365883398689"
: (scl 64) (format (sqrt* 9967.0) *Scl)
-> "99.8348636499294268455686673311236280296661789737252407300182230035"
: (scl 100000) (bench (nil (sqrt* 2.0)))
2.027 sec
# Arbitrary Fixed Precision Long Division
# Compare: bc e.g. "scale=64 1/9967"
: (scl 64) (/* 1.0 9967.0)
-> 1003310926055984749673923949031804956355974716564663389184308
: (scl 64) (format (/* 1.0 9967.0) *Scl)
-> "0.0001003310926055984749673923949031804956355974716564663389184308"
: (scl 64) (format (format (/* 1.0 9967.0) *Scl) *Scl)
-> 1003310926055984749673923949031804956355974716564663389184308
: (scl 32) (/* 22.0 7.0)
-> 314285714285714285714285714285714
: (scl 0) (/* 22.0 7.0)
-> 3
: (scl 1) (/* 22.0 7.0)
-> 31
: (scl 8) (/* 22.0 7.0)
-> 314285714
: (scl 32) (/* 1.0 3.0)
-> 33333333333333333333333333333333
: (scl 32) (/* 10.0 3.0)
-> 333333333333333333333333333333333
: (scl 32) (format (/* 1.0 3.0) *Scl)
-> "0.33333333333333333333333333333333"
: (scl 32) (format (/* 10.0 3.0) *Scl)
-> "3.33333333333333333333333333333333"
: (scl 32) (format (/* 0.22 0.7) *Scl)
-> "0.31428571428571428571428571428571"
: (scl 32) (format (/* 9968.0 32.0) *Scl)
-> "311.50000000000000000000000000000000"
: (scl 10000) (bench (nil (format (/* 1.0 9967.0) *Scl)))
7.685 sec
