star-hitran

shtr.c (5734B)

---

 /* Copyright (C) 2022, 2025, 2026 |Méso|Star> (contact@meso-star.com)
  * Copyright (C) 2025, 2026 Université de Lorraine
  * Copyright (C) 2022 Centre National de la Recherche Scientifique
  * Copyright (C) 2022 Université Paul Sabatier
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #include "shtr.h"
 #include "shtr_c.h"
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static INLINE res_T
 check_shtr_create_args(const struct shtr_create_args* args)
 {
   return args ? RES_OK : RES_BAD_ARG;
 }
 
 static void
 release_shtr(ref_T* ref)
 {
   struct shtr* shtr = CONTAINER_OF(ref, struct shtr, ref);
   ASSERT(ref);
   if(shtr->logger == &shtr->logger__) logger_release(&shtr->logger__);
   MEM_RM(shtr->allocator, shtr);
 }
 
 /*******************************************************************************
  * Local functions
  ******************************************************************************/
 int64_t
 float2fix
   (const double fp,
    const int ipl, /* Length (in bits) of the integer part */
    const int fpl) /* Length (in bits) of the fractional part */
 {
   union { double flt; int64_t i64; } fpart; /* Fractional part */
   int64_t ipart = 0; /* Integer part */
   double fp_abs = 0; /* Absolute value of fp */
   int sign = 0; /* 1 if fp is negative, and 0 otherwise */
 
   ASSERT(ipl + fpl <= 64 && fpl < 51); /* Pre-conditions */
 
   sign = fp < 0;
   fp_abs = fabs(fp);
   ipart = (int64_t)fp_abs & (BIT_I64(ipl)-1);
 
   /* The IEEE-754 encoding of a double precision floating point number `f'
    * whose binary representation is `F' is defined as:
    *    f = (-1)^S * 2^(E-1023)
    *      * (1 + For(i in [0..51]) { M & BIT(51-i) ? 2^i : 0 })
    * with S = F / 2^63; E = F / 2^52 and M = F & (2^52 - 1).
    *
    * By adding 1 to the fractional part of the floating point number submitted
    * as input, its exponent becomes zero (i.e., E = 1023) and its mantissa is
    * then the binary representation of the fractional part as expected by the
    * fixed format. In other words, the (fpl-1-i)^th bit of the mantissa then
    * corresponds to the power of 2 expected in fixed representation. */
   fpart.flt = 1.0 + (fp_abs - (double)ipart);
   fpart.i64 &= BIT_I64(52)-1;
 
   fpart.i64 >>= 52 - fpl; /* Quantification of the fractional part */
 
   return (((sign << ipl) | ipart) << fpl) | fpart.i64;
 }
 
 double
 fix2float
   (const int64_t fix,
    const int ipl, /* Length (in bits) of the integer part */
    const int fpl) /* Length (in bits) of the fractional part */
 {
   union { double flt; int64_t i64; } fpart; /* Fractional part */
   double ipart = 0; /* Integer part */
   double sign = 0; /* -1 if fix is negative, and 1 otherwise */
 
   ASSERT(ipl + fpl <= 64 && fpl < 51); /* Pre-conditions */
 
   sign = (fix & BIT_I64(ipl + fpl)) ? -1.0 : +1.0;
   ipart = (double)((fix >> fpl) & (BIT_I64(ipl)-1));
 
   /* The IEEE-754 encoding of a double precision floating point number `f'
    * whose binary representation is `F' is defined as:
    *    f = (-1)^S * 2^(E-1023)
    *      * (1 + For(i in [0..51]) { M & BIT(51-i) ? 2^i : 0 })
    * with S = F / 2^63; E = F / 2^52 and M = F & (2^52 - 1).
    *
    * Setting a floating-point number to 1 sets its exponent to zero (i.e., E =
    * 1023). Thus, the mantissa bits and the fractional part of a fixed-point
    * representation are identical. All you need to do is subtract 1 from the
    * floating-point representation to retrieve the floating-point encoding of
    * the fractional part */
   fpart.flt  = 1.0;
   fpart.i64 |= (fix & (BIT(fpl) - 1)) << (52 - fpl);
   fpart.flt -= 1.0;
 
   return sign * (ipart + fpart.flt);
 }
 
 /*******************************************************************************
  * Exported functions
  ******************************************************************************/
 res_T
 shtr_create
   (const struct shtr_create_args* args,
    struct shtr** out_shtr)
 {
   struct mem_allocator* allocator = NULL;
   struct shtr* shtr = NULL;
   res_T res = RES_OK;
 
   if(!out_shtr) { res = RES_BAD_ARG; goto error; }
   res = check_shtr_create_args(args);
   if(res != RES_OK) goto error;
 
   allocator = args->allocator ? args->allocator : &mem_default_allocator;
   shtr = MEM_CALLOC(allocator, 1, sizeof(*shtr));
   if(!shtr) {
     #define ERR_STR "Could not allocate the Star-HITRAN data structure.\n"
     if(args->logger) {
       logger_print(args->logger, LOG_ERROR, ERR_STR);
     } else {
       fprintf(stderr, MSG_ERROR_PREFIX ERR_STR);
     }
     #undef ERR_STR
     res = RES_MEM_ERR;
     goto error;
   }
   ref_init(&shtr->ref);
   shtr->allocator = allocator;
   shtr->verbose = args->verbose;
   shtr->logger = args->logger ? args->logger : LOGGER_DEFAULT;
 
 exit:
   if(out_shtr) *out_shtr = shtr;
   return res;
 error:
   if(shtr) { SHTR(ref_put(shtr)); shtr = NULL; }
   goto exit;
 }
 
 res_T
 shtr_ref_get(struct shtr* shtr)
 {
   if(!shtr) return RES_BAD_ARG;
   ref_get(&shtr->ref);
   return RES_OK;
 }
 
 res_T
 shtr_ref_put(struct shtr* shtr)
 {
   if(!shtr) return RES_BAD_ARG;
   ref_put(&shtr->ref, release_shtr);
   return RES_OK;
 }