# HG changeset patch
# User Arun Giridhar <arungiridhar@gmail.com>
# Date 1650288732 14400
# Node ID 84e7222b6b5cf6ea83f6eb037f71a0afa0301fb9
# Parent  272614f05636c236d3f92d18d5bc30593a6cc6e4
lo-specfun.cc: Remove code duplication for airy and biry functions (bug #62321)

Co-author: Zoïs Moitier <zois.moitier@kit.edu>

lo-specfun.cc: The external Fortran library functions ZAIRY, CAIRY, ZBIRY, and
CBIRY can already scale or unscale their respective results depending on
whether an input parameter KODE is set to 1 or to 2. This changeset removes
the duplication of that scaling functionality in lo-specfun.cc, and calls the
library functions with the appropriate value of KODE.

diff -r 272614f05636 -r 84e7222b6b5c liboctave/numeric/lo-specfun.cc
--- a/liboctave/numeric/lo-specfun.cc	Sun Apr 17 17:58:14 2022 -0400
+++ b/liboctave/numeric/lo-specfun.cc	Mon Apr 18 09:32:12 2022 -0400
@@ -139,23 +139,13 @@
       F77_INT id = (deriv ? 1 : 0);
       F77_INT nz, t_ierr;
 
-      F77_FUNC (zairy, ZAIRY) (zr, zi, id, 2, ar, ai, nz, t_ierr);
+      if (scaled)
+        F77_FUNC (zairy, ZAIRY) (zr, zi, id, 2, ar, ai, nz, t_ierr);
+      else
+        F77_FUNC (zairy, ZAIRY) (zr, zi, id, 1, ar, ai, nz, t_ierr);
 
       ierr = t_ierr;
 
-      if (! scaled)
-        {
-          Complex expz = exp (- 2.0 / 3.0 * z * sqrt (z));
-
-          double rexpz = expz.real ();
-          double iexpz = expz.imag ();
-
-          double tmp = ar*rexpz - ai*iexpz;
-
-          ai = ar*iexpz + ai*rexpz;
-          ar = tmp;
-        }
-
       if (zi == 0.0 && (! scaled || zr >= 0.0))
         ai = 0.0;
 
@@ -205,27 +195,18 @@
       F77_INT id = (deriv ? 1 : 0);
       F77_INT nz, t_ierr;
 
-      F77_FUNC (cairy, CAIRY) (F77_CONST_CMPLX_ARG (&z), id, 2,
-                               F77_CMPLX_ARG (&a), nz, t_ierr);
+      if (scaled)
+        F77_FUNC (cairy, CAIRY) (F77_CONST_CMPLX_ARG (&z), id, 2,
+                                 F77_CMPLX_ARG (&a), nz, t_ierr);
+      else
+        F77_FUNC (cairy, CAIRY) (F77_CONST_CMPLX_ARG (&z), id, 1,
+                                 F77_CMPLX_ARG (&a), nz, t_ierr);
 
       ierr = t_ierr;
 
       float ar = a.real ();
       float ai = a.imag ();
 
-      if (! scaled)
-        {
-          FloatComplex expz = exp (- 2.0f / 3.0f * z * sqrt (z));
-
-          float rexpz = expz.real ();
-          float iexpz = expz.imag ();
-
-          float tmp = ar*rexpz - ai*iexpz;
-
-          ai = ar*iexpz + ai*rexpz;
-          ar = tmp;
-        }
-
       if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0))
         ai = 0.0;
 
@@ -1385,23 +1366,13 @@
       F77_INT id = (deriv ? 1 : 0);
       F77_INT t_ierr;
 
-      F77_FUNC (zbiry, ZBIRY) (zr, zi, id, 2, ar, ai, t_ierr);
+      if (scaled)
+        F77_FUNC (zbiry, ZBIRY) (zr, zi, id, 2, ar, ai, t_ierr);
+      else
+        F77_FUNC (zbiry, ZBIRY) (zr, zi, id, 1, ar, ai, t_ierr);
 
       ierr = t_ierr;
 
-      if (! scaled)
-        {
-          Complex expz = exp (std::abs (std::real (2.0 / 3.0 * z * sqrt (z))));
-
-          double rexpz = expz.real ();
-          double iexpz = expz.imag ();
-
-          double tmp = ar*rexpz - ai*iexpz;
-
-          ai = ar*iexpz + ai*rexpz;
-          ar = tmp;
-        }
-
       if (zi == 0.0 && (! scaled || zr >= 0.0))
         ai = 0.0;
 
@@ -1451,28 +1422,18 @@
       F77_INT id = (deriv ? 1 : 0);
       F77_INT t_ierr;
 
-      F77_FUNC (cbiry, CBIRY) (F77_CONST_CMPLX_ARG (&z), id, 2,
-                               F77_CMPLX_ARG (&a), t_ierr);
+      if (scaled)
+        F77_FUNC (cbiry, CBIRY) (F77_CONST_CMPLX_ARG (&z), id, 2,
+                                 F77_CMPLX_ARG (&a), t_ierr);
+      else
+        F77_FUNC (cbiry, CBIRY) (F77_CONST_CMPLX_ARG (&z), id, 1,
+                                 F77_CMPLX_ARG (&a), t_ierr);
 
       ierr = t_ierr;
 
       float ar = a.real ();
       float ai = a.imag ();
 
-      if (! scaled)
-        {
-          FloatComplex expz
-            = exp (std::abs (std::real (2.0f / 3.0f * z * sqrt (z))));
-
-          float rexpz = expz.real ();
-          float iexpz = expz.imag ();
-
-          float tmp = ar*rexpz - ai*iexpz;
-
-          ai = ar*iexpz + ai*rexpz;
-          ar = tmp;
-        }
-
       if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0))
         ai = 0.0;