Algol 68 Genie

rts-real16.c

     
   1  //! @file rts-real16.c
   2  //! @author J. Marcel van der Veer
   3  //
   4  //! @section Copyright
   5  //
   6  // This file is part of VIF - vintage FORTRAN compiler.
   7  // Copyright 2020-2024 J. Marcel van der Veer <algol68g@xs4all.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 support for REAL*16 and COMPLEX*32.
  25  
  26  #include <vif.h>
  27  
  28  #define MAX_DOUBLE_EXPO 4932
  29  
  30  static real_16 pow_10_double[] = {
  31    10.0q, 100.0q, 1.0e4q, 1.0e8q, 1.0e16q, 1.0e32q, 1.0e64q, 1.0e128q, 1.0e256q, 1.0e512q, 1.0e1024q, 1.0e2048q, 1.0e4096q
  32  };
  33  
  34  real_16 _qnint (real_16 x)
  35  {
  36    return (x >= 0 ? floorq (0.5q + x) : - floorq (0.5q - x));
  37  }
  38  
  39  real_16 ten_up_double (int_4 expo)
  40  {
  41  // This way appears sufficiently accurate.
  42    real_16 dbl_expo = 1.0q, *dep;
  43    logical_4 neg_expo;
  44    if (expo == 0) {
  45      return 1.0q;
  46    }
  47    neg_expo = (logical_4) (expo < 0);
  48    if (neg_expo) {
  49      expo = -expo;
  50    }
  51    if (expo > MAX_DOUBLE_EXPO) {
  52      expo = 0;
  53      errno = EDOM;
  54    }
  55  //  ABEND (expo > MAX_DOUBLE_EXPO, ERROR_INVALID_VALUE, __func__);
  56    for (dep = pow_10_double; expo != 0; expo >>= 1, dep++) {
  57      if (expo & 0x1) {
  58        dbl_expo *= *dep;
  59      }
  60    }
  61    return neg_expo ? 1.0q / dbl_expo : dbl_expo;
  62  }
  63  
  64  //! @brief Transform string into real-16.
  65  // From Algol 68 Genie.
  66  
  67  real_16 _strtoquad (char *s, char **end)
  68  {
  69    int_4 i, dot = -1, pos = 0, pow = 0, expo;
  70    real_16 sum, W, y[FLT128_DIG];
  71    errno = 0;
  72    for (i = 0; i < FLT128_DIG; i++) {
  73      y[i] = 0.0q;
  74    }
  75    while (isspace (s[0])) {
  76      s++;
  77    }
  78  // Scan mantissa digits and put them into "y".
  79    if (s[0] == '-') {
  80      W = -1.0q;
  81    } else {
  82      W = 1.0q;
  83    }
  84    if (s[0] == '+' || s[0] == '-') {
  85      s++;
  86    }
  87    while (s[0] == '0') {
  88      s++;
  89    }
  90    while (pow < FLT128_DIG && s[pos] != '\0' && (isdigit (s[pos]) || s[pos] == '.')) {
  91      if (s[pos] == '.') {
  92        dot = pos;
  93      } else {
  94        int_4 val = (int_4) s[pos] - (int_4) '0';
  95        y[pow] = W * val;
  96        W /= 10.0q;
  97        pow++;
  98      }
  99      pos++;
 100    }
 101    if (end != NULL) {
 102      (*end) = &(s[pos]);
 103    }
 104  // Sum from low to high to preserve precision.
 105    sum = 0.0q;
 106    for (i = FLT128_DIG - 1; i >= 0; i--) {
 107      sum = sum + y[i];
 108    }
 109  // See if there is an exponent.
 110    if (_EXPCHAR (s[pos])) {
 111      expo = (int_4) strtol (&(s[++pos]), end, 10);
 112    } else {
 113      expo = 0;
 114    }
 115  // Standardise.
 116    if (dot >= 0) {
 117      expo += dot - 1;
 118    } else {
 119      expo += pow - 1;
 120    }
 121    while (sum != 0.0q && fabsq (sum) < 1.0q) {
 122      sum *= 10.0q;
 123      expo -= 1;
 124    }
 125  //
 126    if (errno == 0) {
 127      return sum * ten_up_double (expo);
 128    } else {
 129      return 0.0q;
 130    }
 131  }
 132
Deze website is gearchiveerd door de KB, nationale bibliotheek.

This website is archived by the National Library of the Netherlands.