[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master d5ab8c9 2/2: Arithmetic: sqrt, log and log10 w
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master d5ab8c9 2/2: Arithmetic: sqrt, log and log10 will always produce floating point |
|
Date: |
Thu, 28 Feb 2019 12:52:53 -0500 (EST) |
branch: master
commit d5ab8c9635f800c7c88e44af1f27fdd05d34ffea
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Arithmetic: sqrt, log and log10 will always produce floating point
Until now, these three operators would save their results in the same type
as the input. But this is not a natural way to use these operators: they
don't map linearly and commonly loosing floating-point precision makes it
hard to deal with them.
With this commit, the Arithmetic libray's internal
`arithmetic_unary_function' will only do in place operations if the input
already has floating point types. If it doesn't, the output will be 32-bit
floating point.
---
doc/gnuastro.texi | 30 +++++++++---------
lib/arithmetic.c | 91 ++++++++++++++++++++++++++++++++++++++++++++-----------
2 files changed, 89 insertions(+), 32 deletions(-)
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 4c90c2a..f9d36b3 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12164,18 +12164,17 @@ non-zero decimals) put an @code{f} after it.
@item sqrt
The square root of the first operand, so address@hidden sqrt}'' is equivalent
-to @mymath{\sqrt{5}}. The output type is determined from the input, so the
-output of this example will be @command{2} (since @command{5} doesn't have
-any non-zero decimal digits). If you want @command{2.23607}, run
address@hidden sqrt} instead, the @command{f} will ensure that a number will
-be read as a floating point number, even if it doesn't have decimal
-digits. If the input image has an integer type, you should explicitly
-convert the image to floating point, for example @command{a.fits float
-sqrt}, see the type conversion operators below.
+to @mymath{\sqrt{5}}. The output will have a floating point type, but its
+precision is determined from the input: if the input is a 64-bit floating
+point, the output will also be 64-bit. Otherwise, the output will be 32-bit
+floating point (see @ref{Numeric data types} for the respective
+precision). Therefore if you require 64-bit precision in estimating the
+square root, convert the input to 64-bit floating point first, for example
+with @code{5 float64 sqrt}.
@item log
Natural logarithm of first operand, so address@hidden log}'' is equivalent to
address@hidden(4)}. The output type is determined from the input, see the
address@hidden(4)}. The output type is determined from the input, see the
explanation under @command{sqrt} for more.
@item log10
@@ -27005,11 +27004,14 @@ polish notation}).
@deffnx Macro GAL_ARITHMETIC_OP_LOG10
Unary operator functions for calculating the square root
(@mymath{\sqrt{i}}), @mymath{ln(i)} and @mymath{log(i)} mathematic
-operators on each element of the input dataset. The output will have the
-same type as the input, so if your inputs are integer types be careful.
-
-If you want your output to be floating point but your input is an integer
-type, you can convert the input to a floating point type with
+operators on each element of the input dataset. The returned dataset will
+have a floating point type, but its precision is determined from the input:
+if the input is a 64-bit floating point, the output will also be
+64-bit. Otherwise, the returned dataset will be 32-bit floating point. See
address@hidden data types} for the respective precision.
+
+If you want your output to be 64-bit floating point but your input is a
+different type, you can convert the input to a floating point type with
@code{gal_data_copy_to_new_type} or
@code{gal_data_copy_to_new_type_free}(see @ref{Copying datasets}).
@end deffn
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index d24e960..aba3df2 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -350,43 +350,82 @@ arithmetic_abs(int flags, gal_data_t *in)
-#define UNIFUNC_RUN_FUNCTION_ON_ELEMENT(IT, OP){ \
- IT *ia=in->array, *oa=o->array, *iaf=ia + in->size; \
+#define UNIFUNC_RUN_FUNCTION_ON_ELEMENT(OT, IT, OP){ \
+ OT *oa=o->array; \
+ IT *ia=in->array, *iaf=ia + in->size; \
do *oa++ = OP(*ia++); while(ia<iaf); \
}
+#define UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(IT, OP) \
+ switch(o->type) \
+ { \
+ case GAL_TYPE_UINT8: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint8_t, IT, OP) \
+ break; \
+ case GAL_TYPE_INT8: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int8_t, IT, OP) \
+ break; \
+ case GAL_TYPE_UINT16: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint16_t, IT, OP) \
+ break; \
+ case GAL_TYPE_INT16: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int16_t, IT, OP) \
+ break; \
+ case GAL_TYPE_UINT32: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint32_t, IT, OP) \
+ break; \
+ case GAL_TYPE_INT32: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int32_t, IT, OP) \
+ break; \
+ case GAL_TYPE_UINT64: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint64_t, IT, OP) \
+ break; \
+ case GAL_TYPE_INT64: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int64_t, IT, OP) \
+ break; \
+ case GAL_TYPE_FLOAT32: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(float, IT, OP) \
+ break; \
+ case GAL_TYPE_FLOAT64: \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT(double, IT, OP) \
+ break; \
+ default: \
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized", \
+ "UNIARY_FUNCTION_ON_ELEMENT", in->type); \
+ }
+
#define UNIARY_FUNCTION_ON_ELEMENT(OP) \
switch(in->type) \
{ \
case GAL_TYPE_UINT8: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint8_t, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(uint8_t, OP) \
break; \
case GAL_TYPE_INT8: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int8_t, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(int8_t, OP) \
break; \
case GAL_TYPE_UINT16: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint16_t, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(uint16_t, OP) \
break; \
case GAL_TYPE_INT16: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int16_t, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(int16_t, OP) \
break; \
case GAL_TYPE_UINT32: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint32_t, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(uint32_t, OP) \
break; \
case GAL_TYPE_INT32: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int32_t, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(int32_t, OP) \
break; \
case GAL_TYPE_UINT64: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint64_t, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(uint64_t, OP) \
break; \
case GAL_TYPE_INT64: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int64_t, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(int64_t, OP) \
break; \
case GAL_TYPE_FLOAT32: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(float, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(float, OP) \
break; \
case GAL_TYPE_FLOAT64: \
- UNIFUNC_RUN_FUNCTION_ON_ELEMENT(double, OP) \
+ UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(double, OP) \
break; \
default: \
error(EXIT_FAILURE, 0, "%s: type code %d not recognized", \
@@ -396,15 +435,31 @@ arithmetic_abs(int flags, gal_data_t *in)
static gal_data_t *
arithmetic_unary_function(int operator, int flags, gal_data_t *in)
{
+ uint8_t otype;
+ int inplace=0;
gal_data_t *o;
- /* If we want inplace output, set the output pointer to the input
- pointer, for every pixel, the operation will be independent. */
- if(flags & GAL_ARITHMETIC_INPLACE)
- o = in;
+ /* See if the operation should be done in place. Note that so far, the
+ output of these operators is defined in the real space (floating
+ point). So even if the user requested inplace opereation, if its not a
+ floating point type, its not useful.*/
+ if( (flags & GAL_ARITHMETIC_INPLACE)
+ && (in->type==GAL_TYPE_FLOAT32 || in->type==GAL_TYPE_FLOAT64) )
+ inplace=1;
+
+ if(inplace)
+ {
+ o = in;
+ otype=in->type;
+ }
else
- o = gal_data_alloc(NULL, in->type, in->ndim, in->dsize, in->wcs,
- 0, in->minmapsize, NULL, NULL, NULL);
+ {
+ otype = ( in->type==GAL_TYPE_FLOAT64
+ ? GAL_TYPE_FLOAT64
+ : GAL_TYPE_FLOAT32 );
+ o = gal_data_alloc(NULL, otype, in->ndim, in->dsize, in->wcs,
+ 0, in->minmapsize, NULL, NULL, NULL);
+ }
/* Start setting the operator and operands. */
switch(operator)