Merge branch 'master' of git://projects.elasticdog.com/git/factor
Slava Pestov
commit
f0552accce
|
|
@ -15,7 +15,7 @@ IN: math.functions
|
|||
PRIVATE>
|
||||
|
||||
: rect> ( x y -- z )
|
||||
over real? over real? and [
|
||||
2dup [ real? ] both? [
|
||||
(rect>)
|
||||
] [
|
||||
"Complex number must have real components" throw
|
||||
|
|
@ -27,10 +27,10 @@ M: real sqrt
|
|||
>float dup 0.0 < [ neg fsqrt 0.0 swap rect> ] [ fsqrt ] if ;
|
||||
|
||||
: each-bit ( n quot: ( ? -- ) -- )
|
||||
over 0 = pick -1 = or [
|
||||
over [ 0 = ] [ -1 = ] bi or [
|
||||
2drop
|
||||
] [
|
||||
2dup >r >r >r odd? r> call r> 2/ r> each-bit
|
||||
2dup { [ odd? ] [ call ] [ 2/ ] [ each-bit ] } spread
|
||||
] if ; inline recursive
|
||||
|
||||
: map-bits ( n quot: ( ? -- obj ) -- seq )
|
||||
|
|
@ -69,8 +69,7 @@ PRIVATE>
|
|||
>rect [ >float ] bi@ ; inline
|
||||
|
||||
: >polar ( z -- abs arg )
|
||||
>float-rect [ [ sq ] bi@ + fsqrt ] [ swap fatan2 ] 2bi ;
|
||||
inline
|
||||
>float-rect [ [ sq ] bi@ + fsqrt ] [ swap fatan2 ] 2bi ; inline
|
||||
|
||||
: cis ( arg -- z ) dup fcos swap fsin rect> ; inline
|
||||
|
||||
|
|
@ -79,11 +78,10 @@ PRIVATE>
|
|||
<PRIVATE
|
||||
|
||||
: ^mag ( w abs arg -- magnitude )
|
||||
>r >r >float-rect swap r> swap fpow r> rot * fexp /f ;
|
||||
inline
|
||||
[ >float-rect swap ] [ swap fpow ] [ rot * fexp /f ] tri* ; inline
|
||||
|
||||
: ^theta ( w abs arg -- theta )
|
||||
>r >r >float-rect r> flog * swap r> * + ; inline
|
||||
[ >float-rect ] [ flog * swap ] [ * + ] tri* ; inline
|
||||
|
||||
: ^complex ( x y -- z )
|
||||
swap >polar [ ^mag ] [ ^theta ] 3bi polar> ; inline
|
||||
|
|
@ -106,18 +104,18 @@ PRIVATE>
|
|||
|
||||
: (^mod) ( n x y -- z )
|
||||
1 swap [
|
||||
[ dupd * pick mod ] when >r sq over mod r>
|
||||
[ dupd * pick mod ] when [ sq over mod ] dip
|
||||
] each-bit 2nip ; inline
|
||||
|
||||
: (gcd) ( b a x y -- a d )
|
||||
over zero? [
|
||||
2nip
|
||||
] [
|
||||
swap [ /mod >r over * swapd - r> ] keep (gcd)
|
||||
swap [ /mod [ over * swapd - ] dip ] keep (gcd)
|
||||
] if ;
|
||||
|
||||
: gcd ( x y -- a d )
|
||||
0 -rot 1 -rot (gcd) dup 0 < [ neg ] when ; foldable
|
||||
[ 0 1 ] 2dip (gcd) dup 0 < [ neg ] when ; foldable
|
||||
|
||||
: lcm ( a b -- c )
|
||||
[ * ] 2keep gcd nip /i ; foldable
|
||||
|
|
@ -131,7 +129,7 @@ PRIVATE>
|
|||
|
||||
: ^mod ( x y n -- z )
|
||||
over 0 < [
|
||||
[ >r neg r> ^mod ] keep mod-inv
|
||||
[ [ neg ] dip ^mod ] keep mod-inv
|
||||
] [
|
||||
-rot (^mod)
|
||||
] if ; foldable
|
||||
|
|
@ -141,14 +139,14 @@ GENERIC: absq ( x -- y ) foldable
|
|||
M: real absq sq ;
|
||||
|
||||
: ~abs ( x y epsilon -- ? )
|
||||
>r - abs r> < ;
|
||||
[ - abs ] dip < ;
|
||||
|
||||
: ~rel ( x y epsilon -- ? )
|
||||
>r [ - abs ] 2keep [ abs ] bi@ + r> * < ;
|
||||
[ [ - abs ] 2keep [ abs ] bi@ + ] dip * < ;
|
||||
|
||||
: ~ ( x y epsilon -- ? )
|
||||
{
|
||||
{ [ pick fp-nan? pick fp-nan? or ] [ 3drop f ] }
|
||||
{ [ 2over [ fp-nan? ] either? ] [ 3drop f ] }
|
||||
{ [ dup zero? ] [ drop number= ] }
|
||||
{ [ dup 0 < ] [ ~rel ] }
|
||||
[ ~abs ]
|
||||
|
|
|
|||
|
|
@ -12,10 +12,10 @@ SYMBOL: full-interval
|
|||
TUPLE: interval { from read-only } { to read-only } ;
|
||||
|
||||
: <interval> ( from to -- int )
|
||||
over first over first {
|
||||
2dup [ first ] bi@ {
|
||||
{ [ 2dup > ] [ 2drop 2drop empty-interval ] }
|
||||
{ [ 2dup = ] [
|
||||
2drop over second over second and
|
||||
2drop 2dup [ second ] both?
|
||||
[ interval boa ] [ 2drop empty-interval ] if
|
||||
] }
|
||||
[ 2drop interval boa ]
|
||||
|
|
@ -26,16 +26,16 @@ TUPLE: interval { from read-only } { to read-only } ;
|
|||
: closed-point ( n -- endpoint ) t 2array ;
|
||||
|
||||
: [a,b] ( a b -- interval )
|
||||
>r closed-point r> closed-point <interval> ; foldable
|
||||
[ closed-point ] dip closed-point <interval> ; foldable
|
||||
|
||||
: (a,b) ( a b -- interval )
|
||||
>r open-point r> open-point <interval> ; foldable
|
||||
[ open-point ] dip open-point <interval> ; foldable
|
||||
|
||||
: [a,b) ( a b -- interval )
|
||||
>r closed-point r> open-point <interval> ; foldable
|
||||
[ closed-point ] dip open-point <interval> ; foldable
|
||||
|
||||
: (a,b] ( a b -- interval )
|
||||
>r open-point r> closed-point <interval> ; foldable
|
||||
[ open-point ] dip closed-point <interval> ; foldable
|
||||
|
||||
: [a,a] ( a -- interval )
|
||||
closed-point dup <interval> ; foldable
|
||||
|
|
@ -51,11 +51,11 @@ TUPLE: interval { from read-only } { to read-only } ;
|
|||
: [-inf,inf] ( -- interval ) full-interval ; inline
|
||||
|
||||
: compare-endpoints ( p1 p2 quot -- ? )
|
||||
>r over first over first r> call [
|
||||
[ 2dup [ first ] bi@ ] dip call [
|
||||
2drop t
|
||||
] [
|
||||
over first over first = [
|
||||
swap second swap second not or
|
||||
2dup [ first ] bi@ = [
|
||||
[ second ] bi@ not or
|
||||
] [
|
||||
2drop f
|
||||
] if
|
||||
|
|
@ -86,7 +86,7 @@ TUPLE: interval { from read-only } { to read-only } ;
|
|||
] if ;
|
||||
|
||||
: (interval-op) ( p1 p2 quot -- p3 )
|
||||
[ [ first ] [ first ] [ ] tri* call ]
|
||||
[ [ first ] [ first ] [ call ] tri* ]
|
||||
[ drop [ second ] both? ]
|
||||
3bi 2array ; inline
|
||||
|
||||
|
|
@ -177,7 +177,7 @@ TUPLE: interval { from read-only } { to read-only } ;
|
|||
drop f
|
||||
] [
|
||||
interval>points
|
||||
2dup [ second ] bi@ and
|
||||
2dup [ second ] both?
|
||||
[ [ first ] bi@ = ]
|
||||
[ 2drop f ] if
|
||||
] if ;
|
||||
|
|
@ -193,9 +193,9 @@ TUPLE: interval { from read-only } { to read-only } ;
|
|||
dup [ interval>points [ first ] bi@ [a,b] ] when ;
|
||||
|
||||
: interval-integer-op ( i1 i2 quot -- i3 )
|
||||
>r 2dup
|
||||
[ interval>points [ first integer? ] both? ] both?
|
||||
r> [ 2drop [-inf,inf] ] if ; inline
|
||||
[
|
||||
2dup [ interval>points [ first integer? ] both? ] both?
|
||||
] dip [ 2drop [-inf,inf] ] if ; inline
|
||||
|
||||
: interval-shift ( i1 i2 -- i3 )
|
||||
#! Inaccurate; could be tighter
|
||||
|
|
@ -302,7 +302,7 @@ SYMBOL: incomparable
|
|||
2tri and and ;
|
||||
|
||||
: (interval<) ( i1 i2 -- i1 i2 ? )
|
||||
over from>> over from>> endpoint< ;
|
||||
2dup [ from>> ] bi@ endpoint< ;
|
||||
|
||||
: interval< ( i1 i2 -- ? )
|
||||
{
|
||||
|
|
@ -314,10 +314,10 @@ SYMBOL: incomparable
|
|||
} cond 2nip ;
|
||||
|
||||
: left-endpoint-<= ( i1 i2 -- ? )
|
||||
>r from>> r> to>> = ;
|
||||
[ from>> ] dip to>> = ;
|
||||
|
||||
: right-endpoint-<= ( i1 i2 -- ? )
|
||||
>r to>> r> from>> = ;
|
||||
[ to>> ] dip from>> = ;
|
||||
|
||||
: interval<= ( i1 i2 -- ? )
|
||||
{
|
||||
|
|
|
|||
|
|
@ -126,7 +126,7 @@ SYMBOL: fast-math-ops
|
|||
|
||||
: math-method* ( word left right -- quot )
|
||||
3dup math-op
|
||||
[ >r 3drop r> 1quotation ] [ drop math-method ] if ;
|
||||
[ [ 3drop ] dip 1quotation ] [ drop math-method ] if ;
|
||||
|
||||
: math-both-known? ( word left right -- ? )
|
||||
3dup math-op
|
||||
|
|
@ -157,13 +157,13 @@ SYMBOL: fast-math-ops
|
|||
] bi@ append ;
|
||||
|
||||
: each-derived-op ( word quot -- )
|
||||
>r derived-ops r> each ; inline
|
||||
[ derived-ops ] dip each ; inline
|
||||
|
||||
: each-fast-derived-op ( word quot -- )
|
||||
>r fast-derived-ops r> each ; inline
|
||||
[ fast-derived-ops ] dip each ; inline
|
||||
|
||||
: each-integer-derived-op ( word quot -- )
|
||||
>r integer-derived-ops r> each ; inline
|
||||
[ integer-derived-ops ] dip each ; inline
|
||||
|
||||
[
|
||||
[
|
||||
|
|
|
|||
|
|
@ -8,7 +8,7 @@ TUPLE: range
|
|||
{ step read-only } ;
|
||||
|
||||
: <range> ( a b step -- range )
|
||||
>r over - r>
|
||||
[ over - ] dip
|
||||
[ / 1+ 0 max >integer ] keep
|
||||
range boa ; inline
|
||||
|
||||
|
|
|
|||
|
|
@ -12,10 +12,10 @@ IN: math.ratios
|
|||
dup 1 number= [ drop ] [ <ratio> ] if ; inline
|
||||
|
||||
: scale ( a/b c/d -- a*d b*c )
|
||||
2>fraction >r * swap r> * swap ; inline
|
||||
2>fraction [ * swap ] dip * swap ; inline
|
||||
|
||||
: ratio+d ( a/b c/d -- b*d )
|
||||
denominator swap denominator * ; inline
|
||||
[ denominator ] bi@ * ; inline
|
||||
|
||||
PRIVATE>
|
||||
|
||||
|
|
@ -24,7 +24,7 @@ M: integer /
|
|||
"Division by zero" throw
|
||||
] [
|
||||
dup 0 < [ [ neg ] bi@ ] when
|
||||
2dup gcd nip tuck /i >r /i r> fraction>
|
||||
2dup gcd nip tuck /i [ /i ] dip fraction>
|
||||
] if ;
|
||||
|
||||
M: ratio hashcode*
|
||||
|
|
@ -52,7 +52,7 @@ M: ratio >= scale >= ;
|
|||
|
||||
M: ratio + 2dup scale + -rot ratio+d / ;
|
||||
M: ratio - 2dup scale - -rot ratio+d / ;
|
||||
M: ratio * 2>fraction * >r * r> / ;
|
||||
M: ratio * 2>fraction * [ * ] dip / ;
|
||||
M: ratio / scale / ;
|
||||
M: ratio /i scale /i ;
|
||||
M: ratio /f scale /f ;
|
||||
|
|
|
|||
|
|
@ -34,7 +34,7 @@ HELP: n*v
|
|||
{ $description "Multiplies each element of " { $snippet "u" } " by " { $snippet "n" } "." } ;
|
||||
|
||||
HELP: v*n
|
||||
{ $values { "n" "a number" } { "u" "a sequence of numbers" } { "v" "a sequence of numbers" } }
|
||||
{ $values { "u" "a sequence of numbers" } { "n" "a number" } { "v" "a sequence of numbers" } }
|
||||
{ $description "Multiplies each element of " { $snippet "u" } " by " { $snippet "n" } "." } ;
|
||||
|
||||
HELP: n/v
|
||||
|
|
|
|||
|
|
@ -25,7 +25,7 @@ IN: math.vectors
|
|||
: normalize ( u -- v ) dup norm v/n ;
|
||||
|
||||
: set-axis ( u v axis -- w )
|
||||
[ >r zero? 2over ? r> swap nth ] map-index 2nip ;
|
||||
[ [ zero? 2over ? ] dip swap nth ] map-index 2nip ;
|
||||
|
||||
HINTS: vneg { array } ;
|
||||
HINTS: norm-sq { array } ;
|
||||
|
|
|
|||
|
|
@ -90,7 +90,6 @@ HELP: derivative-func
|
|||
" [ cos ]"
|
||||
" bi - abs"
|
||||
"] map minmax"
|
||||
|
||||
}
|
||||
}
|
||||
} ;
|
||||
|
|
@ -100,4 +99,5 @@ ARTICLE: "derivatives" "The Derivative Toolkit"
|
|||
{ $subsection derivative }
|
||||
{ $subsection derivative-func }
|
||||
{ $subsection (derivative) } ;
|
||||
|
||||
ABOUT: "derivatives"
|
||||
|
|
|
|||
|
|
@ -0,0 +1,94 @@
|
|||
USING: help.markup help.syntax math sequences ;
|
||||
IN: math.polynomials
|
||||
|
||||
ARTICLE: "polynomials" "Polynomials"
|
||||
"A polynomial is a vector with the highest powers on the right:"
|
||||
{ $code "{ 1 1 0 1 } -> 1 + x + x^3" "{ } -> 0" }
|
||||
"Numerous words are defined to help with polynomial arithmetic:"
|
||||
{ $subsection p= }
|
||||
{ $subsection p+ }
|
||||
{ $subsection p- }
|
||||
{ $subsection p* }
|
||||
{ $subsection p-sq }
|
||||
{ $subsection powers }
|
||||
{ $subsection n*p }
|
||||
{ $subsection p/mod }
|
||||
{ $subsection pgcd }
|
||||
{ $subsection polyval }
|
||||
{ $subsection pdiff }
|
||||
{ $subsection pextend-conv }
|
||||
{ $subsection ptrim }
|
||||
{ $subsection 2ptrim } ;
|
||||
|
||||
ABOUT: "polynomials"
|
||||
|
||||
HELP: powers
|
||||
{ $values { "n" integer } { "x" number } { "seq" sequence } }
|
||||
{ $description "Output a sequence having " { $snippet "n" } " elements in the format: " { $snippet "{ 1 x x^2 x^3 ... }" } "." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "4 2 powers ." "{ 1 2 4 8 }" } } ;
|
||||
|
||||
HELP: p=
|
||||
{ $values { "p" "a polynomial" } { "q" "a polynomial" } { "?" "a boolean" } }
|
||||
{ $description "Tests if two polynomials are equal." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 0 1 } { 0 1 0 } p= ." "t" } } ;
|
||||
|
||||
HELP: ptrim
|
||||
{ $values { "p" "a polynomial" } { "p" "a polynomial" } }
|
||||
{ $description "Trims excess zeros from a polynomial." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 0 1 0 0 } ptrim ." "{ 0 1 }" } } ;
|
||||
|
||||
HELP: 2ptrim
|
||||
{ $values { "p" "a polynomial" } { "q" "a polynomial" } { "p" "a polynomial" } { "q" "a polynomial" } }
|
||||
{ $description "Trims excess zeros from two polynomials." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 0 1 0 0 } { 1 0 0 } 2ptrim swap . ." "{ 0 1 }\n{ 1 }" } } ;
|
||||
|
||||
HELP: p+
|
||||
{ $values { "p" "a polynomial" } { "q" "a polynomial" } { "r" "a polynomial" } }
|
||||
{ $description "Adds " { $snippet "p" } " and " { $snippet "q" } " component-wise." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 1 0 1 } { 0 1 } p+ ." "{ 1 1 1 }" } } ;
|
||||
|
||||
HELP: p-
|
||||
{ $values { "p" "a polynomial" } { "q" "a polynomial" } { "r" "a polynomial" } }
|
||||
{ $description "Subtracts " { $snippet "q" } " from " { $snippet "p" } " component-wise." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 1 1 1 } { 0 1 } p- ." "{ 1 0 1 }" } } ;
|
||||
|
||||
HELP: n*p
|
||||
{ $values { "n" number } { "p" "a polynomial" } { "n*p" "a polynomial" } }
|
||||
{ $description "Multiplies each element of " { $snippet "p" } " by " { $snippet "n" } "." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "4 { 3 0 1 } n*p ." "{ 12 0 4 }" } } ;
|
||||
|
||||
HELP: pextend-conv
|
||||
{ $values { "p" "a polynomial" } { "q" "a polynomial" } { "p" "a polynomial" } { "q" "a polynomial" } }
|
||||
{ $description "Convulution, extending to " { $snippet "p_m + q_n - 1" } "." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 1 0 1 } { 0 1 } pextend-conv swap . ." "V{ 1 0 1 0 }\nV{ 0 1 0 0 }" } } ;
|
||||
|
||||
HELP: p*
|
||||
{ $values { "p" "a polynomial" } { "q" "a polynomial" } { "r" "a polynomial" } }
|
||||
{ $description "Multiplies two polynomials." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 1 2 3 0 0 0 } { 1 2 0 0 } p* ." "{ 1 4 7 6 0 0 0 0 0 }" } } ;
|
||||
|
||||
HELP: p-sq
|
||||
{ $values { "p" "a polynomial" } { "p^2" "a polynomial" } }
|
||||
{ $description "Squares a polynomial." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 1 2 0 } p-sq ." "{ 1 4 4 0 0 }" } } ;
|
||||
|
||||
HELP: p/mod
|
||||
{ $values { "p" "a polynomial" } { "q" "a polynomial" } { "z" "a polynomial" } { "w" "a polynomial" } }
|
||||
{ $description "Computes to quotient " { $snippet "z" } " and remainder " { $snippet "w" } " of dividing " { $snippet "p" } " by " { $snippet "q" } "." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 1 1 1 1 } { 3 1 } p/mod swap . ." "V{ 7 -2 1 }\nV{ -20 0 0 }" } } ;
|
||||
|
||||
HELP: pgcd
|
||||
{ $values { "p" "a polynomial" } { "q" "a polynomial" } { "a" "a polynomial" } { "d" "a polynomial" } }
|
||||
{ $description "Computes the greatest common divisor " { $snippet "d" } " of " { $snippet "p" } " and " { $snippet "q" } ", and another value " { $snippet "a" } " satisfying:" { $code "a*q = d mod p" } }
|
||||
{ $notes "GCD in the case of polynomials is a monic polynomial of the highest possible degree that divides into both " { $snippet "p" } " and " { $snippet "q" } "." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 1 1 1 1} { 1 1 } pgcd swap . ." "{ 0 0 }\n{ 1 1 }" } } ;
|
||||
|
||||
HELP: pdiff
|
||||
{ $values { "p" "a polynomial" } { "p'" "a polynomial" } }
|
||||
{ $description "Finds the derivative of " { $snippet "p" } "." } ;
|
||||
|
||||
HELP: polyval
|
||||
{ $values { "p" "a polynomial" } { "x" number } { "p[x]" number } }
|
||||
{ $description "Evaluate " { $snippet "p" } " with the input " { $snippet "x" } "." }
|
||||
{ $examples { $example "USING: math.polynomials prettyprint ;" "{ 1 0 1 } 2 polyval ." "5" } } ;
|
||||
|
||||
|
|
@ -1,7 +1,6 @@
|
|||
IN: math.polynomials.tests
|
||||
USING: kernel math math.polynomials tools.test ;
|
||||
IN: math.polynomials.tests
|
||||
|
||||
! Tests
|
||||
[ { 0 1 } ] [ { 0 1 0 0 } ptrim ] unit-test
|
||||
[ { 1 } ] [ { 1 0 0 } ptrim ] unit-test
|
||||
[ { 0 } ] [ { 0 } ptrim ] unit-test
|
||||
|
|
|
|||
|
|
@ -4,46 +4,38 @@ USING: arrays kernel make math math.order math.vectors sequences shuffle
|
|||
splitting vectors ;
|
||||
IN: math.polynomials
|
||||
|
||||
! Polynomials are vectors with the highest powers on the right:
|
||||
! { 1 1 0 1 } -> 1 + x + x^3
|
||||
! { } -> 0
|
||||
|
||||
: powers ( n x -- seq )
|
||||
#! Output sequence has n elements, { 1 x x^2 x^3 ... }
|
||||
<array> 1 [ * ] accumulate nip ;
|
||||
|
||||
<PRIVATE
|
||||
|
||||
: 2pad-left ( p p n -- p p ) [ 0 pad-left ] curry bi@ ;
|
||||
: 2pad-right ( p p n -- p p ) [ 0 pad-right ] curry bi@ ;
|
||||
: pextend ( p p -- p p ) 2dup [ length ] bi@ max 2pad-right ;
|
||||
: pextend-left ( p p -- p p ) 2dup [ length ] bi@ max 2pad-left ;
|
||||
: 2pad-left ( p q n -- p q ) [ 0 pad-left ] curry bi@ ;
|
||||
: 2pad-right ( p q n -- p q ) [ 0 pad-right ] curry bi@ ;
|
||||
: pextend ( p q -- p q ) 2dup [ length ] bi@ max 2pad-right ;
|
||||
: pextend-left ( p q -- p q ) 2dup [ length ] bi@ max 2pad-left ;
|
||||
: unempty ( seq -- seq ) [ { 0 } ] when-empty ;
|
||||
: 2unempty ( seq seq -- seq seq ) [ unempty ] bi@ ;
|
||||
|
||||
PRIVATE>
|
||||
|
||||
: p= ( p p -- ? ) pextend = ;
|
||||
: powers ( n x -- seq )
|
||||
<array> 1 [ * ] accumulate nip ;
|
||||
|
||||
: p= ( p q -- ? ) pextend = ;
|
||||
|
||||
: ptrim ( p -- p )
|
||||
dup length 1 = [ [ zero? ] trim-right ] unless ;
|
||||
|
||||
: 2ptrim ( p p -- p p ) [ ptrim ] bi@ ;
|
||||
: p+ ( p p -- p ) pextend v+ ;
|
||||
: p- ( p p -- p ) pextend v- ;
|
||||
: 2ptrim ( p q -- p q ) [ ptrim ] bi@ ;
|
||||
: p+ ( p q -- r ) pextend v+ ;
|
||||
: p- ( p q -- r ) pextend v- ;
|
||||
: n*p ( n p -- n*p ) n*v ;
|
||||
|
||||
! convolution
|
||||
: pextend-conv ( p p -- p p )
|
||||
#! extend to: p_m + p_n - 1
|
||||
: pextend-conv ( p q -- p q )
|
||||
2dup [ length ] bi@ + 1- 2pad-right [ >vector ] bi@ ;
|
||||
|
||||
: p* ( p p -- p )
|
||||
#! Multiply two polynomials.
|
||||
: p* ( p q -- r )
|
||||
2unempty pextend-conv <reversed> dup length
|
||||
[ over length pick <slice> pick [ * ] 2map sum ] map 2nip reverse ;
|
||||
|
||||
: p-sq ( p -- p-sq )
|
||||
: p-sq ( p -- p^2 )
|
||||
dup p* ;
|
||||
|
||||
<PRIVATE
|
||||
|
|
@ -66,10 +58,12 @@ PRIVATE>
|
|||
|
||||
PRIVATE>
|
||||
|
||||
: p/mod ( a b -- / mod )
|
||||
: p/mod ( p q -- z w )
|
||||
p/mod-setup [ [ (p/mod) ] times ] V{ } make
|
||||
reverse nip swap 2ptrim pextend ;
|
||||
|
||||
<PRIVATE
|
||||
|
||||
: (pgcd) ( b a y x -- a d )
|
||||
dup V{ 0 } clone p= [
|
||||
drop nip
|
||||
|
|
@ -77,14 +71,14 @@ PRIVATE>
|
|||
tuck p/mod [ pick p* swap [ swapd p- ] dip ] dip (pgcd)
|
||||
] if ;
|
||||
|
||||
: pgcd ( p p -- p q )
|
||||
PRIVATE>
|
||||
|
||||
: pgcd ( p q -- a d )
|
||||
swap V{ 0 } clone V{ 1 } clone 2swap (pgcd) [ >array ] bi@ ;
|
||||
|
||||
: pdiff ( p -- p' )
|
||||
#! Polynomial derivative.
|
||||
dup length v* { 0 } ?head drop ;
|
||||
|
||||
: polyval ( p x -- p[x] )
|
||||
#! Evaluate a polynomial.
|
||||
[ dup length ] dip powers v. ;
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,46 @@
|
|||
USING: help.markup help.syntax math math.vectors vectors ;
|
||||
IN: math.quaternions
|
||||
|
||||
HELP: q*
|
||||
{ $values { "u" "a quaternion" } { "v" "a quaternion" } { "u*v" "a quaternion" } }
|
||||
{ $description "Multiply quaternions." }
|
||||
{ $examples { $example "USING: math.quaternions prettyprint ;" "{ C{ 0 1 } 0 } { 0 1 } q* ." "{ 0 C{ 0 1 } }" } } ;
|
||||
|
||||
HELP: qconjugate
|
||||
{ $values { "u" "a quaternion" } { "u'" "a quaternion" } }
|
||||
{ $description "Quaternion conjugate." } ;
|
||||
|
||||
HELP: qrecip
|
||||
{ $values { "u" "a quaternion" } { "1/u" "a quaternion" } }
|
||||
{ $description "Quaternion inverse." } ;
|
||||
|
||||
HELP: q/
|
||||
{ $values { "u" "a quaternion" } { "v" "a quaternion" } { "u/v" "a quaternion" } }
|
||||
{ $description "Divide quaternions." }
|
||||
{ $examples { $example "USING: math.quaternions prettyprint ;" "{ 0 C{ 0 1 } } { 0 1 } q/ ." "{ C{ 0 1 } 0 }" } } ;
|
||||
|
||||
HELP: q*n
|
||||
{ $values { "q" "a quaternion" } { "n" number } { "q" "a quaternion" } }
|
||||
{ $description "Multiplies each element of " { $snippet "q" } " by " { $snippet "n" } "." }
|
||||
{ $notes "You will get the wrong result if you try to multiply a quaternion by a complex number on the right using " { $link v*n } ". Use this word instead."
|
||||
$nl "Note that " { $link v*n } " with a quaternion and a real is okay." } ;
|
||||
|
||||
HELP: c>q
|
||||
{ $values { "c" number } { "q" "a quaternion" } }
|
||||
{ $description "Turn a complex number into a quaternion." }
|
||||
{ $examples { $example "USING: math.quaternions prettyprint ;" "C{ 0 1 } c>q ." "{ C{ 0 1 } 0 }" } } ;
|
||||
|
||||
HELP: v>q
|
||||
{ $values { "v" vector } { "q" "a quaternion" } }
|
||||
{ $description "Turn a 3-vector into a quaternion with real part 0." }
|
||||
{ $examples { $example "USING: math.quaternions prettyprint ;" "{ 1 0 0 } v>q ." "{ C{ 0 1 } 0 }" } } ;
|
||||
|
||||
HELP: q>v
|
||||
{ $values { "q" "a quaternion" } { "v" vector } }
|
||||
{ $description "Get the vector part of a quaternion, discarding the real part." }
|
||||
{ $examples { $example "USING: math.quaternions prettyprint ;" "{ C{ 0 1 } 0 } q>v ." "{ 1 0 0 }" } } ;
|
||||
|
||||
HELP: euler
|
||||
{ $values { "phi" number } { "theta" number } { "psi" number } { "q" "a quaternion" } }
|
||||
{ $description "Convert a rotation given by Euler angles (phi, theta, and psi) to a quaternion." } ;
|
||||
|
||||
|
|
@ -1,15 +1,13 @@
|
|||
! Copyright (C) 2005, 2007 Slava Pestov.
|
||||
! See http://factorcode.org/license.txt for BSD license.
|
||||
|
||||
! Everybody's favorite non-commutative skew field, the
|
||||
! quaternions!
|
||||
|
||||
! Quaternions are represented as pairs of complex numbers,
|
||||
! using the identity: (a+bi)+(c+di)j = a+bi+cj+dk.
|
||||
USING: arrays kernel math math.vectors math.functions
|
||||
arrays sequences ;
|
||||
USING: arrays kernel math math.functions math.vectors sequences ;
|
||||
IN: math.quaternions
|
||||
|
||||
! Everybody's favorite non-commutative skew field, the quaternions!
|
||||
|
||||
! Quaternions are represented as pairs of complex numbers, using the
|
||||
! identity: (a+bi)+(c+di)j = a+bi+cj+dk.
|
||||
|
||||
<PRIVATE
|
||||
|
||||
: ** conjugate * ; inline
|
||||
|
|
@ -23,39 +21,27 @@ IN: math.quaternions
|
|||
PRIVATE>
|
||||
|
||||
: q* ( u v -- u*v )
|
||||
#! Multiply quaternions.
|
||||
[ q*a ] [ q*b ] 2bi 2array ;
|
||||
|
||||
: qconjugate ( u -- u' )
|
||||
#! Quaternion conjugate.
|
||||
first2 [ conjugate ] [ neg ] bi* 2array ;
|
||||
|
||||
: qrecip ( u -- 1/u )
|
||||
#! Quaternion inverse.
|
||||
qconjugate dup norm-sq v/n ;
|
||||
|
||||
: q/ ( u v -- u/v )
|
||||
#! Divide quaternions.
|
||||
qrecip q* ;
|
||||
|
||||
: q*n ( q n -- q )
|
||||
#! Note: you will get the wrong result if you try to
|
||||
#! multiply a quaternion by a complex number on the right
|
||||
#! using v*n. Use this word instead. Note that v*n with a
|
||||
#! quaternion and a real is okay.
|
||||
conjugate v*n ;
|
||||
|
||||
: c>q ( c -- q )
|
||||
#! Turn a complex number into a quaternion.
|
||||
0 2array ;
|
||||
|
||||
: v>q ( v -- q )
|
||||
#! Turn a 3-vector into a quaternion with real part 0.
|
||||
first3 rect> [ 0 swap rect> ] dip 2array ;
|
||||
|
||||
: q>v ( q -- v )
|
||||
#! Get the vector part of a quaternion, discarding the real
|
||||
#! part.
|
||||
first2 [ imaginary-part ] dip >rect 3array ;
|
||||
|
||||
! Zero
|
||||
|
|
@ -67,11 +53,14 @@ PRIVATE>
|
|||
: qj { 0 1 } ;
|
||||
: qk { 0 C{ 0 1 } } ;
|
||||
|
||||
! Euler angles -- see
|
||||
! http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/euleranglestoquaternions.html
|
||||
! Euler angles
|
||||
|
||||
<PRIVATE
|
||||
|
||||
: (euler) ( theta unit -- q )
|
||||
[ -0.5 * dup cos c>q swap sin ] dip n*v v- ;
|
||||
[ -0.5 * [ cos c>q ] [ sin ] bi ] dip n*v v- ;
|
||||
|
||||
PRIVATE>
|
||||
|
||||
: euler ( phi theta psi -- q )
|
||||
[ qi (euler) ] [ qj (euler) ] [ qk (euler) ] tri* q* q* ;
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
! Copyright (C) 2008 Doug Coleman, Michael Judge.
|
||||
! See http://factorcode.org/license.txt for BSD license.
|
||||
USING: arrays kernel math math.analysis math.functions sequences sequences.lib
|
||||
sorting ;
|
||||
USING: arrays combinators kernel math math.analysis math.functions sequences
|
||||
sequences.lib sorting ;
|
||||
IN: math.statistics
|
||||
|
||||
: mean ( seq -- n )
|
||||
|
|
@ -63,7 +63,7 @@ IN: math.statistics
|
|||
r sq ;
|
||||
|
||||
: least-squares ( {{x,y}...} -- alpha beta )
|
||||
[r] >r >r >r >r 2dup r> r> r> r>
|
||||
[r] { [ 2dup ] [ ] [ ] [ ] [ ] } spread
|
||||
! stack is mean(x) mean(y) mean(x) mean(y) {x} {y} sx sy
|
||||
[ (r) ] 2keep ! stack is mean(x) mean(y) r sx sy
|
||||
swap / * ! stack is mean(x) mean(y) beta
|
||||
|
|
|
|||
|
|
@ -32,7 +32,7 @@ IN: project-euler.047
|
|||
<PRIVATE
|
||||
|
||||
: (consecutive) ( count goal test -- n )
|
||||
pick pick = [
|
||||
2over = [
|
||||
swap - nip
|
||||
] [
|
||||
dup prime? [ [ drop 0 ] 2dip ] [
|
||||
|
|
|
|||
|
|
@ -0,0 +1,5 @@
|
|||
USING: project-euler.099 project-euler.099.private tools.test ;
|
||||
IN: project-euler.099.tests
|
||||
|
||||
[ 2 ] [ { { 2 11 } { 3 7 } } solve ] unit-test
|
||||
[ 709 ] [ euler099 ] unit-test
|
||||
|
|
@ -0,0 +1,52 @@
|
|||
! Copyright (c) 2008 Aaron Schaefer.
|
||||
! See http://factorcode.org/license.txt for BSD license.
|
||||
USING: io.encodings.ascii io.files kernel math math.functions math.parser
|
||||
math.vectors sequences splitting ;
|
||||
IN: project-euler.099
|
||||
|
||||
! http://projecteuler.net/index.php?section=problems&id=99
|
||||
|
||||
! DESCRIPTION
|
||||
! -----------
|
||||
|
||||
! Comparing two numbers written in index form like 2^11 and 3^7 is not difficult,
|
||||
! as any calculator would confirm that 2^11 = 2048 < 3^7 = 2187.
|
||||
|
||||
! However, confirming that 632382^518061 519432^525806 would be much more
|
||||
! difficult, as both numbers contain over three million digits.
|
||||
|
||||
! Using base_exp.txt (right click and 'Save Link/Target As...'), a 22K text
|
||||
! file containing one thousand lines with a base/exponent pair on each line,
|
||||
! determine which line number has the greatest numerical value.
|
||||
|
||||
! NOTE: The first two lines in the file represent the numbers in the example
|
||||
! given above.
|
||||
|
||||
|
||||
! SOLUTION
|
||||
! --------
|
||||
|
||||
! Use logarithms to make the calculations necessary more manageable.
|
||||
|
||||
<PRIVATE
|
||||
|
||||
: source-099 ( -- seq )
|
||||
"resource:extra/project-euler/099/base_exp.txt"
|
||||
ascii file-lines [ "," split [ string>number ] map ] map ;
|
||||
|
||||
: simplify ( seq -- seq )
|
||||
#! exponent * log(base)
|
||||
flip first2 swap [ log ] map v* ;
|
||||
|
||||
: solve ( seq -- index )
|
||||
simplify [ supremum ] keep index 1+ ;
|
||||
|
||||
PRIVATE>
|
||||
|
||||
: euler099 ( -- answer )
|
||||
source-099 solve ;
|
||||
|
||||
! [ euler099 ] 100 ave-time
|
||||
! 16 ms ave run timen - 1.67 SD (100 trials)
|
||||
|
||||
MAIN: euler099
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,5 +1,5 @@
|
|||
USING: project-euler.203 tools.test ;
|
||||
USING: project-euler.203 project-euler.203.private tools.test ;
|
||||
IN: project-euler.203.tests
|
||||
|
||||
[ 105 ] [ 8 solve ] unit-test
|
||||
[ 34029210557338 ] [ 51 solve ] unit-test
|
||||
[ 34029210557338 ] [ euler203 ] unit-test
|
||||
|
|
|
|||
|
|
@ -1,9 +1,64 @@
|
|||
! Copyright (c) 2008 Eric Mertens.
|
||||
! See http://factorcode.org/license.txt for BSD license.
|
||||
USING: fry kernel math math.primes.factors sequences sets ;
|
||||
IN: project-euler.203
|
||||
|
||||
: iterate ( n initial quot -- results ) swapd '[ @ dup ] replicate nip ; inline
|
||||
: (generate) ( seq -- seq ) [ 0 prefix ] [ 0 suffix ] bi [ + ] 2map ;
|
||||
: generate ( n -- seq ) 1- { 1 } [ (generate) ] iterate concat prune ;
|
||||
: squarefree ( n -- ? ) factors duplicates empty? ;
|
||||
: solve ( n -- n ) generate [ squarefree ] filter sum ;
|
||||
: euler203 ( -- n ) 51 solve ;
|
||||
! http://projecteuler.net/index.php?section=problems&id=203
|
||||
|
||||
! DESCRIPTION
|
||||
! -----------
|
||||
|
||||
! The binomial coefficients nCk can be arranged in triangular form, Pascal's
|
||||
! triangle, like this:
|
||||
|
||||
! 1
|
||||
! 1 1
|
||||
! 1 2 1
|
||||
! 1 3 3 1
|
||||
! 1 4 6 4 1
|
||||
! 1 5 10 10 5 1
|
||||
! 1 6 15 20 15 6 1
|
||||
! 1 7 21 35 35 21 7 1
|
||||
! .........
|
||||
|
||||
! It can be seen that the first eight rows of Pascal's triangle contain twelve
|
||||
! distinct numbers: 1, 2, 3, 4, 5, 6, 7, 10, 15, 20, 21 and 35.
|
||||
|
||||
! A positive integer n is called squarefree if no square of a prime divides n.
|
||||
! Of the twelve distinct numbers in the first eight rows of Pascal's triangle,
|
||||
! all except 4 and 20 are squarefree. The sum of the distinct squarefree numbers
|
||||
! in the first eight rows is 105.
|
||||
|
||||
! Find the sum of the distinct squarefree numbers in the first 51 rows of
|
||||
! Pascal's triangle.
|
||||
|
||||
|
||||
! SOLUTION
|
||||
! --------
|
||||
|
||||
<PRIVATE
|
||||
|
||||
: iterate ( n initial quot -- results )
|
||||
swapd '[ @ dup ] replicate nip ; inline
|
||||
|
||||
: (generate) ( seq -- seq )
|
||||
[ 0 prefix ] [ 0 suffix ] bi [ + ] 2map ;
|
||||
|
||||
: generate ( n -- seq )
|
||||
1- { 1 } [ (generate) ] iterate concat prune ;
|
||||
|
||||
: squarefree ( n -- ? )
|
||||
factors all-unique? ;
|
||||
|
||||
: solve ( n -- n )
|
||||
generate [ squarefree ] filter sum ;
|
||||
|
||||
PRIVATE>
|
||||
|
||||
: euler203 ( -- n )
|
||||
51 solve ;
|
||||
|
||||
! [ euler203 ] 100 ave-time
|
||||
! 12 ms ave run time - 1.6 SD (100 trials)
|
||||
|
||||
MAIN: euler203
|
||||
|
|
|
|||
|
|
@ -9,7 +9,7 @@ IN: project-euler.215
|
|||
! -----------
|
||||
|
||||
! Consider the problem of building a wall out of 2x1 and 3x1 bricks
|
||||
! (horizontalvertical dimensions) such that, for extra strength, the gaps
|
||||
! (horizontal x vertical dimensions) such that, for extra strength, the gaps
|
||||
! between horizontally-adjacent bricks never line up in consecutive layers,
|
||||
! i.e. never form a "running crack".
|
||||
|
||||
|
|
|
|||
|
|
@ -17,10 +17,11 @@ USING: definitions io io.files kernel math math.parser
|
|||
project-euler.052 project-euler.053 project-euler.055 project-euler.056
|
||||
project-euler.059 project-euler.067 project-euler.071 project-euler.073
|
||||
project-euler.075 project-euler.076 project-euler.079 project-euler.092
|
||||
project-euler.097 project-euler.100 project-euler.116 project-euler.117
|
||||
project-euler.134 project-euler.148 project-euler.150 project-euler.151
|
||||
project-euler.164 project-euler.169 project-euler.173 project-euler.175
|
||||
project-euler.186 project-euler.190 project-euler.215 ;
|
||||
project-euler.097 project-euler.099 project-euler.100 project-euler.116
|
||||
project-euler.117 project-euler.134 project-euler.148 project-euler.150
|
||||
project-euler.151 project-euler.164 project-euler.169 project-euler.173
|
||||
project-euler.175 project-euler.186 project-euler.190 project-euler.203
|
||||
project-euler.215 ;
|
||||
IN: project-euler
|
||||
|
||||
<PRIVATE
|
||||
|
|
|
|||
Loading…
Reference in New Issue