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