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