2005-10-23 19:07:16 -04:00
|
|
|
IN: math-contrib
|
2005-11-16 19:39:51 -05:00
|
|
|
USING: errors kernel sequences math sequences-internals namespaces arrays ;
|
2005-10-21 03:42:38 -04:00
|
|
|
|
|
|
|
: 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.
|
|
|
|
2dup gcd nip >r * r> /i ; foldable
|
|
|
|
|
|
|
|
: mod-inv ( x n -- y )
|
|
|
|
#! Compute the multiplicative inverse of x mod n.
|
|
|
|
gcd 1 = [ "Non-trivial divisor found" throw ] unless ;
|
|
|
|
foldable
|
|
|
|
|
2006-08-18 00:54:39 -04:00
|
|
|
: each-bit ( n quot -- )
|
2006-06-15 01:36:48 -04:00
|
|
|
over zero? pick -1 number= or [
|
|
|
|
2drop
|
|
|
|
] [
|
|
|
|
2dup >r >r >r 1 bitand r> call r> -1 shift r> each-bit
|
|
|
|
] if ; inline
|
|
|
|
|
2005-10-21 03:42:38 -04:00
|
|
|
: (^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
|
|
|
|
|
2006-01-17 03:08:47 -05:00
|
|
|
: powers ( n x -- { 1 x x^2 x^3 ... } )
|
2005-12-07 19:43:29 -05:00
|
|
|
#! Output sequence has n elements.
|
2006-01-17 03:08:47 -05:00
|
|
|
<array> 1 [ * ] accumulate ;
|
2005-12-07 19:43:29 -05:00
|
|
|
|
2005-10-21 03:42:38 -04:00
|
|
|
: ** ( u v -- u*v' ) conjugate * ; inline
|
|
|
|
|
|
|
|
: c. ( v v -- x )
|
|
|
|
#! Complex inner product.
|
|
|
|
0 [ ** + ] 2reduce ;
|
2005-11-07 20:26:32 -05:00
|
|
|
|
2006-01-12 00:34:56 -05:00
|
|
|
: sum ( v -- n ) 0 [ + ] reduce ;
|
|
|
|
: product ( v -- n ) 1 [ * ] reduce ;
|
|
|
|
|
2005-12-10 01:02:13 -05:00
|
|
|
: proj ( u v -- w )
|
|
|
|
#! Orthogonal projection of u onto v.
|
|
|
|
[ [ v. ] keep norm-sq v/n ] keep n*v ;
|
|
|
|
|
2005-11-16 19:39:51 -05:00
|
|
|
: minmax ( seq -- min max )
|
|
|
|
#! find the min and max of a seq in one pass
|
2006-01-27 10:17:32 -05:00
|
|
|
1./0. -1./0. rot [ dup pick max -rot nip pick min -rot nip ] each ;
|
2005-11-16 19:39:51 -05:00
|
|
|
|
|
|
|
: absminmax ( seq -- min max )
|
2005-11-17 04:31:36 -05:00
|
|
|
#! find the absolute values of the min and max of a seq in one pass
|
|
|
|
minmax 2dup [ abs ] 2apply > [ swap ] when ;
|
2005-11-16 19:39:51 -05:00
|
|
|
|
|
|
|
SYMBOL: almost=-precision .000001 almost=-precision set
|
|
|
|
: almost= ( a b -- bool )
|
|
|
|
2dup - abs almost=-precision get < [
|
|
|
|
2drop t
|
|
|
|
] [
|
|
|
|
2array absminmax dup almost=-precision get * >r - abs r>
|
|
|
|
dup 0 < [ >= ] [ <= ] if
|
|
|
|
] if ;
|
|
|
|
|
2005-11-07 20:26:32 -05:00
|
|
|
TUPLE: frange from step length ;
|
|
|
|
|
|
|
|
C: frange ( from step to -- seq )
|
|
|
|
#! example: 0 .01 10 <frange> >array
|
2005-11-16 19:39:51 -05:00
|
|
|
>r pick - swap [ / ceiling 1+ ] keep -rot swapd r>
|
2005-11-07 20:26:32 -05:00
|
|
|
[ set-frange-length ] keep
|
|
|
|
[ set-frange-step ] keep
|
|
|
|
[ set-frange-from ] keep ;
|
|
|
|
|
|
|
|
: decrement-length ( frange -- )
|
|
|
|
[ frange-length 1- ] keep set-frange-length ;
|
|
|
|
|
2005-11-16 19:39:51 -05:00
|
|
|
: <frange-no-endpt> ( from step length -- seq )
|
|
|
|
<frange> dup decrement-length ;
|
|
|
|
|
|
|
|
M: frange length ( frange -- n )
|
|
|
|
frange-length ;
|
|
|
|
|
2005-11-07 20:26:32 -05:00
|
|
|
: increment-start ( frange -- )
|
|
|
|
[ [ frange-from ] keep frange-step + ] keep set-frange-from ;
|
|
|
|
|
2005-11-09 17:48:55 -05:00
|
|
|
: frange-range ( frange -- range )
|
|
|
|
[ frange-step ] keep frange-length 1- * ;
|
|
|
|
|
2006-08-18 00:54:39 -04:00
|
|
|
M: frange nth ( n frange -- obj )
|
|
|
|
[ frange-step * ] keep frange-from + ;
|
2005-11-07 20:26:32 -05:00
|
|
|
|
2005-11-16 19:39:51 -05:00
|
|
|
: nseq-swap ( a b seq -- seq )
|
|
|
|
#! swap indices a,b in seq
|
|
|
|
3dup [ nth ] keep swapd [ nth ] keep
|
|
|
|
>r >r rot r> r> swapd set-nth -rot set-nth ;
|
|
|
|
|
|
|
|
! : pivot ( left right index seq -- )
|
|
|
|
! [ nth ] keep [ nseq-swap ] 3keep ;
|
|
|
|
|
|
|
|
SYMBOL: step-size .01 step-size set ! 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 ;
|
|
|
|
|
|
|
|
! take elements n at a time and apply the quotation, forming a new seq
|
|
|
|
: group-map ( seq n quot -- seq )
|
|
|
|
pick length pick /
|
|
|
|
[ [ >r pick pick r> -rot pick over * [ + ] keep swap rot <slice> pick call
|
|
|
|
, ] repeat ] { } make 2nip nip ;
|
|
|
|
|
|
|
|
: nths ( start n seq -- seq )
|
|
|
|
-rot pick length <frange-no-endpt> [ over nth ] map nip ;
|
|
|
|
|
|
|
|
! take a set of every nth element and apply the quotation, forming a new seq
|
|
|
|
! { 1 2 3 4 5 6 } 3 [ sum ] skip-map -> { 1 4 } { 2 5 } { 3 6 } -> { 5 7 9 }
|
|
|
|
: skip-map ( seq n quot -- seq )
|
|
|
|
pick length pick /mod
|
|
|
|
0 = [ "seq length must be a multiple of n" throw ] unless
|
|
|
|
1 <= [ "seq must be 2n or longer" throw ] when
|
|
|
|
over [ [ dup >r >r pick pick r> rot swapd nths over call , r> ] repeat ] { } make 2nip nip ;
|
2005-11-07 20:26:32 -05:00
|
|
|
|
2006-01-23 18:53:58 -05:00
|
|
|
: nth-rand ( seq -- elem ) [ length random-int ] keep nth ;
|
|
|
|
|