add BLI_strcpy_rlen, replace strcat, which was used in misleading way.
[blender.git] / intern / cycles / util / util_math.h
1 /*
2  * Copyright 2011, Blender Foundation.
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
19 #ifndef __UTIL_MATH_H__
20 #define __UTIL_MATH_H__
21
22 /* Math
23  *
24  * Basic math functions on scalar and vector types. This header is used by
25  * both the kernel code when compiled as C++, and other C++ non-kernel code. */
26
27 #ifndef __KERNEL_OPENCL__
28
29 #ifdef _MSC_VER
30 #  define _USE_MATH_DEFINES
31 #endif
32
33 #include <float.h>
34 #include <math.h>
35 #include <stdio.h>
36
37 #endif
38
39 #include "util_types.h"
40
41 CCL_NAMESPACE_BEGIN
42
43 /* Float Pi variations */
44
45 /* Division */
46 #ifndef M_PI_F
47 #define M_PI_F          ((float)3.14159265358979323846264338327950288)          /* pi */
48 #endif
49 #ifndef M_PI_2_F
50 #define M_PI_2_F        ((float)1.57079632679489661923132169163975144)          /* pi/2 */
51 #endif
52 #ifndef M_PI_4_F
53 #define M_PI_4_F        ((float)0.785398163397448309615660845819875721)         /* pi/4 */
54 #endif
55 #ifndef M_1_PI_F
56 #define M_1_PI_F        ((float)0.318309886183790671537767526745028724)         /* 1/pi */
57 #endif
58 #ifndef M_2_PI_F
59 #define M_2_PI_F        ((float)0.636619772367581343075535053490057448)         /* 2/pi */
60 #endif
61
62 /* Multiplication */
63 #ifndef M_2PI_F
64 #define M_2PI_F         ((float)6.283185307179586476925286766559005768)         /* 2*pi */
65 #endif
66 #ifndef M_4PI_F
67 #define M_4PI_F         ((float)12.56637061435917295385057353311801153)         /* 4*pi */
68 #endif
69
70 /* Float sqrt variations */
71
72 #ifndef M_SQRT2_F
73 #define M_SQRT2_F       ((float)1.41421356237309504880)                                         /* sqrt(2) */
74 #endif
75
76
77 /* Scalar */
78
79 #ifdef _WIN32
80
81 #ifndef __KERNEL_GPU__
82
83 #if(!defined(FREE_WINDOWS))
84 #define copysignf(x, y) ((float)_copysign(x, y))
85 #define hypotf(x, y) _hypotf(x, y)
86 #define isnan(x) _isnan(x)
87 #define isfinite(x) _finite(x)
88 #endif
89
90 #endif
91
92 #ifndef __KERNEL_OPENCL__
93
94 __device_inline float fmaxf(float a, float b)
95 {
96         return (a > b)? a: b;
97 }
98
99 __device_inline float fminf(float a, float b)
100 {
101         return (a < b)? a: b;
102 }
103
104 #endif
105
106 #endif
107
108 #ifndef __KERNEL_GPU__
109
110 __device_inline int max(int a, int b)
111 {
112         return (a > b)? a: b;
113 }
114
115 __device_inline int min(int a, int b)
116 {
117         return (a < b)? a: b;
118 }
119
120 __device_inline float max(float a, float b)
121 {
122         return (a > b)? a: b;
123 }
124
125 __device_inline float min(float a, float b)
126 {
127         return (a < b)? a: b;
128 }
129
130 __device_inline double max(double a, double b)
131 {
132         return (a > b)? a: b;
133 }
134
135 __device_inline double min(double a, double b)
136 {
137         return (a < b)? a: b;
138 }
139
140 #endif
141
142 __device_inline float min4(float a, float b, float c, float d)
143 {
144         return min(min(a, b), min(c, d));
145 }
146
147 __device_inline float max4(float a, float b, float c, float d)
148 {
149         return max(max(a, b), max(c, d));
150 }
151
152 #ifndef __KERNEL_OPENCL__
153
154 __device_inline int clamp(int a, int mn, int mx)
155 {
156         return min(max(a, mn), mx);
157 }
158
159 __device_inline float clamp(float a, float mn, float mx)
160 {
161         return min(max(a, mn), mx);
162 }
163
164 #endif
165
166 __device_inline int float_to_int(float f)
167 {
168 #ifdef __KERNEL_SSE2__
169         return _mm_cvtt_ss2si(_mm_load_ss(&f));
170 #else
171         return (int)f;
172 #endif
173 }
174
175 __device_inline int floor_to_int(float f)
176 {
177         return float_to_int(floorf(f));
178 }
179
180 __device_inline int ceil_to_int(float f)
181 {
182         return float_to_int(ceilf(f));
183 }
184
185 __device_inline float signf(float f)
186 {
187         return (f < 0.0f)? -1.0f: 1.0f;
188 }
189
190 __device_inline float nonzerof(float f, float eps)
191 {
192         if(fabsf(f) < eps)
193                 return signf(f)*eps;
194         else
195                 return f;
196 }
197
198 __device_inline float smoothstepf(float f)
199 {
200         float ff = f*f;
201         return (3.0f*ff - 2.0f*ff*f);
202 }
203
204 /* Float2 Vector */
205
206 #ifndef __KERNEL_OPENCL__
207
208 __device_inline bool is_zero(const float2 a)
209 {
210         return (a.x == 0.0f && a.y == 0.0f);
211 }
212
213 #endif
214
215 #ifndef __KERNEL_OPENCL__
216
217 __device_inline float average(const float2 a)
218 {
219         return (a.x + a.y)*(1.0f/2.0f);
220 }
221
222 #endif
223
224 #ifndef __KERNEL_OPENCL__
225
226 __device_inline float2 operator-(const float2 a)
227 {
228         return make_float2(-a.x, -a.y);
229 }
230
231 __device_inline float2 operator*(const float2 a, const float2 b)
232 {
233         return make_float2(a.x*b.x, a.y*b.y);
234 }
235
236 __device_inline float2 operator*(const float2 a, float f)
237 {
238         return make_float2(a.x*f, a.y*f);
239 }
240
241 __device_inline float2 operator*(float f, const float2 a)
242 {
243         return make_float2(a.x*f, a.y*f);
244 }
245
246 __device_inline float2 operator/(float f, const float2 a)
247 {
248         return make_float2(f/a.x, f/a.y);
249 }
250
251 __device_inline float2 operator/(const float2 a, float f)
252 {
253         float invf = 1.0f/f;
254         return make_float2(a.x*invf, a.y*invf);
255 }
256
257 __device_inline float2 operator/(const float2 a, const float2 b)
258 {
259         return make_float2(a.x/b.x, a.y/b.y);
260 }
261
262 __device_inline float2 operator+(const float2 a, const float2 b)
263 {
264         return make_float2(a.x+b.x, a.y+b.y);
265 }
266
267 __device_inline float2 operator-(const float2 a, const float2 b)
268 {
269         return make_float2(a.x-b.x, a.y-b.y);
270 }
271
272 __device_inline float2 operator+=(float2& a, const float2 b)
273 {
274         return a = a + b;
275 }
276
277 __device_inline float2 operator*=(float2& a, const float2 b)
278 {
279         return a = a * b;
280 }
281
282 __device_inline float2 operator*=(float2& a, float f)
283 {
284         return a = a * f;
285 }
286
287 __device_inline float2 operator/=(float2& a, const float2 b)
288 {
289         return a = a / b;
290 }
291
292 __device_inline float2 operator/=(float2& a, float f)
293 {
294         float invf = 1.0f/f;
295         return a = a * invf;
296 }
297
298
299 __device_inline float dot(const float2 a, const float2 b)
300 {
301         return a.x*b.x + a.y*b.y;
302 }
303
304 __device_inline float cross(const float2 a, const float2 b)
305 {
306         return (a.x*b.y - a.y*b.x);
307 }
308
309 #endif
310
311 #ifndef __KERNEL_OPENCL__
312
313 __device_inline bool operator==(const int2 a, const int2 b)
314 {
315         return (a.x == b.x && a.y == b.y);
316 }
317
318 __device_inline float len(const float2 a)
319 {
320         return sqrtf(dot(a, a));
321 }
322
323 __device_inline float2 normalize(const float2 a)
324 {
325         return a/len(a);
326 }
327
328 __device_inline float2 normalize_len(const float2 a, float *t)
329 {
330         *t = len(a);
331         return a/(*t);
332 }
333
334 __device_inline bool operator==(const float2 a, const float2 b)
335 {
336         return (a.x == b.x && a.y == b.y);
337 }
338
339 __device_inline bool operator!=(const float2 a, const float2 b)
340 {
341         return !(a == b);
342 }
343
344 __device_inline float2 min(float2 a, float2 b)
345 {
346         return make_float2(min(a.x, b.x), min(a.y, b.y));
347 }
348
349 __device_inline float2 max(float2 a, float2 b)
350 {
351         return make_float2(max(a.x, b.x), max(a.y, b.y));
352 }
353
354 __device_inline float2 clamp(float2 a, float2 mn, float2 mx)
355 {
356         return min(max(a, mn), mx);
357 }
358
359 __device_inline float2 fabs(float2 a)
360 {
361         return make_float2(fabsf(a.x), fabsf(a.y));
362 }
363
364 __device_inline float2 as_float2(const float4 a)
365 {
366         return make_float2(a.x, a.y);
367 }
368
369 #endif
370
371 #ifndef __KERNEL_GPU__
372
373 __device_inline void print_float2(const char *label, const float2& a)
374 {
375         printf("%s: %.8f %.8f\n", label, (double)a.x, (double)a.y);
376 }
377
378 #endif
379
380 #ifndef __KERNEL_OPENCL__
381
382 __device_inline float2 interp(float2 a, float2 b, float t)
383 {
384         return a + t*(b - a);
385 }
386
387 #endif
388
389 /* Float3 Vector */
390
391 #ifndef __KERNEL_OPENCL__
392
393 __device_inline float3 operator-(const float3 a)
394 {
395         return make_float3(-a.x, -a.y, -a.z);
396 }
397
398 __device_inline float3 operator*(const float3 a, const float3 b)
399 {
400         return make_float3(a.x*b.x, a.y*b.y, a.z*b.z);
401 }
402
403 __device_inline float3 operator*(const float3 a, float f)
404 {
405         return make_float3(a.x*f, a.y*f, a.z*f);
406 }
407
408 __device_inline float3 operator*(float f, const float3 a)
409 {
410         return make_float3(a.x*f, a.y*f, a.z*f);
411 }
412
413 __device_inline float3 operator/(float f, const float3 a)
414 {
415         return make_float3(f/a.x, f/a.y, f/a.z);
416 }
417
418 __device_inline float3 operator/(const float3 a, float f)
419 {
420         float invf = 1.0f/f;
421         return make_float3(a.x*invf, a.y*invf, a.z*invf);
422 }
423
424 __device_inline float3 operator/(const float3 a, const float3 b)
425 {
426         return make_float3(a.x/b.x, a.y/b.y, a.z/b.z);
427 }
428
429 __device_inline float3 operator+(const float3 a, const float3 b)
430 {
431         return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
432 }
433
434 __device_inline float3 operator-(const float3 a, const float3 b)
435 {
436         return make_float3(a.x-b.x, a.y-b.y, a.z-b.z);
437 }
438
439 __device_inline float3 operator+=(float3& a, const float3 b)
440 {
441         return a = a + b;
442 }
443
444 __device_inline float3 operator*=(float3& a, const float3 b)
445 {
446         return a = a * b;
447 }
448
449 __device_inline float3 operator*=(float3& a, float f)
450 {
451         return a = a * f;
452 }
453
454 __device_inline float3 operator/=(float3& a, const float3 b)
455 {
456         return a = a / b;
457 }
458
459 __device_inline float3 operator/=(float3& a, float f)
460 {
461         float invf = 1.0f/f;
462         return a = a * invf;
463 }
464
465 __device_inline float dot(const float3 a, const float3 b)
466 {
467         return a.x*b.x + a.y*b.y + a.z*b.z;
468 }
469
470 __device_inline float3 cross(const float3 a, const float3 b)
471 {
472         float3 r = make_float3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x);
473         return r;
474 }
475
476 #endif
477
478 __device_inline float len(const float3 a)
479 {
480         return sqrtf(dot(a, a));
481 }
482
483 __device_inline float len_squared(const float3 a)
484 {
485         return dot(a, a);
486 }
487
488 #ifndef __KERNEL_OPENCL__
489
490 __device_inline float3 normalize(const float3 a)
491 {
492         return a/len(a);
493 }
494
495 #endif
496
497 __device_inline float3 normalize_len(const float3 a, float *t)
498 {
499         *t = len(a);
500         return a/(*t);
501 }
502
503 #ifndef __KERNEL_OPENCL__
504
505 __device_inline bool operator==(const float3 a, const float3 b)
506 {
507 #ifdef __KERNEL_SSE__
508         return (_mm_movemask_ps(_mm_cmpeq_ps(a.m128, b.m128)) & 7) == 7;
509 #else
510         return (a.x == b.x && a.y == b.y && a.z == b.z);
511 #endif
512 }
513
514 __device_inline bool operator!=(const float3 a, const float3 b)
515 {
516         return !(a == b);
517 }
518
519 __device_inline float3 min(float3 a, float3 b)
520 {
521 #ifdef __KERNEL_SSE__
522         return _mm_min_ps(a.m128, b.m128);
523 #else
524         return make_float3(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z));
525 #endif
526 }
527
528 __device_inline float3 max(float3 a, float3 b)
529 {
530 #ifdef __KERNEL_SSE__
531         return _mm_max_ps(a.m128, b.m128);
532 #else
533         return make_float3(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z));
534 #endif
535 }
536
537 __device_inline float3 clamp(float3 a, float3 mn, float3 mx)
538 {
539         return min(max(a, mn), mx);
540 }
541
542 __device_inline float3 fabs(float3 a)
543 {
544 #ifdef __KERNEL_SSE__
545         __m128 mask = _mm_castsi128_ps(_mm_set1_epi32(0x7fffffff));
546         return _mm_and_ps(a.m128, mask);
547 #else
548         return make_float3(fabsf(a.x), fabsf(a.y), fabsf(a.z));
549 #endif
550 }
551
552 #endif
553
554 __device_inline float3 float2_to_float3(const float2 a)
555 {
556         return make_float3(a.x, a.y, 0.0f);
557 }
558
559 __device_inline float3 float4_to_float3(const float4 a)
560 {
561         return make_float3(a.x, a.y, a.z);
562 }
563
564 __device_inline float4 float3_to_float4(const float3 a)
565 {
566         return make_float4(a.x, a.y, a.z, 1.0f);
567 }
568
569 #ifndef __KERNEL_GPU__
570
571 __device_inline void print_float3(const char *label, const float3& a)
572 {
573         printf("%s: %.8f %.8f %.8f\n", label, (double)a.x, (double)a.y, (double)a.z);
574 }
575
576 __device_inline float3 rcp(const float3& a)
577 {
578 #ifdef __KERNEL_SSE__
579         float4 r = _mm_rcp_ps(a.m128);
580         return _mm_sub_ps(_mm_add_ps(r, r), _mm_mul_ps(_mm_mul_ps(r, r), a));
581 #else
582         return make_float3(1.0f/a.x, 1.0f/a.y, 1.0f/a.z);
583 #endif
584 }
585
586 #endif
587
588 __device_inline float3 interp(float3 a, float3 b, float t)
589 {
590         return a + t*(b - a);
591 }
592
593 __device_inline bool is_zero(const float3 a)
594 {
595 #ifdef __KERNEL_SSE__
596         return a == make_float3(0.0f);
597 #else
598         return (a.x == 0.0f && a.y == 0.0f && a.z == 0.0f);
599 #endif
600 }
601
602 __device_inline float reduce_add(const float3 a)
603 {
604 #ifdef __KERNEL_SSE__
605         return (a.x + a.y + a.z);
606 #else
607         return (a.x + a.y + a.z);
608 #endif
609 }
610
611 __device_inline float average(const float3 a)
612 {
613         return reduce_add(a)*(1.0f/3.0f);
614 }
615
616 /* Float4 Vector */
617
618 #ifdef __KERNEL_SSE__
619
620 template<size_t index_0, size_t index_1, size_t index_2, size_t index_3> __forceinline const float4 shuffle(const float4& b)
621 {
622         return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(b), _MM_SHUFFLE(index_3, index_2, index_1, index_0)));
623 }
624
625 template<> __forceinline const float4 shuffle<0, 0, 2, 2>(const float4& b)
626 {
627         return _mm_moveldup_ps(b);
628 }
629
630 template<> __forceinline const float4 shuffle<1, 1, 3, 3>(const float4& b)
631 {
632         return _mm_movehdup_ps(b);
633 }
634
635 template<> __forceinline const float4 shuffle<0, 1, 0, 1>(const float4& b)
636 {
637         return _mm_castpd_ps(_mm_movedup_pd(_mm_castps_pd(b)));
638 }
639
640 #endif
641
642 #ifndef __KERNEL_OPENCL__
643
644 __device_inline float4 operator-(const float4& a)
645 {
646 #ifdef __KERNEL_SSE__
647         __m128 mask = _mm_castsi128_ps(_mm_set1_epi32(0x80000000));
648         return _mm_xor_ps(a.m128, mask);
649 #else
650         return make_float4(-a.x, -a.y, -a.z, -a.w);
651 #endif
652 }
653
654 __device_inline float4 operator*(const float4& a, const float4& b)
655 {
656 #ifdef __KERNEL_SSE__
657         return _mm_mul_ps(a.m128, b.m128);
658 #else
659         return make_float4(a.x*b.x, a.y*b.y, a.z*b.z, a.w*b.w);
660 #endif
661 }
662
663 __device_inline float4 operator*(const float4& a, float f)
664 {
665 #ifdef __KERNEL_SSE__
666         return a * make_float4(f);
667 #else
668         return make_float4(a.x*f, a.y*f, a.z*f, a.w*f);
669 #endif
670 }
671
672 __device_inline float4 operator*(float f, const float4& a)
673 {
674         return a * f;
675 }
676
677 __device_inline float4 rcp(const float4& a)
678 {
679 #ifdef __KERNEL_SSE__
680         float4 r = _mm_rcp_ps(a.m128);
681         return _mm_sub_ps(_mm_add_ps(r, r), _mm_mul_ps(_mm_mul_ps(r, r), a));
682 #else
683         return make_float4(1.0f/a.x, 1.0f/a.y, 1.0f/a.z, 1.0f/a.w);
684 #endif
685 }
686
687 __device_inline float4 operator/(const float4& a, float f)
688 {
689         return a * (1.0f/f);
690 }
691
692 __device_inline float4 operator/(const float4& a, const float4& b)
693 {
694 #ifdef __KERNEL_SSE__
695         return a * rcp(b);
696 #else
697         return make_float4(a.x/b.x, a.y/b.y, a.z/b.z, a.w/b.w);
698 #endif
699
700 }
701
702 __device_inline float4 operator+(const float4& a, const float4& b)
703 {
704 #ifdef __KERNEL_SSE__
705         return _mm_add_ps(a.m128, b.m128);
706 #else
707         return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
708 #endif
709 }
710
711 __device_inline float4 operator-(const float4& a, const float4& b)
712 {
713 #ifdef __KERNEL_SSE__
714         return _mm_sub_ps(a.m128, b.m128);
715 #else
716         return make_float4(a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w);
717 #endif
718 }
719
720 __device_inline float4 operator+=(float4& a, const float4& b)
721 {
722         return a = a + b;
723 }
724
725 __device_inline float4 operator*=(float4& a, const float4& b)
726 {
727         return a = a * b;
728 }
729
730 __device_inline float4 operator/=(float4& a, float f)
731 {
732         return a = a / f;
733 }
734
735 __device_inline int4 operator<(const float4& a, const float4& b)
736 {
737 #ifdef __KERNEL_SSE__
738         return _mm_cvtps_epi32(_mm_cmplt_ps(a.m128, b.m128)); /* todo: avoid cvt */
739 #else
740         return make_int4(a.x < b.x, a.y < b.y, a.z < b.z, a.w < b.w);
741 #endif
742 }
743
744 __device_inline int4 operator>=(float4 a, float4 b)
745 {
746 #ifdef __KERNEL_SSE__
747         return _mm_cvtps_epi32(_mm_cmpge_ps(a.m128, b.m128)); /* todo: avoid cvt */
748 #else
749         return make_int4(a.x >= b.x, a.y >= b.y, a.z >= b.z, a.w >= b.w);
750 #endif
751 }
752
753 __device_inline int4 operator<=(const float4& a, const float4& b)
754 {
755 #ifdef __KERNEL_SSE__
756         return _mm_cvtps_epi32(_mm_cmple_ps(a.m128, b.m128)); /* todo: avoid cvt */
757 #else
758         return make_int4(a.x <= b.x, a.y <= b.y, a.z <= b.z, a.w <= b.w);
759 #endif
760 }
761
762 __device_inline bool operator==(const float4 a, const float4 b)
763 {
764 #ifdef __KERNEL_SSE__
765         return (_mm_movemask_ps(_mm_cmpeq_ps(a.m128, b.m128)) & 15) == 15;
766 #else
767         return (a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w);
768 #endif
769 }
770
771 __device_inline float4 cross(const float4& a, const float4& b)
772 {
773 #ifdef __KERNEL_SSE__
774         return (shuffle<1,2,0,0>(a)*shuffle<2,0,1,0>(b)) - (shuffle<2,0,1,0>(a)*shuffle<1,2,0,0>(b));
775 #else
776         return make_float4(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x, 0.0f);
777 #endif
778 }
779
780 __device_inline bool is_zero(const float4& a)
781 {
782 #ifdef __KERNEL_SSE__
783         return a == make_float4(0.0f);
784 #else
785         return (a.x == 0.0f && a.y == 0.0f && a.z == 0.0f && a.w == 0.0f);
786 #endif
787 }
788
789 __device_inline float reduce_add(const float4& a)
790 {
791 #ifdef __KERNEL_SSE__
792         float4 h = shuffle<1,0,3,2>(a) + a;
793         return _mm_cvtss_f32(shuffle<2,3,0,1>(h) + h); /* todo: efficiency? */
794 #else
795         return ((a.x + a.y) + (a.z + a.w));
796 #endif
797 }
798
799 __device_inline float average(const float4& a)
800 {
801         return reduce_add(a) * 0.25f;
802 }
803
804 __device_inline float dot(const float4& a, const float4& b)
805 {
806         return reduce_add(a * b);
807 }
808
809 __device_inline float len(const float4 a)
810 {
811         return sqrtf(dot(a, a));
812 }
813
814 __device_inline float4 normalize(const float4 a)
815 {
816         return a/len(a);
817 }
818
819 __device_inline float4 min(float4 a, float4 b)
820 {
821 #ifdef __KERNEL_SSE__
822         return _mm_min_ps(a.m128, b.m128);
823 #else
824         return make_float4(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z), min(a.w, b.w));
825 #endif
826 }
827
828 __device_inline float4 max(float4 a, float4 b)
829 {
830 #ifdef __KERNEL_SSE__
831         return _mm_max_ps(a.m128, b.m128);
832 #else
833         return make_float4(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z), max(a.w, b.w));
834 #endif
835 }
836
837 #endif
838
839 #ifndef __KERNEL_GPU__
840
841 __device_inline float4 select(const int4& mask, const float4& a, const float4& b)
842 {
843 #ifdef __KERNEL_SSE__
844         /* blendv is sse4, and apparently broken on vs2008 */
845         return _mm_or_ps(_mm_and_ps(_mm_cvtepi32_ps(mask), a), _mm_andnot_ps(_mm_cvtepi32_ps(mask), b)); /* todo: avoid cvt */
846 #else
847         return make_float4((mask.x)? a.x: b.x, (mask.y)? a.y: b.y, (mask.z)? a.z: b.z, (mask.w)? a.w: b.w);
848 #endif
849 }
850
851 __device_inline float4 reduce_min(const float4& a)
852 {
853 #ifdef __KERNEL_SSE__
854         float4 h = min(shuffle<1,0,3,2>(a), a);
855         return min(shuffle<2,3,0,1>(h), h);
856 #else
857         return make_float4(min(min(a.x, a.y), min(a.z, a.w)));
858 #endif
859 }
860
861 __device_inline float4 reduce_max(const float4& a)
862 {
863 #ifdef __KERNEL_SSE__
864         float4 h = max(shuffle<1,0,3,2>(a), a);
865         return max(shuffle<2,3,0,1>(h), h);
866 #else
867         return make_float4(max(max(a.x, a.y), max(a.z, a.w)));
868 #endif
869 }
870
871 #if 0
872 __device_inline float4 reduce_add(const float4& a)
873 {
874 #ifdef __KERNEL_SSE__
875         float4 h = shuffle<1,0,3,2>(a) + a;
876         return shuffle<2,3,0,1>(h) + h;
877 #else
878         return make_float4((a.x + a.y) + (a.z + a.w));
879 #endif
880 }
881 #endif
882
883 __device_inline void print_float4(const char *label, const float4& a)
884 {
885         printf("%s: %.8f %.8f %.8f %.8f\n", label, (double)a.x, (double)a.y, (double)a.z, (double)a.w);
886 }
887
888 #endif
889
890 /* Int3 */
891
892 #ifndef __KERNEL_OPENCL__
893
894 __device_inline int3 min(int3 a, int3 b)
895 {
896 #ifdef __KERNEL_SSE__
897         return _mm_min_epi32(a.m128, b.m128);
898 #else
899         return make_int3(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z));
900 #endif
901 }
902
903 __device_inline int3 max(int3 a, int3 b)
904 {
905 #ifdef __KERNEL_SSE__
906         return _mm_max_epi32(a.m128, b.m128);
907 #else
908         return make_int3(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z));
909 #endif
910 }
911
912 __device_inline int3 clamp(const int3& a, int mn, int mx)
913 {
914 #ifdef __KERNEL_SSE__
915         return min(max(a, make_int3(mn)), make_int3(mx));
916 #else
917         return make_int3(clamp(a.x, mn, mx), clamp(a.y, mn, mx), clamp(a.z, mn, mx));
918 #endif
919 }
920
921 __device_inline int3 clamp(const int3& a, int3& mn, int mx)
922 {
923 #ifdef __KERNEL_SSE__
924         return min(max(a, mn), make_int3(mx));
925 #else
926         return make_int3(clamp(a.x, mn.x, mx), clamp(a.y, mn.y, mx), clamp(a.z, mn.z, mx));
927 #endif
928 }
929
930 #endif
931
932 #ifndef __KERNEL_GPU__
933
934 __device_inline void print_int3(const char *label, const int3& a)
935 {
936         printf("%s: %d %d %d\n", label, a.x, a.y, a.z);
937 }
938
939 #endif
940
941 /* Int4 */
942
943 #ifndef __KERNEL_GPU__
944
945 __device_inline int4 operator+(const int4& a, const int4& b)
946 {
947 #ifdef __KERNEL_SSE__
948         return _mm_add_epi32(a.m128, b.m128);
949 #else
950         return make_int4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
951 #endif
952 }
953
954 __device_inline int4 operator+=(int4& a, const int4& b)
955 {
956         return a = a + b;
957 }
958
959 __device_inline int4 operator>>(const int4& a, int i)
960 {
961 #ifdef __KERNEL_SSE__
962         return _mm_srai_epi32(a.m128, i);
963 #else
964         return make_int4(a.x >> i, a.y >> i, a.z >> i, a.w >> i);
965 #endif
966 }
967
968 __device_inline int4 min(int4 a, int4 b)
969 {
970 #ifdef __KERNEL_SSE__
971         return _mm_min_epi32(a.m128, b.m128);
972 #else
973         return make_int4(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z), min(a.w, b.w));
974 #endif
975 }
976
977 __device_inline int4 max(int4 a, int4 b)
978 {
979 #ifdef __KERNEL_SSE__
980         return _mm_max_epi32(a.m128, b.m128);
981 #else
982         return make_int4(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z), max(a.w, b.w));
983 #endif
984 }
985
986 __device_inline int4 clamp(const int4& a, const int4& mn, const int4& mx)
987 {
988         return min(max(a, mn), mx);
989 }
990
991 __device_inline int4 select(const int4& mask, const int4& a, const int4& b)
992 {
993 #ifdef __KERNEL_SSE__
994         __m128 m = _mm_cvtepi32_ps(mask);
995         return _mm_castps_si128(_mm_or_ps(_mm_and_ps(m, _mm_castsi128_ps(a)), _mm_andnot_ps(m, _mm_castsi128_ps(b)))); /* todo: avoid cvt */
996 #else
997         return make_int4((mask.x)? a.x: b.x, (mask.y)? a.y: b.y, (mask.z)? a.z: b.z, (mask.w)? a.w: b.w);
998 #endif
999 }
1000
1001 __device_inline void print_int4(const char *label, const int4& a)
1002 {
1003         printf("%s: %d %d %d %d\n", label, a.x, a.y, a.z, a.w);
1004 }
1005
1006 #endif
1007
1008 /* Int/Float conversion */
1009
1010 #ifndef __KERNEL_OPENCL__
1011
1012 __device_inline int as_int(uint i)
1013 {
1014         union { uint ui; int i; } u;
1015         u.ui = i;
1016         return u.i;
1017 }
1018
1019 __device_inline uint as_uint(int i)
1020 {
1021         union { uint ui; int i; } u;
1022         u.i = i;
1023         return u.ui;
1024 }
1025
1026 __device_inline uint as_uint(float f)
1027 {
1028         union { uint i; float f; } u;
1029         u.f = f;
1030         return u.i;
1031 }
1032
1033 __device_inline int __float_as_int(float f)
1034 {
1035         union { int i; float f; } u;
1036         u.f = f;
1037         return u.i;
1038 }
1039
1040 __device_inline float __int_as_float(int i)
1041 {
1042         union { int i; float f; } u;
1043         u.i = i;
1044         return u.f;
1045 }
1046
1047 __device_inline uint __float_as_uint(float f)
1048 {
1049         union { uint i; float f; } u;
1050         u.f = f;
1051         return u.i;
1052 }
1053
1054 __device_inline float __uint_as_float(uint i)
1055 {
1056         union { uint i; float f; } u;
1057         u.i = i;
1058         return u.f;
1059 }
1060
1061 /* Interpolation */
1062
1063 template<class A, class B> A lerp(const A& a, const A& b, const B& t)
1064 {
1065         return (A)(a * ((B)1 - t) + b * t);
1066 }
1067
1068 /* Triangle */
1069
1070 __device_inline float triangle_area(const float3 v1, const float3 v2, const float3 v3)
1071 {
1072         return len(cross(v3 - v2, v1 - v2))*0.5f;
1073 }
1074
1075 #endif
1076
1077 /* Orthonormal vectors */
1078
1079 __device_inline void make_orthonormals(const float3 N, float3 *a, float3 *b)
1080 {
1081         if(N.x != N.y || N.x != N.z)
1082                 *a = make_float3(N.z-N.y, N.x-N.z, N.y-N.x);  //(1,1,1)x N
1083         else
1084                 *a = make_float3(N.z-N.y, N.x+N.z, -N.y-N.x);  //(-1,1,1)x N
1085
1086         *a = normalize(*a);
1087         *b = cross(N, *a);
1088 }
1089
1090 /* Color division */
1091
1092 __device_inline float3 safe_divide_color(float3 a, float3 b)
1093 {
1094         float x, y, z;
1095
1096         x = (b.x != 0.0f)? a.x/b.x: 0.0f;
1097         y = (b.y != 0.0f)? a.y/b.y: 0.0f;
1098         z = (b.z != 0.0f)? a.z/b.z: 0.0f;
1099
1100         return make_float3(x, y, z);
1101 }
1102
1103 /* Rotation of point around axis and angle */
1104
1105 __device_inline float3 rotate_around_axis(float3 p, float3 axis, float angle)
1106 {
1107         float costheta = cosf(angle);
1108         float sintheta = sinf(angle);
1109         float3 r;
1110
1111         r.x = ((costheta + (1 - costheta) * axis.x * axis.x) * p.x) +
1112                 (((1 - costheta) * axis.x * axis.y - axis.z * sintheta) * p.y) +
1113                 (((1 - costheta) * axis.x * axis.z + axis.y * sintheta) * p.z);
1114
1115         r.y = (((1 - costheta) * axis.x * axis.y + axis.z * sintheta) * p.x) +
1116                 ((costheta + (1 - costheta) * axis.y * axis.y) * p.y) +
1117                 (((1 - costheta) * axis.y * axis.z - axis.x * sintheta) * p.z);
1118
1119         r.z = (((1 - costheta) * axis.x * axis.z - axis.y * sintheta) * p.x) +
1120                 (((1 - costheta) * axis.y * axis.z + axis.x * sintheta) * p.y) +
1121                 ((costheta + (1 - costheta) * axis.z * axis.z) * p.z);
1122
1123         return r;
1124 }
1125
1126 /* NaN-safe math ops */
1127
1128 __device float safe_asinf(float a)
1129 {
1130         if(a <= -1.0f)
1131                 return -M_PI_2_F;
1132         else if(a >= 1.0f)
1133                 return M_PI_2_F;
1134
1135         return asinf(a);
1136 }
1137
1138 __device float safe_acosf(float a)
1139 {
1140         if(a <= -1.0f)
1141                 return M_PI_F;
1142         else if(a >= 1.0f)
1143                 return 0.0f;
1144
1145         return acosf(a);
1146 }
1147
1148 __device float compatible_powf(float x, float y)
1149 {
1150         /* GPU pow doesn't accept negative x, do manual checks here */
1151         if(x < 0.0f) {
1152                 if(fmodf(-y, 2.0f) == 0.0f)
1153                         return powf(-x, y);
1154                 else
1155                         return -powf(-x, y);
1156         }
1157         else if(x == 0.0f)
1158                 return 0.0f;
1159
1160         return powf(x, y);
1161 }
1162
1163 __device float safe_powf(float a, float b)
1164 {
1165         if(b == 0.0f)
1166                 return 1.0f;
1167         if(a == 0.0f)
1168                 return 0.0f;
1169         if(a < 0.0f && b != (int)b)
1170                 return 0.0f;
1171         
1172         return compatible_powf(a, b);
1173 }
1174
1175 __device float safe_logf(float a, float b)
1176 {
1177         if(a < 0.0f || b < 0.0f)
1178                 return 0.0f;
1179
1180         return logf(a)/logf(b);
1181 }
1182
1183 __device float safe_divide(float a, float b)
1184 {
1185         return (b != 0.0f)? a/b: 0.0f;
1186 }
1187
1188 __device float safe_modulo(float a, float b)
1189 {
1190         return (b != 0.0f)? fmodf(a, b): 0.0f;
1191 }
1192
1193 /* Ray Intersection */
1194
1195 __device bool ray_sphere_intersect(
1196         float3 ray_P, float3 ray_D, float ray_t,
1197         float3 sphere_P, float sphere_radius,
1198         float3 *isect_P, float *isect_t)
1199 {
1200         float3 d = sphere_P - ray_P;
1201         float radiussq = sphere_radius*sphere_radius;
1202         float tsq = dot(d, d);
1203
1204         if(tsq > radiussq) { /* ray origin outside sphere */
1205                 float tp = dot(d, ray_D);
1206
1207                 if(tp < 0.0f) /* dir points away from sphere */
1208                         return false;
1209
1210                 float dsq = tsq - tp*tp; /* pythagoras */
1211
1212                 if(dsq > radiussq) /* closest point on ray outside sphere */
1213                         return false;
1214
1215                 float t = tp - sqrtf(radiussq - dsq); /* pythagoras */
1216
1217                 if(t < ray_t) {
1218                         *isect_t = t;
1219                         *isect_P = ray_P + ray_D*t;
1220                         return true;
1221                 }
1222         }
1223
1224         return false;
1225 }
1226
1227 __device bool ray_aligned_disk_intersect(
1228         float3 ray_P, float3 ray_D, float ray_t,
1229         float3 disk_P, float disk_radius,
1230         float3 *isect_P, float *isect_t)
1231 {
1232         /* aligned disk normal */
1233         float disk_t;
1234         float3 disk_N = normalize_len(ray_P - disk_P, &disk_t);
1235         float div = dot(ray_D, disk_N);
1236
1237         if(div == 0.0f)
1238                 return false;
1239
1240         /* compute t to intersection point */
1241         float t = -disk_t/div;
1242         if(t < 0.0f || t > ray_t)
1243                 return false;
1244         
1245         /* test if within radius */
1246         float3 P = ray_P + ray_D*t;
1247         if(len_squared(P - disk_P) > disk_radius*disk_radius)
1248                 return false;
1249
1250         *isect_P = P;
1251         *isect_t = t;
1252
1253         return true;
1254 }
1255
1256 __device bool ray_triangle_intersect(
1257         float3 ray_P, float3 ray_D, float ray_t,
1258         float3 v0, float3 v1, float3 v2,
1259         float3 *isect_P, float *isect_t)
1260 {
1261         /* Calculate intersection */
1262         float3 e1 = v1 - v0;
1263         float3 e2 = v2 - v0;
1264         float3 s1 = cross(ray_D, e2);
1265
1266         const float divisor = dot(s1, e1);
1267         if(divisor == 0.0f)
1268                 return false;
1269
1270         const float invdivisor = 1.0f/divisor;
1271
1272         /* compute first barycentric coordinate */
1273         const float3 d = ray_P - v0;
1274         const float u = dot(d, s1)*invdivisor;
1275         if(u < 0.0f)
1276                 return false;
1277
1278         /* Compute second barycentric coordinate */
1279         const float3 s2 = cross(d, e1);
1280         const float v = dot(ray_D, s2)*invdivisor;
1281         if(v < 0.0f)
1282                 return false;
1283
1284         const float b0 = 1.0f - u - v;
1285         if(b0 < 0.0f)
1286                 return false;
1287
1288         /* compute t to intersection point */
1289         const float t = dot(e2, s2)*invdivisor;
1290         if(t < 0.0f || t > ray_t)
1291                 return false;
1292
1293         *isect_t = t;
1294         *isect_P = ray_P + ray_D*t;
1295
1296         return true;
1297 }
1298
1299 __device bool ray_quad_intersect(
1300         float3 ray_P, float3 ray_D, float ray_t,
1301         float3 quad_P, float3 quad_u, float3 quad_v,
1302         float3 *isect_P, float *isect_t)
1303 {
1304         float3 v0 = quad_P - quad_u*0.5f - quad_v*0.5f;
1305         float3 v1 = quad_P + quad_u*0.5f - quad_v*0.5f;
1306         float3 v2 = quad_P + quad_u*0.5f + quad_v*0.5f;
1307         float3 v3 = quad_P - quad_u*0.5f + quad_v*0.5f;
1308
1309         if(ray_triangle_intersect(ray_P, ray_D, ray_t, v0, v1, v2, isect_P, isect_t))
1310                 return true;
1311         else if(ray_triangle_intersect(ray_P, ray_D, ray_t, v0, v2, v3, isect_P, isect_t))
1312                 return true;
1313         
1314         return false;
1315 }
1316
1317 CCL_NAMESPACE_END
1318
1319 #endif /* __UTIL_MATH_H__ */
1320