1-- "Extended" math module for Lunatic. 2 3local ffi = require("ffi") 4 5local bit = require("bit") 6local math = require("math") 7 8local arshift = bit.arshift 9local abs, sqrt = math.abs, math.sqrt 10 11local assert = assert 12local error = error 13local type = type 14 15local OUR_REQUIRE_STRING = [[ 16 local _xm=require'xmath' 17 local _v,_iv=_xm.vec3,_xm.ivec3 18]] 19local function our_get_require() 20 return OUR_REQUIRE_STRING 21end 22 23 24module(...) 25 26 27---=== TRIGONOMETRY ===--- 28 29local BANG2RAD = math.pi/1024 30local isintab = ffi.new("int16_t [?]", 2048) 31local dsintab = ffi.new("double [?]", 2048) 32 33for a=0,511 do 34 local s = math.sin(a*BANG2RAD) 35 isintab[a] = 16384*s 36 dsintab[a] = s 37end 38 39isintab[512] = 16384 40dsintab[512] = 1 41 42for i=513,1023 do 43 isintab[i] = isintab[1024-i]; 44 dsintab[i] = dsintab[1024-i]; 45end 46 47for i=1024,2047 do 48 isintab[i] = -isintab[i-1024]; 49 dsintab[i] = -dsintab[i-1024]; 50end 51 52 53local band = bit.band 54 55local function ksc_common(ang) 56 ang = band(ang, 2047) 57 assert(ang >= 0 and ang < 2048) -- might have been passed NaN 58 return ang 59end 60 61-- k{sin,cos}: 16384-scaled output, 2048-based angle input 62function ksin(ang) 63 return isintab[ksc_common(ang)] 64end 65 66function kcos(ang) 67 return isintab[ksc_common(ang+512)] 68end 69 70 71local sin, cos = math.sin, math.cos 72 73-- {sin,cos}b: [-1..1] output, 2048-based angle input 74function sinb(ang) 75 return dsintab[ksc_common(ang)] 76end 77 78function cosb(ang) 79 return dsintab[ksc_common(ang+512)] 80end 81 82local cosb, sinb = cosb, sinb 83 84 85---=== Approximations to 2D and 3D Euclidean distances ===--- 86-- (also see common.c) 87 88local function dist_common(pos1, pos2) 89 local x = abs(pos1.x - pos2.x) 90 local y = abs(pos1.y - pos2.y) 91 if (x < y) then 92 x, y = y, x 93 end 94 return x, y 95end 96 97function ldist(pos1, pos2) 98 local x, y = dist_common(pos1, pos2) 99 100 local t = y + arshift(y,1) 101 return x - arshift(x,5) - arshift(x,7) + arshift(t,2) + arshift(t,6) 102end 103 104function dist(pos1, pos2) 105 local x, y = dist_common(pos1, pos2) 106 local z = abs(arshift(pos1.z - pos2.z, 4)) 107 108 if (x < z) then 109 x, z = z, x 110 end 111 112 local t = y + z 113 return x - arshift(x,4) + arshift(t,2) + arshift(t,3) 114end 115 116 117---=== VECTOR TYPES ===--- 118 119 120-- The integer 3-vector can be useful for calculations expecting integer 121-- values, e.g. ivec3(x, y, z) is a reasonable way to round a vec3. It can also 122-- be used as the RHS to the vec2/vec3 arithmetic methods. 123-- NOTE: We must have a typedef with that exact name, because for Lunatic 124-- (i.e. not stand-alone), the type was already declared in defs_common.lua. 125ffi.cdef "typedef struct { int32_t x, y, z; } vec3_t;" 126local ivec3_t = ffi.typeof("vec3_t") 127 128 129local dvec2_t = ffi.typeof("struct { double x, y; }") 130local dvec3_t = ffi.typeof("struct { double x, y, z; }") 131 132local vec2_mt = { 133 __add = function(a, b) return dvec2_t(a.x+b.x, a.y+b.y) end, 134 __sub = function(a, b) return dvec2_t(a.x-b.x, a.y-b.y) end, 135 __unm = function(a) return dvec2_t(-a.x, -a.y) end, 136 137 __mul = function(a,b) 138 if (type(a)=="number") then 139 return dvec2_t(a*b.x, a*b.y) 140 end 141 142 if (type(b)~="number") then 143 error("number expected in vec2 multiplication", 2) 144 end 145 return dvec2_t(a.x*b, a.y*b) 146 end, 147 148 __div = function(a,b) 149 if (type(b)~="number") then 150 error("number expected in vec2 division", 2) 151 end 152 return dvec2_t(a.x/b, a.y/b) 153 end, 154 155 __tostring = function(a) return "vec2("..a.x..", "..a.y..")" end, 156 157 __index = { 158 lensq = function(a) return a.x*a.x + a.y*a.y end, 159 160 mhlen = function(a) return abs(a.x)+abs(a.y) end, 161 }, 162} 163 164local l_rotate -- fwd-decl (XXX: could be the other way around) 165 166-- The vec3 metatable is shared between the integer- and double-based 3-vector 167-- types. However, some operations are slightly different. 168local vec3_mt = { 169 -- Arithmetic operations. Note that they always return a dvec3. 170 __add = function(a, b) return dvec3_t(a.x+b.x, a.y+b.y, a.z+b.z) end, 171 __sub = function(a, b) return dvec3_t(a.x-b.x, a.y-b.y, a.z-b.z) end, 172 __unm = function(a) return dvec3_t(-a.x, -a.y, -a.z) end, 173 174 __mul = function(a,b) 175 if (type(a)=="number") then 176 return dvec3_t(a*b.x, a*b.y, a*b.z) 177 end 178 179 if (type(b)~="number") then 180 error("number expected in vec3 multiplication", 2) 181 end 182 return dvec3_t(a.x*b, a.y*b, a.z*b) 183 end, 184 185 __div = function(a,b) 186 if (type(b)~="number") then 187 error("number expected in vec3 division", 2) 188 end 189 return dvec3_t(a.x/b, a.y/b, a.z/b) 190 end, 191 192 -- '^' is the "translate upwards" operator, returns same-typed vector. 193 __pow = function(v, zofs) 194 return v:_ctor(v.x, v.y, v.z-zofs) 195 end, 196 197 -- Convenience for human-readable display. 198 __tostring = function(a) 199 return (a:_isi() and "i" or "").."vec3("..a.x..", "..a.y..", "..a.z..")" 200 end, 201 202 __index = { 203 -- Euclidean 3D length. 204 len = function(a) return sqrt(a.x*a.x + a.y*a.y + a.z*a.z) end, 205 -- Euclidean 3D squared length. 206 lensq = function(a) return a.x*a.x + a.y*a.y + a.z*a.z end, 207 208 -- Euclidean 2D length. 209 len2 = function(a) return sqrt(a.x*a.x + a.y*a.y) end, 210 -- Euclidean 2D squared length. 211 len2sq = function(a) return a.x*a.x + a.y*a.y end, 212 213 -- Manhattan-distance 3D length: 214 mhlen = function(a) return abs(a.x)+abs(a.y)+abs(a.z) end, 215 216 toivec3 = function(v) return ivec3_t(v.x, v.y, v.z) end, 217 218 -- BUILD-coordinate (z scaled by 16) <-> uniform conversions. 219 touniform = function(v) 220 return v:_isi() 221 and v:_ctor(v.x, v.y, arshift(v.z, 4)) 222 or v:_ctor(v.x, v.y, v.z/16) 223 end, 224 225 tobuild = function(v) return v:_ctor(v.x, v.y, 16*v.z) end, 226 227 rotate = function(v, ang, pivot) return l_rotate(v, ang, pivot) end, 228 229 -- PRIVATE methods -- 230 231 -- Get the type constructor for this vector. 232 _ctor = function(v, ...) 233 return v:_isi() and ivec3_t(...) or dvec3_t(...) 234 end, 235 -- Is <v> integer vec3? INTERNAL. 236 _isi = function(v) 237 return ffi.istype(ivec3_t, v) 238 end, 239 240 --- Serialization --- 241 _get_require = our_get_require, 242 243 _serialize = function(v) 244 return (v:_isi() and "_iv" or "_v").."("..v.x..","..v.y..","..v.z..")" 245 end, 246 }, 247} 248 249ffi.metatype(dvec2_t, vec2_mt) 250ffi.metatype(dvec3_t, vec3_mt) 251ffi.metatype(ivec3_t, vec3_mt) 252 253-- VEC2 user data constructor. 254-- * vec2([x [, y]]), assuming that x and y are numbers. Vacant positions are 255-- assumed to be 0. 256-- * vec2(<compound>), <compound> can be anything indexable with "x" and "y" 257function vec2(...) 258 local x, y = ... 259 if (type(x)=="number" or x==nil) then 260 return dvec2_t(...) 261 else 262 return dvec2_t(x.x, x.y) 263 end 264end 265 266-- VEC3 user data constructor. 267-- Analogous to VEC2. 268function vec3(...) 269 local x, y, z = ... 270 if (type(x)=="number" or x==nil) then 271 return dvec3_t(...) 272 else 273 return dvec3_t(x.x, x.y, x.z) 274 end 275end 276 277-- IVEC3 user data constructor. 278function ivec3(...) 279 local x, y, z = ... 280 if (type(x)=="number" or x==nil) then 281 return ivec3_t(...) 282 else 283 return ivec3_t(x.x, x.y, x.z) 284 end 285end 286 287local vec2, vec3 = vec2, vec3 288 289 290---=== MISCELLANEOUS MATH ===--- 291 292local intarg = ffi.new("int32_t [1]") 293function bangvec(bang) 294 intarg[0] = bang -- round towards zero 295 return dvec3_t(cosb(intarg[0]), sinb(intarg[0])) 296end 297 298function kangvec(bang, z) 299 intarg[0] = bang -- round towards zero 300 return ivec3_t(kcos(intarg[0]), ksin(intarg[0]), z or 0) 301end 302 303function angvec(ang) 304 return dvec3_t(cos(ang), sin(ang)) 305end 306 307 308local zerovec = vec3() 309-- Point rotation. Note the different order of arguments from engine function. 310-- XXX: passing mixed vec2/vec3 is problematic. Get rid of vec2? 311-- <ang>: BUILD angle (0-2047 based) 312function rotate(pos, ang, pivot) 313 pivot = pivot or zerovec 314 local p = vec3(pos)-pivot 315 local c, s = cosb(ang), sinb(ang) 316 local x, y = p.x, p.y 317 p.x = pivot.x + (c*x - s*y) 318 p.y = pivot.y + (c*y + s*x) 319 return p 320end 321 322l_rotate = rotate 323 324 325-- Two-element vector cross product. 326-- Anti-commutative, distributive. 327local function cross2(v, w) 328 return v.y*w.x - v.x*w.y 329end 330 331 332-- Finds the intersection point of two lines given by 333-- point a and vector v 334-- and 335-- point b and vector w 336-- 337-- Returns: 338-- if <TODO>, nil 339-- if retpoint_p evaluates to a non-true value, coefficients cv and cw such that <TODO> 340-- else, the intersection point 341function intersect(a,v, b,w, retpoint_p) 342 local vxw = cross2(v,w) 343 344 if (vxw ~= 0) then 345 local btoa = vec2(a) - vec2(b) 346 local cv, cw = cross2(w, btoa)/vxw, cross2(v, btoa)/vxw 347 348 if (retpoint_p) then 349 return vec2(a) + cv*vec2(v) 350 else 351 return cv, cw 352 end 353 end 354 355 -- return nil if v and w parallel (or either of them is a point), or if 356 -- they contain NaNs 357end 358