1 | /* ------------------ */ |
---|
2 | /* --- nrarith0.c --- */ |
---|
3 | /* ------------------ */ |
---|
4 | |
---|
5 | /* |
---|
6 | * Copyright (c) 2000-2014, Lionel Lacassagne, All rights reserved |
---|
7 | * Univ Paris Sud XI, CNRS |
---|
8 | * |
---|
9 | * Distributed under the Boost Software License, Version 1.0 |
---|
10 | * see accompanying file LICENSE.txt or copy it at |
---|
11 | * http://www.boost.org/LICENSE_1_0.txt |
---|
12 | */ |
---|
13 | |
---|
14 | /* |
---|
15 | * History: |
---|
16 | * modif : log2 -> ilog2 (conflict with math.h on Mac OSX) |
---|
17 | */ |
---|
18 | |
---|
19 | #include <stdio.h> |
---|
20 | #include <stddef.h> |
---|
21 | #include <stdlib.h> |
---|
22 | #include <math.h> |
---|
23 | |
---|
24 | #include "mypredef.h" |
---|
25 | #include "nrtype.h" |
---|
26 | #include "nrdef.h" |
---|
27 | #include "nrmacro.h" |
---|
28 | #include "nrkernel.h" |
---|
29 | |
---|
30 | #include "nrarith0.h" |
---|
31 | |
---|
32 | ROUTINE(void) i8swap(int8 *a, int8 *b) { int8 *t; t=a; a=b; b=t;} |
---|
33 | ROUTINE(void) i16swap(int16 *a, int16 *b) { int16 *t; t=a; a=b; b=t;} |
---|
34 | ROUTINE(void) i32swap(int32 *a, int32 *b) { int32 *t; t=a; a=b; b=t;} |
---|
35 | ROUTINE(void) i64swap(int64 *a, int64 *b) { int64 *t; t=a; a=b; b=t;} |
---|
36 | ROUTINE(void) f32swap(float32 *a, float32 *b) { float32 *t; t=a; a=b; b=t;} |
---|
37 | ROUTINE(void) f64swap(float64 *a, float64 *b) { float64 *t; t=a; a=b; b=t;} |
---|
38 | ROUTINE(void) rgb8swap(rgb8 *a, rgb8 *b) { rgb8 *t; t=a; a=b; b=t;} |
---|
39 | ROUTINE(void) rgbx8swap(rgbx8 *a, rgbx8 *b) { rgbx8 *t; t=a; a=b; b=t;} |
---|
40 | |
---|
41 | /* --------- */ |
---|
42 | /* -- Min -- */ |
---|
43 | /* --------- */ |
---|
44 | |
---|
45 | ROUTINE(float32) f32min (float32 x1, float32 x2) {if (x1<x2) return x1; else return x2;} |
---|
46 | ROUTINE(float32) f32min2 (float32 x1, float32 x2) {if (x1<x2) return x1; else return x2;} |
---|
47 | ROUTINE(float32) f32min3 (float32 x1, float32 x2, float32 x3) {return f32min2(f32min2(x1, x2), x3);} |
---|
48 | ROUTINE(float32) f32min4 (float32 x1, float32 x2, float32 x3, float32 x4) {return f32min2(f32min2(x1, x2), f32min2(x3, x4));} |
---|
49 | ROUTINE(float32) f32min5 (float32 x1, float32 x2, float32 x3, float32 x4, float32 x5) {return f32min3(f32min2(x1, x2), f32min2(x3, x4), x5);} |
---|
50 | |
---|
51 | ROUTINE(float64) f64min (float64 x1, float64 x2) {if (x1<x2) return x1; else return x2;} |
---|
52 | ROUTINE(float64) f64min2 (float64 x1, float64 x2) {if (x1<x2) return x1; else return x2;} |
---|
53 | ROUTINE(float64) f64min3 (float64 x1, float64 x2, float64 x3) {return f64min2(f64min2(x1, x2), x3);} |
---|
54 | ROUTINE(float64) f64min4 (float64 x1, float64 x2, float64 x3, float64 x4) {return f64min2(f64min2(x1, x2), f64min2(x3, x4));} |
---|
55 | ROUTINE(float64) f64min5 (float64 x1, float64 x2, float64 x3, float64 x4, float64 x5) {return f64min3(f64min2(x1, x2), f64min2(x3, x4), x5);} |
---|
56 | |
---|
57 | ROUTINE(uint8) ui8min (uint8 x1, uint8 x2) {if (x1<x2) return x1; else return x2;} |
---|
58 | ROUTINE(uint8) ui8min2(uint8 x1, uint8 x2) {if (x1<x2) return x1; else return x2;} |
---|
59 | ROUTINE(uint8) ui8min3(uint8 x1, uint8 x2, uint8 x3) {return ui8min2(ui8min2(x1, x2), x3);} |
---|
60 | ROUTINE(uint8) ui8min4(uint8 x1, uint8 x2, uint8 x3, uint8 x4) {return ui8min2(ui8min2(x1, x2), ui8min2(x3, x4));} |
---|
61 | ROUTINE(uint8) ui8min5(uint8 x1, uint8 x2, uint8 x3, uint8 x4, uint8 x5) {return ui8min3(ui8min2(x1, x2), ui8min2(x3, x4), x5);} |
---|
62 | |
---|
63 | ROUTINE(uint16) ui16min (uint16 x1, uint16 x2) {if (x1<x2) return x1; else return x2;} |
---|
64 | ROUTINE(uint16) ui16min2(uint16 x1, uint16 x2) {if (x1<x2) return x1; else return x2;} |
---|
65 | ROUTINE(uint16) ui16min3(uint16 x1, uint16 x2, uint16 x3) {return ui16min2(ui16min2(x1, x2), x3);} |
---|
66 | ROUTINE(uint16) ui16min4(uint16 x1, uint16 x2, uint16 x3, uint16 x4) {return ui16min2(ui16min2(x1, x2), ui16min2(x3, x4));} |
---|
67 | ROUTINE(uint16) ui16min5(uint16 x1, uint16 x2, uint16 x3, uint16 x4, uint16 x5) {return ui16min3(ui16min2(x1, x2), ui16min2(x3, x4), x5);} |
---|
68 | |
---|
69 | ROUTINE(int32) ui32min (uint32 x1, uint32 x2) {if (x1<x2) return x1; else return x2;} |
---|
70 | ROUTINE(int32) ui32min2(uint32 x1, uint32 x2) {if (x1<x2) return x1; else return x2;} |
---|
71 | ROUTINE(int32) ui32min3(uint32 x1, uint32 x2, uint32 x3) {return ui32min2(ui32min2(x1, x2), x3);} |
---|
72 | ROUTINE(int32) ui32min4(uint32 x1, uint32 x2, uint32 x3, uint32 x4) {return ui32min2(ui32min2(x1, x2), ui32min2(x3, x4));} |
---|
73 | ROUTINE(int32) ui32min5(uint32 x1, uint32 x2, uint32 x3, uint32 x4, uint32 x5) {return ui32min3(ui32min2(x1, x2), ui32min2(x3, x4), x5);} |
---|
74 | |
---|
75 | ROUTINE(rgb8) rgb8min (rgb8 x1, rgb8 x2) {rgb8 y; y.r = ui8min2(x1.r,x2.r);y.g=ui8min2(x1.g,x2.g);y.b=ui8min2(x1.b,x2.b);return y;} |
---|
76 | ROUTINE(rgb8) rgb8min2(rgb8 x1, rgb8 x2) {rgb8 y; y.r = ui8min2(x1.r,x2.r);y.g=ui8min2(x1.g,x2.g);y.b=ui8min2(x1.b,x2.b);return y;} |
---|
77 | ROUTINE(rgb8) rgb8min3(rgb8 x1, rgb8 x2, rgb8 x3) {return rgb8min2(rgb8min2(x1, x2), x3);} |
---|
78 | ROUTINE(rgb8) rgb8min4(rgb8 x1, rgb8 x2, rgb8 x3, rgb8 x4) {return rgb8min2(rgb8min2(x1, x2), rgb8min2(x3,x4));} |
---|
79 | ROUTINE(rgb8) rgb8min5(rgb8 x1, rgb8 x2, rgb8 x3, rgb8 x4, rgb8 x5) {return rgb8min3(rgb8min2(x1, x2), rgb8min2(x3,x4), x5);} |
---|
80 | |
---|
81 | /* --------- */ |
---|
82 | /* -- Max -- */ |
---|
83 | /* --------- */ |
---|
84 | |
---|
85 | ROUTINE(float32) f32max (float32 x1, float32 x2) {if (x1>x2) return x1; else return x2;} |
---|
86 | ROUTINE(float32) f32max2(float32 x1, float32 x2) {if (x1>x2) return x1; else return x2;} |
---|
87 | ROUTINE(float32) f32max3(float32 x1, float32 x2, float32 x3) {return f32max2(f32max2(x1, x2), x3);} |
---|
88 | ROUTINE(float32) f32max4(float32 x1, float32 x2, float32 x3, float32 x4) {return f32max2(f32max2(x1, x2), f32max2(x3, x4));} |
---|
89 | ROUTINE(float32) f32max5(float32 x1, float32 x2, float32 x3, float32 x4, float32 x5) {return f32max3(f32max2(x1, x2), f32max2(x3, x4), x5);} |
---|
90 | |
---|
91 | ROUTINE(float64) f64max (float64 x1, float64 x2) {if (x1>x2) return x1; else return x2;} |
---|
92 | ROUTINE(float64) f64max2 (float64 x1, float64 x2) {if (x1>x2) return x1; else return x2;} |
---|
93 | ROUTINE(float64) f64max3 (float64 x1, float64 x2, float64 x3) {return f64max2(f64max2(x1, x2), x3);} |
---|
94 | ROUTINE(float64) f64max4 (float64 x1, float64 x2, float64 x3, float64 x4) {return f64max2(f64max2(x1, x2), f64max2(x3, x4));} |
---|
95 | ROUTINE(float64) f64max5 (float64 x1, float64 x2, float64 x3, float64 x4, float64 x5) {return f64max3(f64max2(x1, x2), f64max2(x3, x4), x5);} |
---|
96 | |
---|
97 | ROUTINE(uint8) ui8max (uint8 x1, uint8 x2) {if (x1>x2) return x1; else return x2;} |
---|
98 | ROUTINE(uint8) ui8max2(uint8 x1, uint8 x2) {if (x1>x2) return x1; else return x2;} |
---|
99 | ROUTINE(uint8) ui8max3(uint8 x1, uint8 x2, uint8 x3) {return ui8max2(ui8max2(x1, x2), x3);} |
---|
100 | ROUTINE(uint8) ui8max4(uint8 x1, uint8 x2, uint8 x3, uint8 x4) {return ui8max2(ui8max2(x1, x2), ui8max2(x3, x4));} |
---|
101 | ROUTINE(uint8) ui8max5(uint8 x1, uint8 x2, uint8 x3, uint8 x4, uint8 x5) {return ui8max3(ui8max2(x1, x2), ui8max2(x3, x4), x5);} |
---|
102 | |
---|
103 | ROUTINE(uint16) ui16max (uint16 x1, uint16 x2) {if (x1>x2) return x1; else return x2;} |
---|
104 | ROUTINE(uint16) ui16max2(uint16 x1, uint16 x2) {if (x1>x2) return x1; else return x2;} |
---|
105 | ROUTINE(uint16) ui16max3(uint16 x1, uint16 x2, uint16 x3) {return ui16max2(ui16max2(x1, x2), x3);} |
---|
106 | ROUTINE(uint16) ui16max4(uint16 x1, uint16 x2, uint16 x3, uint16 x4) {return ui16max2(ui16max2(x1, x2), ui16max2(x3, x4));} |
---|
107 | ROUTINE(uint16) ui16max5(uint16 x1, uint16 x2, uint16 x3, uint16 x4, uint16 x5) {return ui16max3(ui16max2(x1, x2), ui16max2(x3, x4), x5);} |
---|
108 | |
---|
109 | ROUTINE(int32) ui32max (uint32 x1, uint32 x2) {if (x1>x2) return x1; else return x2;} |
---|
110 | ROUTINE(int32) ui32max2(uint32 x1, uint32 x2) {if (x1>x2) return x1; else return x2;} |
---|
111 | ROUTINE(int32) ui32max3(uint32 x1, uint32 x2, uint32 x3) {return ui32max2(ui32max2(x1, x2), x3);} |
---|
112 | ROUTINE(int32) ui32max4(uint32 x1, uint32 x2, uint32 x3, uint32 x4) {return ui32max2(ui32max2(x1, x2), ui32max2(x3, x4));} |
---|
113 | ROUTINE(int32) ui32max5(uint32 x1, uint32 x2, uint32 x3, uint32 x4, int32 x5) {return ui32max3(ui32max2(x1, x2), ui32max2(x3, x4), x5);} |
---|
114 | |
---|
115 | ROUTINE(rgb8) rgb8max (rgb8 x1, rgb8 x2) {rgb8 y; y.r = ui8max2(x1.r,x2.r);y.g=ui8max2(x1.g,x2.g);y.b=ui8max2(x1.b,x2.b);return y;} |
---|
116 | ROUTINE(rgb8) rgb8max2(rgb8 x1, rgb8 x2) {rgb8 y; y.r = ui8max2(x1.r,x2.r);y.g=ui8max2(x1.g,x2.g);y.b=ui8max2(x1.b,x2.b);return y;} |
---|
117 | ROUTINE(rgb8) rgb8max3(rgb8 x1, rgb8 x2, rgb8 x3) {return rgb8max2(rgb8max2(x1, x2), x3);} |
---|
118 | ROUTINE(rgb8) rgb8max4(rgb8 x1, rgb8 x2, rgb8 x3, rgb8 x4) {return rgb8max2(rgb8max2(x1, x2), rgb8max2(x3,x4));} |
---|
119 | ROUTINE(rgb8) rgb8max5(rgb8 x1, rgb8 x2, rgb8 x3, rgb8 x4, rgb8 x5) {return rgb8max3(rgb8max2(x1, x2), rgb8max2(x3,x4), x5);} |
---|
120 | |
---|
121 | /* ----------- */ |
---|
122 | /* -- Other -- */ |
---|
123 | /* ----------- */ |
---|
124 | |
---|
125 | /* ------------------------------- */ |
---|
126 | ROUTINE(int32) i32bit(int32 x, int n) |
---|
127 | /* ------------------------------- */ |
---|
128 | { |
---|
129 | return ((x>>n)&1); |
---|
130 | } |
---|
131 | /* --------------------------- */ |
---|
132 | ROUTINE(int32) sym_int32(int32 x) |
---|
133 | /* --------------------------- */ |
---|
134 | { |
---|
135 | int i; |
---|
136 | int32 y = 0; |
---|
137 | for(i=0; i<31; i++) { |
---|
138 | y = y | (x & 1); |
---|
139 | x = x >> 1; |
---|
140 | } |
---|
141 | y = y | x; |
---|
142 | return y; |
---|
143 | } |
---|
144 | |
---|
145 | |
---|
146 | /* ----------------------- */ |
---|
147 | ROUTINE(int) ilog2(int x) |
---|
148 | /* ----------------------- */ |
---|
149 | { |
---|
150 | int s = 0; |
---|
151 | while(x) { |
---|
152 | x >>= 1; |
---|
153 | s++; |
---|
154 | } |
---|
155 | return s - 1; |
---|
156 | } |
---|
157 | /* ----------------------------- */ |
---|
158 | ROUTINE(int) next_power2(int x) |
---|
159 | /* ----------------------------- */ |
---|
160 | { |
---|
161 | int s = ilog2(x); |
---|
162 | int n = 1 << s; |
---|
163 | |
---|
164 | if(x != n) |
---|
165 | return n << 1; |
---|
166 | else |
---|
167 | return n; |
---|
168 | } |
---|
169 | /* ---------------------------- */ |
---|
170 | ROUTINE(int) gcd(int u, int v) |
---|
171 | /* ---------------------------- */ |
---|
172 | { |
---|
173 | int r; |
---|
174 | while(v) { |
---|
175 | r = u % v; |
---|
176 | u = v; |
---|
177 | v = r; |
---|
178 | } |
---|
179 | return u; |
---|
180 | } |
---|
181 | /* ---------------------------- */ |
---|
182 | ROUTINE(int) lcm(int u, int v) |
---|
183 | /* ---------------------------- */ |
---|
184 | { |
---|
185 | return (u*v)/gcd(u,v); |
---|
186 | } |
---|