single-laplace.c
1 //! @file single-laplace.c
2 //! @author J. Marcel van der Veer
3
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
8
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! REAL GSL laplace routines.
25
26 #include "a68g.h"
27 #include "a68g-genie.h"
28 #include "a68g-prelude.h"
29
30 #if defined (HAVE_GSL)
31
32 //! @brief Map GSL error handler onto a68g error handler.
33
34 static void laplace_error_handler (const char *reason, const char *file, int line, int gsl_errno)
35 {
36 if (line != 0) {
37 ASSERT (snprintf (A68 (edit_line), SNPRINTF_SIZE, "%s in line %d of file %s", reason, line, file) >= 0);
38 } else {
39 ASSERT (snprintf (A68 (edit_line), SNPRINTF_SIZE, "%s", reason) >= 0);
40 }
41 diagnostic (A68_RUNTIME_ERROR, A68 (f_entry), ERROR_LAPLACE, A68 (edit_line), gsl_strerror (gsl_errno));
42 exit_genie (A68 (f_entry), A68_RUNTIME_ERROR);
43 }
44
45 //! @brief Detect math errors.
46
47 static void laplace_test_error (int ret)
48 {
49 if (ret != 0) {
50 laplace_error_handler ("math error", "", 0, ret);
51 }
52 }
53
54 //! @brief PROC (PROC (REAL) REAL, REAL, REF REAL) REAL laplace
55
56 #define LAPLACE_DIVISIONS 1024
57
58 typedef struct A68_LAPLACE A68_LAPLACE;
59
60 struct A68_LAPLACE
61 {
62 NODE_T *p;
63 A68_PROCEDURE f;
64 REAL_T s;
65 };
66
67 //! @brief Evaluate function for Laplace transform.
68
69 static REAL_T laplace_f (REAL_T t, void *z)
70 {
71 A68_LAPLACE *l = (A68_LAPLACE *) z;
72 ADDR_T pop_sp = A68_SP, pop_fp = A68_FP;
73 MOID_T *u = M_PROC_REAL_REAL;
74 A68_REAL *ft = (A68_REAL *) STACK_TOP;
75 PUSH_VALUE (P (l), t, A68_REAL);
76 genie_call_procedure (P (l), MOID (&(F (l))), u, u, &(F (l)), pop_sp, pop_fp);
77 A68_SP = pop_sp;
78 return VALUE (ft) * a68_exp_real (-(S (l)) * t);
79 }
80
81 //! @brief PROC laplace = (PROC (REAL) REAL, REAL, REF REAL) REAL
82
83 void genie_laplace (NODE_T * p)
84 {
85 gsl_error_handler_t *save_handler = gsl_set_error_handler (laplace_error_handler);
86 A68_REF ref_err;
87 POP_REF (p, &ref_err);
88 CHECK_REF (p, ref_err, M_REF_REAL);
89 A68_REAL *err = (A68_REAL *) ADDRESS (&ref_err);
90 A68_REAL s;
91 POP_OBJECT (p, &s, A68_REAL);
92 A68_PROCEDURE f;
93 POP_PROCEDURE (p, &f);
94 A68_LAPLACE l;
95 P (&l) = p;
96 F (&l) = f;
97 S (&l) = VALUE (&s);
98 gsl_function g;
99 FUNCTION (&g) = &laplace_f;
100 GSL_PARAMS (&g) = &l;
101 gsl_integration_workspace *w = gsl_integration_workspace_alloc (LAPLACE_DIVISIONS);
102 REAL_T res, estimated_err; int ret;
103 if (VALUE (err) >= 0.0) {
104 ret = gsl_integration_qagiu (&g, 0.0, VALUE (err), 0.0, LAPLACE_DIVISIONS, w, &res, &estimated_err);
105 } else {
106 ret = gsl_integration_qagiu (&g, 0.0, 0.0, -VALUE (err), LAPLACE_DIVISIONS, w, &res, &estimated_err);
107 }
108 laplace_test_error (ret);
109 VALUE (err) = estimated_err;
110 PUSH_VALUE (p, res, A68_REAL);
111 gsl_integration_workspace_free (w);
112 (void) gsl_set_error_handler (save_handler);
113 }
114
115 #endif
