Mercurial > hg > CbC > CbC_gcc
comparison libgfortran/m4/bessel.m4 @ 111:04ced10e8804
gcc 7
author | kono |
---|---|
date | Fri, 27 Oct 2017 22:46:09 +0900 |
parents | |
children | 84e7813d76e9 |
comparison
equal
deleted
inserted
replaced
68:561a7518be6b | 111:04ced10e8804 |
---|---|
1 `/* Implementation of the BESSEL_JN and BESSEL_YN transformational | |
2 function using a recurrence algorithm. | |
3 Copyright (C) 2010-2017 Free Software Foundation, Inc. | |
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 | |
29 include(iparm.m4)dnl | |
30 include(`mtype.m4')dnl | |
31 | |
32 mathfunc_macro | |
33 | |
34 `#if defined (HAVE_'rtype_name`) | |
35 | |
36 | |
37 | |
38 #if 'hasmathfunc(jn)` | |
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 | |
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)` | |
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 | |
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 |