For the math library

[Jeremy Dunn]

Thought I would drop some math functions on the board to hopefully drop into a math library that we can all build up over time.



The first group are bounds checking functions. I have added alternate names to the math operator symbols to provide an alternate naming that is consistent with the naming pattern of the other boolean functions in NewLISP for those who are picky about consistency.

Code Select Expand


;Command name alternates for naming consistency
(define-macro (less?)(apply < (args)))
(define-macro (greater?)(apply > (args)))
(define-macro (ge?)(apply >= (args)))
(define-macro (le?)(apply <= (args)))
(define-macro (e?)(apply = (args)))
(define-macro (ne?)(apply != (args)))

;These functions are shorthands for testing whether the number a lies within
; or without of various ranges of x and y
(define (lle?  x a y)(and (<  x a)(<= a y)))
(define (lel?  x a y)(and (<= x a)(<  a y)))
(define (gge?  x a y)(and (>  x a)(>= a y)))
(define (geg?  x a y)(and (>= x a)(>  a y)))
(define (lg?   x a y)(or  (<  a x)(>  a y)))
(define (lege? x a y)(or  (<= a x)(>= a y)))


And here are some functions for dealing with signs

Code Select Expand


;Strictly positive
(define (spos? x)(> x 0))

(define (pos? x)(>= x 0))

(define (neg? x)(<= x 0))

;Strictly negative
(define (sneg? x)(< x 0))


;Similar to the sgn function except that it returns 1 if x=0 rather than 0
(define (sgn2 x)(if (zero? x) 1 (sgn x)))

;Returns true if x and y have the same sign
(define (sgn=? x y)(= 1 (* (sgn2 x)(sgn2 y))))

;Returns y with the same sign as x
(define (sgncast x y)(if (sgn=? x y) y (sub y)))


And finally we have some functions for approximation

Code Select Expand


;; Sets the fuzz factor for functions that use it. (fuzzy) with no arguments resets *fuzz*
;; to its default value of 0.000001. If there is an integer argument then *fuzz* is set to
;; 10^n. If there is a real number as an argument then *fuzz* is set to that value.
(define (fuzzy (x 0.000001))
  (setq *fuzz* (if (integer? x) (pow 10 x) x)))

;Return true if x=y to within +-fz
(define (approx? x (y 0)(fz *fuzz*))
  (< (sub y fz) x (add y fz)))

;If x=y to within the fuzz factor then return y else return nil
(define (approx x (y 0)(fz *fuzz*))
  (if (approx? x y fz) y))


The fuzzy function must be called somewhere at the start of your program so that all functions that use the *fuzz* value will know what level of approximation to use. I find that this is very important to have fuzziness in geometric algorithms.



I will soon post my own rational arithmetic set of functions that is more complete than the current library.

[Lutz]

For the first group of these you don't need to write a macro, i.e. you could just say:

Code Select Expand

