contrib/math/utils.factor

   1 IN: math-contrib
   2 USING: errors kernel sequences math sequences-internals namespaces arrays ;
   3
   4 : deg>rad pi * 180 / ; inline
   5 : rad>deg 180 * pi / ; inline
   6
   7 : lcm ( a b -- c )
   8     #! Smallest integer such that c/a and c/b are both integers.
   9     2dup gcd nip >r * r> /i ; foldable
  10
  11 : mod-inv ( x n -- y )
  12     #! Compute the multiplicative inverse of x mod n.
  13     gcd 1 = [ "Non-trivial divisor found" throw ] unless ;
  14     foldable
  15
  16 : each-bit ( n quot -- )
  17     over zero? pick -1 number= or [
  18         2drop
  19     ] [
  20         2dup >r >r >r 1 bitand r> call r> -1 shift r> each-bit
  21     ] if ; inline
  22
  23 : (^mod) ( n z w -- z^w )
  24     1 swap [
  25         1 number= [ dupd * pick mod ] when >r sq over mod r>
  26     ] each-bit 2nip ; inline
  27
  28 : ^mod ( z w n -- z^w )
  29     #! Compute z^w mod n.
  30     over 0 < [
  31         [ >r neg r> ^mod ] keep mod-inv
  32     ] [
  33         -rot (^mod)
  34     ] if ; foldable
  35
  36 : powers ( n x -- { 1 x x^2 x^3 ... } )
  37     #! Output sequence has n elements.
  38     <array> 1 [ * ] accumulate ;
  39
  40 : ** ( u v -- u*v' ) conjugate * ; inline
  41
  42 : c. ( v v -- x )
  43     #! Complex inner product.
  44     0 [ ** + ] 2reduce ;
  45
  46 : proj ( u v -- w )
  47     #! Orthogonal projection of u onto v.
  48     [ [ v. ] keep norm-sq v/n ] keep n*v ;
  49
  50 : minmax ( seq -- min max )
  51     #! find the min and max of a seq in one pass
  52     1./0. -1./0. rot [ dup pick max -rot nip pick min -rot nip ] each ;
  53
  54 : absminmax ( seq -- min max )
  55     #! find the absolute values of the min and max of a seq in one pass
  56     minmax 2dup [ abs ] 2apply > [ swap ] when ;
  57
  58 SYMBOL: almost=-precision .000001 almost=-precision set-global
  59 : almost= ( a b -- bool )
  60     - abs almost=-precision get < ;
  61
  62 TUPLE: frange from step length ;
  63
  64 C: frange ( from step to -- seq )
  65     #! example: 0 .01 10 <frange> >array
  66     >r pick - swap [ / ceiling 1+ ] keep -rot swapd r>
  67     [ set-frange-length ] keep
  68     [ set-frange-step ] keep
  69     [ set-frange-from ] keep ;
  70
  71 : decrement-length ( frange -- )
  72     [ frange-length 1- ] keep set-frange-length ;
  73
  74 : <frange-no-endpt> ( from step length -- seq )
  75     <frange> dup decrement-length ;
  76
  77 M: frange length ( frange -- n )
  78     frange-length ;
  79
  80 : increment-start ( frange -- )
  81     [ [ frange-from ] keep frange-step + ] keep set-frange-from ;
  82
  83 : frange-range ( frange -- range )
  84     [ frange-step ] keep frange-length 1- * ;
  85
  86 M: frange nth ( n frange -- obj )
  87     [ frange-step * ] keep frange-from + ;
  88
  89 ! : pivot ( left right index seq -- )
  90     ! [ nth ] keep [ exchange ] 3keep ;
  91
  92 SYMBOL: step-size .01 step-size set  ! base on arguments
  93 : (limit) ( count diff quot -- x quot )
  94     pick 10 > [ "Not converging fast enough" throw ] when
  95     [ call ] keep >r 2dup swap - 0 < [ "not converging" throw ] when
  96     2dup almost= rot drop r>
  97     swap [ step-size [ 2 / ] change rot 1+ -rot (limit) ] unless ;
  98
  99 : limit ( quot -- x )
 100     .1 step-size set [ call ] keep step-size [ 2 / ] change 0 -rot (limit) 2drop ;
 101
 102 : nth-rand ( seq -- elem ) [ length random-int ] keep nth ;
 103
 104 : count-end ( seq quot -- count )
 105     >r [ length ] keep r> find-last drop dup -1 = [ 2drop 0 ] [ - 1- ] if ;
 106