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