]>
Commit | Line | Data |
---|---|---|
faa9fea4 | 1 | /* Implementation of the BESSEL_JN and BESSEL_YN transformational |
2 | function using a recurrence algorithm. | |
f1717362 | 3 | Copyright (C) 2010-2016 Free Software Foundation, Inc. |
faa9fea4 | 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 | ||
9dc59fdf | 32 | |
33 | #define MATHFUNC(funcname) funcname ## l | |
34 | ||
faa9fea4 | 35 | #if defined (HAVE_GFC_REAL_10) |
36 | ||
37 | ||
38 | ||
39 | #if defined (HAVE_JNL) | |
40 | extern void bessel_jn_r10 (gfc_array_r10 * const restrict ret, int n1, | |
41 | int n2, GFC_REAL_10 x); | |
42 | export_proto(bessel_jn_r10); | |
43 | ||
44 | void | |
45 | bessel_jn_r10 (gfc_array_r10 * const restrict ret, int n1, int n2, GFC_REAL_10 x) | |
46 | { | |
47 | int i; | |
48 | index_type stride; | |
49 | ||
50 | GFC_REAL_10 last1, last2, x2rev; | |
51 | ||
52 | stride = GFC_DESCRIPTOR_STRIDE(ret,0); | |
53 | ||
553877d9 | 54 | if (ret->base_addr == NULL) |
faa9fea4 | 55 | { |
56 | size_t size = n2 < n1 ? 0 : n2-n1+1; | |
57 | GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); | |
af1e9051 | 58 | ret->base_addr = xmallocarray (size, sizeof (GFC_REAL_10)); |
faa9fea4 | 59 | ret->offset = 0; |
60 | } | |
61 | ||
62 | if (unlikely (n2 < n1)) | |
63 | return; | |
64 | ||
65 | if (unlikely (compile_options.bounds_check) | |
66 | && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) | |
67 | runtime_error("Incorrect extent in return value of BESSEL_JN " | |
68 | "(%ld vs. %ld)", (long int) n2-n1, | |
96138327 | 69 | (long int) GFC_DESCRIPTOR_EXTENT(ret,0)); |
faa9fea4 | 70 | |
71 | stride = GFC_DESCRIPTOR_STRIDE(ret,0); | |
72 | ||
9dc59fdf | 73 | if (unlikely (x == 0)) |
faa9fea4 | 74 | { |
553877d9 | 75 | ret->base_addr[0] = 1; |
faa9fea4 | 76 | for (i = 1; i <= n2-n1; i++) |
553877d9 | 77 | ret->base_addr[i*stride] = 0; |
faa9fea4 | 78 | return; |
79 | } | |
80 | ||
9dc59fdf | 81 | last1 = MATHFUNC(jn) (n2, x); |
553877d9 | 82 | ret->base_addr[(n2-n1)*stride] = last1; |
faa9fea4 | 83 | |
84 | if (n1 == n2) | |
85 | return; | |
86 | ||
9dc59fdf | 87 | last2 = MATHFUNC(jn) (n2 - 1, x); |
553877d9 | 88 | ret->base_addr[(n2-n1-1)*stride] = last2; |
faa9fea4 | 89 | |
90 | if (n1 + 1 == n2) | |
91 | return; | |
92 | ||
9dc59fdf | 93 | x2rev = GFC_REAL_10_LITERAL(2.)/x; |
faa9fea4 | 94 | |
95 | for (i = n2-n1-2; i >= 0; i--) | |
96 | { | |
553877d9 | 97 | ret->base_addr[i*stride] = x2rev * (i+1+n1) * last2 - last1; |
faa9fea4 | 98 | last1 = last2; |
553877d9 | 99 | last2 = ret->base_addr[i*stride]; |
faa9fea4 | 100 | } |
101 | } | |
102 | ||
103 | #endif | |
104 | ||
105 | #if defined (HAVE_YNL) | |
106 | extern void bessel_yn_r10 (gfc_array_r10 * const restrict ret, | |
107 | int n1, int n2, GFC_REAL_10 x); | |
108 | export_proto(bessel_yn_r10); | |
109 | ||
110 | void | |
111 | bessel_yn_r10 (gfc_array_r10 * const restrict ret, int n1, int n2, | |
112 | GFC_REAL_10 x) | |
113 | { | |
114 | int i; | |
115 | index_type stride; | |
116 | ||
117 | GFC_REAL_10 last1, last2, x2rev; | |
118 | ||
119 | stride = GFC_DESCRIPTOR_STRIDE(ret,0); | |
120 | ||
553877d9 | 121 | if (ret->base_addr == NULL) |
faa9fea4 | 122 | { |
123 | size_t size = n2 < n1 ? 0 : n2-n1+1; | |
124 | GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); | |
af1e9051 | 125 | ret->base_addr = xmallocarray (size, sizeof (GFC_REAL_10)); |
faa9fea4 | 126 | ret->offset = 0; |
127 | } | |
128 | ||
129 | if (unlikely (n2 < n1)) | |
130 | return; | |
131 | ||
132 | if (unlikely (compile_options.bounds_check) | |
133 | && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) | |
134 | runtime_error("Incorrect extent in return value of BESSEL_JN " | |
135 | "(%ld vs. %ld)", (long int) n2-n1, | |
96138327 | 136 | (long int) GFC_DESCRIPTOR_EXTENT(ret,0)); |
faa9fea4 | 137 | |
138 | stride = GFC_DESCRIPTOR_STRIDE(ret,0); | |
139 | ||
9dc59fdf | 140 | if (unlikely (x == 0)) |
faa9fea4 | 141 | { |
142 | for (i = 0; i <= n2-n1; i++) | |
143 | #if defined(GFC_REAL_10_INFINITY) | |
553877d9 | 144 | ret->base_addr[i*stride] = -GFC_REAL_10_INFINITY; |
faa9fea4 | 145 | #else |
553877d9 | 146 | ret->base_addr[i*stride] = -GFC_REAL_10_HUGE; |
faa9fea4 | 147 | #endif |
148 | return; | |
149 | } | |
150 | ||
9dc59fdf | 151 | last1 = MATHFUNC(yn) (n1, x); |
553877d9 | 152 | ret->base_addr[0] = last1; |
faa9fea4 | 153 | |
154 | if (n1 == n2) | |
155 | return; | |
156 | ||
9dc59fdf | 157 | last2 = MATHFUNC(yn) (n1 + 1, x); |
553877d9 | 158 | ret->base_addr[1*stride] = last2; |
faa9fea4 | 159 | |
160 | if (n1 + 1 == n2) | |
161 | return; | |
162 | ||
9dc59fdf | 163 | x2rev = GFC_REAL_10_LITERAL(2.)/x; |
faa9fea4 | 164 | |
9872e76f | 165 | for (i = 2; i <= n2 - n1; i++) |
faa9fea4 | 166 | { |
167 | #if defined(GFC_REAL_10_INFINITY) | |
168 | if (unlikely (last2 == -GFC_REAL_10_INFINITY)) | |
169 | { | |
553877d9 | 170 | ret->base_addr[i*stride] = -GFC_REAL_10_INFINITY; |
faa9fea4 | 171 | } |
172 | else | |
173 | #endif | |
174 | { | |
553877d9 | 175 | ret->base_addr[i*stride] = x2rev * (i-1+n1) * last2 - last1; |
faa9fea4 | 176 | last1 = last2; |
553877d9 | 177 | last2 = ret->base_addr[i*stride]; |
faa9fea4 | 178 | } |
179 | } | |
180 | } | |
181 | #endif | |
182 | ||
183 | #endif | |
184 |