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