Algol 68 Genie

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-2026 J. Marcel van der Veer [algol68g@algol68genie.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 (a68g_bufprt (A68G (edit_line), SNPRINTF_SIZE, "%s in line %d of file %s", reason, line, file) >= 0);
  38    } else {
  39      ASSERT (a68g_bufprt (A68G (edit_line), SNPRINTF_SIZE, "%s", reason) >= 0);
  40    }
  41    diagnostic (A68G_RUNTIME_ERROR, A68G (f_entry), ERROR_LAPLACE, A68G (edit_line), gsl_strerror (gsl_errno));
  42    exit_genie (A68G (f_entry), A68G_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 A68G_LAPLACE A68G_LAPLACE;
  59  
  60  struct A68G_LAPLACE
  61  {
  62    NODE_T *p;
  63    A68G_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    A68G_LAPLACE *l = (A68G_LAPLACE *) z;
  72    ADDR_T pop_sp = A68G_SP, pop_fp = A68G_FP;
  73    MOID_T *u = M_PROC_REAL_REAL;
  74    A68G_REAL *ft = (A68G_REAL *) STACK_TOP;
  75    PUSH_VALUE (P (l), t, A68G_REAL);
  76    genie_call_procedure (P (l), MOID (&(F (l))), u, u, &(F (l)), pop_sp, pop_fp);
  77    A68G_SP = pop_sp;
  78    return VALUE (ft) * a68g_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    A68G_REF ref_err;
  87    POP_REF (p, &ref_err);
  88    CHECK_REF (p, ref_err, M_REF_REAL);
  89    A68G_REAL *err = (A68G_REAL *) ADDRESS (&ref_err);
  90    A68G_REAL s;
  91    POP_OBJECT (p, &s, A68G_REAL);
  92    A68G_PROCEDURE f;
  93    POP_PROCEDURE (p, &f);
  94    A68G_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, A68G_REAL);
 111    gsl_integration_workspace_free (w);
 112    (void) gsl_set_error_handler (save_handler);
 113  }
 114  
 115  #endif