factor/libs/math/utils.factor

IN: math-contrib
USING: errors kernel sequences math sequences-internals namespaces arrays ;

: deg>rad pi * 180 / ; inline
: rad>deg 180 * pi / ; inline

: lcm ( a b -- c )
    #! Smallest integer such that c/a and c/b are both integers.
    [ * ] 2keep gcd nip /i ; foldable

: mod-inv ( x n -- y )
    #! Compute the multiplicative inverse of x mod n.
    gcd 1 = [ "Non-trivial divisor found" throw ] unless ;
    foldable

: each-bit ( n quot -- )
    over zero? pick -1 number= or [
        2drop
    ] [
        2dup >r >r >r 1 bitand r> call r> -1 shift r> each-bit
    ] if ; inline

: (^mod) ( n z w -- z^w )
    1 swap [
        1 number= [ dupd * pick mod ] when >r sq over mod r>
    ] each-bit 2nip ; inline

: ^mod ( z w n -- z^w )
    #! Compute z^w mod n.
    over 0 < [
        [ >r neg r> ^mod ] keep mod-inv
    ] [
        -rot (^mod)
    ] if ; foldable

: powers ( n x -- seq )
    #! Output sequence has n elements, { 1 x x^2 x^3 ... }
    <array> 1 [ * ] accumulate nip ;

: ** ( u v -- u*v' ) conjugate * ; inline

: c. ( v v -- x )
    #! Complex inner product.
    0 [ ** + ] 2reduce ;

: proj ( u v -- w )
    #! Orthogonal projection of u onto v.
    [ [ v. ] keep norm-sq v/n ] keep n*v ;

: minmax ( seq -- min max )
    #! find the min and max of a seq in one pass
    1./0. -1./0. rot [ dup pick max -rot nip pick min -rot nip ] each ;

: absminmax ( seq -- min max )
    #! find the absolute values of the min and max of a seq in one pass
    minmax 2dup [ abs ] 2apply > [ swap ] when ;

SYMBOL: almost=-precision .000001 almost=-precision set-global
: almost= ( a b -- bool )
    - abs almost=-precision get < ;

TUPLE: frange from step length ;

C: frange ( from step to -- seq )
    #! example: 0 .01 10 <frange> >array
    >r pick - swap [ / ceiling 1+ >bignum ] keep r> 
    [ set-frange-step ] keep
    [ set-frange-length ] keep
    [ set-frange-from ] keep ;

: <frange-no-endpt> ( from step to -- seq )
    over - <frange> ;

M: frange length ( frange -- n ) frange-length ;
M: frange nth ( n frange -- obj ) [ frange-step * ] keep frange-from + ;

: frange-range ( frange -- range )
    [ frange-step ] keep frange-length 1- * ;

SYMBOL: step-size .01 step-size set-global  ! TODO: base on arguments
: (limit) ( count diff quot -- x quot )
    pick 10 > [ "Not converging fast enough" throw ] when
    [ call ] keep >r 2dup swap - 0 < [ "not converging" throw ] when
    2dup almost= rot drop r>
    swap [ step-size [ 2 / ] change rot 1+ -rot (limit) ] unless ;

: limit ( quot -- x )
    .1 step-size set [ call ] keep
    step-size [ 2 / ] change 0 -rot (limit) 2drop ;

: nth-rand ( seq -- elem ) [ length random-int ] keep nth ;