Re: Scaled division and sqrt

2017-03-01 Thread Lindsay John Lawrence 
Thanks Eric.

This was also an opportunity to start dabbling a bit with the app ui
framework. Very interesting.

At least initially, it reminds me a bit of the way you could write
PHP... but this feels more elegant.
Even lisp looks less like line noise than php ;)

The javascript canvas pipeline I utilized there is quite powerful.
There is a lot more I can do with it.
For the denser images, the server side code is generating 1M+ points
which are then  piped back to the client (~3Mb) for rendering.

Once I am more familiar with it I'd like to exploit the session state
a bit more to dribble data streams in smaller chunks.

/Lindsay
-- 
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-03-01 Thread Erik Gustafson 
Hi Lindsay,

I practiced a bit more with the built-in functions...

https://github.com/thinknlive/picolisp-gosper.git


Very cool! Looks like you've got '*/' all figured out now :)

Thanks for sharing

- Erik

Re: Scaled division and sqrt

2017-02-28 Thread Lindsay John Lawrence 
allbase   -->  tankfeeder / exercism-io / a-f.l

/Lindsay


On Tue, Feb 28, 2017 at 3:41 AM, Danilo Kordic  wrote:
>> If you are going to use a list, could also just put the 'scale' in a property
>> of the list?
>
>   I wouldn't call it a propert to avoid confusion with PicoLisp `prop'erties.
>
>   I couldn't find `allbase'.
>
>   I had [[https://en.wikipedia.org/wiki/Algebraic_data_type][Algebraic
> DataTypes]]
> in mind, just like fractions.
>
> (= (Float B S E)
>(* S (** B E)) )
>
>
>   ``(314 E -2)'' looks like
> [[https://en.wikipedia.org/wiki/Scientific_notation][Scientific
> notation]].
> And wouldn't cause an error in PicoLisp.
>
>   Current implementation of PicoLisp numbers is also a linked list:
> [path "@doc64/structures"].  Too bad they can not be manipulated like any 
> other
> list.  Base is (** 2 64).
>
>   [[https://en.wikipedia.org/wiki/Quote_notation][Quote notation]]: If we 
> allow
> circular reference in the `Num' structure it will result in a subset of
> rational numbers.  Arithmetic Shift, `AS', will give the rest.
> ``I mean, how hard could it be?  ''
> --
> UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe
-- 
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-02-28 Thread Lindsay John Lawrence 
Hi,

I practiced a bit more with the built-in functions...

https://github.com/thinknlive/picolisp-gosper.git

/Lindsay
-- 
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-02-28 Thread Danilo Kordic 
> If you are going to use a list, could also just put the 'scale' in a property
> of the list?

  I wouldn't call it a propert to avoid confusion with PicoLisp `prop'erties.

  I couldn't find `allbase'.

  I had [[https://en.wikipedia.org/wiki/Algebraic_data_type][Algebraic
DataTypes]]
in mind, just like fractions.

(= (Float B S E)
   (* S (** B E)) )


  ``(314 E -2)'' looks like
[[https://en.wikipedia.org/wiki/Scientific_notation][Scientific
notation]].
And wouldn't cause an error in PicoLisp.

  Current implementation of PicoLisp numbers is also a linked list:
[path "@doc64/structures"].  Too bad they can not be manipulated like any other
list.  Base is (** 2 64).

  [[https://en.wikipedia.org/wiki/Quote_notation][Quote notation]]: If we allow
circular reference in the `Num' structure it will result in a subset of
rational numbers.  Arithmetic Shift, `AS', will give the rest.
``I mean, how hard could it be?  ''
-- 
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-02-27 Thread Lindsay John Lawrence 
Hi Danilo,

Ah. I get it. Interesting idea. If you are going to use a list, could also
just put the 'scale' in a property of the list?

When I came across the 'allbase function in Mike Penchkin's code... I
thought of the following...
It would be a bit wasteful of storage, but could represent arbitrarily
sized fixed or floating point numbers in any base as a list, each list atom
being a digit in that base. If you keep the digits ordered (least...most)
significant,  the performance may be quite reasonable as well.

The arithmetic functions to manipulate those 'lists' would be an
interesting exercise in basic concepts... and with picolisp numbers, you
could work with not just bignums, but big number bases. For example..
simple positive addition, base10, might look something like this (code at
end)...

: (num2L 1234)
-> (4 3 2 1)
: (L2Num (num2L 1234))
-> "1234"
: (format (L2Num (num2L 1234)))
-> 1234
: (L+ (num2L 1234) (num2L 654321))
-> (5 5 5 5 5 6)
: (num2L (L+ (num2L 1234) (num2L 654321)))
-> (6 5 5 5 5 5)
: (format (L2Num (L+ (num2L 1234) (num2L 654321
-> 65


/Lindsay

# TODO: Generalize Num2L, L2Num, or add similar functions to accept any
number base and scale.
(de num2L (N)
   (mapcar format (reverse (chop (format N
)

(de L2Num (L)
   (pack (reverse L))
)

# Positive number addition (the easiest one)
(de L+ (L1 L2)
   (let (R NIL  S 0  D 0  C 0)
  (default L1 '(0) L2 '(0))
  (if (> (length L2) (length L1)) (swap 'L1 (swap 'L2 L1)))
  (mapcar
 '((D1 D2)
(setq S
   (+ (or D1 0) (or D2 0) C) )
(until (< S 10)
   (setq S (- S 10))
   (inc 'C) )
(setq R (cons S R)) )
 L1
 L2 )
  (if (> 0 C) (setq R (cons C R)))
  (reverse R)) )


On Mon, Feb 27, 2017 at 8:25 AM, Danilo Kordic 
wrote:

> By constructing base 10 floating point numbers with
> ``[list Significand 'E Exponent]''.
> As if `E' is a binary operator.  Of course that makes no sense in Lisp.
> ``(Decimal Significand Exponent)'' would make sense.
>
> For example:
>   : (sqrt 2)
>   -> (141421356237309505 E -17)
> --
> UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe
>

Re: Scaled division and sqrt

2017-02-27 Thread Lindsay John Lawrence 
Hi Alex, Erik,

I definitely took the 'long way around' on this one... and got some good
exercise along the way! =)

Thank you for the detailed explanation. I get it now... although I need a
bit more practice and care in using the functions comfortably.

My initial impetus was to be able to easily, and repeatably, generate long
sequences of digits. Some with repeating patterns. Others not.
Rational and irrational numbers  are one simple way to do that and the
support for arbitrary size numbers in picolisp make it easy to do so.
Hence the examples with sqrt(2) and 1/[prime], etc.

Kind Regards,
/Lindsay

Re: Scaled division and sqrt

2017-02-27 Thread Danilo Kordic 
I forgot to say: ``As if `E' is an _infix_ binary operator.  ''
-- 
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-02-27 Thread Danilo Kordic 
By constructing base 10 floating point numbers with
``[list Significand 'E Exponent]''.
As if `E' is a binary operator.  Of course that makes no sense in Lisp.
``(Decimal Significand Exponent)'' would make sense.

For example:
  : (sqrt 2)
  -> (141421356237309505 E -17)
-- 
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-02-27 Thread Alexander Burger 
Hi Lindsay,

> Having just tried it... I am surprised that the 'sqrt*  I implemented is
> quite a bit faster! than the built-in 'sqrt. I'll try the long division
> later.

Wow, that's indeed a bit surprising. However, it depends on the scale, and I
think a can explain (see below).

> : (scl 10) (bench (nil (sqrt 2.0)))   # Built-in
> 30.917 sec
> -> NIL
> : (scl 10) (bench (nil (sqrt* 2.0)))   #
> 2.906 sec

Which scale did you use here? If I try 2

   : (scl 2)
   -> 2

   : (bench (do 100 (sqrt 2.0 1.0)))  # Built-in
   0.106 sec
   -> 141

   : (bench (do 100 (sqrt* 2.0)))
   5.052 sec
   -> 141

and then bigger ones,

   : (scl 20)
   -> 20

   : (bench (do 10 (sqrt 2.0 1.0)))  # Built-in
   1.435 sec
   -> 141421356237309504880

   : (bench (do 10 (sqrt* 2.0)))
   6.103 sec
   -> 141421356237309504880

the built-in is a quite faster.


But for bigger scales

   : (scl 64)
   -> 64

   : (bench (do 1 (sqrt 2.0 1.0)))  # Built-in
   1.269 sec
   -> 14142135623730950488016887242096980785696718753769480731766797380

   : (bench (do 1 (sqrt* 2.0)))
   0.762 sec
   -> 14142135623730950488016887242096980785696718753769480731766797379

the built-in is slower.


I think this is because the built-in 'sqrt' first multiplies with the scale
(just as we do before a division), and thus calculates with bignums of double
length.

I think this is necessary for a precise result.

♪♫ Alex
-- 
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-02-27 Thread Alexander Burger 
Hi Lindsay,

On Mon, Feb 27, 2017 at 06:28:15AM -0600, Erik Gustafson wrote:
> With '*/', I like to think of it like this: for multiplication, the scale
> factor is last, e.g.
> 
>(*/ X Y Scl)
> 
> For division, the scale factor comes first, e.g.
> 
>(*/ Scl X Y)

Exactly. The reason is simply that if you multiply two scaled integers, you have
to divide by the scale factor afterwards.

For example, assume your integers denote centimeters. You set *Scl to 2, so that

   100 cm = 1.0 m

Now, one square meter is

   1.0 m * 1.0 m  =  100 cm * 100 cm  =  1.0 m^2  =  1 cm^2

As we want the result not in cm but in meters, we must divide by 100.
Example:

   :(scl 2)
   -> 2

   : (*/ 3.0 4.0 1.0)  # This is equal to (*/ 300 400 100)
   -> 1200

   : (round @)
   -> "12.00"


For division you need the opposite, multiply with 100 instead of dividing.
So you first multiply thi dividend by the scale factor before you divide:

   (*/ 12.0 1.0 3.0)  # this is 12.0 / 3.0


If you multiply more than two numbers, you have to use the corresponding scale
factor:

   : (*/ 3.0 4.0 5.0 `(* 1.0 1.0))  # 3.0 * 4.0 * 5.0
   -> 6000

   : (round @)
   -> "60.00"

♪♫ Alex
-- 
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-02-27 Thread Lindsay John Lawrence 
Thanks Eric,

I am going to have to play with those two particular built-in functions a
bit more. Their use, especially '*/ was not intuitive to me. At first, or
second, reading.

Having just tried it... I am surprised that the 'sqrt*  I implemented is
quite a bit faster! than the built-in 'sqrt. I'll try the long division
later.

: (scl 10) (bench (nil (sqrt 2.0)))   # Built-in
30.917 sec
-> NIL
: (scl 10) (bench (nil (sqrt* 2.0)))   #
2.906 sec
-> NIL
: (scl 1) (bench (nil (sqrt 2.0))) # Built-in
0.253 sec
-> NIL
: (scl 1) (bench (nil (sqrt* 2.0)))
0.044 sec
-> NIL


1 Million digits of sqrt(2):

The 'sqrt* function implemented does need more than 13 iterations to
converge accurately when the number of digits are > 300K or so. The
corrected version is below.  It now checks for convergence. 1M digits of
the sqrt(2) takes 20 iterations to converge. The results, dumped to a file,
are the same as the first 1M digits here...
https://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil

I'll say it again.. I am really impressed by picolisp. I have not had this
much fun expressing ideas in code in some time.

/Lindsay


: (scl 100) (bench (nil (sqrt* 2.0)))
323.805 sec
: (scl 100) (bench (out "sqrt2-1M.txt" (nil (prin (sqrt* 2.0)
932.620 sec


# Scaled fixed point sqrt()
# ... check for convergence added.
(de sqrt* (N S)
   (let (P (/ N 2)  M 63  C (0 0 .))
  (default S *Scl)
  (set C P  (cdr C) 0)
  (setq N (* N (** 10 S)))
  (while
 (and
(ge0 (dec 'M))
(<> (car C) (cadr C) P) )
 (setq P (/ (+ P (/ N P)) 2))
 (pop 'C)
 (set C P))
 (set C 0  (cdr C) 0)
  P ) )

Re: Scaled division and sqrt

2017-02-27 Thread Erik Gustafson 
Hi Lindsay,


With the functions I implemented I can write something like...

: (scl 64) (format (sqrt* 2.0) *Scl)
-> "1.4142135623730950488016887242096980785696718753769480731766797379"


'sqrt' also accepts a 'scl' argument

: (format (sqrt 2.0 1.0) *Scl)
->
"1.4142135623730950488016887242096980785696718753769480731766797380"


: (scl 64) (format (/* 1.0 9967.0) *Scl)
-> "0.0001003310926055984749673923949031804956355974716564663389184308"


: (format (*/ 1.0 1.0 9967.0) *Scl)
-> "0.0001003310926055984749673923949031804956355974716564663389184308"

With '*/', I like to think of it like this: for multiplication, the scale
factor is last, e.g.

   (*/ X Y Scl)

For division, the scale factor comes first, e.g.

   (*/ Scl X Y)

Pretty weird at first, indeed!

- Erik

Re: Scaled division and sqrt

2017-02-27 Thread Lindsay John Lawrence 
Hi Danilo,

Sorry, I do not follow this. Can you explain a bit more?

Thanks!
/Lindsay


On Mon, Feb 27, 2017 at 3:32 AM, Danilo Kordic 
wrote:

> It seems You reimplemented `*/'.
>
> How about something like:
>   # [name 'sqrt "isqrt"]  # Somthing like that.
>   [scl [if native 17 8]]
>   [de sqrt [N]
> [list [isqrt (** 100 *Scl) N] 'E [- *Scl]] ]
>
>
>

Re: Scaled division and sqrt

2017-02-27 Thread Lindsay John Lawrence 
Hi Alex,

I am missing something basic here that article is how I ended up
writing the functions.

With the functions I implemented I can write something like...

: (scl 64) (format (sqrt* 2.0) *Scl)
-> "1.4142135623730950488016887242096980785696718753769480731766797379"

: (scl 64) (format (/* 1.0 9967.0) *Scl)
-> "0.0001003310926055984749673923949031804956355974716564663389184308"

How do I get similar results with the built-in functions?


/Lindsay

Re: Scaled division and sqrt

2017-02-27 Thread Danilo Kordic 
It seems You reimplemented `*/'.

How about something like:
  # [name 'sqrt "isqrt"]  # Somthing like that.
  [scl [if native 17 8]]
  [de sqrt [N]
[list [isqrt (** 100 *Scl) N] 'E [- *Scl]] ]



On 2/27/17, Joh-Tob Schäg  wrote:
> Now there are before there were not.
> Am 27.02.2017 10:07 schrieb "Lindsay John Lawrence" <
> lawrence.lindsayj...@gmail.com>:
>
>> 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 10) (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)
>> -> 
>> : (scl 32) (/* 10.0 3.0)
>> -> 3
>> : (scl 32) (format (/* 1.0 3.0) *Scl)
>> -> "0."
>> : (scl 32) (format (/* 10.0 3.0) *Scl)
>> -> "3."
>> : (scl 32) (format (/* 0.22 0.7) *Scl)
>> -> "0.31428571428571428571428571428571"
>> : (scl 32) (format (/* 9968.0 32.0) *Scl)
>> -> "311.5000"
>> : (scl 1) (bench (nil (format (/* 1.0 9967.0) *Scl)))
>> 7.685 sec
>>
>
--
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-02-27 Thread Alexander Burger 
Hi Lindsay,

> Are there scaled (fixed point) versions of division and sqrt in picolisp? I
> may have missed them...

The built-in '*/' function is both for multiplication *and* division. And 'sqrt'
accepts an optional scale factor.

For details, take a look at Rick's excellent article on the matter:

   http://the-m6.net/blog/fixed-point-arithmetic-in-picolisp.html


In any case, congratulation for writing your own functions! That is not trivial
at all! :)

♪♫ Alex
-- 
UNSUBSCRIBE: mailto:picolisp@software-lab.de?subject=Unsubscribe

Re: Scaled division and sqrt

2017-02-27 Thread Joh-Tob Schäg 
Now there are before there were not.
Am 27.02.2017 10:07 schrieb "Lindsay John Lawrence" <
lawrence.lindsayj...@gmail.com>:

> 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 10) (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)
> -> 
> : (scl 32) (/* 10.0 3.0)
> -> 3
> : (scl 32) (format (/* 1.0 3.0) *Scl)
> -> "0."
> : (scl 32) (format (/* 10.0 3.0) *Scl)
> -> "3."
> : (scl 32) (format (/* 0.22 0.7) *Scl)
> -> "0.31428571428571428571428571428571"
> : (scl 32) (format (/* 9968.0 32.0) *Scl)
> -> "311.5000"
> : (scl 1) (bench (nil (format (/* 1.0 9967.0) *Scl)))
> 7.685 sec
>

Scaled division and sqrt

2017-02-27 Thread Lindsay John Lawrence 
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 10) (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)
-> 
: (scl 32) (/* 10.0 3.0)
-> 3
: (scl 32) (format (/* 1.0 3.0) *Scl)
-> "0."
: (scl 32) (format (/* 10.0 3.0) *Scl)
-> "3."
: (scl 32) (format (/* 0.22 0.7) *Scl)
-> "0.31428571428571428571428571428571"
: (scl 32) (format (/* 9968.0 32.0) *Scl)
-> "311.5000"
: (scl 1) (bench (nil (format (/* 1.0 9967.0) *Scl)))
7.685 sec