[772] | 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 | } |
---|