home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fs.zip / octave / octave-2.1.23 / liboctave / lo-specfun.cc < prev    next >
C/C++ Source or Header  |  2000-01-15  |  19KB  |  926 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. /* Modified by Klaus Gebhardt, 1997-1999 */
  24.  
  25. #ifdef HAVE_CONFIG_H
  26. #include <config.h>
  27. #endif
  28.  
  29. #ifdef __EMX__
  30. #include <float.h>
  31. #endif
  32.  
  33. #include "Range.h"
  34. #include "CColVector.h"
  35. #include "CMatrix.h"
  36. #include "dRowVector.h"
  37. #include "dMatrix.h"
  38. #include "f77-fcn.h"
  39. #include "lo-error.h"
  40. #include "lo-ieee.h"
  41. #include "lo-specfun.h"
  42. #include "mx-inlines.cc"
  43.  
  44. extern "C"
  45. {
  46.   int F77_FCN (zbesj, ZBESJ) (const double&, const double&,
  47.                   const double&, const int&, const int&,
  48.                   double*, double*, int&, int&);
  49.  
  50.   int F77_FCN (zbesy, ZBESY) (const double&, const double&,
  51.                   const double&, const int&, const int&,
  52.                   double*, double*, int&,
  53.                   double*, double*, int&);
  54.  
  55.   int F77_FCN (zbesi, ZBESI) (const double&, const double&,
  56.                   const double&, const int&, const int&,
  57.                   double*, double*, int&, int&);
  58.  
  59.   int F77_FCN (zbesk, ZBESK) (const double&, const double&,
  60.                   const double&, const int&, const int&,
  61.                   double*, double*, int&, int&);
  62.  
  63.   int F77_FCN (zbesh, ZBESH) (const double&, const double&,
  64.                   const double&, const int&, const int&,
  65.                   const int&, double*, double*, int&, int&);
  66.  
  67.   int F77_FCN (zairy, ZAIRY) (const double&, const double&,
  68.                   const int&, const int&,
  69.                   double&, double&, int&, int&);
  70.  
  71.   int F77_FCN (zbiry, ZBIRY) (const double&, const double&,
  72.                   const int&, const int&,
  73.                   double&, double&, int&);
  74.  
  75.   int F77_FCN (xdacosh, XDACOSH) (const double&, double&);
  76.  
  77.   int F77_FCN (xdasinh, XDASINH) (const double&, double&);
  78.  
  79.   int F77_FCN (xdatanh, XDATANH) (const double&, double&);
  80.  
  81.   int F77_FCN (xderf, XDERF) (const double&, double&);
  82.  
  83.   int F77_FCN (xderfc, XDERFC) (const double&, double&);
  84.  
  85.   int F77_FCN (xdbetai, XDBETAI) (const double&, const double&,
  86.                   const double&, double&);
  87.  
  88.   int F77_FCN (xdgamma, XDGAMMA) (const double&, double&);
  89.  
  90.   int F77_FCN (xdgamit, XDGAMIT) (const double&, const double&, double&);
  91.  
  92.   int F77_FCN (dlgams, DLGAMS) (const double&, double&, double&);
  93. }
  94.  
  95. #if !defined (HAVE_ACOSH)
  96. double
  97. acosh (double x)
  98. {
  99.   double retval;
  100.   F77_XFCN (xdacosh, DACOSH, (x, retval));
  101.   return retval;
  102. }
  103. #endif
  104.  
  105. #if !defined (HAVE_ASINH)
  106. double
  107. asinh (double x)
  108. {
  109.   double retval;
  110.   F77_XFCN (xdasinh, DASINH, (x, retval));
  111.   return retval;
  112. }
  113. #endif
  114.  
  115. #if !defined (HAVE_ATANH)
  116. double
  117. atanh (double x)
  118. {
  119.   double retval;
  120.   F77_XFCN (xdatanh, DATANH, (x, retval));
  121.   return retval;
  122. }
  123. #endif
  124.  
  125. #if !defined (HAVE_ERF)
  126. double
  127. erf (double x)
  128. {
  129.   double retval;
  130.   F77_XFCN (xderf, DERF, (x, retval));
  131.   return retval;
  132. }
  133. #endif
  134.  
  135. #if !defined (HAVE_ERFC)
  136. double
  137. erfc (double x)
  138. {
  139.   double retval;
  140.   F77_XFCN (xderfc, DERFC, (x, retval));
  141.   return retval;
  142. }
  143. #endif
  144.  
  145. double
  146. xgamma (double x)
  147. {
  148.   double result;
  149.   F77_XFCN (xdgamma, XDGAMMA, (x, result));
  150.   return result;
  151. }
  152.  
  153. double
  154. xlgamma (double x)
  155. {
  156.   double result;
  157.   double sgngam;
  158.   F77_XFCN (dlgams, DLGAMS, (x, result, sgngam));
  159.   return result;
  160. }
  161.  
  162. static inline Complex
  163. zbesj (const Complex& z, double alpha, int kode, int& ierr);
  164.  
  165. static inline Complex
  166. zbesy (const Complex& z, double alpha, int kode, int& ierr);
  167.  
  168. static inline Complex
  169. zbesi (const Complex& z, double alpha, int kode, int& ierr);
  170.  
  171. static inline Complex
  172. zbesk (const Complex& z, double alpha, int kode, int& ierr);
  173.  
  174. static inline Complex
  175. zbesh1 (const Complex& z, double alpha, int kode, int& ierr);
  176.  
  177. static inline Complex
  178. zbesh2 (const Complex& z, double alpha, int kode, int& ierr);
  179.  
  180. static inline Complex
  181. bessel_return_value (const Complex& val, int ierr)
  182. {
  183.   static const Complex inf_val = Complex (octave_Inf, octave_Inf);
  184.   static const Complex nan_val = Complex (octave_NaN, octave_NaN);
  185.  
  186.   Complex retval;
  187.  
  188.   switch (ierr)
  189.     {
  190.     case 0:
  191.     case 3:
  192.       retval = val;
  193.       break;
  194.  
  195.     case 2:
  196.       retval = inf_val;
  197.       break;
  198.  
  199.     default:
  200.       retval = nan_val;
  201.       break;
  202.     }
  203.  
  204.   return retval;
  205. }
  206.  
  207. static inline Complex
  208. zbesj (const Complex& z, double alpha, int kode, int& ierr)
  209. {
  210.   Complex retval;
  211.  
  212.   if (alpha >= 0.0)
  213.     {
  214.       double yr = 0.0;
  215.       double yi = 0.0;
  216.  
  217.       int nz;
  218.  
  219.       double zr = z.real ();
  220.       double zi = z.imag ();
  221.  
  222.       F77_FCN (zbesj, ZBESJ) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr);
  223.  
  224.       if (zi == 0.0 && zr > 0.0)
  225.     yi = 0.0;
  226.  
  227.       retval = bessel_return_value (Complex (yr, yi), ierr);
  228.     }
  229.   else
  230.     {
  231.       alpha = -alpha;
  232.  
  233.       Complex tmp = cos (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
  234.  
  235.       if (ierr == 0 || ierr == 3)
  236.     {
  237.       tmp -= sin (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
  238.  
  239.       retval = bessel_return_value (tmp, ierr);
  240.     }
  241.       else
  242.     retval = Complex (octave_NaN, octave_NaN);
  243.     }
  244.  
  245.   return retval;
  246. }
  247.  
  248. static inline Complex
  249. zbesy (const Complex& z, double alpha, int kode, int& ierr)
  250. {
  251.   Complex retval;
  252.  
  253.   if (alpha >= 0.0)
  254.     {
  255.       double yr = 0.0;
  256.       double yi = 0.0;
  257.  
  258.       int nz;
  259.  
  260.       double wr, wi;
  261.  
  262.       double zr = z.real ();
  263.       double zi = z.imag ();
  264.  
  265.       ierr = 0;
  266.  
  267.       if (zr == 0.0 && zi == 0.0)
  268.     {
  269.       yr = -octave_Inf;
  270.       yi = 0.0;
  271.     }
  272.       else
  273.     {
  274.       F77_FCN (zbesy, ZBESY) (zr, zi, alpha, kode, 1, &yr, &yi, nz,
  275.                   &wr, &wi, ierr);
  276.  
  277.       if (zi == 0.0 && zr > 0.0)
  278.         yi = 0.0;
  279.     }
  280.  
  281.       return bessel_return_value (Complex (yr, yi), ierr);
  282.     }
  283.   else
  284.     {
  285.       alpha = -alpha;
  286.  
  287.       Complex tmp = cos (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
  288.  
  289.       if (ierr == 0 || ierr == 3)
  290.     {
  291.       tmp += sin (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
  292.  
  293.       retval = bessel_return_value (tmp, ierr);
  294.     }
  295.       else
  296.     retval = Complex (octave_NaN, octave_NaN);
  297.     }
  298.  
  299.   return retval;
  300. }
  301.  
  302. static inline Complex
  303. zbesi (const Complex& z, double alpha, int kode, int& ierr)
  304. {
  305.   Complex retval;
  306.  
  307.   if (alpha >= 0.0)
  308.     {
  309.       double yr = 0.0;
  310.       double yi = 0.0;
  311.  
  312.       int nz;
  313.  
  314.       double zr = z.real ();
  315.       double zi = z.imag ();
  316.  
  317.       F77_FCN (zbesi, ZBESI) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr);
  318.  
  319.       if (zi == 0.0 && zr > 0.0)
  320.     yi = 0.0;
  321.  
  322.       retval = bessel_return_value (Complex (yr, yi), ierr);
  323.     }
  324.   else
  325.     {
  326.       alpha = -alpha;
  327.  
  328.       Complex tmp = zbesi (z, alpha, kode, ierr);
  329.  
  330.       if (ierr == 0 || ierr == 3)
  331.     {
  332.       tmp += (2.0 / M_PI) * sin (M_PI * alpha)
  333.         * zbesk (z, alpha, kode, ierr);
  334.  
  335.       retval = bessel_return_value (tmp, ierr);
  336.     }
  337.       else
  338.     retval = Complex (octave_NaN, octave_NaN);
  339.     }
  340.  
  341.   return retval;
  342. }
  343.  
  344. static inline Complex
  345. zbesk (const Complex& z, double alpha, int kode, int& ierr)
  346. {
  347.   Complex retval;
  348.  
  349.   if (alpha >= 0.0)
  350.     {
  351.       double yr = 0.0;
  352.       double yi = 0.0;
  353.  
  354.       int nz;
  355.  
  356.       double zr = z.real ();
  357.       double zi = z.imag ();
  358.  
  359.       ierr = 0;
  360.  
  361.       if (zr == 0.0 && zi == 0.0)
  362.     {
  363.       yr = octave_Inf;
  364.       yi = 0.0;
  365.     }
  366.       else
  367.     {
  368.       F77_FCN (zbesk, ZBESK) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr);
  369.  
  370.       if (zi == 0.0 && zr > 0.0)
  371.         yi = 0.0;
  372.     }
  373.  
  374.       retval = bessel_return_value (Complex (yr, yi), ierr);
  375.     }
  376.   else
  377.     {
  378.       Complex tmp = zbesk (z, -alpha, kode, ierr);
  379.  
  380.       retval = bessel_return_value (tmp, ierr);
  381.     }
  382.  
  383.   return retval;
  384. }
  385.  
  386. static inline Complex
  387. zbesh1 (const Complex& z, double alpha, int kode, int& ierr)
  388. {
  389.   Complex retval;
  390.  
  391.   if (alpha >= 0.0)
  392.     {
  393.       double yr = 0.0;
  394.       double yi = 0.0;
  395.  
  396.       int nz;
  397.  
  398.       double zr = z.real ();
  399.       double zi = z.imag ();
  400.  
  401.       F77_FCN (zbesh, ZBESH) (zr, zi, alpha, kode, 1, 1, &yr, &yi, nz, ierr);
  402.  
  403.       retval = bessel_return_value (Complex (yr, yi), ierr);
  404.     }
  405.   else
  406.     {
  407.       alpha = -alpha;
  408.  
  409.       static const Complex eye = Complex (0.0, 1.0);
  410.  
  411.       Complex tmp = exp (M_PI * alpha * eye) * zbesh1 (z, alpha, kode, ierr);
  412.  
  413.       retval = bessel_return_value (tmp, ierr);
  414.     }
  415.  
  416.   return retval;
  417. }
  418.  
  419. static inline Complex
  420. zbesh2 (const Complex& z, double alpha, int kode, int& ierr)
  421. {
  422.   Complex retval;
  423.  
  424.   if (alpha >= 0.0)
  425.     {
  426.       double yr = 0.0;
  427.       double yi = 0.0;
  428.  
  429.       int nz;
  430.  
  431.       double zr = z.real ();
  432.       double zi = z.imag ();
  433.  
  434.       F77_FCN (zbesh, ZBESH) (zr, zi, alpha, kode, 2, 1, &yr, &yi, nz, ierr);
  435.  
  436.       retval = bessel_return_value (Complex (yr, yi), ierr);
  437.     }
  438.   else
  439.     {
  440.       alpha = -alpha;
  441.  
  442.       static const Complex eye = Complex (0.0, 1.0);
  443.  
  444.       Complex tmp = exp (-M_PI * alpha * eye) * zbesh2 (z, alpha, kode, ierr);
  445.  
  446.       retval = bessel_return_value (tmp, ierr);
  447.     }
  448.  
  449.   return retval;
  450. }
  451.  
  452. typedef Complex (*fptr) (const Complex&, double, int, int&);
  453.  
  454. static inline Complex
  455. do_bessel (fptr f, const char *, double alpha, const Complex& x,
  456.        bool scaled, int& ierr)
  457. {
  458.   Complex retval;
  459.  
  460.   retval = f (x, alpha, (scaled ? 2 : 1), ierr);
  461.  
  462.   return retval;
  463. }
  464.  
  465. static inline ComplexMatrix
  466. do_bessel (fptr f, const char *, double alpha, const ComplexMatrix& x,
  467.        bool scaled, Array2<int>& ierr)
  468. {
  469.   int nr = x.rows ();
  470.   int nc = x.cols ();
  471.  
  472.   ComplexMatrix retval (nr, nc);
  473.  
  474.   ierr.resize (nr, nc);
  475.  
  476.   for (int j = 0; j < nc; j++)
  477.     for (int i = 0; i < nr; i++)
  478.       retval(i,j) = f (x(i,j), alpha, (scaled ? 2 : 1), ierr(i,j));
  479.  
  480.   return retval;
  481. }
  482.  
  483. static inline ComplexMatrix
  484. do_bessel (fptr f, const char *, const Matrix& alpha, const Complex& x,
  485.        bool scaled, Array2<int>& ierr)
  486. {
  487.   int nr = alpha.rows ();
  488.   int nc = alpha.cols ();
  489.  
  490.   ComplexMatrix retval (nr, nc);
  491.  
  492.   ierr.resize (nr, nc);
  493.  
  494.   for (int j = 0; j < nc; j++)
  495.     for (int i = 0; i < nr; i++)
  496.       retval(i,j) = f (x, alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
  497.  
  498.   return retval;
  499. }
  500.  
  501. static inline ComplexMatrix
  502. do_bessel (fptr f, const char *fn, const Matrix& alpha,
  503.        const ComplexMatrix& x, bool scaled, Array2<int>& ierr)
  504. {
  505.   ComplexMatrix retval;
  506.  
  507.   int x_nr = x.rows ();
  508.   int x_nc = x.cols ();
  509.  
  510.   int alpha_nr = alpha.rows ();
  511.   int alpha_nc = alpha.cols ();
  512.  
  513.   if (x_nr == alpha_nr && x_nc == alpha_nc)
  514.     {
  515.       int nr = x_nr;
  516.       int nc = x_nc;
  517.  
  518.       retval.resize (nr, nc);
  519.  
  520.       ierr.resize (nr, nc);
  521.  
  522.       for (int j = 0; j < nc; j++)
  523.     for (int i = 0; i < nr; i++)
  524.       retval(i,j) = f (x(i,j), alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
  525.     }
  526.   else
  527.     (*current_liboctave_error_handler)
  528.       ("%s: the sizes of alpha and x must conform", fn);
  529.  
  530.   return retval;
  531. }
  532.  
  533. static inline ComplexMatrix
  534. do_bessel (fptr f, const char *, const RowVector& alpha,
  535.        const ComplexColumnVector& x, bool scaled, Array2<int>& ierr)
  536. {
  537.   int nr = x.length ();
  538.   int nc = alpha.length ();
  539.  
  540.   ComplexMatrix retval (nr, nc);
  541.  
  542.   ierr.resize (nr, nc);
  543.  
  544.   for (int j = 0; j < nc; j++)
  545.     for (int i = 0; i < nr; i++)
  546.       retval(i,j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i,j));
  547.  
  548.   return retval;
  549. }
  550.  
  551. #define SS_BESSEL(name, fcn) \
  552.   Complex \
  553.   name (double alpha, const Complex& x, bool scaled, int& ierr) \
  554.   { \
  555.     return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
  556.   }
  557.  
  558. #define SM_BESSEL(name, fcn) \
  559.   ComplexMatrix \
  560.   name (double alpha, const ComplexMatrix& x, bool scaled, \
  561.     Array2<int>& ierr) \
  562.   { \
  563.     return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
  564.   }
  565.  
  566. #define MS_BESSEL(name, fcn) \
  567.   ComplexMatrix \
  568.   name (const Matrix& alpha, const Complex& x, bool scaled, \
  569.     Array2<int>& ierr) \
  570.   { \
  571.     return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
  572.   }
  573.  
  574. #define MM_BESSEL(name, fcn) \
  575.   ComplexMatrix \
  576.   name (const Matrix& alpha, const ComplexMatrix& x, bool scaled, \
  577.     Array2<int>& ierr) \
  578.   { \
  579.     return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
  580.   }
  581.  
  582. #define RC_BESSEL(name, fcn) \
  583.   ComplexMatrix \
  584.   name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
  585.         Array2<int>& ierr) \
  586.   { \
  587.     return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
  588.   }
  589.  
  590. #define ALL_BESSEL(name, fcn) \
  591.   SS_BESSEL (name, fcn) \
  592.   SM_BESSEL (name, fcn) \
  593.   MS_BESSEL (name, fcn) \
  594.   MM_BESSEL (name, fcn) \
  595.   RC_BESSEL (name, fcn)
  596.  
  597. ALL_BESSEL (besselj, zbesj)
  598. ALL_BESSEL (bessely, zbesy)
  599. ALL_BESSEL (besseli, zbesi)
  600. ALL_BESSEL (besselk, zbesk)
  601. ALL_BESSEL (besselh1, zbesh1)
  602. ALL_BESSEL (besselh2, zbesh2)
  603.  
  604. Complex
  605. airy (const Complex& z, bool deriv, bool scaled, int& ierr)
  606. {
  607.   double ar = 0.0;
  608.   double ai = 0.0;
  609.  
  610.   int nz;
  611.  
  612.   double zr = z.real ();
  613.   double zi = z.imag ();
  614.  
  615.   int id = deriv ? 1 : 0;
  616.  
  617.   int kode = scaled ? 2 : 1;
  618.  
  619.   F77_FCN (zairy, ZAIRY) (zr, zi, id, kode, ar, ai, nz, ierr);
  620.  
  621.   if (zi == 0.0 && (! scaled || zr > 0.0))
  622.     ai = 0.0;
  623.  
  624.   return bessel_return_value (Complex (ar, ai), ierr);
  625. }
  626.  
  627. Complex
  628. biry (const Complex& z, bool deriv, bool scaled, int& ierr)
  629. {
  630.   double ar = 0.0;
  631.   double ai = 0.0;
  632.  
  633.   double zr = z.real ();
  634.   double zi = z.imag ();
  635.  
  636.   int id = deriv ? 1 : 0;
  637.  
  638.   int kode = scaled ? 2 : 1;
  639.  
  640.   F77_FCN (zbiry, ZBIRY) (zr, zi, id, kode, ar, ai, ierr);
  641.  
  642.   if (zi == 0.0 && (! scaled || zr > 0.0))
  643.     ai = 0.0;
  644.  
  645.   return bessel_return_value (Complex (ar, ai), ierr);
  646. }
  647.  
  648. ComplexMatrix
  649. airy (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr)
  650. {
  651.   int nr = z.rows ();
  652.   int nc = z.cols ();
  653.  
  654.   ComplexMatrix retval (nr, nc);
  655.  
  656.   ierr.resize (nr, nc);
  657.  
  658.   for (int j = 0; j < nc; j++)
  659.     for (int i = 0; i < nr; i++)
  660.       retval(i,j) = airy (z(i,j), deriv, scaled, ierr(i,j));
  661.  
  662.   return retval;
  663. }
  664.  
  665. ComplexMatrix
  666. biry (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr)
  667. {
  668.   int nr = z.rows ();
  669.   int nc = z.cols ();
  670.  
  671.   ComplexMatrix retval (nr, nc);
  672.  
  673.   ierr.resize (nr, nc);
  674.  
  675.   for (int j = 0; j < nc; j++)
  676.     for (int i = 0; i < nr; i++)
  677.       retval(i,j) = biry (z(i,j), deriv, scaled, ierr(i,j));
  678.  
  679.   return retval;
  680. }
  681.  
  682. static void
  683. gripe_betainc_nonconformant (int r1, int c1, int r2, int c2, int r3,
  684.                  int c3)
  685. {
  686.   (*current_liboctave_error_handler)
  687.    ("betainc: nonconformant arguments (x is %dx%d, a is %dx%d, b is %dx%d)",
  688.      r1, c1, r2, c2, r3, c3);
  689. }
  690.  
  691. double
  692. betainc (double x, double a, double b)
  693. {
  694.   double retval;
  695.   F77_XFCN (xdbetai, XDBETAI, (x, a, b, retval));
  696.   return retval;
  697. }
  698.  
  699. Matrix
  700. betainc (double x, double a, const Matrix& b)
  701. {
  702.   int nr = b.rows ();
  703.   int nc = b.cols ();
  704.  
  705.   Matrix retval (nr, nc);
  706.  
  707.   for (int j = 0; j < nc; j++)
  708.     for (int i = 0; i < nr; i++)
  709.       retval(i,j) = betainc (x, a, b(i,j));
  710.  
  711.   return retval;
  712. }
  713.  
  714. Matrix
  715. betainc (double x, const Matrix& a, double b)
  716. {
  717.   int nr = a.rows ();
  718.   int nc = a.cols ();
  719.  
  720.   Matrix retval (nr, nc);
  721.  
  722.   for (int j = 0; j < nc; j++)
  723.     for (int i = 0; i < nr; i++)
  724.       retval(i,j) = betainc (x, a(i,j), b);
  725.  
  726.   return retval;
  727. }
  728.  
  729. Matrix
  730. betainc (double x, const Matrix& a, const Matrix& b)
  731. {
  732.   Matrix retval;
  733.  
  734.   int a_nr = a.rows ();
  735.   int a_nc = a.cols ();
  736.  
  737.   int b_nr = b.rows ();
  738.   int b_nc = b.cols ();
  739.  
  740.   if (a_nr == b_nr && a_nc == b_nc)
  741.     {
  742.       retval.resize (a_nr, a_nc);
  743.  
  744.       for (int j = 0; j < a_nc; j++)
  745.     for (int i = 0; i < a_nr; i++)
  746.       retval(i,j) = betainc (x, a(i,j), b(i,j));
  747.     }
  748.   else
  749.     gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc);
  750.  
  751.   return retval;
  752. }
  753.  
  754. Matrix
  755. betainc (const Matrix& x, double a, double b)
  756. {
  757.   int nr = x.rows ();
  758.   int nc = x.cols ();
  759.  
  760.   Matrix retval (nr, nc);
  761.  
  762.   for (int j = 0; j < nc; j++)
  763.     for (int i = 0; i < nr; i++)
  764.       retval(i,j) = betainc (x(i,j), a, b);
  765.  
  766.   return retval;
  767. }
  768.  
  769. Matrix
  770. betainc (const Matrix& x, double a, const Matrix& b)
  771. {
  772.   Matrix retval;
  773.  
  774.   int nr = x.rows ();
  775.   int nc = x.cols ();
  776.  
  777.   int b_nr = b.rows ();
  778.   int b_nc = b.cols ();
  779.  
  780.   if (nr == b_nr && nc == b_nc)
  781.     {
  782.       retval.resize (nr, nc);
  783.  
  784.       for (int j = 0; j < nc; j++)
  785.     for (int i = 0; i < nr; i++)
  786.       retval(i,j) = betainc (x(i,j), a, b(i,j));
  787.     }
  788.   else
  789.     gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc);
  790.  
  791.   return retval;
  792. }
  793.  
  794. Matrix
  795. betainc (const Matrix& x, const Matrix& a, double b)
  796. {
  797.   Matrix retval;
  798.  
  799.   int nr = x.rows ();
  800.   int nc = x.cols ();
  801.  
  802.   int a_nr = a.rows ();
  803.   int a_nc = a.cols ();
  804.  
  805.   if (nr == a_nr && nc == a_nc)
  806.     {
  807.       retval.resize (nr, nc);
  808.  
  809.       for (int j = 0; j < nc; j++)
  810.     for (int i = 0; i < nr; i++)
  811.       retval(i,j) = betainc (x(i,j), a(i,j), b);
  812.     }
  813.   else
  814.     gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1);
  815.  
  816.   return retval;
  817. }
  818.  
  819. Matrix
  820. betainc (const Matrix& x, const Matrix& a, const Matrix& b)
  821. {
  822.   Matrix retval;
  823.  
  824.   int nr = x.rows ();
  825.   int nc = x.cols ();
  826.  
  827.   int a_nr = a.rows ();
  828.   int a_nc = a.cols ();
  829.  
  830.   int b_nr = b.rows ();
  831.   int b_nc = b.cols ();
  832.  
  833.   if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc)
  834.     {
  835.       retval.resize (nr, nc);
  836.  
  837.       for (int j = 0; j < nc; j++)
  838.     for (int i = 0; i < nr; i++)
  839.       retval(i,j) = betainc (x(i,j), a(i,j), b(i,j));
  840.     }
  841.   else
  842.     gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc);
  843.  
  844.   return retval;
  845. }
  846.  
  847. // XXX FIXME XXX -- there is still room for improvement here...
  848.  
  849. double
  850. gammainc (double x, double a)
  851. {
  852.   double retval;
  853.  
  854.   F77_XFCN (xdgamit, XDGAMIT, (a, x, retval));
  855.  
  856.   if (x == 0.0)
  857.     retval = 0.0;
  858.   else if (x > 0.0)
  859.     retval = exp (a * log (x) + log (retval));
  860.  
  861.   return retval;
  862. }
  863.  
  864. Matrix
  865. gammainc (double x, const Matrix& a)
  866. {
  867.   int nr = a.rows ();
  868.   int nc = a.cols ();
  869.  
  870.   Matrix retval (nr, nc);
  871.  
  872.   for (int j = 0; j < nc; j++)
  873.     for (int i = 0; i < nr; i++)
  874.       retval(i,j) = gammainc (x, a(i,j));
  875.  
  876.   return retval;
  877. }
  878.  
  879. Matrix
  880. gammainc (const Matrix& x, double a)
  881. {
  882.   int nr = x.rows ();
  883.   int nc = x.cols ();
  884.  
  885.   Matrix retval (nr, nc);
  886.  
  887.   for (int j = 0; j < nc; j++)
  888.     for (int i = 0; i < nr; i++)
  889.       retval(i,j) = gammainc (x(i,j), a);
  890.  
  891.   return retval;
  892. }
  893.  
  894. Matrix
  895. gammainc (const Matrix& x, const Matrix& a)
  896. {
  897.   Matrix retval;
  898.  
  899.   int nr = x.rows ();
  900.   int nc = x.cols ();
  901.  
  902.   int a_nr = a.rows ();
  903.   int a_nc = a.cols ();
  904.  
  905.   if (nr == a_nr && nc == a_nc)
  906.     {
  907.       retval.resize (nr, nc);
  908.  
  909.       for (int j = 0; j < nc; j++)
  910.     for (int i = 0; i < nr; i++)
  911.       retval(i,j) = gammainc (x(i,j), a(i,j));
  912.     }
  913.   else
  914.     (*current_liboctave_error_handler)
  915.       ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
  916.        nr, nc, a_nr, a_nc);
  917.  
  918.   return retval;
  919. }
  920.  
  921. /*
  922. ;;; Local Variables: ***
  923. ;;; mode: C++ ***
  924. ;;; End: ***
  925. */
  926.