add an influence slider to mesh cache.
[blender.git] / source / blender / blenlib / intern / math_vector.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: some of this file.
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  * */
25
26 /** \file blender/blenlib/intern/math_vector.c
27  *  \ingroup bli
28  */
29
30
31
32 #include "BLI_math.h"
33
34 //******************************* Interpolation *******************************/
35
36 void interp_v2_v2v2(float target[2], const float a[2], const float b[2], const float t)
37 {
38         float s = 1.0f - t;
39
40         target[0] = s * a[0] + t * b[0];
41         target[1] = s * a[1] + t * b[1];
42 }
43
44 /* weight 3 2D vectors,
45  * 'w' must be unit length but is not a vector, just 3 weights */
46 void interp_v2_v2v2v2(float p[2], const float v1[2], const float v2[2], const float v3[2], const float w[3])
47 {
48         p[0] = v1[0] * w[0] + v2[0] * w[1] + v3[0] * w[2];
49         p[1] = v1[1] * w[0] + v2[1] * w[1] + v3[1] * w[2];
50 }
51
52 void interp_v3_v3v3(float target[3], const float a[3], const float b[3], const float t)
53 {
54         float s = 1.0f - t;
55
56         target[0] = s * a[0] + t * b[0];
57         target[1] = s * a[1] + t * b[1];
58         target[2] = s * a[2] + t * b[2];
59 }
60
61 void interp_v4_v4v4(float target[4], const float a[4], const float b[4], const float t)
62 {
63         float s = 1.0f - t;
64
65         target[0] = s * a[0] + t * b[0];
66         target[1] = s * a[1] + t * b[1];
67         target[2] = s * a[2] + t * b[2];
68         target[3] = s * a[3] + t * b[3];
69 }
70
71 /* weight 3 vectors,
72  * 'w' must be unit length but is not a vector, just 3 weights */
73 void interp_v3_v3v3v3(float p[3], const float v1[3], const float v2[3], const float v3[3], const float w[3])
74 {
75         p[0] = v1[0] * w[0] + v2[0] * w[1] + v3[0] * w[2];
76         p[1] = v1[1] * w[0] + v2[1] * w[1] + v3[1] * w[2];
77         p[2] = v1[2] * w[0] + v2[2] * w[1] + v3[2] * w[2];
78 }
79
80 /* weight 3 vectors,
81  * 'w' must be unit length but is not a vector, just 4 weights */
82 void interp_v3_v3v3v3v3(float p[3], const float v1[3], const float v2[3], const float v3[3], const float v4[3], const float w[4])
83 {
84         p[0] = v1[0] * w[0] + v2[0] * w[1] + v3[0] * w[2] + v4[0] * w[3];
85         p[1] = v1[1] * w[0] + v2[1] * w[1] + v3[1] * w[2] + v4[1] * w[3];
86         p[2] = v1[2] * w[0] + v2[2] * w[1] + v3[2] * w[2] + v4[2] * w[3];
87 }
88
89 void interp_v4_v4v4v4(float p[4], const float v1[4], const float v2[4], const float v3[4], const float w[3])
90 {
91         p[0] = v1[0] * w[0] + v2[0] * w[1] + v3[0] * w[2];
92         p[1] = v1[1] * w[0] + v2[1] * w[1] + v3[1] * w[2];
93         p[2] = v1[2] * w[0] + v2[2] * w[1] + v3[2] * w[2];
94         p[3] = v1[3] * w[0] + v2[3] * w[1] + v3[3] * w[2];
95 }
96
97 void interp_v4_v4v4v4v4(float p[4], const float v1[4], const float v2[4], const float v3[4], const float v4[4], const float w[4])
98 {
99         p[0] = v1[0] * w[0] + v2[0] * w[1] + v3[0] * w[2] + v4[0] * w[3];
100         p[1] = v1[1] * w[0] + v2[1] * w[1] + v3[1] * w[2] + v4[1] * w[3];
101         p[2] = v1[2] * w[0] + v2[2] * w[1] + v3[2] * w[2] + v4[2] * w[3];
102         p[3] = v1[3] * w[0] + v2[3] * w[1] + v3[3] * w[2] + v4[3] * w[3];
103 }
104
105 void mid_v3_v3v3(float v[3], const float v1[3], const float v2[3])
106 {
107         v[0] = 0.5f * (v1[0] + v2[0]);
108         v[1] = 0.5f * (v1[1] + v2[1]);
109         v[2] = 0.5f * (v1[2] + v2[2]);
110 }
111
112 void mid_v2_v2v2(float v[2], const float v1[2], const float v2[2])
113 {
114         v[0] = 0.5f * (v1[0] + v2[0]);
115         v[1] = 0.5f * (v1[1] + v2[1]);
116 }
117
118 void mid_v3_v3v3v3(float v[3], const float v1[3], const float v2[3], const float v3[3])
119 {
120         v[0] = (v1[0] + v2[0] + v3[0]) / 3.0f;
121         v[1] = (v1[1] + v2[1] + v3[1]) / 3.0f;
122         v[2] = (v1[2] + v2[2] + v3[2]) / 3.0f;
123 }
124
125 /**
126  * Equivalent to:
127  * interp_v3_v3v3(v, v1, v2, -1.0f);
128  */
129
130 void flip_v4_v4v4(float v[4], const float v1[4], const float v2[4])
131 {
132         v[0] = v1[0] + (v1[0] - v2[0]);
133         v[1] = v1[1] + (v1[1] - v2[1]);
134         v[2] = v1[2] + (v1[2] - v2[2]);
135         v[3] = v1[3] + (v1[3] - v2[3]);
136 }
137
138 void flip_v3_v3v3(float v[3], const float v1[3], const float v2[3])
139 {
140         v[0] = v1[0] + (v1[0] - v2[0]);
141         v[1] = v1[1] + (v1[1] - v2[1]);
142         v[2] = v1[2] + (v1[2] - v2[2]);
143 }
144
145 void flip_v2_v2v2(float v[2], const float v1[2], const float v2[2])
146 {
147         v[0] = v1[0] + (v1[0] - v2[0]);
148         v[1] = v1[1] + (v1[1] - v2[1]);
149 }
150
151 /********************************** Angles ***********************************/
152
153 /* Return the angle in radians between vecs 1-2 and 2-3 in radians
154  * If v1 is a shoulder, v2 is the elbow and v3 is the hand,
155  * this would return the angle at the elbow.
156  *
157  * note that when v1/v2/v3 represent 3 points along a straight line
158  * that the angle returned will be pi (180deg), rather then 0.0
159  */
160 float angle_v3v3v3(const float v1[3], const float v2[3], const float v3[3])
161 {
162         float vec1[3], vec2[3];
163
164         sub_v3_v3v3(vec1, v2, v1);
165         sub_v3_v3v3(vec2, v2, v3);
166         normalize_v3(vec1);
167         normalize_v3(vec2);
168
169         return angle_normalized_v3v3(vec1, vec2);
170 }
171
172 /* Quicker than full angle computation */
173 float cos_v3v3v3(const float p1[3], const float p2[3], const float p3[3])
174 {
175         float vec1[3], vec2[3];
176
177         sub_v3_v3v3(vec1, p2, p1);
178         sub_v3_v3v3(vec2, p2, p3);
179         normalize_v3(vec1);
180         normalize_v3(vec2);
181
182         return dot_v3v3(vec1, vec2);
183 }
184
185 /* Return the shortest angle in radians between the 2 vectors */
186 float angle_v3v3(const float v1[3], const float v2[3])
187 {
188         float vec1[3], vec2[3];
189
190         normalize_v3_v3(vec1, v1);
191         normalize_v3_v3(vec2, v2);
192
193         return angle_normalized_v3v3(vec1, vec2);
194 }
195
196 float angle_v2v2v2(const float v1[2], const float v2[2], const float v3[2])
197 {
198         float vec1[2], vec2[2];
199
200         vec1[0] = v2[0] - v1[0];
201         vec1[1] = v2[1] - v1[1];
202
203         vec2[0] = v2[0] - v3[0];
204         vec2[1] = v2[1] - v3[1];
205
206         normalize_v2(vec1);
207         normalize_v2(vec2);
208
209         return angle_normalized_v2v2(vec1, vec2);
210 }
211
212 /* Return the shortest angle in radians between the 2 vectors */
213 float angle_v2v2(const float v1[2], const float v2[2])
214 {
215         float vec1[2], vec2[2];
216
217         vec1[0] = v1[0];
218         vec1[1] = v1[1];
219
220         vec2[0] = v2[0];
221         vec2[1] = v2[1];
222
223         normalize_v2(vec1);
224         normalize_v2(vec2);
225
226         return angle_normalized_v2v2(vec1, vec2);
227 }
228
229 float angle_signed_v2v2(const float v1[2], const float v2[2])
230 {
231         const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]);
232         return atan2f(perp_dot, dot_v2v2(v1, v2));
233 }
234
235 float angle_normalized_v3v3(const float v1[3], const float v2[3])
236 {
237         /* double check they are normalized */
238 #ifdef DEBUG
239         float test;
240         BLI_assert(fabsf((test = len_squared_v3(v1)) - 1.0f) < 0.0001f || fabsf(test) < 0.0001f);
241         BLI_assert(fabsf((test = len_squared_v3(v2)) - 1.0f) < 0.0001f || fabsf(test) < 0.0001f);
242 #endif
243
244         /* this is the same as acos(dot_v3v3(v1, v2)), but more accurate */
245         if (dot_v3v3(v1, v2) < 0.0f) {
246                 float vec[3];
247
248                 vec[0] = -v2[0];
249                 vec[1] = -v2[1];
250                 vec[2] = -v2[2];
251
252                 return (float)M_PI - 2.0f * (float)saasin(len_v3v3(vec, v1) / 2.0f);
253         }
254         else
255                 return 2.0f * (float)saasin(len_v3v3(v2, v1) / 2.0f);
256 }
257
258 float angle_normalized_v2v2(const float v1[2], const float v2[2])
259 {
260         /* double check they are normalized */
261 #ifdef DEBUG
262         float test;
263         BLI_assert(fabsf((test = len_squared_v2(v1)) - 1.0f) < 0.0001f || fabsf(test) < 0.0001f);
264         BLI_assert(fabsf((test = len_squared_v2(v2)) - 1.0f) < 0.0001f || fabsf(test) < 0.0001f);
265 #endif
266
267         /* this is the same as acos(dot_v3v3(v1, v2)), but more accurate */
268         if (dot_v2v2(v1, v2) < 0.0f) {
269                 float vec[2];
270
271                 vec[0] = -v2[0];
272                 vec[1] = -v2[1];
273
274                 return (float)M_PI - 2.0f * saasin(len_v2v2(vec, v1) / 2.0f);
275         }
276         else
277                 return 2.0f * (float)saasin(len_v2v2(v2, v1) / 2.0f);
278 }
279
280 /**
281  * angle between 2 vectors defined by 3 coords, about an axis. */
282 float angle_on_axis_v3v3v3_v3(const float v1[3], const float v2[3], const float v3[3], const float axis[3])
283 {
284         float v1_proj[3], v2_proj[3], tproj[3];
285
286         sub_v3_v3v3(v1_proj, v1, v2);
287         sub_v3_v3v3(v2_proj, v3, v2);
288
289         /* project the vectors onto the axis */
290         project_v3_v3v3(tproj, v1_proj, axis);
291         sub_v3_v3(v1_proj, tproj);
292
293         project_v3_v3v3(tproj, v2_proj, axis);
294         sub_v3_v3(v2_proj, tproj);
295
296         return angle_v3v3(v1_proj, v2_proj);
297 }
298
299 void angle_tri_v3(float angles[3], const float v1[3], const float v2[3], const float v3[3])
300 {
301         float ed1[3], ed2[3], ed3[3];
302
303         sub_v3_v3v3(ed1, v3, v1);
304         sub_v3_v3v3(ed2, v1, v2);
305         sub_v3_v3v3(ed3, v2, v3);
306
307         normalize_v3(ed1);
308         normalize_v3(ed2);
309         normalize_v3(ed3);
310
311         angles[0] = (float)M_PI - angle_normalized_v3v3(ed1, ed2);
312         angles[1] = (float)M_PI - angle_normalized_v3v3(ed2, ed3);
313         // face_angles[2] = M_PI - angle_normalized_v3v3(ed3, ed1);
314         angles[2] = (float)M_PI - (angles[0] + angles[1]);
315 }
316
317 void angle_quad_v3(float angles[4], const float v1[3], const float v2[3], const float v3[3], const float v4[3])
318 {
319         float ed1[3], ed2[3], ed3[3], ed4[3];
320
321         sub_v3_v3v3(ed1, v4, v1);
322         sub_v3_v3v3(ed2, v1, v2);
323         sub_v3_v3v3(ed3, v2, v3);
324         sub_v3_v3v3(ed4, v3, v4);
325
326         normalize_v3(ed1);
327         normalize_v3(ed2);
328         normalize_v3(ed3);
329         normalize_v3(ed4);
330
331         angles[0] = (float)M_PI - angle_normalized_v3v3(ed1, ed2);
332         angles[1] = (float)M_PI - angle_normalized_v3v3(ed2, ed3);
333         angles[2] = (float)M_PI - angle_normalized_v3v3(ed3, ed4);
334         angles[3] = (float)M_PI - angle_normalized_v3v3(ed4, ed1);
335 }
336
337 void angle_poly_v3(float *angles, const float *verts[3], int len)
338 {
339         int i;
340         float vec[3][3];
341
342         sub_v3_v3v3(vec[2], verts[len - 1], verts[0]);
343         normalize_v3(vec[2]);
344         for (i = 0; i < len; i++) {
345                 sub_v3_v3v3(vec[i % 3], verts[i % len], verts[(i + 1) % len]);
346                 normalize_v3(vec[i % 3]);
347                 angles[i] = (float)M_PI - angle_normalized_v3v3(vec[(i + 2) % 3], vec[i % 3]);
348         }
349 }
350
351 /********************************* Geometry **********************************/
352
353 /* Project v1 on v2 */
354 void project_v2_v2v2(float c[2], const float v1[2], const float v2[2])
355 {
356         float mul;
357         mul = dot_v2v2(v1, v2) / dot_v2v2(v2, v2);
358
359         c[0] = mul * v2[0];
360         c[1] = mul * v2[1];
361 }
362
363 /* Project v1 on v2 */
364 void project_v3_v3v3(float c[3], const float v1[3], const float v2[3])
365 {
366         float mul;
367         mul = dot_v3v3(v1, v2) / dot_v3v3(v2, v2);
368
369         c[0] = mul * v2[0];
370         c[1] = mul * v2[1];
371         c[2] = mul * v2[2];
372 }
373
374 /* project a vector on a plane defined by normal and a plane point p */
375 void project_v3_plane(float v[3], const float n[3], const float p[3])
376 {
377         float vector[3];
378         float mul;
379
380         sub_v3_v3v3(vector, v, p);
381         mul = dot_v3v3(vector, n) / len_squared_v3(n);
382
383         mul_v3_v3fl(vector, n, mul);
384
385         sub_v3_v3(v, vector);
386 }
387
388 /* Returns a vector bisecting the angle at v2 formed by v1, v2 and v3 */
389 void bisect_v3_v3v3v3(float out[3], const float v1[3], const float v2[3], const float v3[3])
390 {
391         float d_12[3], d_23[3];
392         sub_v3_v3v3(d_12, v2, v1);
393         sub_v3_v3v3(d_23, v3, v2);
394         normalize_v3(d_12);
395         normalize_v3(d_23);
396         add_v3_v3v3(out, d_12, d_23);
397         normalize_v3(out);
398 }
399
400 /* Returns a reflection vector from a vector and a normal vector
401  * reflect = vec - ((2 * DotVecs(vec, mirror)) * mirror)
402  */
403 void reflect_v3_v3v3(float out[3], const float v1[3], const float v2[3])
404 {
405         float vec[3], normal[3];
406         float reflect[3] = {0.0f, 0.0f, 0.0f};
407         float dot2;
408
409         copy_v3_v3(vec, v1);
410         copy_v3_v3(normal, v2);
411
412         dot2 = 2 * dot_v3v3(vec, normal);
413
414         reflect[0] = vec[0] - (dot2 * normal[0]);
415         reflect[1] = vec[1] - (dot2 * normal[1]);
416         reflect[2] = vec[2] - (dot2 * normal[2]);
417
418         copy_v3_v3(out, reflect);
419 }
420
421 void ortho_basis_v3v3_v3(float v1[3], float v2[3], const float v[3])
422 {
423         const float f = (float)sqrt(v[0] * v[0] + v[1] * v[1]);
424
425         if (f < 1e-35f) {
426                 // degenerate case
427                 v1[0] = (v[2] < 0.0f) ? -1.0f : 1.0f;
428                 v1[1] = v1[2] = v2[0] = v2[2] = 0.0f;
429                 v2[1] = 1.0f;
430         }
431         else {
432                 const float d = 1.0f / f;
433
434                 v1[0] = v[1] * d;
435                 v1[1] = -v[0] * d;
436                 v1[2] = 0.0f;
437                 v2[0] = -v[2] * v1[1];
438                 v2[1] = v[2] * v1[0];
439                 v2[2] = v[0] * v1[1] - v[1] * v1[0];
440         }
441 }
442
443 /* Rotate a point p by angle theta around an arbitrary axis r
444  * http://local.wasp.uwa.edu.au/~pbourke/geometry/
445  */
446 void rotate_normalized_v3_v3v3fl(float r[3], const float p[3], const float axis[3], const float angle)
447 {
448         const float costheta = cos(angle);
449         const float sintheta = sin(angle);
450
451         /* double check they are normalized */
452 #ifdef DEBUG
453         float test;
454         BLI_assert(fabsf((test = len_squared_v3(axis)) - 1.0f) < 0.0001f || fabsf(test) < 0.0001f);
455 #endif
456
457         r[0] = ((costheta + (1 - costheta) * axis[0] * axis[0]) * p[0]) +
458                (((1 - costheta) * axis[0] * axis[1] - axis[2] * sintheta) * p[1]) +
459                (((1 - costheta) * axis[0] * axis[2] + axis[1] * sintheta) * p[2]);
460
461         r[1] = (((1 - costheta) * axis[0] * axis[1] + axis[2] * sintheta) * p[0]) +
462                ((costheta + (1 - costheta) * axis[1] * axis[1]) * p[1]) +
463                (((1 - costheta) * axis[1] * axis[2] - axis[0] * sintheta) * p[2]);
464
465         r[2] = (((1 - costheta) * axis[0] * axis[2] - axis[1] * sintheta) * p[0]) +
466                (((1 - costheta) * axis[1] * axis[2] + axis[0] * sintheta) * p[1]) +
467                ((costheta + (1 - costheta) * axis[2] * axis[2]) * p[2]);
468 }
469
470 void rotate_v3_v3v3fl(float r[3], const float p[3], const float axis[3], const float angle)
471 {
472         float axis_n[3];
473
474         normalize_v3_v3(axis_n, axis);
475
476         rotate_normalized_v3_v3v3fl(r, p, axis_n, angle);
477 }
478
479 /*********************************** Other ***********************************/
480
481 void print_v2(const char *str, const float v[2])
482 {
483         printf("%s: %.3f %.3f\n", str, v[0], v[1]);
484 }
485
486 void print_v3(const char *str, const float v[3])
487 {
488         printf("%s: %.3f %.3f %.3f\n", str, v[0], v[1], v[2]);
489 }
490
491 void print_v4(const char *str, const float v[4])
492 {
493         printf("%s: %.3f %.3f %.3f %.3f\n", str, v[0], v[1], v[2], v[3]);
494 }
495
496 void minmax_v3v3_v3(float min[3], float max[3], const float vec[3])
497 {
498         if (min[0] > vec[0]) min[0] = vec[0];
499         if (min[1] > vec[1]) min[1] = vec[1];
500         if (min[2] > vec[2]) min[2] = vec[2];
501
502         if (max[0] < vec[0]) max[0] = vec[0];
503         if (max[1] < vec[1]) max[1] = vec[1];
504         if (max[2] < vec[2]) max[2] = vec[2];
505 }
506
507 void minmax_v2v2_v2(float min[2], float max[2], const float vec[2])
508 {
509         if (min[0] > vec[0]) min[0] = vec[0];
510         if (min[1] > vec[1]) min[1] = vec[1];
511
512         if (max[0] < vec[0]) max[0] = vec[0];
513         if (max[1] < vec[1]) max[1] = vec[1];
514 }
515
516 /** ensure \a v1 is \a dist from \a v2 */
517 void dist_ensure_v3_v3fl(float v1[3], const float v2[3], const float dist)
518 {
519         if (!equals_v3v3(v2, v1)) {
520                 float nor[3];
521
522                 sub_v3_v3v3(nor, v1, v2);
523                 normalize_v3(nor);
524                 madd_v3_v3v3fl(v1, v2, nor, dist);
525         }
526 }
527
528 void dist_ensure_v2_v2fl(float v1[2], const float v2[2], const float dist)
529 {
530         if (!equals_v2v2(v2, v1)) {
531                 float nor[2];
532
533                 sub_v2_v2v2(nor, v1, v2);
534                 normalize_v2(nor);
535                 madd_v2_v2v2fl(v1, v2, nor, dist);
536         }
537 }
538
539 /***************************** Array Functions *******************************/
540
541 double dot_vn_vn(const float *array_src_a, const float *array_src_b, const int size)
542 {
543         double d = 0.0f;
544         const float *array_pt_a = array_src_a + (size - 1);
545         const float *array_pt_b = array_src_b + (size - 1);
546         int i = size;
547         while (i--) {
548                 d += (double)(*(array_pt_a--) * *(array_pt_b--));
549         }
550         return d;
551 }
552
553 float normalize_vn_vn(float *array_tar, const float *array_src, const int size)
554 {
555         double d = dot_vn_vn(array_tar, array_src, size);
556         float d_sqrt;
557         if (d > 1.0e-35) {
558                 d_sqrt = (float)sqrt(d);
559                 mul_vn_vn_fl(array_tar, array_src, size, 1.0f / d_sqrt);
560         }
561         else {
562                 fill_vn_fl(array_tar, size, 0.0f);
563                 d_sqrt = 0.0f;
564         }
565         return d_sqrt;
566 }
567
568 float normalize_vn(float *array_tar, const int size)
569 {
570         return normalize_vn_vn(array_tar, array_tar, size);
571 }
572
573 void range_vn_i(int *array_tar, const int size, const int start)
574 {
575         int *array_pt = array_tar + (size - 1);
576         int j = start + (size - 1);
577         int i = size;
578         while (i--) {
579                 *(array_pt--) = j--;
580         }
581 }
582
583 void range_vn_fl(float *array_tar, const int size, const float start, const float step)
584 {
585         float *array_pt = array_tar + (size - 1);
586         int i = size;
587         while (i--) {
588                 *(array_pt--) = start + step * (float)(i);
589         }
590 }
591
592 void negate_vn(float *array_tar, const int size)
593 {
594         float *array_pt = array_tar + (size - 1);
595         int i = size;
596         while (i--) {
597                 *(array_pt--) *= -1.0f;
598         }
599 }
600
601 void negate_vn_vn(float *array_tar, const float *array_src, const int size)
602 {
603         float *tar = array_tar + (size - 1);
604         const float *src = array_src + (size - 1);
605         int i = size;
606         while (i--) {
607                 *(tar--) = -*(src--);
608         }
609 }
610
611 void mul_vn_fl(float *array_tar, const int size, const float f)
612 {
613         float *array_pt = array_tar + (size - 1);
614         int i = size;
615         while (i--) {
616                 *(array_pt--) *= f;
617         }
618 }
619
620 void mul_vn_vn_fl(float *array_tar, const float *array_src, const int size, const float f)
621 {
622         float *tar = array_tar + (size - 1);
623         const float *src = array_src + (size - 1);
624         int i = size;
625         while (i--) {
626                 *(tar--) = *(src--) * f;
627         }
628 }
629
630 void add_vn_vn(float *array_tar, const float *array_src, const int size)
631 {
632         float *tar = array_tar + (size - 1);
633         const float *src = array_src + (size - 1);
634         int i = size;
635         while (i--) {
636                 *(tar--) += *(src--);
637         }
638 }
639
640 void add_vn_vnvn(float *array_tar, const float *array_src_a, const float *array_src_b, const int size)
641 {
642         float *tar = array_tar + (size - 1);
643         const float *src_a = array_src_a + (size - 1);
644         const float *src_b = array_src_b + (size - 1);
645         int i = size;
646         while (i--) {
647                 *(tar--) = *(src_a--) + *(src_b--);
648         }
649 }
650
651 void madd_vn_vn(float *array_tar, const float *array_src, const float f, const int size)
652 {
653         float *tar = array_tar + (size - 1);
654         const float *src = array_src + (size - 1);
655         int i = size;
656         while (i--) {
657                 *(tar--) += *(src--) * f;
658         }
659 }
660
661 void madd_vn_vnvn(float *array_tar, const float *array_src_a, const float *array_src_b, const float f, const int size)
662 {
663         float *tar = array_tar + (size - 1);
664         const float *src_a = array_src_a + (size - 1);
665         const float *src_b = array_src_b + (size - 1);
666         int i = size;
667         while (i--) {
668                 *(tar--) = *(src_a--) + (*(src_b--) * f);
669         }
670 }
671
672 void sub_vn_vn(float *array_tar, const float *array_src, const int size)
673 {
674         float *tar = array_tar + (size - 1);
675         const float *src = array_src + (size - 1);
676         int i = size;
677         while (i--) {
678                 *(tar--) -= *(src--);
679         }
680 }
681
682 void sub_vn_vnvn(float *array_tar, const float *array_src_a, const float *array_src_b, const int size)
683 {
684         float *tar = array_tar + (size - 1);
685         const float *src_a = array_src_a + (size - 1);
686         const float *src_b = array_src_b + (size - 1);
687         int i = size;
688         while (i--) {
689                 *(tar--) = *(src_a--) - *(src_b--);
690         }
691 }
692
693 void msub_vn_vn(float *array_tar, const float *array_src, const float f, const int size)
694 {
695         float *tar = array_tar + (size - 1);
696         const float *src = array_src + (size - 1);
697         int i = size;
698         while (i--) {
699                 *(tar--) -= *(src--) * f;
700         }
701 }
702
703 void msub_vn_vnvn(float *array_tar, const float *array_src_a, const float *array_src_b, const float f, const int size)
704 {
705         float *tar = array_tar + (size - 1);
706         const float *src_a = array_src_a + (size - 1);
707         const float *src_b = array_src_b + (size - 1);
708         int i = size;
709         while (i--) {
710                 *(tar--) = *(src_a--) - (*(src_b--) * f);
711         }
712 }
713
714 void interp_vn_vn(float *array_tar, const float *array_src, const float t, const int size)
715 {
716         const float s = 1.0f - t;
717         float *tar = array_tar + (size - 1);
718         const float *src = array_src + (size - 1);
719         int i = size;
720         while (i--) {
721                 *(tar) = (s * *(tar)) + (t * *(src));
722                 tar--;
723                 src--;
724         }
725 }
726
727 void fill_vn_i(int *array_tar, const int size, const int val)
728 {
729         int *tar = array_tar + (size - 1);
730         int i = size;
731         while (i--) {
732                 *(tar--) = val;
733         }
734 }
735
736 void fill_vn_ushort(unsigned short *array_tar, const int size, const unsigned short val)
737 {
738         unsigned short *tar = array_tar + (size - 1);
739         int i = size;
740         while (i--) {
741                 *(tar--) = val;
742         }
743 }
744
745 void fill_vn_fl(float *array_tar, const int size, const float val)
746 {
747         float *tar = array_tar + (size - 1);
748         int i = size;
749         while (i--) {
750                 *(tar--) = val;
751         }
752 }