]>
Commit | Line | Data |
---|---|---|
47b99694 TB |
1 | `/* Implementation of the BESSEL_JN and BESSEL_YN transformational |
2 | function using a recurrence algorithm. | |
8d9254fc | 3 | Copyright (C) 2010-2020 Free Software Foundation, Inc. |
47b99694 TB |
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 | ||
887d9b8b | 27 | #include "libgfortran.h"' |
47b99694 TB |
28 | |
29 | include(iparm.m4)dnl | |
30 | include(`mtype.m4')dnl | |
31 | ||
08fd13d4 FXC |
32 | mathfunc_macro |
33 | ||
47b99694 TB |
34 | `#if defined (HAVE_'rtype_name`) |
35 | ||
36 | ||
37 | ||
08fd13d4 | 38 | #if 'hasmathfunc(jn)` |
47b99694 TB |
39 | extern void bessel_jn_r'rtype_kind` ('rtype` * const restrict ret, int n1, |
40 | int n2, 'rtype_name` x); | |
41 | export_proto(bessel_jn_r'rtype_kind`); | |
42 | ||
43 | void | |
44 | bessel_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 | ||
21d1335b | 53 | if (ret->base_addr == NULL) |
47b99694 TB |
54 | { |
55 | size_t size = n2 < n1 ? 0 : n2-n1+1; | |
56 | GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); | |
92e6f3a4 | 57 | ret->base_addr = xmallocarray (size, sizeof ('rtype_name`)); |
47b99694 TB |
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, | |
28cc1e9a | 68 | (long int) GFC_DESCRIPTOR_EXTENT(ret,0)); |
47b99694 TB |
69 | |
70 | stride = GFC_DESCRIPTOR_STRIDE(ret,0); | |
71 | ||
08fd13d4 | 72 | if (unlikely (x == 0)) |
47b99694 | 73 | { |
21d1335b | 74 | ret->base_addr[0] = 1; |
47b99694 | 75 | for (i = 1; i <= n2-n1; i++) |
21d1335b | 76 | ret->base_addr[i*stride] = 0; |
47b99694 TB |
77 | return; |
78 | } | |
79 | ||
08fd13d4 | 80 | last1 = MATHFUNC(jn) (n2, x); |
21d1335b | 81 | ret->base_addr[(n2-n1)*stride] = last1; |
47b99694 TB |
82 | |
83 | if (n1 == n2) | |
84 | return; | |
85 | ||
08fd13d4 | 86 | last2 = MATHFUNC(jn) (n2 - 1, x); |
21d1335b | 87 | ret->base_addr[(n2-n1-1)*stride] = last2; |
47b99694 TB |
88 | |
89 | if (n1 + 1 == n2) | |
90 | return; | |
91 | ||
08fd13d4 | 92 | x2rev = GFC_REAL_'rtype_kind`_LITERAL(2.)/x; |
47b99694 TB |
93 | |
94 | for (i = n2-n1-2; i >= 0; i--) | |
95 | { | |
21d1335b | 96 | ret->base_addr[i*stride] = x2rev * (i+1+n1) * last2 - last1; |
47b99694 | 97 | last1 = last2; |
21d1335b | 98 | last2 = ret->base_addr[i*stride]; |
47b99694 TB |
99 | } |
100 | } | |
101 | ||
102 | #endif | |
103 | ||
08fd13d4 | 104 | #if 'hasmathfunc(yn)` |
47b99694 TB |
105 | extern void bessel_yn_r'rtype_kind` ('rtype` * const restrict ret, |
106 | int n1, int n2, 'rtype_name` x); | |
107 | export_proto(bessel_yn_r'rtype_kind`); | |
108 | ||
109 | void | |
110 | bessel_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 | ||
21d1335b | 120 | if (ret->base_addr == NULL) |
47b99694 TB |
121 | { |
122 | size_t size = n2 < n1 ? 0 : n2-n1+1; | |
123 | GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); | |
92e6f3a4 | 124 | ret->base_addr = xmallocarray (size, sizeof ('rtype_name`)); |
47b99694 TB |
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, | |
28cc1e9a | 135 | (long int) GFC_DESCRIPTOR_EXTENT(ret,0)); |
47b99694 TB |
136 | |
137 | stride = GFC_DESCRIPTOR_STRIDE(ret,0); | |
138 | ||
08fd13d4 | 139 | if (unlikely (x == 0)) |
47b99694 TB |
140 | { |
141 | for (i = 0; i <= n2-n1; i++) | |
142 | #if defined('rtype_name`_INFINITY) | |
21d1335b | 143 | ret->base_addr[i*stride] = -'rtype_name`_INFINITY; |
47b99694 | 144 | #else |
21d1335b | 145 | ret->base_addr[i*stride] = -'rtype_name`_HUGE; |
47b99694 TB |
146 | #endif |
147 | return; | |
148 | } | |
149 | ||
08fd13d4 | 150 | last1 = MATHFUNC(yn) (n1, x); |
21d1335b | 151 | ret->base_addr[0] = last1; |
47b99694 TB |
152 | |
153 | if (n1 == n2) | |
154 | return; | |
155 | ||
08fd13d4 | 156 | last2 = MATHFUNC(yn) (n1 + 1, x); |
21d1335b | 157 | ret->base_addr[1*stride] = last2; |
47b99694 TB |
158 | |
159 | if (n1 + 1 == n2) | |
160 | return; | |
161 | ||
08fd13d4 | 162 | x2rev = GFC_REAL_'rtype_kind`_LITERAL(2.)/x; |
47b99694 | 163 | |
0d43a91d | 164 | for (i = 2; i <= n2 - n1; i++) |
47b99694 TB |
165 | { |
166 | #if defined('rtype_name`_INFINITY) | |
167 | if (unlikely (last2 == -'rtype_name`_INFINITY)) | |
168 | { | |
21d1335b | 169 | ret->base_addr[i*stride] = -'rtype_name`_INFINITY; |
47b99694 TB |
170 | } |
171 | else | |
172 | #endif | |
173 | { | |
21d1335b | 174 | ret->base_addr[i*stride] = x2rev * (i-1+n1) * last2 - last1; |
47b99694 | 175 | last1 = last2; |
21d1335b | 176 | last2 = ret->base_addr[i*stride]; |
47b99694 TB |
177 | } |
178 | } | |
179 | } | |
180 | #endif | |
181 | ||
182 | #endif' | |
183 |