From 672cdba2383536ee3e71f83782ba11204d8653ba Mon Sep 17 00:00:00 2001 From: erg Date: Thu, 7 Sep 2006 23:29:13 +0000 Subject: [PATCH] contrib/math cleanup, SIGFPE bug on intel mac (other platforms?) --- contrib/crypto/base64.factor | 10 +--- contrib/crypto/common.factor | 5 +- contrib/crypto/montgomery.factor | 5 +- contrib/math/analysis.factor | 4 +- contrib/math/polynomials.factor | 84 ++++++++------------------------ contrib/math/statistics.factor | 9 ++-- contrib/math/test.factor | 26 +++++----- contrib/math/utils.factor | 33 ++++++------- 8 files changed, 61 insertions(+), 115 deletions(-) diff --git a/contrib/crypto/base64.factor b/contrib/crypto/base64.factor index fd3cdce583..01a7516112 100644 --- a/contrib/crypto/base64.factor +++ b/contrib/crypto/base64.factor @@ -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 ; diff --git a/contrib/crypto/common.factor b/contrib/crypto/common.factor index 7179a01049..30a2af4958 100644 --- a/contrib/crypto/common.factor +++ b/contrib/crypto/common.factor @@ -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 ; diff --git a/contrib/crypto/montgomery.factor b/contrib/crypto/montgomery.factor index 46e9cb8bef..1151c0e524 100644 --- a/contrib/crypto/montgomery.factor +++ b/contrib/crypto/montgomery.factor @@ -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) diff --git a/contrib/math/analysis.factor b/contrib/math/analysis.factor index 75cc732565..b9905b0d74 100644 --- a/contrib/math/analysis.factor +++ b/contrib/math/analysis.factor @@ -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 diff --git a/contrib/math/polynomials.factor b/contrib/math/polynomials.factor index 89f1eec7c8..8e9299d24b 100644 --- a/contrib/math/polynomials.factor +++ b/contrib/math/polynomials.factor @@ -5,74 +5,34 @@ 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 >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 dup length + [ over length pick 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 ) @@ -104,7 +60,7 @@ IN: math-contrib ] [ tuck p/mod >r pick p* swap >r swapd p- r> r> (pgcd) ] if ; - + : pgcd ( p p -- p ) swap V{ 0 } clone V{ 1 } clone 2swap (pgcd) ; diff --git a/contrib/math/statistics.factor b/contrib/math/statistics.factor index 8e1c6ef5cd..cfb00c015b 100644 --- a/contrib/math/statistics.factor +++ b/contrib/math/statistics.factor @@ -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 / diff --git a/contrib/math/test.factor b/contrib/math/test.factor index e4f7430cc6..ebf23e9bcb 100644 --- a/contrib/math/test.factor +++ b/contrib/math/test.factor @@ -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 diff --git a/contrib/math/utils.factor b/contrib/math/utils.factor index f61aa0b845..ffa74db832 100644 --- a/contrib/math/utils.factor +++ b/contrib/math/utils.factor @@ -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 pick call - , ] repeat ] { } make 2nip nip ; + >r group r> map ; : nths ( start n seq -- seq ) -rot pick length [ 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 ;