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-2025 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 ten_up_double (int_4 expo)
    35  {
    36  // This way appears sufficiently accurate.
    37    real_16 dbl_expo = 1.0q, *dep;
    38    logical_4 neg_expo;
    39    if (expo == 0) {
    40      return 1.0q;
    41    }
    42    neg_expo = (logical_4) (expo < 0);
    43    if (neg_expo) {
    44      expo = -expo;
    45    }
    46    if (expo > MAX_DOUBLE_EXPO) {
    47      expo = 0;
    48      errno = EDOM;
    49    }
    50  //  ABEND (expo > MAX_DOUBLE_EXPO, ERROR_INVALID_VALUE, __func__);
    51    for (dep = pow_10_double; expo != 0; expo >>= 1, dep++) {
    52      if (expo & 0x1) {
    53        dbl_expo *= *dep;
    54      }
    55    }
    56    return neg_expo ? 1.0q / dbl_expo : dbl_expo;
    57  }
    58  
    59  //! @brief Transform string into real-16.
    60  // From Algol 68 Genie.
    61  
    62  real_16 _strtoquad (char *s, char **end)
    63  {
    64    int_4 i, dot = -1, pos = 0, pow = 0, expo;
    65    real_16 sum, W, y[FLT128_DIG];
    66    errno = 0;
    67    for (i = 0; i < FLT128_DIG; i++) {
    68      y[i] = 0.0q;
    69    }
    70    while (isspace (s[0])) {
    71      s++;
    72    }
    73  // Scan mantissa digits and put them into "y".
    74    if (s[0] == '-') {
    75      W = -1.0q;
    76    } else {
    77      W = 1.0q;
    78    }
    79    if (s[0] == '+' || s[0] == '-') {
    80      s++;
    81    }
    82    while (s[0] == '0') {
    83      s++;
    84    }
    85    while (pow < FLT128_DIG && s[pos] != '\0' && (isdigit (s[pos]) || s[pos] == '.')) {
    86      if (s[pos] == '.') {
    87        dot = pos;
    88      } else {
    89        int_4 val = (int_4) s[pos] - (int_4) '0';
    90        y[pow] = W * val;
    91        W /= 10.0q;
    92        pow++;
    93      }
    94      pos++;
    95    }
    96    if (end != NO_REF_TEXT) {
    97      (*end) = &(s[pos]);
    98    }
    99  // Sum from low to high to preserve precision.
   100    sum = 0.0q;
   101    for (i = FLT128_DIG - 1; i >= 0; i--) {
   102      sum = sum + y[i];
   103    }
   104  // See if there is an exponent.
   105    if (_EXPCHAR (s[pos])) {
   106      expo = (int_4) strtol (&(s[++pos]), end, 10);
   107    } else {
   108      expo = 0;
   109    }
   110  // Standardise.
   111    if (dot >= 0) {
   112      expo += dot - 1;
   113    } else {
   114      expo += pow - 1;
   115    }
   116    while (sum != 0.0q && fabsq (sum) < 1.0q) {
   117      sum *= 10.0q;
   118      expo -= 1;
   119    }
   120  //
   121    if (errno == 0) {
   122      return sum * ten_up_double (expo);
   123    } else {
   124      return 0.0q;
   125    }
   126  }
   127
Deze website is gearchiveerd door de KB, nationale bibliotheek.

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