(constant (global 'less?) <)

> (less? 1 2)
true
>


less? will be just as fast as the original < . Whenever you use a new name you can use the above form. 'constant' protects your new name from accidental change and 'global' makes it work in contexts too, just like a built-in primitive.



Lutz

[Jeremy Dunn]

Thought I would revisit the GCD (greatest common divisor) algorithm again. In the code snippets on the website Eddie Rucker contributed this GCD algorithm:

Code Select Expand


(define (gcd_ a b)
  (let (r (% b a))
    (if (= r 0) a (gcd_ r a))))

(define-macro (gcd)(apply gcd_ (args) 2))


This is a simple algorithm not optimised for speed. A better one is

Code Select Expand


(define (gcd_ x y , r)
  (while (> y 0)
    (setq r (% x y))
    (if (> r (>> y 1))(setq r (- y r)))
    (setq x y y r)
  )
  x
)


An even better one is the right-shift binary algorithm. I tried to transcribe it from a C example out of my number theory book and came up with this

Code Select Expand


(define (gcd_ x y , g t)

  (setq g 0)

  (while (~ (& (| x y) 1))
    (setq x (>> x 1)
          y (>> y 1)
          g (+ g 1)
    )
  )

  (while (~ (& x 1))(setq x (>> x 1)))
  (while (~ (& y 1))(setq y (>> y 1)))

  (while (!= x y)
    (if (< x y)
        (setq t x x y y t))
    (setq x (- x y))
    (while (~ (& x 1))
        (setq x (>> x 1)))
  )

  (<< x g)
)


However, when I run this program it seems to get stuck in an infinite loop. The C-code I was trying to translate from was as follows:



NAT RSBGCD(NAT a, NAT b)

{

   NAT c, g, t;

/* find the even part of the gcd */



g = 0

while (!((a | b) & 1))

{

   a >>= 1;

   b >>= 1;

   g++;

}



/* remove extra factors of 2 */



while (!(a & 1)) a >>= 1;

while (!(b & 1)) b >>= 1;



/* find the odd part of the gcd */



while (a != b)

{

   if (a < b) { t = a; a = b; b = t; }

   a -= b;

   while (!(a & 1)) a >>= 1;

}

return(a << g);

}



Any eagle eyes out there spot where I screwed up?

[Lutz]

The infinite loop is probably here:

Code Select Expand

(while (~ (& x 1))(setq x (>> x 1)))
(while (~ (& y 1))(setq y (>> y 1)))
; or here
(while (~ (& (| x y) 1))


In the C-language 0 (zero) is taken as logical false and the C-while loop exists on 0, but in newLISP the number 0 (zero) is still a logical true, and ~ (squivel) is a binary negation (bit flip) not a logical. You could change to:

Code Select Expand

(while (= 0 (& x 1)) ...)
; or
(while (zero? (& x 1)) ...)


Lutz



ps: note that zero? is broken for floats in 8.9.0 (fixed after), but you are dealing with integers



ps: I edited the squivel out again and reversed logic after looking at the C code. The C code is shifting out all 0s.

[Jeremy Dunn]

Good call, this one seems to work

Code Select Expand


(define (gcd_ x y , g t)

  (setq g 0)

  (while (zero? (& (| x y) 1))
    (setq x (>> x 1)
          y (>> y 1)
          g (+ g 1)
    ))

  (while (zero? (& x 1))(setq x (>> x 1)))
  (while (zero? (& y 1))(setq y (>> y 1)))

  (while (!= x y)
    (if (< x y)
        (setq t x x y y t))
    (setq x (- x y))
    (while (zero? (& x 1))
        (setq x (>> x 1))))

  (<< x g)
)

(define-macro (gcd)(apply gcd_ (args) 2))


The author of my number theory book tested this algorithm on thousand digit numbers and found that it was at least twice as fast as the 2nd algorithm  I showed. And for good measure the least common multiple function is easily derived as

Code Select Expand


;Returns the least common multiple of two integers x and y.
(define (lcm_ x y)(* x (/ y (gcd x y))))
(define-macro (lcm)(apply lcm_ (args) 2))


It occurs to me that the << and >> operators could be improved by defaulting to 1 if they only have one argument i.e. (<< x) is the same as (<< x 1). Many algorithms could be further trimmed of fat since shifting by 1 bit is very common.

[Lutz]

occurs to me that the << and >> operators could be improved by defaulting to 1 if they only have one argument i.e. (<< x) is the same as (<< x 1).

Look for it in 8.9.4:

Code Select Expand


(>> num) ; (>> num 1)
(<< num) ; (<< num 1)
(div num) ; (div 1.0 num)


Lutz

[Jeremy Dunn]

Here we continue with math functions for performing vector operations of various kinds. The QUAD function is introduced because it is often useful to avoid the square root that a distance calculation entails. Often we can convert an algorithm using distance to one using QUAD instead. The MAG function uses QUAD to define the conventional magnitude of a vector. And the UNITV function gives an example of how to avoid the square root.

Code Select Expand


;Calculate the square of a number x
(define (sq x)(mul x x))

;Calculate the quadrance of a vector
(define (quad u)(if (vector? u)(apply add (map sq u))))

;Calculate the magnitude of a vector
(define (mag u)(sqrt (quad u)))

;Return true if u is a vector i.e a list of numbers
(define (vector? u)
  (and (list? u)
       (apply and (map number? u))))

;The unit vector of a vector u, the distance from origin to u divided into the coordinates of u
(define (unitv u)
  (if (vector? u)
      (scalar/ (vec* u u)(quad u))))

;Add two or more vectors together. It is assumed that all vectors are the same length.
(define (v+ u v)(map add u v))
(define-macro (vec+)(apply v+ (map eval (args)) 2))

;Vector subtraction. If there is only a single
;vector then the negation of all its coordinates is returned.
(define (v- u v)(map sub u v))
(define-macro (vec-)
  (if (= 1 (length (args)))
      (map sub (eval (args 0)))
      (apply v- (map eval (args)) 2)))

(define (v* u v)(map mul u v))
(define-macro (vec*)
  (apply v* (map eval (args)) 2))

(define-macro (vec/)
  (map div (eval (args 0))
           (if (> (length (args)) 2)
               (apply vec* (map eval (rest (args)))) ;we test to perform only one division
               (eval (args 1)))))

;Multiply a vector by a scalar. Vector first, scalar second
(define (scalar* v n)
  (map mul v (dup n (length v))))

;Divide a vector by a scalar. Vector first, scalar second
(define (scalar/ v n)
  (map div v (dup n (length v))))

;Vector Dot product of two vectors - produces a scalar
(define (dot u v)(apply add (vec* u v)))

;Vector Cross product of two vectors.
(define (cross u v)
  (setq u1 (u 0) u2 (u 1) u3 (u 2)
        v1 (v 0) v2 (v 1) v3 (v 2))
  (list (sub (mul u2 v3)(mul u3 v2))
        (sub (mul u3 v1)(mul u1 v3))
        (sub (mul u1 v2)(mul u2 v1))))

;Tensor product of three vectors
(define (tensor u v w)(vec* (dot u v) w))

;Scalar triple product
(define (striple u v w)(dot u (cross v w)))

;Vector triple product of three vectors
(define (vtriple u v w)(cross u (cross v w)))


[Jeremy Dunn]

Lutz,



Another couple of simple suggestions. Change the POW function so that it defaults to squaring its input i.e. (pow x) = (pow x 2). This would balance out the SQRT function.



Also, we have INC and DEC for x+1 and x-1 but we also need 1-x. I suggest calling this CED since it is the inverse of DEC.

[Jeremy Dunn]

Here is a variation of the FACTOR function that returns the prime factors paired with their exponent. I would like to suggest adding a flg argument onto the FACTOR function that when true returns the list in this fashion.

Code Select Expand


;;; Return a list of lists of the prime factors with the number of times that they appear
;;; (factor 280) -> (2 2 2 5 7)
;;; (primepowers 280) -> ((2 3)(5 1)(7 1))
(define (primepowers n , f pf cnts)
  (setq f    (factor n)
        pf   (unique f)
        cnts (count pf f)
  )
  (map (fn (x y)(list x y)) pf cnts)
)

[Fanda]

I have a wish: Could we speed up (pow x) ?

Code Select Expand

> (set 'N 1000000)
1000000
> (time (pow 5) N)
234
> (time (mul 5 5) N)
94


In nl-math.c, we could just change:

Code Select Expand

result = pow(result, 2.0);

to

Code Select Expand

result = result * result;


Thanks, Fanda

[Lutz]

These things are very hardware, compiler and OS dependent. My benchmarks show the opposite on OS X and the difference before and after the change is minimal:

Code Select Expand

BEFORE THE CHANGE
newLISP v.8.9.8 on BSD, execute 'newlisp -h' for more info.

>  (time (pow 5) 1000000)
188
>  (time (pow 5) 1000000)
191
> (time (mul 5 5) 1000000)
176
> (time (mul 5 5) 1000000)
177
>

newLISP v.8.9.8 on OSX UTF-8, execute 'newlisp -h' for more info.

> (time (pow 5) 1000000)
307
> (time (pow 5) 1000000)
310
> (time (mul 5 5) 1000000)
426
> (time (mul 5 5) 1000000)
423

AFTER CHANGING TO result = result * result;

newLISP v.8.9.8 on BSD, execute 'newlisp -h' for more info.

>  (time (pow 5) 1000000)
153
>  (time (pow 5) 1000000)
143
> (time (mul 5 5) 1000000)
192
> (time (mul 5 5) 1000000)
190
>

newLISP v.8.9.8 on OSX UTF-8, execute 'newlisp -h' for more info.

> (time (pow 5) 1000000)
328
> (time (pow 5) 1000000)
332
> (time (mul 5 5) 1000000)
443
> (time (mul 5 5) 1000000)
438
>


OSX on a Mac Mini PPC 1.42 Ghz and FreeBSD 4.9 on a 2 Ghz PC. What OS and hardware have you tested on?



Lutz

[Fanda]

Oh, wow! I wouldn't expect such differences between OSes and HW.



I am running Win XP on AMD 64 3.2 GHz.



I guess, it makes difference only for PCs.



Fanda

[Fanda]

PS: I also tested speed ratio:

Code Select Expand

> (set 'N 1000000)
1000000
> (div (time (pow 5) N) (time (mul 5 5) N))
2.329787234
> (div (time (pow 5) N) (time (mul 5 5) N))
2.319148936
> (div (time (pow 5) N) (time (mul 5 5) N))
2.35483871


Decision is yours, Lutz ;-)



Fanda

[Jeremy Dunn]

Here are some handy functions for the math toolbox

Code Select Expand


(define (divadd a b (c 2))(div (add a b) c))
(define (divsub a b (c 2))(div (sub a b) c))


These are functions for (a+b)/c and (a-b)/c. If c is not supplied it defaults to 2. I have found these quite useful.