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#if defined(POWER_IEEE128) 32#define MATHFUNC(funcname) __ ## funcname ## ieee128 33#else 34#define MATHFUNC(funcname) funcname ## q 35#endif 36 37#if defined (HAVE_GFC_REAL_17) 38 39 40 41#if 1 /* FIXME: figure this out later. */ 42extern void bessel_jn_r17 (gfc_array_r17 * const restrict ret, int n1, 43 int n2, GFC_REAL_17 x); 44export_proto(bessel_jn_r17); 45 46void 47bessel_jn_r17 (gfc_array_r17 * const restrict ret, int n1, int n2, GFC_REAL_17 x) 48{ 49 int i; 50 index_type stride; 51 52 GFC_REAL_17 last1, last2, x2rev; 53 54 stride = GFC_DESCRIPTOR_STRIDE(ret,0); 55 56 if (ret->base_addr == NULL) 57 { 58 size_t size = n2 < n1 ? 0 : n2-n1+1; 59 GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); 60 ret->base_addr = xmallocarray (size, sizeof (GFC_REAL_17)); 61 ret->offset = 0; 62 } 63 64 if (unlikely (n2 < n1)) 65 return; 66 67 if (unlikely (compile_options.bounds_check) 68 && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) 69 runtime_error("Incorrect extent in return value of BESSEL_JN " 70 "(%ld vs. %ld)", (long int) n2-n1, 71 (long int) GFC_DESCRIPTOR_EXTENT(ret,0)); 72 73 stride = GFC_DESCRIPTOR_STRIDE(ret,0); 74 75 if (unlikely (x == 0)) 76 { 77 ret->base_addr[0] = 1; 78 for (i = 1; i <= n2-n1; i++) 79 ret->base_addr[i*stride] = 0; 80 return; 81 } 82 83 last1 = MATHFUNC(jn) (n2, x); 84 ret->base_addr[(n2-n1)*stride] = last1; 85 86 if (n1 == n2) 87 return; 88 89 last2 = MATHFUNC(jn) (n2 - 1, x); 90 ret->base_addr[(n2-n1-1)*stride] = last2; 91 92 if (n1 + 1 == n2) 93 return; 94 95 x2rev = GFC_REAL_17_LITERAL(2.)/x; 96 97 for (i = n2-n1-2; i >= 0; i--) 98 { 99 ret->base_addr[i*stride] = x2rev * (i+1+n1) * last2 - last1; 100 last1 = last2; 101 last2 = ret->base_addr[i*stride]; 102 } 103} 104 105#endif 106 107#if 1 /* FIXME: figure this out later. */ 108extern void bessel_yn_r17 (gfc_array_r17 * const restrict ret, 109 int n1, int n2, GFC_REAL_17 x); 110export_proto(bessel_yn_r17); 111 112void 113bessel_yn_r17 (gfc_array_r17 * const restrict ret, int n1, int n2, 114 GFC_REAL_17 x) 115{ 116 int i; 117 index_type stride; 118 119 GFC_REAL_17 last1, last2, x2rev; 120 121 stride = GFC_DESCRIPTOR_STRIDE(ret,0); 122 123 if (ret->base_addr == NULL) 124 { 125 size_t size = n2 < n1 ? 0 : n2-n1+1; 126 GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); 127 ret->base_addr = xmallocarray (size, sizeof (GFC_REAL_17)); 128 ret->offset = 0; 129 } 130 131 if (unlikely (n2 < n1)) 132 return; 133 134 if (unlikely (compile_options.bounds_check) 135 && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) 136 runtime_error("Incorrect extent in return value of BESSEL_JN " 137 "(%ld vs. %ld)", (long int) n2-n1, 138 (long int) GFC_DESCRIPTOR_EXTENT(ret,0)); 139 140 stride = GFC_DESCRIPTOR_STRIDE(ret,0); 141 142 if (unlikely (x == 0)) 143 { 144 for (i = 0; i <= n2-n1; i++) 145#if defined(GFC_REAL_17_INFINITY) 146 ret->base_addr[i*stride] = -GFC_REAL_17_INFINITY; 147#else 148 ret->base_addr[i*stride] = -GFC_REAL_17_HUGE; 149#endif 150 return; 151 } 152 153 last1 = MATHFUNC(yn) (n1, x); 154 ret->base_addr[0] = last1; 155 156 if (n1 == n2) 157 return; 158 159 last2 = MATHFUNC(yn) (n1 + 1, x); 160 ret->base_addr[1*stride] = last2; 161 162 if (n1 + 1 == n2) 163 return; 164 165 x2rev = GFC_REAL_17_LITERAL(2.)/x; 166 167 for (i = 2; i <= n2 - n1; i++) 168 { 169#if defined(GFC_REAL_17_INFINITY) 170 if (unlikely (last2 == -GFC_REAL_17_INFINITY)) 171 { 172 ret->base_addr[i*stride] = -GFC_REAL_17_INFINITY; 173 } 174 else 175#endif 176 { 177 ret->base_addr[i*stride] = x2rev * (i-1+n1) * last2 - last1; 178 last1 = last2; 179 last2 = ret->base_addr[i*stride]; 180 } 181 } 182} 183#endif 184 185#endif 186 187