libgfortran/generated/bessel_r16.c

   1 /* Implementation of the BESSEL_JN and BESSEL_YN transformational
   2    function using a recurrence algorithm.
   3    Copyright 2010 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 #include <stdlib.h>
  29 #include <assert.h>
  30
  31
  32
  33 #if defined(GFC_REAL_16_IS_FLOAT128)
  34 #define MATHFUNC(funcname) funcname ## q
  35 #else
  36 #define MATHFUNC(funcname) funcname ## l
  37 #endif
  38
  39 #if defined (HAVE_GFC_REAL_16)
  40
  41
  42
  43 #if (defined(GFC_REAL_16_IS_FLOAT128) || defined(HAVE_JNL))
  44 extern void bessel_jn_r16 (gfc_array_r16 * const restrict ret, int n1,
  45                                      int n2, GFC_REAL_16 x);
  46 export_proto(bessel_jn_r16);
  47
  48 void
  49 bessel_jn_r16 (gfc_array_r16 * const restrict ret, int n1, int n2, GFC_REAL_16 x)
  50 {
  51   int i;
  52   index_type stride;
  53
  54   GFC_REAL_16 last1, last2, x2rev;
  55
  56   stride = GFC_DESCRIPTOR_STRIDE(ret,0);
  57
  58   if (ret->data == NULL)
  59     {
  60       size_t size = n2 < n1 ? 0 : n2-n1+1;
  61       GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1);
  62       ret->data = internal_malloc_size (sizeof (GFC_REAL_16) * size);
  63       ret->offset = 0;
  64     }
  65
  66   if (unlikely (n2 < n1))
  67     return;
  68
  69   if (unlikely (compile_options.bounds_check)
  70       && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1))
  71     runtime_error("Incorrect extent in return value of BESSEL_JN "
  72                   "(%ld vs. %ld)", (long int) n2-n1,
  73                   (long int) GFC_DESCRIPTOR_EXTENT(ret,0));
  74
  75   stride = GFC_DESCRIPTOR_STRIDE(ret,0);
  76
  77   if (unlikely (x == 0))
  78     {
  79       ret->data[0] = 1;
  80       for (i = 1; i <= n2-n1; i++)
  81         ret->data[i*stride] = 0;
  82       return;
  83     }
  84
  85   ret->data = ret->data;
  86   last1 = MATHFUNC(jn) (n2, x);
  87   ret->data[(n2-n1)*stride] = last1;
  88
  89   if (n1 == n2)
  90     return;
  91
  92   last2 = MATHFUNC(jn) (n2 - 1, x);
  93   ret->data[(n2-n1-1)*stride] = last2;
  94
  95   if (n1 + 1 == n2)
  96     return;
  97
  98   x2rev = GFC_REAL_16_LITERAL(2.)/x;
  99
 100   for (i = n2-n1-2; i >= 0; i--)
 101     {
 102       ret->data[i*stride] = x2rev * (i+1+n1) * last2 - last1;
 103       last1 = last2;
 104       last2 = ret->data[i*stride];
 105     }
 106 }
 107
 108 #endif
 109
 110 #if (defined(GFC_REAL_16_IS_FLOAT128) || defined(HAVE_YNL))
 111 extern void bessel_yn_r16 (gfc_array_r16 * const restrict ret,
 112                                      int n1, int n2, GFC_REAL_16 x);
 113 export_proto(bessel_yn_r16);
 114
 115 void
 116 bessel_yn_r16 (gfc_array_r16 * const restrict ret, int n1, int n2,
 117                          GFC_REAL_16 x)
 118 {
 119   int i;
 120   index_type stride;
 121
 122   GFC_REAL_16 last1, last2, x2rev;
 123
 124   stride = GFC_DESCRIPTOR_STRIDE(ret,0);
 125
 126   if (ret->data == NULL)
 127     {
 128       size_t size = n2 < n1 ? 0 : n2-n1+1;
 129       GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1);
 130       ret->data = internal_malloc_size (sizeof (GFC_REAL_16) * size);
 131       ret->offset = 0;
 132     }
 133
 134   if (unlikely (n2 < n1))
 135     return;
 136
 137   if (unlikely (compile_options.bounds_check)
 138       && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1))
 139     runtime_error("Incorrect extent in return value of BESSEL_JN "
 140                   "(%ld vs. %ld)", (long int) n2-n1,
 141                   (long int) GFC_DESCRIPTOR_EXTENT(ret,0));
 142
 143   stride = GFC_DESCRIPTOR_STRIDE(ret,0);
 144
 145   if (unlikely (x == 0))
 146     {
 147       for (i = 0; i <= n2-n1; i++)
 148 #if defined(GFC_REAL_16_INFINITY)
 149         ret->data[i*stride] = -GFC_REAL_16_INFINITY;
 150 #else
 151         ret->data[i*stride] = -GFC_REAL_16_HUGE;
 152 #endif
 153       return;
 154     }
 155
 156   ret->data = ret->data;
 157   last1 = MATHFUNC(yn) (n1, x);
 158   ret->data[0] = last1;
 159
 160   if (n1 == n2)
 161     return;
 162
 163   last2 = MATHFUNC(yn) (n1 + 1, x);
 164   ret->data[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 <= n1+n2; i++)
 172     {
 173 #if defined(GFC_REAL_16_INFINITY)
 174       if (unlikely (last2 == -GFC_REAL_16_INFINITY))
 175         {
 176           ret->data[i*stride] = -GFC_REAL_16_INFINITY;
 177         }
 178       else
 179 #endif
 180         {
 181           ret->data[i*stride] = x2rev * (i-1+n1) * last2 - last1;
 182           last1 = last2;
 183           last2 = ret->data[i*stride];
 184         }
 185     }
 186 }
 187 #endif
 188
 189 #endif
 190