Move factor-2s from miller-rabin to math.functions; use it to make ^ more efficient
parent
f1a1a4c1cb
commit
4c0f2cc3f5
|
|
@ -134,3 +134,6 @@ IN: math.functions.tests
|
||||||
[ -4.0 ] [ -4.4 round ] unit-test
|
[ -4.0 ] [ -4.4 round ] unit-test
|
||||||
[ 5.0 ] [ 4.5 round ] unit-test
|
[ 5.0 ] [ 4.5 round ] unit-test
|
||||||
[ 4.0 ] [ 4.4 round ] unit-test
|
[ 4.0 ] [ 4.4 round ] unit-test
|
||||||
|
|
||||||
|
[ 6 59967 ] [ 3837888 factor-2s ] unit-test
|
||||||
|
[ 6 -59967 ] [ -3837888 factor-2s ] unit-test
|
||||||
|
|
|
||||||
|
|
@ -1,7 +1,7 @@
|
||||||
! Copyright (C) 2004, 2007 Slava Pestov.
|
! Copyright (C) 2004, 2008 Slava Pestov.
|
||||||
! See http://factorcode.org/license.txt for BSD license.
|
! See http://factorcode.org/license.txt for BSD license.
|
||||||
USING: math kernel math.constants math.private
|
USING: math kernel math.constants math.private
|
||||||
math.libm combinators math.order ;
|
math.libm combinators math.order math.ratios sequences ;
|
||||||
IN: math.functions
|
IN: math.functions
|
||||||
|
|
||||||
<PRIVATE
|
<PRIVATE
|
||||||
|
|
@ -30,14 +30,35 @@ M: real sqrt
|
||||||
2dup >r >r >r odd? r> call r> 2/ r> each-bit
|
2dup >r >r >r odd? r> call r> 2/ r> each-bit
|
||||||
] if ; inline recursive
|
] if ; inline recursive
|
||||||
|
|
||||||
: ^n ( z w -- z^w )
|
: map-bits ( n quot: ( ? -- obj ) -- seq )
|
||||||
1 swap [
|
accumulator [ each-bit ] dip ; inline
|
||||||
[ dupd * ] when >r sq r>
|
|
||||||
] each-bit nip ; inline
|
: factor-2s ( n -- r s )
|
||||||
|
#! factor an integer into 2^r * s
|
||||||
|
dup 0 = [ 1 ] [
|
||||||
|
0 swap [ dup even? ] [ [ 1+ ] [ 2/ ] bi* ] [ ] while
|
||||||
|
] if ; inline
|
||||||
|
|
||||||
|
<PRIVATE
|
||||||
|
|
||||||
|
GENERIC# ^n 1 ( z w -- z^w )
|
||||||
|
|
||||||
|
: (^n) 1 swap [ [ dupd * ] when [ sq ] dip ] each-bit nip ; inline
|
||||||
|
|
||||||
|
M: integer ^n
|
||||||
|
[ factor-2s ] dip [ (^n) ] keep rot * shift ;
|
||||||
|
|
||||||
|
M: ratio ^n
|
||||||
|
[ >fraction ] dip tuck [ ^n ] 2bi@ / ;
|
||||||
|
|
||||||
|
M: float ^n
|
||||||
|
(^n) ;
|
||||||
|
|
||||||
: integer^ ( x y -- z )
|
: integer^ ( x y -- z )
|
||||||
dup 0 > [ ^n ] [ neg ^n recip ] if ; inline
|
dup 0 > [ ^n ] [ neg ^n recip ] if ; inline
|
||||||
|
|
||||||
|
PRIVATE>
|
||||||
|
|
||||||
: >rect ( z -- x y )
|
: >rect ( z -- x y )
|
||||||
[ real-part ] [ imaginary-part ] bi ; inline
|
[ real-part ] [ imaginary-part ] bi ; inline
|
||||||
|
|
||||||
|
|
@ -52,6 +73,8 @@ M: real sqrt
|
||||||
|
|
||||||
: polar> ( abs arg -- z ) cis * ; inline
|
: polar> ( abs arg -- z ) cis * ; inline
|
||||||
|
|
||||||
|
<PRIVATE
|
||||||
|
|
||||||
: ^mag ( w abs arg -- magnitude )
|
: ^mag ( w abs arg -- magnitude )
|
||||||
>r >r >float-rect swap r> swap fpow r> rot * fexp /f ;
|
>r >r >float-rect swap r> swap fpow r> rot * fexp /f ;
|
||||||
inline
|
inline
|
||||||
|
|
@ -68,6 +91,8 @@ M: real sqrt
|
||||||
: 0^ ( x -- z )
|
: 0^ ( x -- z )
|
||||||
dup zero? [ drop 0./0. ] [ 0 < 1./0. 0 ? ] if ; inline
|
dup zero? [ drop 0./0. ] [ 0 < 1./0. 0 ? ] if ; inline
|
||||||
|
|
||||||
|
PRIVATE>
|
||||||
|
|
||||||
: ^ ( x y -- z )
|
: ^ ( x y -- z )
|
||||||
{
|
{
|
||||||
{ [ over zero? ] [ nip 0^ ] }
|
{ [ over zero? ] [ nip 0^ ] }
|
||||||
|
|
|
||||||
|
|
@ -11,13 +11,6 @@ IN: math.miller-rabin
|
||||||
|
|
||||||
TUPLE: positive-even-expected n ;
|
TUPLE: positive-even-expected n ;
|
||||||
|
|
||||||
: (factor-2s) ( r s -- r s )
|
|
||||||
dup even? [ -1 shift [ 1+ ] dip (factor-2s) ] when ;
|
|
||||||
|
|
||||||
: factor-2s ( n -- r s )
|
|
||||||
#! factor an integer into s * 2^r
|
|
||||||
0 swap (factor-2s) ;
|
|
||||||
|
|
||||||
:: (miller-rabin) ( n trials -- ? )
|
:: (miller-rabin) ( n trials -- ? )
|
||||||
[let | r [ n 1- factor-2s drop ]
|
[let | r [ n 1- factor-2s drop ]
|
||||||
s [ n 1- factor-2s nip ]
|
s [ n 1- factor-2s nip ]
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue