
| =============================================== |
." | \\/ MATH \\/ | "
 
| use Pi defined in ` ./CoSy/Furniture.f

: fpi* fpi f* ;		| this is the monadic circle fn in APL

: Rad fpi 180. f/ _f ; 

: Tau Pi 2. _f *f ;

: Phi ( the golden ratio ) .5 5. f^ 1. f+ 2. f/  _f ; | at the x87 level then  _f 


| \/ | Matrix & Complex | \/ | 
 
: Im ( n -- n*n_IdentMat ) 1. 0. 2_f --abca rep take cL swap 2 _take take ;

: 1i 0. 1. 2_f cL  .. -1. _f *f reverse ,L ; | matrix form of imaginary unit 

| \/ | ===== |  Most basic pythagorean ( aka euclidian ) computations | ======== | \/ |

: dot *f +/ ; 		 : doti *i +/ ;	| in case you need integer
| f( -1 0 1 )f f( 1 2 3 )f dot  |>| 2.00 
: norm^2 .. dot ;
|   f( -1 0 1 )f norm^2 		|>| 2.00 
: norm norm^2 sqrtf ;


| \/ | misc  | \/ |
: gcd	( a b -- c) | Jack Browns recursive  greatest common divisor 
    dup if swap over mod gcd else drop then ;
 | From http://ronware.org/reva/wiki/index.php/Intermediate_Tutorial

: Fact >_ : _Fact _iota i1 +i ['] *i ./ ; | Factorial | 20191211
: !*! iota 1+i  ' *i ./ ; 	| 20211210 |

|  PascalTri ( n-1 -- n row Pascal triangle ) | 20190101 
: PascalTri i1 enc >a { a@ a@ -1 _at ['] + ': i1 braket cL /\a } nxtimes a> ;
|  33 PascalTri 	| That's the largest triangle before integer overflow .

| /\ MATH /\ |
| =============================================== |
| \/ | Stats | \/ |
| cf Arthur Whitney's K math functions | 
|  http://cosy.com/K/Math_AW.txt

: avgf  1p> +/ R@ rho i>f %f 1P> ; 

: var 1p> ^2f avgf R@ avgf ^2f -f 1P> ;
: dev var sqrtf ; 

: cov 2p> *f avgf LR@ ['] avgf on2 *f -f 2P> ;
: cor 2p> cov LR@ ['] dev on2 *f %f 2P> ; 

| 20 _iota i>f .. reverse cor 10 _i fmtnF
| need to deal w rounding . 

| Geometric mean | 20190106 
: gavg  1p> ['] *f ./ R@ rho i>f 1%f ^f 1P> ; 

| /\ | Stats | /\ |
| =============================================== |

| \ 32-bit PRNG XorShift algorithm by Michel Jean 
| https://www.facebook.com/groups/PROGRAMMINGFORTH/permalink/1991574754475882/
0 [IF]
code um*
mov ecx, edx 
pop eax 
mul ebx 
push eax 
mov ebx, edx 
mov edx, ecx
next,
end-code

Maybe useful? In a few lines, a pseudo-random generator very reliable :
\ 32-bit PRNG XorShift algorithm
variable seed
seed seed ! \ initialize seed with its address
: random32 ( n1 -- n2 ) \ return a number < n1
seed @
dup 13 lshift xor
dup 17 rshift xor
dup 5 lshift xor
seed !
seed @  ( swap mod ) UM* NIP 	| improvement by Johan Kotlinski 
 ;

[THEN]

needs random/gm
 
: _rand_ ( n -- a ) dup intVecInit >aux
   0 ?do rand aux@ i ii! loop aux> ; 
 
: _rand ( i n -- iv ) | n rands in  i iota raw 
   dup intVecInit >aux
   0 ?do dup rand um* nip aux@ i ii! loop drop aux> ; 

: rand ( i n -- iv )  2p L@ >_ R@ >_ _rand 2P> ; 	| n rands in  i iota . 

: randf $7FFFFFFF _i >a> swap rand a> ' i>f on2 %f ; | 20201027

0 [IF] | 20201027
: perm >_ : _perm ( n -- [ random permutation of n items ] )
  dup _iota >aux dup dup _rand 
  swap 0 ?do dup i ii@ aux@ i ii@
  swap aux@ swap ix xchg aux@ i ii! loop free aux> ; 
| not sure who in the Reva mail group offered this algo . Neat because it only
| requires 1 rand , but it's not complete and uniform . see
| http://math.stackexchange.com/questions/1003779/show-whether-this-algorithm-produces-a-uniform-random-permutation
[THEN]

: factors ( n _i -- factors ) .. i>f sqrtf f>i iota i1 +i .. --bac modi 0=i & at ;
| Returns smaller factors of n . See  Thu.Nov,20171130 
| eg : | 60 _i factors     |>| 1 2 3 4 5 6 

| Exponentiation by squaring
| See http://ronware.org/reva/viewtopic.php?id=341 
|  from | http://en.wikipedia.org/wiki/Exponentiating_by_squaring

0 [IF]    | Ruby version  
def power(x,n)
  result = 1
  while n.nonzero?
    if n.modulo(2).nonzero?
      result = result * x
      n = n-1
    else
      x = x*x
      n = n/2
    end
  end
  return result
end

| using  aux  stack 
: i^n ( i n -- i^n ) | converted from Ruby on
    | http://en.wikipedia.org/wiki/Exponentiating_by_squaring
  dup 0if 2drop 1 ;; then
  1 >aux  repeat 
   dup 2 mod
   if over aux> * >aux  1-
    else swap dup * swap 2/ 
   then 
   dup while   
   2drop aux> ;

[THEN]

 |  using  put  ( and --abcab )
 
: i^n ( i n -- i^n )
  dup 0if 2drop 1 ;; then
  1 -rot repeat 
   dup 2 mod
   if --abcab * 2 put 1-
   else swap dup * swap 2/ 
   then 
   dup while   
   2drop ;

." | /\\ MATH /\\ | " cr
