a6bc478ee64d79dbddafde7ae8ce0927be461a94
[blender-staging.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 #define _USE_MATH_DEFINES
30
31 #include <float.h>
32 #include <math.h>
33 #include <stdio.h>
34
35 #endif
36
37 #include "util_types.h"
38
39 CCL_NAMESPACE_BEGIN
40
41 /* Float Pi variations */
42
43 #ifndef M_PI_F
44 #define M_PI_F          ((float)3.14159265358979323846264338327950288)
45 #endif
46 #ifndef M_PI_2_F
47 #define M_PI_2_F        ((float)1.57079632679489661923132169163975144)
48 #endif
49 #ifndef M_PI_4_F
50 #define M_PI_4_F        ((float)0.785398163397448309615660845819875721)
51 #endif
52 #ifndef M_1_PI_F
53 #define M_1_PI_F        ((float)0.318309886183790671537767526745028724)
54 #endif
55 #ifndef M_2_PI_F
56 #define M_2_PI_F        ((float)0.636619772367581343075535053490057448)
57 #endif
58 #ifndef M_SQRT2_F
59 #define M_SQRT2_F       ((float)1.41421356237309504880)
60 #endif
61
62
63 /* Scalar */
64
65 #ifdef _WIN32
66
67 #ifndef __KERNEL_GPU__
68
69 #if(!defined(FREE_WINDOWS))
70 #define copysignf(x, y) ((float)_copysign(x, y))
71 #define hypotf(x, y) _hypotf(x, y)
72 #define isnan(x) _isnan(x)
73 #define isfinite(x) _finite(x)
74 #endif
75
76 #endif
77
78 #ifndef __KERNEL_OPENCL__
79
80 __device_inline float fmaxf(float a, float b)
81 {
82         return (a > b)? a: b;
83 }
84
85 __device_inline float fminf(float a, float b)
86 {
87         return (a < b)? a: b;
88 }
89
90 #endif
91
92 #endif
93
94 #ifndef __KERNEL_GPU__
95
96 __device_inline int max(int a, int b)
97 {
98         return (a > b)? a: b;
99 }
100
101 __device_inline int min(int a, int b)
102 {
103         return (a < b)? a: b;
104 }
105
106 __device_inline float max(float a, float b)
107 {
108         return (a > b)? a: b;
109 }
110
111 __device_inline float min(float a, float b)
112 {
113         return (a < b)? a: b;
114 }
115
116 __device_inline double max(double a, double b)
117 {
118         return (a > b)? a: b;
119 }
120
121 __device_inline double min(double a, double b)
122 {
123         return (a < b)? a: b;
124 }
125
126 #endif
127
128 __device_inline float min4(float a, float b, float c, float d)
129 {
130         return min(min(a, b), min(c, d));
131 }
132
133 __device_inline float max4(float a, float b, float c, float d)
134 {
135         return max(max(a, b), max(c, d));
136 }
137
138 #ifndef __KERNEL_OPENCL__
139
140 __device_inline int clamp(int a, int mn, int mx)
141 {
142         return min(max(a, mn), mx);
143 }
144
145 __device_inline float clamp(float a, float mn, float mx)
146 {
147         return min(max(a, mn), mx);
148 }
149
150 #endif
151
152 __device_inline float signf(float f)
153 {
154         return (f < 0.0f)? -1.0f: 1.0f;
155 }
156
157 __device_inline float nonzerof(float f, float eps)
158 {
159         if(fabsf(f) < eps)
160                 return signf(f)*eps;
161         else
162                 return f;
163 }
164
165 __device_inline float smoothstepf(float f)
166 {
167         float ff = f*f;
168         return (3.0f*ff - 2.0f*ff*f);
169 }
170
171 /* Float2 Vector */
172
173 #ifndef __KERNEL_OPENCL__
174
175 __device_inline bool is_zero(const float2 a)
176 {
177         return (a.x == 0.0f && a.y == 0.0f);
178 }
179
180 #endif
181
182 #ifndef __KERNEL_OPENCL__
183
184 __device_inline float average(const float2 a)
185 {
186         return (a.x + a.y)*(1.0f/2.0f);
187 }
188
189 #endif
190
191 #ifndef __KERNEL_OPENCL__
192
193 __device_inline float2 operator-(const float2 a)
194 {
195         return make_float2(-a.x, -a.y);
196 }
197
198 __device_inline float2 operator*(const float2 a, const float2 b)
199 {
200         return make_float2(a.x*b.x, a.y*b.y);
201 }
202
203 __device_inline float2 operator*(const float2 a, float f)
204 {
205         return make_float2(a.x*f, a.y*f);
206 }
207
208 __device_inline float2 operator*(float f, const float2 a)
209 {
210         return make_float2(a.x*f, a.y*f);
211 }
212
213 __device_inline float2 operator/(float f, const float2 a)
214 {
215         return make_float2(f/a.x, f/a.y);
216 }
217
218 __device_inline float2 operator/(const float2 a, float f)
219 {
220         float invf = 1.0f/f;
221         return make_float2(a.x*invf, a.y*invf);
222 }
223
224 __device_inline float2 operator/(const float2 a, const float2 b)
225 {
226         return make_float2(a.x/b.x, a.y/b.y);
227 }
228
229 __device_inline float2 operator+(const float2 a, const float2 b)
230 {
231         return make_float2(a.x+b.x, a.y+b.y);
232 }
233
234 __device_inline float2 operator-(const float2 a, const float2 b)
235 {
236         return make_float2(a.x-b.x, a.y-b.y);
237 }
238
239 __device_inline float2 operator+=(float2& a, const float2 b)
240 {
241         return a = a + b;
242 }
243
244 __device_inline float2 operator*=(float2& a, const float2 b)
245 {
246         return a = a * b;
247 }
248
249 __device_inline float2 operator*=(float2& a, float f)
250 {
251         return a = a * f;
252 }
253
254 __device_inline float2 operator/=(float2& a, const float2 b)
255 {
256         return a = a / b;
257 }
258
259 __device_inline float2 operator/=(float2& a, float f)
260 {
261         float invf = 1.0f/f;
262         return a = a * invf;
263 }
264
265
266 __device_inline float dot(const float2 a, const float2 b)
267 {
268         return a.x*b.x + a.y*b.y;
269 }
270
271 __device_inline float cross(const float2 a, const float2 b)
272 {
273         return (a.x*b.y - a.y*b.x);
274 }
275
276 #endif
277
278 #ifndef __KERNEL_OPENCL__
279
280 __device_inline float len(const float2 a)
281 {
282         return sqrtf(dot(a, a));
283 }
284
285 __device_inline float2 normalize(const float2 a)
286 {
287         return a/len(a);
288 }
289
290 __device_inline float2 normalize_len(const float2 a, float *t)
291 {
292         *t = len(a);
293         return a/(*t);
294 }
295
296 __device_inline bool operator==(const float2 a, const float2 b)
297 {
298         return (a.x == b.x && a.y == b.y);
299 }
300
301 __device_inline bool operator!=(const float2 a, const float2 b)
302 {
303         return !(a == b);
304 }
305
306 __device_inline float2 min(float2 a, float2 b)
307 {
308         return make_float2(min(a.x, b.x), min(a.y, b.y));
309 }
310
311 __device_inline float2 max(float2 a, float2 b)
312 {
313         return make_float2(max(a.x, b.x), max(a.y, b.y));
314 }
315
316 __device_inline float2 clamp(float2 a, float2 mn, float2 mx)
317 {
318         return min(max(a, mn), mx);
319 }
320
321 __device_inline float2 fabs(float2 a)
322 {
323         return make_float2(fabsf(a.x), fabsf(a.y));
324 }
325
326 __device_inline float2 as_float2(const float4 a)
327 {
328         return make_float2(a.x, a.y);
329 }
330
331 #endif
332
333 #ifndef __KERNEL_GPU__
334
335 __device_inline void print_float2(const char *label, const float2& a)
336 {
337         printf("%s: %.8f %.8f\n", label, (double)a.x, (double)a.y);
338 }
339
340 #endif
341
342 #ifndef __KERNEL_OPENCL__
343
344 __device_inline float2 interp(float2 a, float2 b, float t)
345 {
346         return a + t*(b - a);
347 }
348
349 #endif
350
351 /* Float3 Vector */
352
353 #ifndef __KERNEL_OPENCL__
354
355 __device_inline float3 operator-(const float3 a)
356 {
357         return make_float3(-a.x, -a.y, -a.z);
358 }
359
360 __device_inline float3 operator*(const float3 a, const float3 b)
361 {
362         return make_float3(a.x*b.x, a.y*b.y, a.z*b.z);
363 }
364
365 __device_inline float3 operator*(const float3 a, float f)
366 {
367         return make_float3(a.x*f, a.y*f, a.z*f);
368 }
369
370 __device_inline float3 operator*(float f, const float3 a)
371 {
372         return make_float3(a.x*f, a.y*f, a.z*f);
373 }
374
375 __device_inline float3 operator/(float f, const float3 a)
376 {
377         return make_float3(f/a.x, f/a.y, f/a.z);
378 }
379
380 __device_inline float3 operator/(const float3 a, float f)
381 {
382         float invf = 1.0f/f;
383         return make_float3(a.x*invf, a.y*invf, a.z*invf);
384 }
385
386 __device_inline float3 operator/(const float3 a, const float3 b)
387 {
388         return make_float3(a.x/b.x, a.y/b.y, a.z/b.z);
389 }
390
391 __device_inline float3 operator+(const float3 a, const float3 b)
392 {
393         return make_float3(a.x+b.x, a.y+b.y, a.z+b.z);
394 }
395
396 __device_inline float3 operator-(const float3 a, const float3 b)
397 {
398         return make_float3(a.x-b.x, a.y-b.y, a.z-b.z);
399 }
400
401 __device_inline float3 operator+=(float3& a, const float3 b)
402 {
403         return a = a + b;
404 }
405
406 __device_inline float3 operator*=(float3& a, const float3 b)
407 {
408         return a = a * b;
409 }
410
411 __device_inline float3 operator*=(float3& a, float f)
412 {
413         return a = a * f;
414 }
415
416 __device_inline float3 operator/=(float3& a, const float3 b)
417 {
418         return a = a / b;
419 }
420
421 __device_inline float3 operator/=(float3& a, float f)
422 {
423         float invf = 1.0f/f;
424         return a = a * invf;
425 }
426
427 __device_inline float dot(const float3 a, const float3 b)
428 {
429         return a.x*b.x + a.y*b.y + a.z*b.z;
430 }
431
432 __device_inline float3 cross(const float3 a, const float3 b)
433 {
434         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);
435         return r;
436 }
437
438 #endif
439
440 __device_inline float len(const float3 a)
441 {
442         return sqrtf(dot(a, a));
443 }
444
445 __device_inline float len_squared(const float3 a)
446 {
447         return dot(a, a);
448 }
449
450 #ifndef __KERNEL_OPENCL__
451
452 __device_inline float3 normalize(const float3 a)
453 {
454         return a/len(a);
455 }
456
457 #endif
458
459 __device_inline float3 normalize_len(const float3 a, float *t)
460 {
461         *t = len(a);
462         return a/(*t);
463 }
464
465 #ifndef __KERNEL_OPENCL__
466
467 __device_inline bool operator==(const float3 a, const float3 b)
468 {
469 #ifdef __KERNEL_SSE__
470         return (_mm_movemask_ps(_mm_cmpeq_ps(a.m128, b.m128)) & 7) == 7;
471 #else
472         return (a.x == b.x && a.y == b.y && a.z == b.z);
473 #endif
474 }
475
476 __device_inline bool operator!=(const float3 a, const float3 b)
477 {
478         return !(a == b);
479 }
480
481 __device_inline float3 min(float3 a, float3 b)
482 {
483 #ifdef __KERNEL_SSE__
484         return _mm_min_ps(a.m128, b.m128);
485 #else
486         return make_float3(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z));
487 #endif
488 }
489
490 __device_inline float3 max(float3 a, float3 b)
491 {
492 #ifdef __KERNEL_SSE__
493         return _mm_max_ps(a.m128, b.m128);
494 #else
495         return make_float3(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z));
496 #endif
497 }
498
499 __device_inline float3 clamp(float3 a, float3 mn, float3 mx)
500 {
501         return min(max(a, mn), mx);
502 }
503
504 __device_inline float3 fabs(float3 a)
505 {
506 #ifdef __KERNEL_SSE__
507         __m128 mask = _mm_castsi128_ps(_mm_set1_epi32(0x7fffffff));
508         return _mm_and_ps(a.m128, mask);
509 #else
510         return make_float3(fabsf(a.x), fabsf(a.y), fabsf(a.z));
511 #endif
512 }
513
514 #endif
515
516 __device_inline float3 float2_to_float3(const float2 a)
517 {
518         return make_float3(a.x, a.y, 0.0f);
519 }
520
521 __device_inline float3 float4_to_float3(const float4 a)
522 {
523         return make_float3(a.x, a.y, a.z);
524 }
525
526 __device_inline float4 float3_to_float4(const float3 a)
527 {
528         return make_float4(a.x, a.y, a.z, 1.0f);
529 }
530
531 #ifndef __KERNEL_GPU__
532
533 __device_inline void print_float3(const char *label, const float3& a)
534 {
535         printf("%s: %.8f %.8f %.8f\n", label, (double)a.x, (double)a.y, (double)a.z);
536 }
537
538 __device_inline float3 rcp(const float3& a)
539 {
540 #ifdef __KERNEL_SSE__
541         float4 r = _mm_rcp_ps(a.m128);
542         return _mm_sub_ps(_mm_add_ps(r, r), _mm_mul_ps(_mm_mul_ps(r, r), a));
543 #else
544         return make_float3(1.0f/a.x, 1.0f/a.y, 1.0f/a.z);
545 #endif
546 }
547
548 #endif
549
550 __device_inline float3 interp(float3 a, float3 b, float t)
551 {
552         return a + t*(b - a);
553 }
554
555 __device_inline bool is_zero(const float3 a)
556 {
557 #ifdef __KERNEL_SSE__
558         return a == make_float3(0.0f);
559 #else
560         return (a.x == 0.0f && a.y == 0.0f && a.z == 0.0f);
561 #endif
562 }
563
564 __device_inline float reduce_add(const float3 a)
565 {
566 #ifdef __KERNEL_SSE__
567         return (a.x + a.y + a.z);
568 #else
569         return (a.x + a.y + a.z);
570 #endif
571 }
572
573 __device_inline float average(const float3 a)
574 {
575         return reduce_add(a)*(1.0f/3.0f);
576 }
577
578 /* Float4 Vector */
579
580 #ifdef __KERNEL_SSE__
581
582 template<size_t index_0, size_t index_1, size_t index_2, size_t index_3> __forceinline const float4 shuffle(const float4& b)
583 {
584         return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(b), _MM_SHUFFLE(index_3, index_2, index_1, index_0)));
585 }
586
587 template<> __forceinline const float4 shuffle<0, 0, 2, 2>(const float4& b)
588 {
589         return _mm_moveldup_ps(b);
590 }
591
592 template<> __forceinline const float4 shuffle<1, 1, 3, 3>(const float4& b)
593 {
594         return _mm_movehdup_ps(b);
595 }
596
597 template<> __forceinline const float4 shuffle<0, 1, 0, 1>(const float4& b)
598 {
599         return _mm_castpd_ps(_mm_movedup_pd(_mm_castps_pd(b)));
600 }
601
602 #endif
603
604 #ifndef __KERNEL_OPENCL__
605
606 __device_inline float4 operator-(const float4& a)
607 {
608 #ifdef __KERNEL_SSE__
609         __m128 mask = _mm_castsi128_ps(_mm_set1_epi32(0x80000000));
610         return _mm_xor_ps(a.m128, mask);
611 #else
612         return make_float4(-a.x, -a.y, -a.z, -a.w);
613 #endif
614 }
615
616 __device_inline float4 operator*(const float4& a, const float4& b)
617 {
618 #ifdef __KERNEL_SSE__
619         return _mm_mul_ps(a.m128, b.m128);
620 #else
621         return make_float4(a.x*b.x, a.y*b.y, a.z*b.z, a.w*b.w);
622 #endif
623 }
624
625 __device_inline float4 operator*(const float4& a, float f)
626 {
627 #ifdef __KERNEL_SSE__
628         return a * make_float4(f);
629 #else
630         return make_float4(a.x*f, a.y*f, a.z*f, a.w*f);
631 #endif
632 }
633
634 __device_inline float4 operator*(float f, const float4& a)
635 {
636         return a * f;
637 }
638
639 __device_inline float4 rcp(const float4& a)
640 {
641 #ifdef __KERNEL_SSE__
642         float4 r = _mm_rcp_ps(a.m128);
643         return _mm_sub_ps(_mm_add_ps(r, r), _mm_mul_ps(_mm_mul_ps(r, r), a));
644 #else
645         return make_float4(1.0f/a.x, 1.0f/a.y, 1.0f/a.z, 1.0f/a.w);
646 #endif
647 }
648
649 __device_inline float4 operator/(const float4& a, float f)
650 {
651         return a * (1.0f/f);
652 }
653
654 __device_inline float4 operator/(const float4& a, const float4& b)
655 {
656 #ifdef __KERNEL_SSE__
657         return a * rcp(b);
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
664 __device_inline float4 operator+(const float4& a, const float4& b)
665 {
666 #ifdef __KERNEL_SSE__
667         return _mm_add_ps(a.m128, b.m128);
668 #else
669         return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
670 #endif
671 }
672
673 __device_inline float4 operator-(const float4& a, const float4& b)
674 {
675 #ifdef __KERNEL_SSE__
676         return _mm_sub_ps(a.m128, b.m128);
677 #else
678         return make_float4(a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w);
679 #endif
680 }
681
682 __device_inline float4 operator+=(float4& a, const float4& b)
683 {
684         return a = a + b;
685 }
686
687 __device_inline float4 operator*=(float4& a, const float4& b)
688 {
689         return a = a * b;
690 }
691
692 __device_inline float4 operator/=(float4& a, float f)
693 {
694         return a = a / f;
695 }
696
697 __device_inline int4 operator<(const float4& a, const float4& b)
698 {
699 #ifdef __KERNEL_SSE__
700         return _mm_cvtps_epi32(_mm_cmplt_ps(a.m128, b.m128)); /* todo: avoid cvt */
701 #else
702         return make_int4(a.x < b.x, a.y < b.y, a.z < b.z, a.w < b.w);
703 #endif
704 }
705
706 __device_inline int4 operator>=(float4 a, float4 b)
707 {
708 #ifdef __KERNEL_SSE__
709         return _mm_cvtps_epi32(_mm_cmpge_ps(a.m128, b.m128)); /* todo: avoid cvt */
710 #else
711         return make_int4(a.x >= b.x, a.y >= b.y, a.z >= b.z, a.w >= b.w);
712 #endif
713 }
714
715 __device_inline int4 operator<=(const float4& a, const float4& b)
716 {
717 #ifdef __KERNEL_SSE__
718         return _mm_cvtps_epi32(_mm_cmple_ps(a.m128, b.m128)); /* todo: avoid cvt */
719 #else
720         return make_int4(a.x <= b.x, a.y <= b.y, a.z <= b.z, a.w <= b.w);
721 #endif
722 }
723
724 __device_inline bool operator==(const float4 a, const float4 b)
725 {
726 #ifdef __KERNEL_SSE__
727         return (_mm_movemask_ps(_mm_cmpeq_ps(a.m128, b.m128)) & 15) == 15;
728 #else
729         return (a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w);
730 #endif
731 }
732
733 __device_inline float4 cross(const float4& a, const float4& b)
734 {
735 #ifdef __KERNEL_SSE__
736         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));
737 #else
738         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);
739 #endif
740 }
741
742 __device_inline bool is_zero(const float4& a)
743 {
744 #ifdef __KERNEL_SSE__
745         return a == make_float4(0.0f);
746 #else
747         return (a.x == 0.0f && a.y == 0.0f && a.z == 0.0f && a.w == 0.0f);
748 #endif
749 }
750
751 __device_inline float reduce_add(const float4& a)
752 {
753 #ifdef __KERNEL_SSE__
754         float4 h = shuffle<1,0,3,2>(a) + a;
755         return _mm_cvtss_f32(shuffle<2,3,0,1>(h) + h); /* todo: efficiency? */
756 #else
757         return ((a.x + a.y) + (a.z + a.w));
758 #endif
759 }
760
761 __device_inline float average(const float4& a)
762 {
763         return reduce_add(a) * 0.25f;
764 }
765
766 __device_inline float dot(const float4& a, const float4& b)
767 {
768         return reduce_add(a * b);
769 }
770
771 __device_inline float len(const float4 a)
772 {
773         return sqrtf(dot(a, a));
774 }
775
776 __device_inline float4 normalize(const float4 a)
777 {
778         return a/len(a);
779 }
780
781 __device_inline float4 min(float4 a, float4 b)
782 {
783 #ifdef __KERNEL_SSE__
784         return _mm_min_ps(a.m128, b.m128);
785 #else
786         return make_float4(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z), min(a.w, b.w));
787 #endif
788 }
789
790 __device_inline float4 max(float4 a, float4 b)
791 {
792 #ifdef __KERNEL_SSE__
793         return _mm_max_ps(a.m128, b.m128);
794 #else
795         return make_float4(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z), max(a.w, b.w));
796 #endif
797 }
798
799 #endif
800
801 #ifndef __KERNEL_GPU__
802
803 __device_inline float4 select(const int4& mask, const float4& a, const float4& b)
804 {
805 #ifdef __KERNEL_SSE__
806         /* blendv is sse4, and apparently broken on vs2008 */
807         return _mm_or_ps(_mm_and_ps(_mm_cvtepi32_ps(mask), a), _mm_andnot_ps(_mm_cvtepi32_ps(mask), b)); /* todo: avoid cvt */
808 #else
809         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);
810 #endif
811 }
812
813 __device_inline float4 reduce_min(const float4& a)
814 {
815 #ifdef __KERNEL_SSE__
816         float4 h = min(shuffle<1,0,3,2>(a), a);
817         return min(shuffle<2,3,0,1>(h), h);
818 #else
819         return make_float4(min(min(a.x, a.y), min(a.z, a.w)));
820 #endif
821 }
822
823 __device_inline float4 reduce_max(const float4& a)
824 {
825 #ifdef __KERNEL_SSE__
826         float4 h = max(shuffle<1,0,3,2>(a), a);
827         return max(shuffle<2,3,0,1>(h), h);
828 #else
829         return make_float4(max(max(a.x, a.y), max(a.z, a.w)));
830 #endif
831 }
832
833 #if 0
834 __device_inline float4 reduce_add(const float4& a)
835 {
836 #ifdef __KERNEL_SSE__
837         float4 h = shuffle<1,0,3,2>(a) + a;
838         return shuffle<2,3,0,1>(h) + h;
839 #else
840         return make_float4((a.x + a.y) + (a.z + a.w));
841 #endif
842 }
843 #endif
844
845 __device_inline void print_float4(const char *label, const float4& a)
846 {
847         printf("%s: %.8f %.8f %.8f %.8f\n", label, (double)a.x, (double)a.y, (double)a.z, (double)a.w);
848 }
849
850 #endif
851
852 /* Int3 */
853
854 #ifndef __KERNEL_OPENCL__
855
856 __device_inline int3 min(int3 a, int3 b)
857 {
858 #ifdef __KERNEL_SSE__
859         return _mm_min_epi32(a.m128, b.m128);
860 #else
861         return make_int3(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z));
862 #endif
863 }
864
865 __device_inline int3 max(int3 a, int3 b)
866 {
867 #ifdef __KERNEL_SSE__
868         return _mm_max_epi32(a.m128, b.m128);
869 #else
870         return make_int3(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z));
871 #endif
872 }
873
874 __device_inline int3 clamp(const int3& a, int mn, int mx)
875 {
876 #ifdef __KERNEL_SSE__
877         return min(max(a, make_int3(mn)), make_int3(mx));
878 #else
879         return make_int3(clamp(a.x, mn, mx), clamp(a.y, mn, mx), clamp(a.z, mn, mx));
880 #endif
881 }
882
883 __device_inline int3 clamp(const int3& a, int3& mn, int mx)
884 {
885 #ifdef __KERNEL_SSE__
886         return min(max(a, mn), make_int3(mx));
887 #else
888         return make_int3(clamp(a.x, mn.x, mx), clamp(a.y, mn.y, mx), clamp(a.z, mn.z, mx));
889 #endif
890 }
891
892 #endif
893
894 #ifndef __KERNEL_GPU__
895
896 __device_inline void print_int3(const char *label, const int3& a)
897 {
898         printf("%s: %d %d %d\n", label, a.x, a.y, a.z);
899 }
900
901 #endif
902
903 /* Int4 */
904
905 #ifndef __KERNEL_GPU__
906
907 __device_inline int4 operator+(const int4& a, const int4& b)
908 {
909 #ifdef __KERNEL_SSE__
910         return _mm_add_epi32(a.m128, b.m128);
911 #else
912         return make_int4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
913 #endif
914 }
915
916 __device_inline int4 operator+=(int4& a, const int4& b)
917 {
918         return a = a + b;
919 }
920
921 __device_inline int4 operator>>(const int4& a, int i)
922 {
923 #ifdef __KERNEL_SSE__
924         return _mm_srai_epi32(a.m128, i);
925 #else
926         return make_int4(a.x >> i, a.y >> i, a.z >> i, a.w >> i);
927 #endif
928 }
929
930 __device_inline int4 min(int4 a, int4 b)
931 {
932 #ifdef __KERNEL_SSE__
933         return _mm_min_epi32(a.m128, b.m128);
934 #else
935         return make_int4(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z), min(a.w, b.w));
936 #endif
937 }
938
939 __device_inline int4 max(int4 a, int4 b)
940 {
941 #ifdef __KERNEL_SSE__
942         return _mm_max_epi32(a.m128, b.m128);
943 #else
944         return make_int4(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z), max(a.w, b.w));
945 #endif
946 }
947
948 __device_inline int4 clamp(const int4& a, const int4& mn, const int4& mx)
949 {
950         return min(max(a, mn), mx);
951 }
952
953 __device_inline int4 select(const int4& mask, const int4& a, const int4& b)
954 {
955 #ifdef __KERNEL_SSE__
956         __m128 m = _mm_cvtepi32_ps(mask);
957         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 */
958 #else
959         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);
960 #endif
961 }
962
963 __device_inline void print_int4(const char *label, const int4& a)
964 {
965         printf("%s: %d %d %d %d\n", label, a.x, a.y, a.z, a.w);
966 }
967
968 #endif
969
970 /* Int/Float conversion */
971
972 #ifndef __KERNEL_OPENCL__
973
974 __device_inline unsigned int as_int(uint i)
975 {
976         union { unsigned int ui; int i; } u;
977         u.ui = i;
978         return u.i;
979 }
980
981 __device_inline unsigned int as_uint(int i)
982 {
983         union { unsigned int ui; int i; } u;
984         u.i = i;
985         return u.ui;
986 }
987
988 __device_inline unsigned int as_uint(float f)
989 {
990         union { unsigned int i; float f; } u;
991         u.f = f;
992         return u.i;
993 }
994
995 __device_inline int __float_as_int(float f)
996 {
997         union { int i; float f; } u;
998         u.f = f;
999         return u.i;
1000 }
1001
1002 __device_inline float __int_as_float(int i)
1003 {
1004         union { int i; float f; } u;
1005         u.i = i;
1006         return u.f;
1007 }
1008
1009 __device_inline uint __float_as_uint(float f)
1010 {
1011         union { uint i; float f; } u;
1012         u.f = f;
1013         return u.i;
1014 }
1015
1016 __device_inline float __uint_as_float(uint i)
1017 {
1018         union { uint i; float f; } u;
1019         u.i = i;
1020         return u.f;
1021 }
1022
1023 /* Interpolation */
1024
1025 template<class A, class B> A lerp(const A& a, const A& b, const B& t)
1026 {
1027         return (A)(a * ((B)1 - t) + b * t);
1028 }
1029
1030 /* Triangle */
1031
1032 __device_inline float triangle_area(const float3 v1, const float3 v2, const float3 v3)
1033 {
1034         return len(cross(v3 - v2, v1 - v2))*0.5f;
1035 }
1036
1037 #endif
1038
1039 /* Orthonormal vectors */
1040
1041 __device_inline void make_orthonormals(const float3 N, float3 *a, float3 *b)
1042 {
1043         if(N.x != N.y || N.x != N.z)
1044                 *a = make_float3(N.z-N.y, N.x-N.z, N.y-N.x);  //(1,1,1)x N
1045         else
1046                 *a = make_float3(N.z-N.y, N.x+N.z, -N.y-N.x);  //(-1,1,1)x N
1047
1048         *a = normalize(*a);
1049         *b = cross(N, *a);
1050 }
1051
1052 /* Color division */
1053
1054 __device_inline float3 safe_divide_color(float3 a, float3 b)
1055 {
1056         float x, y, z;
1057
1058         x = (b.x != 0.0f)? a.x/b.x: 0.0f;
1059         y = (b.y != 0.0f)? a.y/b.y: 0.0f;
1060         z = (b.z != 0.0f)? a.z/b.z: 0.0f;
1061
1062         return make_float3(x, y, z);
1063 }
1064
1065 CCL_NAMESPACE_END
1066
1067 #endif /* __UTIL_MATH_H__ */
1068