basis/math/vectors/vectors.factor

   1 ! Copyright (C) 2005, 2010 Slava Pestov, Joe Groff.
   2 ! See https://factorcode.org/license.txt for BSD license.
   3 USING: arrays assocs combinators grouping kernel math
   4 math.functions math.libm math.order sequences ;
   5 QUALIFIED-WITH: alien.c-types c
   6 IN: math.vectors
   7
   8 GENERIC: vneg ( v -- w )
   9 M: object vneg [ neg ] map ; inline
  10
  11 GENERIC#: v+n 1 ( v n -- w )
  12 M: object v+n [ + ] curry map ; inline
  13
  14 GENERIC: n+v ( n v -- w )
  15 M: object n+v [ + ] with map ; inline
  16
  17 GENERIC#: v-n 1 ( v n -- w )
  18 M: object v-n [ - ] curry map ; inline
  19
  20 GENERIC: n-v ( n v -- w )
  21 M: object n-v [ - ] with map ; inline
  22
  23 GENERIC#: v*n 1 ( v n -- w )
  24 M: object v*n [ * ] curry map ; inline
  25
  26 GENERIC: n*v ( n v -- w )
  27 M: object n*v [ * ] with map ; inline
  28
  29 GENERIC#: v/n 1 ( v n -- w )
  30 M: object v/n [ / ] curry map ; inline
  31
  32 GENERIC: n/v ( n v -- w )
  33 M: object n/v [ / ] with map ; inline
  34
  35 GENERIC: v+ ( u v -- w )
  36 M: object v+ [ + ] 2map ; inline
  37
  38 GENERIC: v- ( u v -- w )
  39 M: object v- [ - ] 2map ; inline
  40
  41 GENERIC: [v-] ( u v -- w )
  42 M: object [v-] [ [-] ] 2map ; inline
  43
  44 GENERIC: v* ( u v -- w )
  45 M: object v* [ * ] 2map ; inline
  46
  47 GENERIC: v*high ( u v -- w )
  48
  49 <PRIVATE
  50 : (h+) ( v -- w ) 2 <groups> [ first2 + ] map ;
  51 PRIVATE>
  52
  53 GENERIC: v*hs+ ( u v -- w )
  54 M: object v*hs+ [ * ] 2map (h+) ; inline
  55
  56 GENERIC: v/ ( u v -- w )
  57 M: object v/ [ / ] 2map ; inline
  58
  59 GENERIC: v^ ( u v -- w )
  60 M: object v^ [ ^ ] 2map ; inline
  61
  62 GENERIC: v^n ( v n -- w )
  63 M: object v^n [ ^ ] curry map ; inline
  64
  65 GENERIC: n^v ( n v -- w )
  66 M: object n^v [ ^ ] with map ; inline
  67
  68 GENERIC: vavg ( u v -- w )
  69 M: object vavg [ + 2 / ] 2map ; inline
  70
  71 GENERIC: vmax ( u v -- w )
  72 M: object vmax [ max ] 2map ; inline
  73
  74 GENERIC: vmin ( u v -- w )
  75 M: object vmin [ min ] 2map ; inline
  76
  77 GENERIC: v+- ( u v -- w )
  78 M: object v+-
  79     [ t ] 2dip
  80     [ [ not ] 2dip pick [ + ] [ - ] if ] 2map
  81     nip ; inline
  82
  83 GENERIC: vs+ ( u v -- w )
  84 M: object vs+ [ + ] 2map ; inline
  85
  86 GENERIC: vs- ( u v -- w )
  87 M: object vs- [ - ] 2map ; inline
  88
  89 GENERIC: vs* ( u v -- w )
  90 M: object vs* [ * ] 2map ; inline
  91
  92 GENERIC: vabs ( v -- w )
  93 M: object vabs [ abs ] map ; inline
  94
  95 GENERIC: vsqrt ( v -- w )
  96 M: object vsqrt [ >float fsqrt ] map ; inline
  97
  98 GENERIC: vsad ( u v -- n )
  99 M: object vsad [ - abs ] [ + ] 2map-reduce ; inline
 100
 101 <PRIVATE
 102
 103 : bitandn ( x y -- z ) [ bitnot ] dip bitand ; inline
 104
 105 PRIVATE>
 106
 107 GENERIC: vbitand ( u v -- w )
 108 M: object vbitand [ bitand ] 2map ; inline
 109 GENERIC: vbitandn ( u v -- w )
 110 M: object vbitandn [ bitandn ] 2map ; inline
 111 GENERIC: vbitor ( u v -- w )
 112 M: object vbitor [ bitor ] 2map ; inline
 113 GENERIC: vbitxor ( u v -- w )
 114 M: object vbitxor [ bitxor ] 2map ; inline
 115 GENERIC: vbitnot ( v -- w )
 116 M: object vbitnot [ bitnot ] map ; inline
 117
 118 GENERIC#: vbroadcast 1 ( u n -- v )
 119 M:: object vbroadcast ( u n -- v ) u length n u nth <repetition> u like ; inline
 120
 121 GENERIC#: vshuffle-elements 1 ( v perm -- w )
 122 M: object vshuffle-elements
 123     over length 0 pad-tail
 124     swap [ '[ _ nth ] ] keep map-as ; inline
 125
 126 GENERIC#: vshuffle2-elements 1 ( u v perm -- w )
 127 M: object vshuffle2-elements
 128     [ append ] dip vshuffle-elements ; inline
 129
 130 GENERIC#: vshuffle-bytes 1 ( v perm -- w )
 131
 132 GENERIC: vshuffle ( v perm -- w )
 133 M: array vshuffle
 134     vshuffle-elements ; inline
 135
 136 GENERIC#: vlshift 1 ( v n -- w )
 137 M: object vlshift '[ _ shift ] map ; inline
 138 GENERIC#: vrshift 1 ( v n -- w )
 139 M: object vrshift neg '[ _ shift ] map ; inline
 140
 141 GENERIC#: hlshift 1 ( v n -- w )
 142 GENERIC#: hrshift 1 ( v n -- w )
 143
 144 GENERIC: (vmerge-head) ( u v -- h )
 145 M: object (vmerge-head) over length 2 /i '[ _ head-slice ] bi@ [ zip ] keep concat-as ; inline
 146 GENERIC: (vmerge-tail) ( u v -- t )
 147 M: object (vmerge-tail) over length 2 /i '[ _ tail-slice ] bi@ [ zip ] keep concat-as ; inline
 148
 149 : (vmerge) ( u v -- h t )
 150     [ (vmerge-head) ] [ (vmerge-tail) ] 2bi ; inline
 151
 152 GENERIC: vmerge ( u v -- w )
 153 M: object vmerge [ zip ] keep concat-as ; inline
 154
 155 GENERIC: vand ( u v -- w )
 156 M: object vand [ and ] 2map ; inline
 157
 158 GENERIC: vandn ( u v -- w )
 159 M: object vandn [ [ not ] dip and ] 2map ; inline
 160
 161 GENERIC: vor ( u v -- w )
 162 M: object vor [ or ] 2map ; inline
 163
 164 GENERIC: vxor ( u v -- w )
 165 M: object vxor [ xor ] 2map ; inline
 166
 167 GENERIC: vnot ( v -- w )
 168 M: object vnot [ not ] map ; inline
 169
 170 GENERIC: vall? ( v -- ? )
 171 M: object vall? [ ] all? ; inline
 172
 173 GENERIC: vcount ( v -- count )
 174 M: object vcount [ ] count ; inline
 175
 176 GENERIC: vany? ( v -- ? )
 177 M: object vany? [ ] any? ; inline
 178
 179 GENERIC: vnone? ( v -- ? )
 180 M: object vnone? [ not ] all? ; inline
 181
 182 GENERIC: v< ( u v -- w )
 183 M: object v< [ < ] 2map ; inline
 184
 185 GENERIC: v<= ( u v -- w )
 186 M: object v<= [ <= ] 2map ; inline
 187
 188 GENERIC: v>= ( u v -- w )
 189 M: object v>= [ >= ] 2map ; inline
 190
 191 GENERIC: v> ( u v -- w )
 192 M: object v> [ > ] 2map ; inline
 193
 194 GENERIC: vunordered? ( u v -- w )
 195 M: object vunordered? [ unordered? ] 2map ; inline
 196
 197 GENERIC: v= ( u v -- w )
 198 M: object v= [ = ] 2map ; inline
 199
 200 GENERIC: v? ( mask true false -- result )
 201 M: object v?
 202     [ vand ] [ vandn ] bi-curry* bi vor ; inline
 203
 204 : vif ( mask true-quot: ( -- vector ) false-quot: ( -- vector ) -- result )
 205     {
 206         { [ pick vall? ] [ drop nip call ] }
 207         { [ pick vnone? ] [ 2nip call ] }
 208         [ [ call ] dip call v? ]
 209     } cond ; inline
 210
 211 : v>integer ( v -- w ) [ >integer ] map ;
 212 : vfloor ( v -- w ) [ floor ] map ;
 213 : vceiling ( v -- w ) [ ceiling ] map ;
 214 : vtruncate ( v -- w ) [ truncate ] map ;
 215
 216 : vsupremum ( seq -- vmax ) [ ] [ vmax ] map-reduce ; inline
 217 : vinfimum ( seq -- vmin ) [ ] [ vmin ] map-reduce ; inline
 218
 219 GENERIC: vdot ( u v -- x )
 220 M: object vdot [ * ] [ + ] 2map-reduce ; inline
 221
 222 GENERIC: hdot ( u v -- x )
 223 M: object hdot [ conjugate * ] [ + ] 2map-reduce ; inline
 224
 225 GENERIC: norm-sq ( v -- x )
 226 M: object norm-sq [ absq ] [ + ] map-reduce ; inline
 227
 228 : l1-norm ( k -- x ) [ abs ] map-sum ; inline
 229
 230 : l2-norm ( k -- x ) norm-sq sqrt ; inline
 231
 232 ALIAS: norm l2-norm
 233
 234 : l-infinity-norm ( k -- x ) supremum ; inline
 235
 236 : p-norm-default ( k p -- x )
 237     [ [ [ abs ] dip ^ ] curry map-sum ] keep recip ^ ; inline
 238
 239 : p-norm ( k p -- x )
 240     {
 241         { [ dup 1 = ] [ drop l1-norm ] }
 242         { [ dup 2 = ] [ drop l2-norm ] }
 243         { [ dup fp-infinity? ] [ drop supremum ] }
 244         [ p-norm-default ]
 245     } cond ;
 246
 247 : normalize ( v -- w ) dup l2-norm v/n ; inline
 248
 249 GENERIC: distance ( u v -- x )
 250 M: object distance [ - absq ] [ + ] 2map-reduce sqrt ; inline
 251
 252 : set-axis ( u v axis -- w )
 253     [ [ zero? 2over ? ] dip swap nth ] map-index 2nip ;
 254
 255 <PRIVATE
 256
 257 : 2tetra@ ( p q r s t u v w quot -- )
 258     dup [ [ 2bi@ ] curry 4dip ] dip 2bi@ ; inline
 259
 260 PRIVATE>
 261
 262 : trilerp ( aaa baa aba bba aab bab abb bbb {t,u,v} -- a_tuv )
 263     [ first lerp ] [ second lerp ] [ third lerp ] tri-curry
 264     [ 2tetra@ ] [ 2bi@ ] [ call ] tri* ; inline
 265
 266 : bilerp ( aa ba ab bb {t,u} -- a_tu )
 267     [ first lerp ] [ second lerp ] bi-curry
 268     [ 2bi@ ] [ call ] bi* ; inline
 269
 270 : vlerp ( a b t -- a_t )
 271     [ over v- ] dip v* v+ ; inline
 272
 273 : vnlerp ( a b t -- a_t )
 274     [ over v- ] dip v*n v+ ; inline
 275
 276 : vbilerp ( aa ba ab bb {t,u} -- a_tu )
 277     [ first vnlerp ] [ second vnlerp ] bi-curry
 278     [ 2bi@ ] [ call ] bi* ; inline
 279
 280 : v~ ( a b epsilon -- ? )
 281     [ ~ ] curry 2all? ; inline
 282
 283 : vclamp ( v min max -- w )
 284     rot vmin vmax ; inline
 285
 286 : cross ( vec1 vec2 -- vec3 )
 287     [ [ { 1 2 0 } vshuffle ] [ { 2 0 1 } vshuffle ] bi* v* ]
 288     [ [ { 2 0 1 } vshuffle ] [ { 1 2 0 } vshuffle ] bi* v* ] 2bi v- ; inline
 289
 290 :: normal ( vec1 vec2 vec3 -- vec4 )
 291     vec2 vec1 v- vec3 vec1 v- cross normalize ; inline
 292
 293 : proj ( v u -- w )
 294     [ [ vdot ] [ norm-sq ] bi / ] keep n*v ;
 295
 296 : perp ( v u -- w )
 297     dupd proj v- ;
 298
 299 : angle-between ( v u -- a )
 300     [ normalize ] bi@ hdot acos ;