libgfortran/generated/bessel_r16.c

   1 /* Implementation of the BESSEL_JN and BESSEL_YN transformational
   2    function using a recurrence algorithm.
   3    Copyright (C) 2010-2023 Free Software Foundation, Inc.
   4    Contributed by Tobias Burnus <burnus@net-b.de>
   5
   6 This file is part of the GNU Fortran runtime library (libgfortran).
   7
   8 Libgfortran is free software; you can redistribute it and/or
   9 modify it under the terms of the GNU General Public
  10 License as published by the Free Software Foundation; either
  11 version 3 of the License, or (at your option) any later version.
  12
  13 Libgfortran is distributed in the hope that it will be useful,
  14 but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16 GNU General Public License for more details.
  17
  18 Under Section 7 of GPL version 3, you are granted additional
  19 permissions described in the GCC Runtime Library Exception, version
  20 3.1, as published by the Free Software Foundation.
  21
  22 You should have received a copy of the GNU General Public License and
  23 a copy of the GCC Runtime Library Exception along with this program;
  24 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  25 <http://www.gnu.org/licenses/>.  */
  26
  27 #include "libgfortran.h"
  28
  29
  30
  31 #if defined(GFC_REAL_16_IS_FLOAT128)
  32 #if defined(GFC_REAL_16_USE_IEC_60559)
  33 #define MATHFUNC(funcname) funcname ## f128
  34 #else
  35 #define MATHFUNC(funcname) funcname ## q
  36 #endif
  37 #else
  38 #define MATHFUNC(funcname) funcname ## l
  39 #endif
  40
  41 #if defined (HAVE_GFC_REAL_16)
  42
  43
  44
  45 #if (defined(GFC_REAL_16_IS_FLOAT128) || defined(HAVE_JNL))
  46 extern void bessel_jn_r16 (gfc_array_r16 * const restrict ret, int n1,
  47                                      int n2, GFC_REAL_16 x);
  48 export_proto(bessel_jn_r16);
  49
  50 void
  51 bessel_jn_r16 (gfc_array_r16 * const restrict ret, int n1, int n2, GFC_REAL_16 x)
  52 {
  53   int i;
  54   index_type stride;
  55
  56   GFC_REAL_16 last1, last2, x2rev;
  57
  58   stride = GFC_DESCRIPTOR_STRIDE(ret,0);
  59
  60   if (ret->base_addr == NULL)
  61     {
  62       size_t size = n2 < n1 ? 0 : n2-n1+1;
  63       GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1);
  64       ret->base_addr = xmallocarray (size, sizeof (GFC_REAL_16));
  65       ret->offset = 0;
  66     }
  67
  68   if (unlikely (n2 < n1))
  69     return;
  70
  71   if (unlikely (compile_options.bounds_check)
  72       && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1))
  73     runtime_error("Incorrect extent in return value of BESSEL_JN "
  74                   "(%ld vs. %ld)", (long int) n2-n1,
  75                   (long int) GFC_DESCRIPTOR_EXTENT(ret,0));
  76
  77   stride = GFC_DESCRIPTOR_STRIDE(ret,0);
  78
  79   if (unlikely (x == 0))
  80     {
  81       ret->base_addr[0] = 1;
  82       for (i = 1; i <= n2-n1; i++)
  83         ret->base_addr[i*stride] = 0;
  84       return;
  85     }
  86
  87   last1 = MATHFUNC(jn) (n2, x);
  88   ret->base_addr[(n2-n1)*stride] = last1;
  89
  90   if (n1 == n2)
  91     return;
  92
  93   last2 = MATHFUNC(jn) (n2 - 1, x);
  94   ret->base_addr[(n2-n1-1)*stride] = last2;
  95
  96   if (n1 + 1 == n2)
  97     return;
  98
  99   x2rev = GFC_REAL_16_LITERAL(2.)/x;
 100
 101   for (i = n2-n1-2; i >= 0; i--)
 102     {
 103       ret->base_addr[i*stride] = x2rev * (i+1+n1) * last2 - last1;
 104       last1 = last2;
 105       last2 = ret->base_addr[i*stride];
 106     }
 107 }
 108
 109 #endif
 110
 111 #if (defined(GFC_REAL_16_IS_FLOAT128) || defined(HAVE_YNL))
 112 extern void bessel_yn_r16 (gfc_array_r16 * const restrict ret,
 113                                      int n1, int n2, GFC_REAL_16 x);
 114 export_proto(bessel_yn_r16);
 115
 116 void
 117 bessel_yn_r16 (gfc_array_r16 * const restrict ret, int n1, int n2,
 118                          GFC_REAL_16 x)
 119 {
 120   int i;
 121   index_type stride;
 122
 123   GFC_REAL_16 last1, last2, x2rev;
 124
 125   stride = GFC_DESCRIPTOR_STRIDE(ret,0);
 126
 127   if (ret->base_addr == NULL)
 128     {
 129       size_t size = n2 < n1 ? 0 : n2-n1+1;
 130       GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1);
 131       ret->base_addr = xmallocarray (size, sizeof (GFC_REAL_16));
 132       ret->offset = 0;
 133     }
 134
 135   if (unlikely (n2 < n1))
 136     return;
 137
 138   if (unlikely (compile_options.bounds_check)
 139       && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1))
 140     runtime_error("Incorrect extent in return value of BESSEL_JN "
 141                   "(%ld vs. %ld)", (long int) n2-n1,
 142                   (long int) GFC_DESCRIPTOR_EXTENT(ret,0));
 143
 144   stride = GFC_DESCRIPTOR_STRIDE(ret,0);
 145
 146   if (unlikely (x == 0))
 147     {
 148       for (i = 0; i <= n2-n1; i++)
 149 #if defined(GFC_REAL_16_INFINITY)
 150         ret->base_addr[i*stride] = -GFC_REAL_16_INFINITY;
 151 #else
 152         ret->base_addr[i*stride] = -GFC_REAL_16_HUGE;
 153 #endif
 154       return;
 155     }
 156
 157   last1 = MATHFUNC(yn) (n1, x);
 158   ret->base_addr[0] = last1;
 159
 160   if (n1 == n2)
 161     return;
 162
 163   last2 = MATHFUNC(yn) (n1 + 1, x);
 164   ret->base_addr[1*stride] = last2;
 165
 166   if (n1 + 1 == n2)
 167     return;
 168
 169   x2rev = GFC_REAL_16_LITERAL(2.)/x;
 170
 171   for (i = 2; i <= n2 - n1; i++)
 172     {
 173 #if defined(GFC_REAL_16_INFINITY)
 174       if (unlikely (last2 == -GFC_REAL_16_INFINITY))
 175         {
 176           ret->base_addr[i*stride] = -GFC_REAL_16_INFINITY;
 177         }
 178       else
 179 #endif
 180         {
 181           ret->base_addr[i*stride] = x2rev * (i-1+n1) * last2 - last1;
 182           last1 = last2;
 183           last2 = ret->base_addr[i*stride];
 184         }
 185     }
 186 }
 187 #endif
 188
 189 #endif
 190