initial math cleanup

db4
Doug Coleman 2008-10-03 02:19:03 -05:00
parent e7a785c58f
commit 56a0af9628
17 changed files with 101 additions and 80 deletions

View File

@ -1,8 +1,10 @@
! Copyright (c) 2007 Samuel Tardieu ! Copyright (c) 2007 Samuel Tardieu
! See http://factorcode.org/license.txt for BSD license. ! See http://factorcode.org/license.txt for BSD license.
USING: kernel math math.functions sequences ; USING: kernel math math.functions sequences fry ;
IN: math.algebra IN: math.algebra
: chinese-remainder ( aseq nseq -- x ) : chinese-remainder ( aseq nseq -- x )
dup product dup product
[ [ over / [ swap gcd drop ] keep * * ] curry 2map sum ] keep rem ; foldable [
'[ _ over / [ swap gcd drop ] keep * * ] 2map sum
] keep rem ; foldable

View File

@ -1,5 +1,7 @@
! Copyright (C) 2008 Doug Coleman, Slava Pestov.
! See http://factorcode.org/license.txt for BSD license.
USING: kernel math math.constants math.functions math.intervals USING: kernel math math.constants math.functions math.intervals
math.vectors namespaces sequences ; math.vectors namespaces sequences combinators.short-circuit ;
IN: math.analysis IN: math.analysis
<PRIVATE <PRIVATE
@ -20,8 +22,8 @@ IN: math.analysis
: (gamma-lanczos6) ( x -- log[gamma[x+1]] ) : (gamma-lanczos6) ( x -- log[gamma[x+1]] )
#! log(gamma(x+1) #! log(gamma(x+1)
dup 0.5 + dup gamma-g6 + dup >r log * r> - [ 0.5 + dup gamma-g6 + dup [ log * ] dip - ]
swap 6 gamma-z gamma-p6 v. log + ; [ 6 gamma-z gamma-p6 v. log ] bi + ;
: gamma-lanczos6 ( x -- gamma[x] ) : gamma-lanczos6 ( x -- gamma[x] )
#! gamma(x) = gamma(x+1) / x #! gamma(x) = gamma(x+1) / x
@ -39,7 +41,7 @@ PRIVATE>
: gamma ( x -- y ) : gamma ( x -- y )
#! gamma(x) = integral 0..inf [ t^(x-1) exp(-t) ] dt #! gamma(x) = integral 0..inf [ t^(x-1) exp(-t) ] dt
#! gamma(n+1) = n! for n > 0 #! gamma(n+1) = n! for n > 0
dup 0.0 <= over 1.0 mod zero? and [ dup { [ 0.0 <= ] [ 1.0 mod zero? ] } 1&& [
drop 1./0. drop 1./0.
] [ ] [
dup abs gamma-lanczos6 swap dup 0 > [ drop ] [ gamma-neg ] if dup abs gamma-lanczos6 swap dup 0 > [ drop ] [ gamma-neg ] if
@ -55,7 +57,7 @@ PRIVATE>
] if ; ] if ;
: nth-root ( n x -- y ) : nth-root ( n x -- y )
over 0 = [ "0th root is undefined" throw ] when >r recip r> swap ^ ; [ recip ] dip swap ^ ;
! Forth Scientific Library Algorithm #1 ! Forth Scientific Library Algorithm #1
! !

View File

