From 47b996944dcb50a831c1332b8ea667ff6f95fa95 Mon Sep 17 00:00:00 2001 From: Tobias Burnus Date: Sat, 21 Aug 2010 12:12:53 +0200 Subject: re PR fortran/36158 (Transformational function BESSEL_YN(n1,n2,x) and BESSEL_JN missing) 2010-08-21 Tobias Burnus PR fortran/36158 PR fortran/33197 * intrinsic.c (add_sym): Init value attribute. (set_attr_value): New function. (add_functions) Use it and add JN/YN resolvers. * symbol.c (gfc_copy_formal_args_intr): Copy value attr. * intrinsic.h (gfc_resolve_bessel_n2): New prototype. * gfortran.h (gfc_intrinsic_arg): Add value attribute. * iresolve.c (gfc_resolve_bessel_n2): New function. * trans-intrinsic.c (gfc_get_symbol_for_expr): Create formal arg list. (gfc_conv_intrinsic_function,gfc_is_intrinsic_libcall): Add GFC_ISYM_JN2/GFC_ISYM_YN2 as case value. * simplify.c (): For YN set to -INF if previous values was -INF. * trans-expr.c (gfc_conv_procedure_call): Don't crash if sym->as is NULL. * iresolve.c (gfc_resolve_extends_type_of): Set the type of the dummy argument to the one of the actual. 2010-08-21 Tobias Burnus PR fortran/36158 PR fortran/33197 * m4/bessel.m4: Implement bessel_jn and bessel_yn. * gfortran.map: Add the generated bessel_jn_r{4,8,10,16} and bessel_yn_r{4,8,10,16}. * Makefile.am: Add bessel.m4. * Makefile.in: Regenerated. * generated/bessel_r4.c: Generated. * generated/bessel_r16.c: Generated. * generated/bessel_r8.c: Generated. * generated/bessel_r10.c: Generated. 2010-08-21 Tobias Burnus PR fortran/36158 PR fortran/33197 * gfortran.dg/bessel_6.f90: New. * gfortran.dg/bessel_7.f90: New. From-SVN: r163440 --- libgfortran/generated/bessel_r10.c | 183 +++++++++++++++++++++++++++++++++++++ libgfortran/generated/bessel_r16.c | 183 +++++++++++++++++++++++++++++++++++++ libgfortran/generated/bessel_r4.c | 183 +++++++++++++++++++++++++++++++++++++ libgfortran/generated/bessel_r8.c | 183 +++++++++++++++++++++++++++++++++++++ 4 files changed, 732 insertions(+) create mode 100644 libgfortran/generated/bessel_r10.c create mode 100644 libgfortran/generated/bessel_r16.c create mode 100644 libgfortran/generated/bessel_r4.c create mode 100644 libgfortran/generated/bessel_r8.c (limited to 'libgfortran/generated') diff --git a/libgfortran/generated/bessel_r10.c b/libgfortran/generated/bessel_r10.c new file mode 100644 index 0000000..579de84 --- /dev/null +++ b/libgfortran/generated/bessel_r10.c @@ -0,0 +1,183 @@ +/* Implementation of the BESSEL_JN and BESSEL_YN transformational + function using a recurrence algorithm. + Copyright 2010 Free Software Foundation, Inc. + Contributed by Tobias Burnus + +This file is part of the GNU Fortran runtime library (libgfortran). + +Libgfortran is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public +License as published by the Free Software Foundation; either +version 3 of the License, or (at your option) any later version. + +Libgfortran is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +Under Section 7 of GPL version 3, you are granted additional +permissions described in the GCC Runtime Library Exception, version +3.1, as published by the Free Software Foundation. + +You should have received a copy of the GNU General Public License and +a copy of the GCC Runtime Library Exception along with this program; +see the files COPYING3 and COPYING.RUNTIME respectively. If not, see +. */ + +#include "libgfortran.h" +#include +#include + + +#if defined (HAVE_GFC_REAL_10) + + + +#if defined (HAVE_JNL) +extern void bessel_jn_r10 (gfc_array_r10 * const restrict ret, int n1, + int n2, GFC_REAL_10 x); +export_proto(bessel_jn_r10); + +void +bessel_jn_r10 (gfc_array_r10 * const restrict ret, int n1, int n2, GFC_REAL_10 x) +{ + int i; + index_type stride; + + GFC_REAL_10 last1, last2, x2rev; + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (ret->data == NULL) + { + size_t size = n2 < n1 ? 0 : n2-n1+1; + GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); + ret->data = internal_malloc_size (sizeof (GFC_REAL_10) * size); + ret->offset = 0; + } + + if (unlikely (n2 < n1)) + return; + + if (unlikely (compile_options.bounds_check) + && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) + runtime_error("Incorrect extent in return value of BESSEL_JN " + "(%ld vs. %ld)", (long int) n2-n1, + GFC_DESCRIPTOR_EXTENT(ret,0)); + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (unlikely (x == 0.0L)) + { + ret->data[0] = 1.0L; + for (i = 1; i <= n2-n1; i++) + ret->data[i*stride] = 0.0L; + return; + } + + ret->data = ret->data; + last1 = jnl (n2, x); + ret->data[(n2-n1)*stride] = last1; + + if (n1 == n2) + return; + + last2 = jnl (n2 - 1, x); + ret->data[(n2-n1-1)*stride] = last2; + + if (n1 + 1 == n2) + return; + + x2rev = 2.0L/x; + + for (i = n2-n1-2; i >= 0; i--) + { + ret->data[i*stride] = x2rev * (i+1+n1) * last2 - last1; + last1 = last2; + last2 = ret->data[i*stride]; + } +} + +#endif + +#if defined (HAVE_YNL) +extern void bessel_yn_r10 (gfc_array_r10 * const restrict ret, + int n1, int n2, GFC_REAL_10 x); +export_proto(bessel_yn_r10); + +void +bessel_yn_r10 (gfc_array_r10 * const restrict ret, int n1, int n2, + GFC_REAL_10 x) +{ + int i; + index_type stride; + + GFC_REAL_10 last1, last2, x2rev; + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (ret->data == NULL) + { + size_t size = n2 < n1 ? 0 : n2-n1+1; + GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); + ret->data = internal_malloc_size (sizeof (GFC_REAL_10) * size); + ret->offset = 0; + } + + if (unlikely (n2 < n1)) + return; + + if (unlikely (compile_options.bounds_check) + && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) + runtime_error("Incorrect extent in return value of BESSEL_JN " + "(%ld vs. %ld)", (long int) n2-n1, + GFC_DESCRIPTOR_EXTENT(ret,0)); + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (unlikely (x == 0.0L)) + { + for (i = 0; i <= n2-n1; i++) +#if defined(GFC_REAL_10_INFINITY) + ret->data[i*stride] = -GFC_REAL_10_INFINITY; +#else + ret->data[i*stride] = -GFC_REAL_10_HUGE; +#endif + return; + } + + ret->data = ret->data; + last1 = ynl (n1, x); + ret->data[0] = last1; + + if (n1 == n2) + return; + + last2 = ynl (n1 + 1, x); + ret->data[1*stride] = last2; + + if (n1 + 1 == n2) + return; + + x2rev = 2.0L/x; + + for (i = 2; i <= n1+n2; i++) + { +#if defined(GFC_REAL_10_INFINITY) + if (unlikely (last2 == -GFC_REAL_10_INFINITY)) + { + ret->data[i*stride] = -GFC_REAL_10_INFINITY; + } + else +#endif + { + ret->data[i*stride] = x2rev * (i-1+n1) * last2 - last1; + last1 = last2; + last2 = ret->data[i*stride]; + } + } +} +#endif + +#endif + diff --git a/libgfortran/generated/bessel_r16.c b/libgfortran/generated/bessel_r16.c new file mode 100644 index 0000000..a3ef181 --- /dev/null +++ b/libgfortran/generated/bessel_r16.c @@ -0,0 +1,183 @@ +/* Implementation of the BESSEL_JN and BESSEL_YN transformational + function using a recurrence algorithm. + Copyright 2010 Free Software Foundation, Inc. + Contributed by Tobias Burnus + +This file is part of the GNU Fortran runtime library (libgfortran). + +Libgfortran is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public +License as published by the Free Software Foundation; either +version 3 of the License, or (at your option) any later version. + +Libgfortran is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +Under Section 7 of GPL version 3, you are granted additional +permissions described in the GCC Runtime Library Exception, version +3.1, as published by the Free Software Foundation. + +You should have received a copy of the GNU General Public License and +a copy of the GCC Runtime Library Exception along with this program; +see the files COPYING3 and COPYING.RUNTIME respectively. If not, see +. */ + +#include "libgfortran.h" +#include +#include + + +#if defined (HAVE_GFC_REAL_16) + + + +#if defined (HAVE_JNL) +extern void bessel_jn_r16 (gfc_array_r16 * const restrict ret, int n1, + int n2, GFC_REAL_16 x); +export_proto(bessel_jn_r16); + +void +bessel_jn_r16 (gfc_array_r16 * const restrict ret, int n1, int n2, GFC_REAL_16 x) +{ + int i; + index_type stride; + + GFC_REAL_16 last1, last2, x2rev; + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (ret->data == NULL) + { + size_t size = n2 < n1 ? 0 : n2-n1+1; + GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); + ret->data = internal_malloc_size (sizeof (GFC_REAL_16) * size); + ret->offset = 0; + } + + if (unlikely (n2 < n1)) + return; + + if (unlikely (compile_options.bounds_check) + && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) + runtime_error("Incorrect extent in return value of BESSEL_JN " + "(%ld vs. %ld)", (long int) n2-n1, + GFC_DESCRIPTOR_EXTENT(ret,0)); + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (unlikely (x == 0.0L)) + { + ret->data[0] = 1.0L; + for (i = 1; i <= n2-n1; i++) + ret->data[i*stride] = 0.0L; + return; + } + + ret->data = ret->data; + last1 = jnl (n2, x); + ret->data[(n2-n1)*stride] = last1; + + if (n1 == n2) + return; + + last2 = jnl (n2 - 1, x); + ret->data[(n2-n1-1)*stride] = last2; + + if (n1 + 1 == n2) + return; + + x2rev = 2.0L/x; + + for (i = n2-n1-2; i >= 0; i--) + { + ret->data[i*stride] = x2rev * (i+1+n1) * last2 - last1; + last1 = last2; + last2 = ret->data[i*stride]; + } +} + +#endif + +#if defined (HAVE_YNL) +extern void bessel_yn_r16 (gfc_array_r16 * const restrict ret, + int n1, int n2, GFC_REAL_16 x); +export_proto(bessel_yn_r16); + +void +bessel_yn_r16 (gfc_array_r16 * const restrict ret, int n1, int n2, + GFC_REAL_16 x) +{ + int i; + index_type stride; + + GFC_REAL_16 last1, last2, x2rev; + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (ret->data == NULL) + { + size_t size = n2 < n1 ? 0 : n2-n1+1; + GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); + ret->data = internal_malloc_size (sizeof (GFC_REAL_16) * size); + ret->offset = 0; + } + + if (unlikely (n2 < n1)) + return; + + if (unlikely (compile_options.bounds_check) + && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) + runtime_error("Incorrect extent in return value of BESSEL_JN " + "(%ld vs. %ld)", (long int) n2-n1, + GFC_DESCRIPTOR_EXTENT(ret,0)); + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (unlikely (x == 0.0L)) + { + for (i = 0; i <= n2-n1; i++) +#if defined(GFC_REAL_16_INFINITY) + ret->data[i*stride] = -GFC_REAL_16_INFINITY; +#else + ret->data[i*stride] = -GFC_REAL_16_HUGE; +#endif + return; + } + + ret->data = ret->data; + last1 = ynl (n1, x); + ret->data[0] = last1; + + if (n1 == n2) + return; + + last2 = ynl (n1 + 1, x); + ret->data[1*stride] = last2; + + if (n1 + 1 == n2) + return; + + x2rev = 2.0L/x; + + for (i = 2; i <= n1+n2; i++) + { +#if defined(GFC_REAL_16_INFINITY) + if (unlikely (last2 == -GFC_REAL_16_INFINITY)) + { + ret->data[i*stride] = -GFC_REAL_16_INFINITY; + } + else +#endif + { + ret->data[i*stride] = x2rev * (i-1+n1) * last2 - last1; + last1 = last2; + last2 = ret->data[i*stride]; + } + } +} +#endif + +#endif + diff --git a/libgfortran/generated/bessel_r4.c b/libgfortran/generated/bessel_r4.c new file mode 100644 index 0000000..961c67c --- /dev/null +++ b/libgfortran/generated/bessel_r4.c @@ -0,0 +1,183 @@ +/* Implementation of the BESSEL_JN and BESSEL_YN transformational + function using a recurrence algorithm. + Copyright 2010 Free Software Foundation, Inc. + Contributed by Tobias Burnus + +This file is part of the GNU Fortran runtime library (libgfortran). + +Libgfortran is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public +License as published by the Free Software Foundation; either +version 3 of the License, or (at your option) any later version. + +Libgfortran is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +Under Section 7 of GPL version 3, you are granted additional +permissions described in the GCC Runtime Library Exception, version +3.1, as published by the Free Software Foundation. + +You should have received a copy of the GNU General Public License and +a copy of the GCC Runtime Library Exception along with this program; +see the files COPYING3 and COPYING.RUNTIME respectively. If not, see +. */ + +#include "libgfortran.h" +#include +#include + + +#if defined (HAVE_GFC_REAL_4) + + + +#if defined (HAVE_JNF) +extern void bessel_jn_r4 (gfc_array_r4 * const restrict ret, int n1, + int n2, GFC_REAL_4 x); +export_proto(bessel_jn_r4); + +void +bessel_jn_r4 (gfc_array_r4 * const restrict ret, int n1, int n2, GFC_REAL_4 x) +{ + int i; + index_type stride; + + GFC_REAL_4 last1, last2, x2rev; + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (ret->data == NULL) + { + size_t size = n2 < n1 ? 0 : n2-n1+1; + GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); + ret->data = internal_malloc_size (sizeof (GFC_REAL_4) * size); + ret->offset = 0; + } + + if (unlikely (n2 < n1)) + return; + + if (unlikely (compile_options.bounds_check) + && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) + runtime_error("Incorrect extent in return value of BESSEL_JN " + "(%ld vs. %ld)", (long int) n2-n1, + GFC_DESCRIPTOR_EXTENT(ret,0)); + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (unlikely (x == 0.0F)) + { + ret->data[0] = 1.0F; + for (i = 1; i <= n2-n1; i++) + ret->data[i*stride] = 0.0F; + return; + } + + ret->data = ret->data; + last1 = jnf (n2, x); + ret->data[(n2-n1)*stride] = last1; + + if (n1 == n2) + return; + + last2 = jnf (n2 - 1, x); + ret->data[(n2-n1-1)*stride] = last2; + + if (n1 + 1 == n2) + return; + + x2rev = 2.0F/x; + + for (i = n2-n1-2; i >= 0; i--) + { + ret->data[i*stride] = x2rev * (i+1+n1) * last2 - last1; + last1 = last2; + last2 = ret->data[i*stride]; + } +} + +#endif + +#if defined (HAVE_YNF) +extern void bessel_yn_r4 (gfc_array_r4 * const restrict ret, + int n1, int n2, GFC_REAL_4 x); +export_proto(bessel_yn_r4); + +void +bessel_yn_r4 (gfc_array_r4 * const restrict ret, int n1, int n2, + GFC_REAL_4 x) +{ + int i; + index_type stride; + + GFC_REAL_4 last1, last2, x2rev; + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (ret->data == NULL) + { + size_t size = n2 < n1 ? 0 : n2-n1+1; + GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); + ret->data = internal_malloc_size (sizeof (GFC_REAL_4) * size); + ret->offset = 0; + } + + if (unlikely (n2 < n1)) + return; + + if (unlikely (compile_options.bounds_check) + && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) + runtime_error("Incorrect extent in return value of BESSEL_JN " + "(%ld vs. %ld)", (long int) n2-n1, + GFC_DESCRIPTOR_EXTENT(ret,0)); + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (unlikely (x == 0.0F)) + { + for (i = 0; i <= n2-n1; i++) +#if defined(GFC_REAL_4_INFINITY) + ret->data[i*stride] = -GFC_REAL_4_INFINITY; +#else + ret->data[i*stride] = -GFC_REAL_4_HUGE; +#endif + return; + } + + ret->data = ret->data; + last1 = ynf (n1, x); + ret->data[0] = last1; + + if (n1 == n2) + return; + + last2 = ynf (n1 + 1, x); + ret->data[1*stride] = last2; + + if (n1 + 1 == n2) + return; + + x2rev = 2.0F/x; + + for (i = 2; i <= n1+n2; i++) + { +#if defined(GFC_REAL_4_INFINITY) + if (unlikely (last2 == -GFC_REAL_4_INFINITY)) + { + ret->data[i*stride] = -GFC_REAL_4_INFINITY; + } + else +#endif + { + ret->data[i*stride] = x2rev * (i-1+n1) * last2 - last1; + last1 = last2; + last2 = ret->data[i*stride]; + } + } +} +#endif + +#endif + diff --git a/libgfortran/generated/bessel_r8.c b/libgfortran/generated/bessel_r8.c new file mode 100644 index 0000000..c896ae4 --- /dev/null +++ b/libgfortran/generated/bessel_r8.c @@ -0,0 +1,183 @@ +/* Implementation of the BESSEL_JN and BESSEL_YN transformational + function using a recurrence algorithm. + Copyright 2010 Free Software Foundation, Inc. + Contributed by Tobias Burnus + +This file is part of the GNU Fortran runtime library (libgfortran). + +Libgfortran is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public +License as published by the Free Software Foundation; either +version 3 of the License, or (at your option) any later version. + +Libgfortran is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +Under Section 7 of GPL version 3, you are granted additional +permissions described in the GCC Runtime Library Exception, version +3.1, as published by the Free Software Foundation. + +You should have received a copy of the GNU General Public License and +a copy of the GCC Runtime Library Exception along with this program; +see the files COPYING3 and COPYING.RUNTIME respectively. If not, see +. */ + +#include "libgfortran.h" +#include +#include + + +#if defined (HAVE_GFC_REAL_8) + + + +#if defined (HAVE_JN) +extern void bessel_jn_r8 (gfc_array_r8 * const restrict ret, int n1, + int n2, GFC_REAL_8 x); +export_proto(bessel_jn_r8); + +void +bessel_jn_r8 (gfc_array_r8 * const restrict ret, int n1, int n2, GFC_REAL_8 x) +{ + int i; + index_type stride; + + GFC_REAL_8 last1, last2, x2rev; + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (ret->data == NULL) + { + size_t size = n2 < n1 ? 0 : n2-n1+1; + GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); + ret->data = internal_malloc_size (sizeof (GFC_REAL_8) * size); + ret->offset = 0; + } + + if (unlikely (n2 < n1)) + return; + + if (unlikely (compile_options.bounds_check) + && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) + runtime_error("Incorrect extent in return value of BESSEL_JN " + "(%ld vs. %ld)", (long int) n2-n1, + GFC_DESCRIPTOR_EXTENT(ret,0)); + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (unlikely (x == 0.0)) + { + ret->data[0] = 1.0; + for (i = 1; i <= n2-n1; i++) + ret->data[i*stride] = 0.0; + return; + } + + ret->data = ret->data; + last1 = jn (n2, x); + ret->data[(n2-n1)*stride] = last1; + + if (n1 == n2) + return; + + last2 = jn (n2 - 1, x); + ret->data[(n2-n1-1)*stride] = last2; + + if (n1 + 1 == n2) + return; + + x2rev = 2.0/x; + + for (i = n2-n1-2; i >= 0; i--) + { + ret->data[i*stride] = x2rev * (i+1+n1) * last2 - last1; + last1 = last2; + last2 = ret->data[i*stride]; + } +} + +#endif + +#if defined (HAVE_YN) +extern void bessel_yn_r8 (gfc_array_r8 * const restrict ret, + int n1, int n2, GFC_REAL_8 x); +export_proto(bessel_yn_r8); + +void +bessel_yn_r8 (gfc_array_r8 * const restrict ret, int n1, int n2, + GFC_REAL_8 x) +{ + int i; + index_type stride; + + GFC_REAL_8 last1, last2, x2rev; + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (ret->data == NULL) + { + size_t size = n2 < n1 ? 0 : n2-n1+1; + GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1); + ret->data = internal_malloc_size (sizeof (GFC_REAL_8) * size); + ret->offset = 0; + } + + if (unlikely (n2 < n1)) + return; + + if (unlikely (compile_options.bounds_check) + && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1)) + runtime_error("Incorrect extent in return value of BESSEL_JN " + "(%ld vs. %ld)", (long int) n2-n1, + GFC_DESCRIPTOR_EXTENT(ret,0)); + + stride = GFC_DESCRIPTOR_STRIDE(ret,0); + + if (unlikely (x == 0.0)) + { + for (i = 0; i <= n2-n1; i++) +#if defined(GFC_REAL_8_INFINITY) + ret->data[i*stride] = -GFC_REAL_8_INFINITY; +#else + ret->data[i*stride] = -GFC_REAL_8_HUGE; +#endif + return; + } + + ret->data = ret->data; + last1 = yn (n1, x); + ret->data[0] = last1; + + if (n1 == n2) + return; + + last2 = yn (n1 + 1, x); + ret->data[1*stride] = last2; + + if (n1 + 1 == n2) + return; + + x2rev = 2.0/x; + + for (i = 2; i <= n1+n2; i++) + { +#if defined(GFC_REAL_8_INFINITY) + if (unlikely (last2 == -GFC_REAL_8_INFINITY)) + { + ret->data[i*stride] = -GFC_REAL_8_INFINITY; + } + else +#endif + { + ret->data[i*stride] = x2rev * (i-1+n1) * last2 - last1; + last1 = last2; + last2 = ret->data[i*stride]; + } + } +} +#endif + +#endif + -- cgit v1.1