contrib/math cleanup, SIGFPE bug on intel mac (other platforms?)
erg
parent
2067553f3e
commit
672cdba238
|
|
@ -1,14 +1,6 @@
|
|||
USING: kernel math sequences namespaces io strings hashtables ;
|
||||
USING: kernel math math-contrib sequences namespaces io strings hashtables ;
|
||||
IN: crypto-internals
|
||||
|
||||
: (count-end) ( elt count seq -- elt count seq )
|
||||
2dup length < [
|
||||
3dup [ length swap - 1- ] keep nth = [ >r 1+ r> (count-end) ] when
|
||||
] when ;
|
||||
|
||||
: count-end ( elt seq -- n )
|
||||
#! count the number of elem at the end of the seq
|
||||
0 swap (count-end) drop nip ;
|
||||
|
||||
: ch>base64 ( ch -- ch )
|
||||
"ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/" nth ;
|
||||
|
|
|
|||
|
|
@ -22,7 +22,9 @@ USING: kernel io strings sequences namespaces math parser ;
|
|||
|
||||
: shift-mod ( n s w -- n ) >r shift r> 1 swap shift 1 - bitand ; inline
|
||||
|
||||
|
||||
IN: crypto
|
||||
|
||||
: bitroll ( n s w -- n )
|
||||
#! Roll n by s bits to the left, wrapping around after
|
||||
#! w bits.
|
||||
|
|
@ -31,4 +33,5 @@ IN: crypto
|
|||
[ shift-mod ] 3keep
|
||||
[ - ] keep shift-mod bitor ; inline
|
||||
|
||||
: hex-string ( str -- str ) [ [ >hex 2 48 pad-left % ] each ] "" make ;
|
||||
: hex-string ( str -- str )
|
||||
[ [ >hex 2 48 pad-left % ] each ] "" make ;
|
||||
|
|
|
|||
|
|
@ -6,9 +6,8 @@ IN: crypto
|
|||
#! a' = a * nextpowerof2(a) mod n
|
||||
>r dup next-power-of-2 * r> mod ;
|
||||
|
||||
: montgomery* ( a b -- a*b )
|
||||
"todo" throw
|
||||
;
|
||||
! : montgomery* ( a b -- a*b )
|
||||
! "todo" throw ;
|
||||
|
||||
: montgomery-r^2 ( n -- a )
|
||||
#! ans = r^2 mod n, where r = nextpowerof2(n)
|
||||
|
|
|
|||
|
|
@ -17,7 +17,7 @@ USING: kernel sequences errors namespaces math ;
|
|||
} ; inline
|
||||
|
||||
: gamma-z ( x n -- seq )
|
||||
[ [ over + 1.0 swap / , ] each ] { } make 1.0 0 pick set-nth nip ;
|
||||
[ + recip ] map-with 1.0 0 pick set-nth ;
|
||||
|
||||
: (gamma-lanczos6) ( x -- log[gamma[x+1]] )
|
||||
#! log(gamma(x+1)
|
||||
|
|
@ -55,7 +55,7 @@ IN: math-contrib
|
|||
dup abs gammaln-lanczos6 swap dup 0 > [ drop ] [ gamma-neg ] if
|
||||
] if ;
|
||||
|
||||
: nth-root ( n x -- )
|
||||
: nth-root ( n x -- y )
|
||||
over 0 = [ "0th root is undefined" throw ] when >r recip r> swap ^ ;
|
||||
|
||||
! Forth Scientific Library Algorithm #1
|
||||
|
|
|
|||
|
|
@ -5,73 +5,33 @@ USING: kernel sequences vectors math math-internals namespaces arrays ;
|
|||
! { 1 1 0 1 } -> 1 + x + x^3
|
||||
! { } -> 0
|
||||
|
||||
: 2length ( seq seq -- ) [ length ] 2apply ;
|
||||
|
||||
: zero-vector ( n -- vector ) 0 <array> >vector ;
|
||||
|
||||
: zero-pad ( n seq -- seq )
|
||||
#! extend seq by n zeros
|
||||
>r zero-vector r> swap append ;
|
||||
|
||||
: zero-pad-front ( n seq -- seq )
|
||||
>r zero-vector r> append ;
|
||||
|
||||
: nzero-pad ( n seq -- )
|
||||
#! extend seq by n zeros
|
||||
>r zero-vector r> swap nappend ;
|
||||
|
||||
: zero-extend ( n seq -- )
|
||||
#! extend seq to max(n,length) with 0s
|
||||
[ length ] keep -rot - swap nzero-pad ;
|
||||
|
||||
: 2zero-extend ( seq seq -- )
|
||||
2dup max-length [ swap zero-extend ] keep swap zero-extend ;
|
||||
|
||||
: pextend ( p p -- p p )
|
||||
#! make two polynomials the same length, if empty, make length 1
|
||||
[ >vector ] 2apply 2dup 2zero-extend ;
|
||||
|
||||
: unempty ( seq seq -- )
|
||||
#! turn { } into { 0 }
|
||||
dup empty? [ 1 swap zero-pad ] when ;
|
||||
|
||||
: 2unempty ( seq seq -- )
|
||||
[ unempty ] 2apply ;
|
||||
: pextend* ( p p n -- p p ) 0 [ pad-right swap ] 2keep pad-right swap ;
|
||||
: pextend ( p p -- p p ) 2dup [ length ] 2apply max pextend* ;
|
||||
: unempty ( seq -- seq ) dup empty? [ drop { 0 } ] when ;
|
||||
: 2unempty ( seq seq -- seq seq ) [ unempty ] 2apply ;
|
||||
|
||||
IN: math-contrib
|
||||
|
||||
: p= ( p p -- )
|
||||
pextend = ;
|
||||
: p= ( p p -- ? ) pextend = ;
|
||||
|
||||
: ptrim ( p -- p )
|
||||
>vector
|
||||
dup length 1 > [ dup peek 0 = [ dup pop drop ptrim ] when ] when ;
|
||||
dup length 1 > [
|
||||
0 over count-end >r dup length r> - dup zero? [ drop ] [ head ] if
|
||||
] when ;
|
||||
|
||||
: 2ptrim ( p -- p )
|
||||
[ ptrim ] 2apply ;
|
||||
|
||||
: p+ ( p p -- p )
|
||||
pextend v+ ;
|
||||
|
||||
: p- ( p p -- p )
|
||||
pextend v- ;
|
||||
|
||||
: n*p ( n p -- n*p )
|
||||
n*v ;
|
||||
: 2ptrim ( p p -- p p ) [ ptrim ] 2apply ;
|
||||
: p+ ( p p -- p ) pextend v+ ;
|
||||
: p- ( p p -- p ) pextend v- ;
|
||||
: n*p ( n p -- n*p ) n*v ;
|
||||
|
||||
! convolution
|
||||
: (conv*a)
|
||||
2dup swap length - rot zero-pad-front ;
|
||||
|
||||
: conv*a ( seq seq -- seq seq )
|
||||
2dup 2length + 1- (conv*a) reverse -rot (conv*a) swap ;
|
||||
|
||||
: conv*b ( seq -- seq )
|
||||
rot dup pop drop 1 zero-vector swap append -rot ;
|
||||
: pextend-conv ( p p -- p p )
|
||||
#! extend to: p_m + p_n - 1
|
||||
2dup [ length ] 2apply + 1- pextend* [ >vector ] 2apply ;
|
||||
|
||||
: p* ( p p -- p )
|
||||
#! Multiply two polynomials.
|
||||
2unempty conv*a [ 3dup -rot v* sum >r pick r> -rot set-nth conv*b ] repeat nip ;
|
||||
2unempty pextend-conv <reversed> dup length
|
||||
[ over length pick <slice> pick [ * ] 2map sum ] map 2nip reverse ;
|
||||
|
||||
: p-sq ( p -- p-sq )
|
||||
dup p* ;
|
||||
|
|
@ -85,17 +45,13 @@ IN: polynomials-internals
|
|||
#! divide the last two numbers in the sequences
|
||||
[ peek ] 2apply /i ;
|
||||
|
||||
: p/mod-setup
|
||||
2ptrim 2dup 2length - dup 1 < [ drop 1 ] when
|
||||
dup >r swap zero-pad-front pextend r> 1+ ;
|
||||
|
||||
: (p/mod)
|
||||
2dup /-last 2dup , n*p swapd p- dup pop drop swap pop-front ;
|
||||
|
||||
IN: math-contrib
|
||||
|
||||
: p/mod
|
||||
p/mod-setup [ [ (p/mod) ] times ] V{ } make
|
||||
pextend dup length 1- [ [ (p/mod) ] times ] V{ } make
|
||||
reverse nip swap 2ptrim pextend ;
|
||||
|
||||
: (pgcd) ( b a y x -- a d )
|
||||
|
|
|
|||
|
|
@ -14,13 +14,12 @@ USING: kernel math sequences ;
|
|||
#! positive reals only
|
||||
0 [ recip + ] reduce recip ;
|
||||
|
||||
! : number-sort [ - ] sort ;
|
||||
: median ( seq -- n )
|
||||
#! middle number if odd, avg of two middle numbers if even
|
||||
[ - ] sort dup length dup even? [
|
||||
1+ 2 /i dup 1- rot [ nth ] keep swapd nth + 2 /
|
||||
natural-sort dup length dup even? [
|
||||
1- 2 / swap [ nth ] 2keep >r 1+ r> nth + 2 /
|
||||
] [
|
||||
2 /i swap nth
|
||||
2 / swap nth
|
||||
] if ;
|
||||
|
||||
: range ( seq -- n )
|
||||
|
|
@ -29,7 +28,7 @@ USING: kernel math sequences ;
|
|||
|
||||
: var ( seq -- x )
|
||||
#! variance, normalize by N-1
|
||||
dup length 1- dup 0 = [
|
||||
dup length 1- dup zero? [
|
||||
0 2nip
|
||||
] [
|
||||
swap [ mean ] keep 0 [ pick - sq + ] reduce nip swap /
|
||||
|
|
|
|||
|
|
@ -2,21 +2,21 @@ IN: temporary
|
|||
USING: kernel math test sequences math-contrib ;
|
||||
|
||||
! Tests
|
||||
[ V{ 0 1 } ] [ { 0 1 0 0 } ptrim ] unit-test
|
||||
[ V{ 1 } ] [ { 1 0 0 } ptrim ] unit-test
|
||||
[ V{ 0 } ] [ { 0 } ptrim ] unit-test
|
||||
[ V{ 3 10 8 } ] [ { 1 2 } { 3 4 } p* ] unit-test
|
||||
[ V{ 3 10 8 } ] [ { 3 4 } { 1 2 } p* ] unit-test
|
||||
[ V{ 0 0 0 0 0 0 0 0 0 0 } ] [ { 0 0 0 } { 0 0 0 0 0 0 0 0 } p* ] unit-test
|
||||
[ V{ 0 1 } ] [ { 0 1 } { 1 } p* ] unit-test
|
||||
[ V{ 0 } ] [ { } { } p* ] unit-test
|
||||
[ V{ 0 } ] [ { 0 } { } p* ] unit-test
|
||||
[ V{ 0 } ] [ { } { 0 } p* ] unit-test
|
||||
[ V{ 0 0 0 } ] [ { 0 0 0 } { 0 0 0 } p+ ] unit-test
|
||||
[ V{ 0 0 0 } ] [ { 0 0 0 } { 0 0 0 } p- ] unit-test
|
||||
[ { 0 1 } ] [ { 0 1 0 0 } ptrim ] unit-test
|
||||
[ { 1 } ] [ { 1 0 0 } ptrim ] unit-test
|
||||
[ { 0 } ] [ { 0 } ptrim ] unit-test
|
||||
[ { 3 10 8 } ] [ { 1 2 } { 3 4 } p* ] unit-test
|
||||
[ { 3 10 8 } ] [ { 3 4 } { 1 2 } p* ] unit-test
|
||||
[ { 0 0 0 0 0 0 0 0 0 0 } ] [ { 0 0 0 } { 0 0 0 0 0 0 0 0 } p* ] unit-test
|
||||
[ { 0 1 } ] [ { 0 1 } { 1 } p* ] unit-test
|
||||
[ { 0 } ] [ { } { } p* ] unit-test
|
||||
[ { 0 } ] [ { 0 } { } p* ] unit-test
|
||||
[ { 0 } ] [ { } { 0 } p* ] unit-test
|
||||
[ { 0 0 0 } ] [ { 0 0 0 } { 0 0 0 } p+ ] unit-test
|
||||
[ { 0 0 0 } ] [ { 0 0 0 } { 0 0 0 } p- ] unit-test
|
||||
[ { 0 0 0 } ] [ 4 { 0 0 0 } n*p ] unit-test
|
||||
[ { 4 8 0 12 } ] [ 4 { 1 2 0 3 } n*p ] unit-test
|
||||
[ V{ 1 4 7 6 0 0 0 0 0 } ] [ { 1 2 3 0 0 0 } { 1 2 0 0 } p* ] unit-test
|
||||
[ { 1 4 7 6 0 0 0 0 0 } ] [ { 1 2 3 0 0 0 } { 1 2 0 0 } p* ] unit-test
|
||||
[ V{ 7 -2 1 } V{ -20 0 0 } ] [ { 1 1 1 1 } { 3 1 } p/mod ] unit-test
|
||||
[ V{ 0 0 } V{ 1 1 } ] [ { 1 1 } { 1 1 1 1 } p/mod ] unit-test
|
||||
[ V{ 1 0 1 } V{ 0 0 0 } ] [ { 1 1 1 1 } { 1 1 } p/mod ] unit-test
|
||||
|
|
|
|||
|
|
@ -4,6 +4,15 @@ USING: errors kernel sequences math sequences-internals namespaces arrays ;
|
|||
: deg>rad pi * 180 / ; inline
|
||||
: rad>deg 180 * pi / ; inline
|
||||
|
||||
: (count-end) ( elt count seq -- elt count seq )
|
||||
2dup length < [
|
||||
3dup [ length swap - 1- ] keep nth = [ >r 1+ r> (count-end) ] when
|
||||
] when ;
|
||||
|
||||
: count-end ( elt seq -- n )
|
||||
#! count the number of elem at the end of the seq
|
||||
0 swap (count-end) drop nip ;
|
||||
|
||||
: lcm ( a b -- c )
|
||||
#! Smallest integer such that c/a and c/b are both integers.
|
||||
2dup gcd nip >r * r> /i ; foldable
|
||||
|
|
@ -43,9 +52,6 @@ USING: errors kernel sequences math sequences-internals namespaces arrays ;
|
|||
#! Complex inner product.
|
||||
0 [ ** + ] 2reduce ;
|
||||
|
||||
: sum ( v -- n ) 0 [ + ] reduce ;
|
||||
: product ( v -- n ) 1 [ * ] reduce ;
|
||||
|
||||
: proj ( u v -- w )
|
||||
#! Orthogonal projection of u onto v.
|
||||
[ [ v. ] keep norm-sq v/n ] keep n*v ;
|
||||
|
|
@ -89,13 +95,8 @@ M: frange length ( frange -- n )
|
|||
M: frange nth ( n frange -- obj )
|
||||
[ frange-step * ] keep frange-from + ;
|
||||
|
||||
: 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 ;
|
||||
! [ nth ] keep [ exchange ] 3keep ;
|
||||
|
||||
SYMBOL: step-size .01 step-size set ! base on arguments
|
||||
: (limit) ( count diff quot -- x quot )
|
||||
|
|
@ -109,20 +110,16 @@ SYMBOL: step-size .01 step-size set ! base on arguments
|
|||
|
||||
! 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 ;
|
||||
>r group r> map ;
|
||||
|
||||
: 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 ;
|
||||
! : skip-map ( seq n quot -- seq )
|
||||
! >r
|
||||
|
||||
: nth-rand ( seq -- elem ) [ length random-int ] keep nth ;
|
||||
: nth-rand ( seq -- elem )
|
||||
[ length random-int ] keep nth ;
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue