; *********************************************************************
;     Math Utilities
; *********************************************************************


; *********************************************************************
;  3-D vector arithmatic. Where a vector is a list of three numbers,
;                         and a basis is a list of three 3-D vectors.

make "-v [
   procedure [ [ :a ] ]
   output vscale :a -1 ]

make "vadd [
   procedure [ [ :a :b ] ]
   output ( list  + first :a first :b
                  + item 2 :a item 2 :b
                  + item 3 :a item 3 :b ) ]

make "vsub [
   procedure [ [ :a :b ] ]
   output ( list  - first :a first :b
                  - item 2 :a item 2 :b
                  - item 3 :a item 3 :b ) ]

make "vscale [
   procedure [ [ :v :s ] ]
   output ( list  * first :v :s
                  * item 2 :v :s
                  * item 3 :v :s ) ]

make "vdot [
   procedure [ [ :a :b ] ]
   output ( +  * first :a first :b
               * item 2 :a item 2 :b
               * item 3 :a item 3 :b ) ]

make "vcross [
   procedure [ [ :a :b ] ]
   output ( list  -  * item 2 :a item 3 :b
                     * item 3 :a item 2 :b
                  -  * item 3 :a item 1 :b
                     * item 1 :a item 3 :b
                  -  * item 1 :a item 2 :b
                     * item 2 :a item 1 :b ) ]

make "vmag [
   procedure [ [ :v ] ]
   output sqrt vdot :v :v ]

make "vrotate [
   procedure [ [ :v :pv :a ] ] 
   output vadd vscale :v cos :a vscale :pv sin :a ]

make "vbasis [
   procedure [ [ :vector :basis ] ]
   output ( list  vdot :vector item 1 :basis
                  vdot :vector item 2 :basis
                  vdot :vector item 3 :basis ) ]

make "rotate-xy [
   procedure [ [ :basis :angle ] ]
   output ( list  vrotate item 1 :basis item 2 :basis :angle
                  vrotate item 2 :basis -v item 1 :basis :angle
                  item 3 :basis ) ]

make "rotate-yz [
   procedure [ [ :basis :angle ] ]
   output ( list  item 1 :basis
                  vrotate item 2 :basis item 3 :basis :angle
                  vrotate item 3 :basis -v item 2 :basis :angle ) ]

make "rotate-zx [
   procedure [ [ :basis :angle ] ]
   output ( list  vrotate item 1 :basis -v item 3 :basis :angle
                  item 2 :basis
                  vrotate item 3 :basis item 1 :basis :angle ) ]


; *********************************************************************
;  Complex Numbers. Where a complex number is a list of two numbers.

make "cadd [
   procedure [ [ :a :b ] ]
   output   list  + first :a first :b
                  + item 2 :a item 2 :b ]

make "csub [
   procedure [ [ :a :b ] ]
   output   list  - first :a first :b
                  - item 2 :a item 2 :b ]

make "cmult [
   procedure [ [ :a :b ] ]
   output   list  -  * first :a first :b
                     * item 2 :a item 2 :b
                  +  * item 2 :a first :b
                     * first :a item 2 :b ]

make "cdiv [
   procedure [ [ :a :b ] [ ] [ :num :den :con ] ]
   make "con conjg :b
   make "den first cmult :b :con
   make "num cmult :a :con
   output list / first :num :den
               / item 2 :num :den ]

make "cconj [
   procedure [ [ :a ] ]
   output   list  first :a
                  +- item 2 :a ]

make "cabs [
   procedure [ [ :a ] ]
   output sqrt + * first :a first :a * item 2 :a item 2 :a ]

make "csqrt [
   procedure [ [ :a ] [ ] [ :x :y :w :r ] ]
   make "x abs first :a
   make "y abs item 2 :a
   if ( =0 :x :y )
   [  output [ 0 0 ] ]
   [  if >= :x :y
      [  make "r / :y :x
         make "w * sqrt :x sqrt * 0.5 + 1 sqrt + 1 * :r :r ]
      [  make "r / :x :y
         make "w * sqrt :y sqrt * 0.5 + :r sqrt + 1 * :r :r ]
      output if >= first :a 0
      [  list :w / item 2 :a * 2 :w ]
      [  make "r if >= item 2 :a 0 [ :w ] [ +- :w ]
         list  :r
               / item 2 :a * 2 :r ] ] ]

make "cscale [
   procedure [ [ :a :b ] ]
   output   list  * first :a :b
                  * item 2 :a :b ]


; *********************************************************************

; FACTORIAL ************************************************************
;    Factorial function.
;  ! num

make "! [
  procedure [ [ :n ] ]
  op if <= :n 1  [ 1 ]  [ * :n ! - :n 1 ] ]


; GCD ******************************************************************
;    Greatest common divisor of two numbers.
;  gcd n1 n2

make "gcd [
  procedure [ [ :a :b ] ]
  if =0 :b  [ op :a ]  [ op gcd :b remainder :a :b ] ]


; LCM ******************************************************************
;    Least common multiple of two numbers.
;  lcm n1 n2

make "lcm [
  procedure [ [ :a :b ] ]
  op * / :a gcd :a :b :b ]


; Useful numerical constants *******************************************

make "e  2.71828182845904523536
make "pi 3.14159265358979323846
make "phi / + 1 sqrt 5 2        ; golden ratio