Algol 68 Genie

rts-math.c

     
   1  //! @file rts-math.c
   2  //! @author J. Marcel van der Veer
   3  
   4  //! @section Copyright
   5  //!
   6  //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
   7  //! Copyright 2001-2026 J. Marcel van der Veer [algol68g@algol68genie.nl].
   8  
   9  //! @section License
  10  //!
  11  //! This program is free software; you can redistribute it and/or modify it 
  12  //! under the terms of the GNU General Public License as published by the 
  13  //! Free Software Foundation; either version 3 of the License, or 
  14  //! (at your option) any later version.
  15  //!
  16  //! This program is distributed in the hope that it will be useful, but 
  17  //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
  18  //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 
  19  //! more details. You should have received a copy of the GNU General Public 
  20  //! License along with this program. If not, see [http://www.gnu.org/licenses/].
  21  
  22  //! @section Synopsis
  23  //!
  24  //! Runtime routines.
  25  
  26  #include "a68g.h"
  27  #include "a68g-prelude.h"
  28  #include "a68g-genie.h"
  29  #include "a68g-numbers.h"
  30  #include "a68g-double.h"
  31  
  32  #define OVF_MASK (0x7fff000000000000ULL)
  33  
  34  DOUBLE_T nudge_double (DOUBLE_T u, int n)
  35  {
  36  // Return the nth next, or previous, DOUBLE_T number representation.
  37  // n > 0: move away from 0; n < 0 means move towards zero.
  38  //
  39  #if (A68G_LEVEL <= 2)
  40    (void) n;
  41    return u; // TODO
  42  #else 
  43    if (!a68g_finite_double (u)) {
  44      errno = EDOM;
  45      return u;
  46    } else if (n == 0) {
  47      return u;
  48    } else {
  49      DOUBLE_NUM_T w;
  50      w.f = ABS (u);
  51      for (int k = 0; k < ABS (n); k++) {
  52        if (n > 0) { // x += ulp
  53          LW (w)++;
  54          if (LW (w) == 0) {
  55            HW (w)++;
  56          }
  57        } else { // x -= ulp
  58          if (LW (w) == 0) {
  59             HW (w)--;
  60          }
  61          LW (w)--;
  62        }
  63        if ((HW (w) & OVF_MASK) == OVF_MASK) {
  64          // Overflow
  65          return u;
  66        } else if ((HW (w) & OVF_MASK) == 0) {
  67          // Underflow
  68          return u;
  69        } 
  70      }
  71      return (u >= 0.0q ? w.f : -w.f);
  72    }
  73  #endif
  74  }