basis/math/functions/functions.factor

   1 ! Copyright (C) 2004, 2007 Slava Pestov.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: math kernel math.constants math.private
   4 math.libm combinators math.order ;
   5 IN: math.functions
   6
   7 <PRIVATE
   8
   9 : (rect>) ( x y -- z )
  10     dup zero? [ drop ] [ <complex> ] if ; inline
  11
  12 PRIVATE>
  13
  14 : rect> ( x y -- z )
  15     over real? over real? and [
  16         (rect>)
  17     ] [
  18         "Complex number must have real components" throw
  19     ] if ; inline
  20
  21 GENERIC: sqrt ( x -- y ) foldable
  22
  23 M: real sqrt
  24     >float dup 0.0 < [ neg fsqrt 0.0 swap rect> ] [ fsqrt ] if ;
  25
  26 : each-bit ( n quot: ( ? -- ) -- )
  27     over 0 number= pick -1 number= or [
  28         2drop
  29     ] [
  30         2dup >r >r >r odd? r> call r> 2/ r> each-bit
  31     ] if ; inline recursive
  32
  33 GENERIC: (^) ( x y -- z ) foldable
  34
  35 : ^n ( z w -- z^w )
  36     1 swap [
  37         [ dupd * ] when >r sq r>
  38     ] each-bit nip ; inline
  39
  40 M: integer (^)
  41     dup 0 < [ neg ^n recip ] [ ^n ] if ;
  42
  43 : ^ ( x y -- z )
  44     over zero? [
  45         dup zero?
  46         [ 2drop 0.0 0.0 / ] [ 0 < [ drop 1.0 0.0 / ] when ] if
  47     ] [
  48         (^)
  49     ] if ; inline
  50
  51 : (^mod) ( n x y -- z )
  52     1 swap [
  53         [ dupd * pick mod ] when >r sq over mod r>
  54     ] each-bit 2nip ; inline
  55
  56 : (gcd) ( b a x y -- a d )
  57     over zero? [
  58         2nip
  59     ] [
  60         swap [ /mod >r over * swapd - r> ] keep (gcd)
  61     ] if ;
  62
  63 : gcd ( x y -- a d )
  64     0 -rot 1 -rot (gcd) dup 0 < [ neg ] when ; foldable
  65
  66 : lcm ( a b -- c )
  67     [ * ] 2keep gcd nip /i ; foldable
  68
  69 : mod-inv ( x n -- y )
  70     tuck gcd 1 = [
  71         dup 0 < [ + ] [ nip ] if
  72     ] [
  73         "Non-trivial divisor found" throw
  74     ] if ; foldable
  75
  76 : ^mod ( x y n -- z )
  77     over 0 < [
  78         [ >r neg r> ^mod ] keep mod-inv
  79     ] [
  80         -rot (^mod)
  81     ] if ; foldable
  82
  83 GENERIC: absq ( x -- y ) foldable
  84
  85 M: real absq sq ;
  86
  87 : ~abs ( x y epsilon -- ? )
  88     >r - abs r> < ;
  89
  90 : ~rel ( x y epsilon -- ? )
  91     >r [ - abs ] 2keep [ abs ] bi@ + r> * < ;
  92
  93 : ~ ( x y epsilon -- ? )
  94     {
  95         { [ pick fp-nan? pick fp-nan? or ] [ 3drop f ] }
  96         { [ dup zero? ] [ drop number= ] }
  97         { [ dup 0 < ] [ ~rel ] }
  98         [ ~abs ]
  99     } cond ;
 100
 101 : >rect ( z -- x y ) dup real-part swap imaginary-part ; inline
 102
 103 : conjugate ( z -- z* ) >rect neg rect> ; inline
 104
 105 : >float-rect ( z -- x y )
 106     >rect swap >float swap >float ; inline
 107
 108 : arg ( z -- arg ) >float-rect swap fatan2 ; inline
 109
 110 : >polar ( z -- abs arg )
 111     >float-rect [ [ sq ] bi@ + fsqrt ] 2keep swap fatan2 ;
 112     inline
 113
 114 : cis ( arg -- z ) dup fcos swap fsin rect> ; inline
 115
 116 : polar> ( abs arg -- z ) cis * ; inline
 117
 118 : ^mag ( w abs arg -- magnitude )
 119     >r >r >float-rect swap r> swap fpow r> rot * fexp /f ;
 120     inline
 121
 122 : ^theta ( w abs arg -- theta )
 123     >r >r >float-rect r> flog * swap r> * + ; inline
 124
 125 M: number (^)
 126     swap >polar 3dup ^theta >r ^mag r> polar> ;
 127
 128 : [-1,1]? ( x -- ? )
 129     dup complex? [ drop f ] [ abs 1 <= ] if ; inline
 130
 131 : >=1? ( x -- ? )
 132     dup complex? [ drop f ] [ 1 >= ] if ; inline
 133
 134 : exp ( x -- y ) >rect swap fexp swap polar> ; inline
 135
 136 : log ( x -- y ) >polar swap flog swap rect> ; inline
 137
 138 : cos ( x -- y )
 139     dup complex? [
 140         >float-rect 2dup
 141         fcosh swap fcos * -rot
 142         fsinh swap fsin neg * rect>
 143     ] [ fcos ] if ; foldable
 144
 145 : sec ( x -- y ) cos recip ; inline
 146
 147 : cosh ( x -- y )
 148     dup complex? [
 149         >float-rect 2dup
 150         fcos swap fcosh * -rot
 151         fsin swap fsinh * rect>
 152     ] [ fcosh ] if ; foldable
 153
 154 : sech ( x -- y ) cosh recip ; inline
 155
 156 : sin ( x -- y )
 157     dup complex? [
 158         >float-rect 2dup
 159         fcosh swap fsin * -rot
 160         fsinh swap fcos * rect>
 161     ] [ fsin ] if ; foldable
 162
 163 : cosec ( x -- y ) sin recip ; inline
 164
 165 : sinh ( x -- y )
 166     dup complex? [
 167         >float-rect 2dup
 168         fcos swap fsinh * -rot
 169         fsin swap fcosh * rect>
 170     ] [ fsinh ] if ; foldable
 171
 172 : cosech ( x -- y ) sinh recip ; inline
 173
 174 : tan ( x -- y )
 175     dup complex? [ dup sin swap cos / ] [ ftan ] if ; inline
 176
 177 : tanh ( x -- y )
 178     dup complex? [ dup sinh swap cosh / ] [ ftanh ] if ; inline
 179
 180 : cot ( x -- y ) tan recip ; inline
 181
 182 : coth ( x -- y ) tanh recip ; inline
 183
 184 : acosh ( x -- y )
 185     dup sq 1- sqrt + log ; inline
 186
 187 : asech ( x -- y ) recip acosh ; inline
 188
 189 : asinh ( x -- y )
 190     dup sq 1+ sqrt + log ; inline
 191
 192 : acosech ( x -- y ) recip asinh ; inline
 193
 194 : atanh ( x -- y )
 195     dup 1+ swap 1- neg / log 2 / ; inline
 196
 197 : acoth ( x -- y ) recip atanh ; inline
 198
 199 : i* ( x -- y ) >rect neg swap rect> ;
 200
 201 : -i* ( x -- y ) >rect swap neg rect> ;
 202
 203 : asin ( x -- y )
 204     dup [-1,1]? [ fasin ] [ i* asinh -i* ] if ; inline
 205
 206 : acos ( x -- y )
 207     dup [-1,1]? [ facos ] [ asin pi 2 / swap - ] if ;
 208     inline
 209
 210 : atan ( x -- y )
 211     dup complex? [ i* atanh i* ] [ fatan ] if ; inline
 212
 213 : asec ( x -- y ) recip acos ; inline
 214
 215 : acosec ( x -- y ) recip asin ; inline
 216
 217 : acot ( x -- y ) recip atan ; inline
 218
 219 : truncate ( x -- y ) dup 1 mod - ; inline
 220
 221 : round ( x -- y ) dup sgn 2 / + truncate ; inline
 222
 223 : floor ( x -- y )
 224     dup 1 mod dup zero?
 225     [ drop ] [ dup 0 < [ - 1- ] [ - ] if ] if ; foldable
 226
 227 : ceiling ( x -- y ) neg floor neg ; foldable