@ -1,7 +1,7 @@
! Copyright (c) 2007, 2008 Slava Pestov, Doug Coleman, Aaron Schaefer. ! Copyright (c) 2007, 2008 Slava Pestov, Doug Coleman, Aaron Schaefer.
! See http://factorcode.org/license.txt for BSD license. ! See http://factorcode.org/license.txt for BSD license.
USING: assocs kernel math math.order math.ranges mirrors USING: assocs kernel math math.order math.ranges mirrors
namespaces make sequences sequences.lib sorting ; namespaces sequences sorting fry ;
IN: math.combinatorics IN: math.combinatorics
<PRIVATE <PRIVATE
@ -13,7 +13,7 @@ IN: math.combinatorics
2dup - dupd > [ dupd - ] when ; inline 2dup - dupd > [ dupd - ] when ; inline
! See this article for explanation of the factoradic-based permutation methodology: ! See this article for explanation of the factoradic-based permutation methodology:
! http://msdn2.microsoft.com/en-us/library/aa302371.aspx ! http://msdn2.microsoft.com/en-us/library/aa302371.aspx
: factoradic ( n -- factoradic ) : factoradic ( n -- factoradic )
0 [ over 0 > ] [ 1+ [ /mod ] keep swap ] [ ] produce reverse 2nip ; 0 [ over 0 > ] [ 1+ [ /mod ] keep swap ] [ ] produce reverse 2nip ;
@ -39,13 +39,10 @@ PRIVATE>
twiddle [ nPk ] keep factorial / ; twiddle [ nPk ] keep factorial / ;
: permutation ( n seq -- seq ) : permutation ( n seq -- seq )
tuck permutation-indices swap nths ; [ permutation-indices ] keep nths ;
: all-permutations ( seq -- seq ) : all-permutations ( seq -- seq )
[ [ length factorial ] keep '[ _ permutation ] map ;
[ length factorial ] keep [ permutation , ] curry each
] { } make ;
: inverse-permutation ( seq -- permutation ) : inverse-permutation ( seq -- permutation )
<enum> >alist sort-values keys ; <enum> >alist sort-values keys ;

View File

@ -19,4 +19,3 @@ IN: math.compare
: clamp ( a value b -- x ) : clamp ( a value b -- x )
min max ; min max ;

View File

@ -1,4 +1,3 @@
USING: kernel continuations combinators sequences math USING: kernel continuations combinators sequences math
math.order math.ranges accessors float-arrays ; math.order math.ranges accessors float-arrays ;
@ -7,11 +6,11 @@ IN: math.derivatives
TUPLE: state x func h err i j errt fac hh ans a done ; TUPLE: state x func h err i j errt fac hh ans a done ;
: largest-float ( -- x ) HEX: 7fefffffffffffff bits>double ; foldable : largest-float ( -- x ) HEX: 7fefffffffffffff bits>double ; foldable
: ntab ( -- val ) 8 ; : ntab ( -- val ) 8 ; inline
: con ( -- val ) 1.6 ; : con ( -- val ) 1.6 ; inline
: con2 ( -- val ) con con * ; : con2 ( -- val ) con con * ; inline
: big ( -- val ) largest-float ; : big ( -- val ) largest-float ; inline
: safe ( -- val ) 2.0 ; : safe ( -- val ) 2.0 ; inline
! Yes, this was ported from C code. ! Yes, this was ported from C code.
: a[i][i] ( state -- elt ) [ i>> ] [ i>> ] [ a>> ] tri nth nth ; : a[i][i] ( state -- elt ) [ i>> ] [ i>> ] [ a>> ] tri nth nth ;
@ -120,4 +119,4 @@ TUPLE: state x func h err i j errt fac hh ans a done ;
bi ; bi ;
: derivative ( x func -- m ) 0.01 2.0 (derivative) drop ; : derivative ( x func -- m ) 0.01 2.0 (derivative) drop ;
: derivative-func ( func -- der ) [ derivative ] curry ; : derivative-func ( func -- der ) [ derivative ] curry ;

View File

@ -11,8 +11,8 @@ TUPLE: erato limit bits latest ;
: ind ( n -- i ) : ind ( n -- i )
2/ 1- ; inline 2/ 1- ; inline
: is-prime ( n erato -- bool ) : is-prime ( n limit -- bool )
>r ind r> bits>> nth ; inline [ ind ] [ bits>> ] bi* nth ; inline
: indices ( n erato -- range ) : indices ( n erato -- range )
limit>> ind over 3 * ind swap rot <range> ; limit>> ind over 3 * ind swap rot <range> ;

View File

@ -9,7 +9,7 @@ IN: math.fft
: odd ( seq -- seq ) 2 group 1 <column> ; : odd ( seq -- seq ) 2 group 1 <column> ;
DEFER: fft DEFER: fft
: two ( seq -- seq ) fft 2 v/n dup append ; : two ( seq -- seq ) fft 2 v/n dup append ;
: omega ( n -- n ) recip -2 pi i* * * exp ; : omega ( n -- n' ) recip -2 pi i* * * exp ;
: twiddle ( seq -- seq ) dup length dup omega swap n^v v* ; : twiddle ( seq -- seq ) dup length dup omega swap n^v v* ;
: (fft) ( seq -- seq ) dup odd two twiddle swap even two v+ ; : (fft) ( seq -- seq ) dup odd two twiddle swap even two v+ ;
: fft ( seq -- seq ) dup length 1 = [ (fft) ] unless ; : fft ( seq -- seq ) dup length 1 = [ (fft) ] unless ;

View File

@ -1,3 +1,5 @@
! Copyright (C) 2008 Doug Coleman.
! See http://factorcode.org/license.txt for BSD license.
USING: combinators combinators.lib io locals kernel math USING: combinators combinators.lib io locals kernel math
math.functions math.ranges namespaces random sequences math.functions math.ranges namespaces random sequences
hashtables sets ; hashtables sets ;

View File

@ -1,11 +1,17 @@
! Copyright © 2008 Reginald Keith Ford II ! Copyright © 2008 Reginald Keith Ford II
! See http://factorcode.org/license.txt for BSD license.
! Newton's Method of approximating roots ! Newton's Method of approximating roots
USING: kernel math math.derivatives ; USING: kernel math math.derivatives ;
IN: math.newtons-method IN: math.newtons-method
<PRIVATE <PRIVATE
: newton-step ( x function -- x2 ) dupd [ call ] [ derivative ] 2bi / - ;
: newton-precision ( -- n ) 13 ; : newton-step ( x function -- x2 )
dupd [ call ] [ derivative ] 2bi / - ; inline
: newton-precision ( -- n ) 13 ; inline
PRIVATE> PRIVATE>
: newtons-method ( guess function -- x ) newton-precision [ [ newton-step ] keep ] times drop ;
: newtons-method ( guess function -- x )
newton-precision [ [ newton-step ] keep ] times drop ;

View File

@ -1,18 +1,20 @@
! Copyright (C) 2008 Doug Coleman.
! See http://factorcode.org/license.txt for BSD license.
USING: arrays kernel sequences namespaces make math math.ranges USING: arrays kernel sequences namespaces make math math.ranges
math.vectors vectors ; math.vectors vectors ;
IN: math.numerical-integration IN: math.numerical-integration
SYMBOL: num-steps 180 num-steps set-global SYMBOL: num-steps 180 num-steps set-global
: setup-simpson-range ( from to -- frange ) : setup-simpson-range ( from to -- frange )
2dup swap - num-steps get / <range> ; 2dup swap - num-steps get / <range> ;
: generate-simpson-weights ( seq -- seq ) : generate-simpson-weights ( seq -- seq )
[ { 1 4 }
{ 1 4 } % length 2 / 2 - { 2 4 } <repetition> concat % 1 , swap length 2 / 2 - { 2 4 } <repetition> concat
] { } make ; { 1 } 3append ;
: integrate-simpson ( from to f -- x ) : integrate-simpson ( from to f -- x )
>r setup-simpson-range r> [ setup-simpson-range dup ] dip
dupd map dup generate-simpson-weights map dup generate-simpson-weights
v. swap [ third ] keep first - 6 / * ; v. swap [ third ] keep first - 6 / * ;

View File

@ -1,3 +1,5 @@
! Copyright (C) 2008 Doug Coleman.
! See http://factorcode.org/license.txt for BSD license.
USING: arrays kernel sequences vectors math math.vectors USING: arrays kernel sequences vectors math math.vectors
namespaces make shuffle splitting sequences.lib math.order ; namespaces make shuffle splitting sequences.lib math.order ;
IN: math.polynomials IN: math.polynomials
@ -82,5 +84,5 @@ PRIVATE>
: polyval ( p x -- p[x] ) : polyval ( p x -- p[x] )
#! Evaluate a polynomial. #! Evaluate a polynomial.
>r dup length r> powers v. ; [ dup length ] dip powers v. ;

View File

@ -8,44 +8,45 @@ IN: math.primes
<PRIVATE <PRIVATE
: find-prime-miller-rabin ( n -- p ) : find-prime-miller-rabin ( n -- p )
dup miller-rabin [ 2 + find-prime-miller-rabin ] unless ; foldable dup miller-rabin [ 2 + find-prime-miller-rabin ] unless ; foldable
PRIVATE> PRIVATE>
: next-prime ( n -- p ) : next-prime ( n -- p )
dup 999983 < [ dup 999983 < [
primes-under-million [ natural-search drop 1+ ] keep nth primes-under-million [ natural-search drop 1+ ] keep nth
] [ ] [
next-odd find-prime-miller-rabin next-odd find-prime-miller-rabin
] if ; foldable ] if ; foldable
: prime? ( n -- ? ) : prime? ( n -- ? )
dup 1000000 < [ dup 1000000 < [
dup primes-under-million natural-search nip = dup primes-under-million natural-search nip =
] [ ] [
miller-rabin miller-rabin
] if ; foldable ] if ; foldable
: lprimes ( -- list ) : lprimes ( -- list )
0 primes-under-million seq>list 0 primes-under-million seq>list
1000003 [ 2 + find-prime-miller-rabin ] lfrom-by 1000003 [ 2 + find-prime-miller-rabin ] lfrom-by
lappend ; lappend ;
: lprimes-from ( n -- list ) : lprimes-from ( n -- list )
dup 3 < [ drop lprimes ] [ 1- next-prime [ next-prime ] lfrom-by ] if ; dup 3 < [ drop lprimes ] [ 1- next-prime [ next-prime ] lfrom-by ] if ;
: primes-upto ( n -- seq ) : primes-upto ( n -- seq )
{ {
{ [ dup 2 < ] [ drop { } ] } { [ dup 2 < ] [ drop { } ] }
{ [ dup 1000003 < ] { [ dup 1000003 < ] [
[ primes-under-million [ natural-search drop 1+ 0 swap ] keep <slice> ] } primes-under-million [ natural-search drop 1+ 0 swap ] keep <slice>
[ primes-under-million 1000003 lprimes-from ] }
rot [ <= ] curry lwhile list>array append ] [ primes-under-million 1000003 lprimes-from
} cond ; foldable rot [ <= ] curry lwhile list>array append ]
} cond ; foldable
: primes-between ( low high -- seq ) : primes-between ( low high -- seq )
primes-upto primes-upto
[ 1- next-prime ] dip [ 1- next-prime ] dip
[ natural-search drop ] keep [ length ] keep <slice> ; foldable [ natural-search drop ] keep [ length ] keep <slice> ; foldable
: coprime? ( a b -- ? ) gcd nip 1 = ; foldable : coprime? ( a b -- ? ) gcd nip 1 = ; foldable

View File

@ -28,7 +28,7 @@ PRIVATE>
: qconjugate ( u -- u' ) : qconjugate ( u -- u' )
#! Quaternion conjugate. #! Quaternion conjugate.
first2 neg >r conjugate r> 2array ; first2 [ conjugate ] [ neg ] bi* 2array ;
: qrecip ( u -- 1/u ) : qrecip ( u -- 1/u )
#! Quaternion inverse. #! Quaternion inverse.

View File

@ -1,14 +1,26 @@
! Copyright © 2008 Reginald Keith Ford II ! Copyright © 2008 Reginald Keith Ford II
! See http://factorcode.org/license.txt for BSD license.
! Secant Method of approximating roots ! Secant Method of approximating roots
USING: kernel math math.function-tools math.points math.vectors ; USING: kernel math math.function-tools math.points math.vectors ;
IN: math.secant-method IN: math.secant-method
<PRIVATE <PRIVATE
: secant-solution ( x1 x2 function -- solution ) [ eval ] curry bi@ linear-solution ;
: secant-step ( x1 x2 func -- x2 x3 func ) 2dup [ secant-solution ] 2dip swapd ; : secant-solution ( x1 x2 function -- solution )
: secant-precision ( -- n ) 15 ; [ eval ] curry bi@ linear-solution ;
: secant-step ( x1 x2 func -- x2 x3 func )
[ secant-solution ] 2keep swapd ;
: secant-precision ( -- n ) 15 ; inline
PRIVATE> PRIVATE>
: secant-method ( left right function -- x ) secant-precision [ secant-step ] times drop + 2 / ;
: secant-method ( left right function -- x )
secant-precision [ secant-step ] times drop + 2 / ;
! : close-enough? ( a b -- t/f ) - abs tiny-amount < ; ! : close-enough? ( a b -- t/f ) - abs tiny-amount < ;
! : secant-method2 ( left right function -- x ) 2over close-enough? [ drop average ] [ secant-step secant-method ] if ;
! : secant-method2 ( left right function -- x )
! 2over close-enough?
! [ drop average ] [ secant-step secant-method ] if ;

View File

@ -1,5 +1,7 @@
! Copyright (C) 2008 Doug Coleman, Michael Judge.
! See http://factorcode.org/license.txt for BSD license.
USING: kernel math math.analysis math.functions math.vectors sequences USING: kernel math math.analysis math.functions math.vectors sequences
sequences.lib sorting ; sequences.lib sorting ;
IN: math.statistics IN: math.statistics
: mean ( seq -- n ) : mean ( seq -- n )
@ -18,7 +20,7 @@ IN: math.statistics
: median ( seq -- n ) : median ( seq -- n )
#! middle number if odd, avg of two middle numbers if even #! middle number if odd, avg of two middle numbers if even
natural-sort dup length dup even? [ natural-sort dup length dup even? [
1- 2 / swap [ nth ] [ >r 1+ r> nth ] 2bi + 2 / 1- 2 / swap [ nth ] [ [ 1+ ] dip nth ] 2bi + 2 /
] [ ] [
2 / swap nth 2 / swap nth
] if ; ] if ;

View File

@ -1,8 +1,7 @@
! Copyright (c) 2007 Aaron Schaefer. ! Copyright (c) 2007 Aaron Schaefer.
! See http://factorcode.org/license.txt for BSD license. ! See http://factorcode.org/license.txt for BSD license.
USING: combinators.lib kernel math math.functions math.parser namespaces USING: combinators.lib kernel math math.functions math.parser namespaces
sequences splitting grouping sequences.lib sequences splitting grouping combinators.short-circuit ;
combinators.short-circuit ;
IN: math.text.english IN: math.text.english
<PRIVATE <PRIVATE
@ -86,14 +85,10 @@ SYMBOL: and-needed?
] if ; ] if ;
: (number>text) ( n -- str ) : (number>text) ( n -- str )
dup negative-text swap abs 3digit-groups recombine append ; [ negative-text ] [ abs 3digit-groups recombine ] bi append ;
PRIVATE> PRIVATE>
: number>text ( n -- str ) : number>text ( n -- str )
dup zero? [ dup zero? [ small-numbers ] [ [ (number>text) ] with-scope ] if ;
small-numbers
] [
[ (number>text) ] with-scope
] if ;

View File

@ -1,6 +1,6 @@
! Copyright (C) 2008 Eduardo Cavazos.
! See http://factorcode.org/license.txt for BSD license.
USING: math math.constants ; USING: math math.constants ;
IN: math.trig IN: math.trig
: deg>rad pi * 180 / ; inline : deg>rad pi * 180 / ; inline