home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-base.tgz / octave-1.1.1p1-base.tar / fsf / octave / src / xpow.cc < prev    next >
C/C++ Source or Header  |  1995-01-03  |  15KB  |  704 lines
  1. // xpow.cc                                               -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  
  22. */
  23.  
  24. #ifdef HAVE_CONFIG_H
  25. #include "config.h"
  26. #endif
  27.  
  28. #include <assert.h>
  29. #include <Complex.h>
  30.  
  31. #include "xpow.h"
  32. #include "dMatrix.h"
  33. #include "CMatrix.h"
  34. #include "dDiagMatrix.h"
  35. #include "CDiagMatrix.h"
  36. #include "CColVector.h"
  37. #include "EIG.h"
  38. #include "tree-const.h"
  39. #include "error.h"
  40.  
  41. // This function also appears in tree-const.cc.  Maybe it should be a
  42. // member function of the Matrix class.
  43.  
  44. static int
  45. any_element_is_negative (const Matrix& a)
  46. {
  47.   int nr = a.rows ();
  48.   int nc = a.columns ();
  49.   for (int j = 0; j < nc; j++)
  50.     for (int i = 0; i < nr; i++)
  51.       if (a.elem (i, j) < 0.0)
  52.     return 1;
  53.   return 0;
  54. }
  55.  
  56. // Safer pow functions.
  57. //
  58. //       op2 \ op1:   s   m   cs   cm
  59. //            +--   +---+---+----+----+
  60. //   scalar   |     | 1 | 5 |  7 | 11 |
  61. //                  +---+---+----+----+
  62. //   matrix         | 2 | E |  8 |  E |
  63. //                  +---+---+----+----+
  64. //   complex_scalar | 3 | 6 |  9 | 12 |
  65. //                  +---+---+----+----+
  66. //   complex_matrix | 4 | E | 10 |  E |
  67. //                  +---+---+----+----+
  68. //
  69. //   E -> error, trapped in arith-ops.cc.
  70.  
  71. // -*- 1 -*-
  72. tree_constant
  73. xpow (double a, double b)
  74. {
  75.   if (a < 0.0 && (int) b != b)
  76.     {
  77.       Complex atmp (a);
  78.       return tree_constant (pow (atmp, b));
  79.     }
  80.   else
  81.     return tree_constant (pow (a, b));
  82. }
  83.  
  84. // -*- 2 -*-
  85. tree_constant
  86. xpow (double a, const Matrix& b)
  87. {
  88.   tree_constant retval;
  89.  
  90.   int nr = b.rows ();
  91.   int nc = b.columns ();
  92.  
  93.   if (nr == 0 || nc == 0 || nr != nc)
  94.     error ("for x^A, A must be square");
  95.   else
  96.     {
  97.       EIG b_eig (b);
  98.       ComplexColumnVector lambda (b_eig.eigenvalues ());
  99.       ComplexMatrix Q (b_eig.eigenvectors ());
  100.  
  101.       for (int i = 0; i < nr; i++)
  102.     {
  103.       Complex elt = lambda.elem (i);
  104.       if (imag (elt) == 0.0)
  105.         lambda.elem (i) = pow (a, real (elt));
  106.       else
  107.         lambda.elem (i) = pow (a, elt);
  108.     }
  109.       ComplexDiagMatrix D (lambda);
  110.  
  111.       ComplexMatrix result = Q * D * Q.inverse ();
  112.       retval = tree_constant (result);
  113.     }
  114.  
  115.   return retval;
  116. }
  117.  
  118. // -*- 3 -*-
  119. tree_constant
  120. xpow (double a, const Complex& b)
  121. {
  122.   Complex result;
  123.   Complex atmp (a);
  124.   result = pow (atmp, b);
  125.   return tree_constant (result);
  126. }
  127.  
  128. // -*- 4 -*-
  129. tree_constant
  130. xpow (double a, const ComplexMatrix& b)
  131. {
  132.   tree_constant retval;
  133.  
  134.   int nr = b.rows ();
  135.   int nc = b.columns ();
  136.  
  137.   if (nr == 0 || nc == 0 || nr != nc)
  138.     error ("for x^A, A must be square");
  139.   else
  140.     {
  141.       EIG b_eig (b);
  142.       ComplexColumnVector lambda (b_eig.eigenvalues ());
  143.       ComplexMatrix Q (b_eig.eigenvectors ());
  144.  
  145.       for (int i = 0; i < nr; i++)
  146.     {
  147.       Complex elt = lambda.elem (i);
  148.       if (imag (elt) == 0.0)
  149.         lambda.elem (i) = pow (a, real (elt));
  150.       else
  151.         lambda.elem (i) = pow (a, elt);
  152.     }
  153.       ComplexDiagMatrix D (lambda);
  154.  
  155.       ComplexMatrix result = Q * D * Q.inverse ();
  156.       retval = tree_constant (result);
  157.     }
  158.  
  159.   return retval;
  160. }
  161.  
  162. // -*- 5 -*-
  163. tree_constant
  164. xpow (const Matrix& a, double b)
  165. {
  166.   tree_constant retval;
  167.  
  168.   int nr = a.rows ();
  169.   int nc = a.columns ();
  170.  
  171.   if (nr == 0 || nc == 0 || nr != nc)
  172.     {
  173.       error ("for A^b, A must be square");
  174.       return retval;
  175.     }
  176.  
  177.   if ((int) b == b)
  178.     {
  179.       int btmp = (int) b;
  180.       if (btmp == 0)
  181.     {
  182.       DiagMatrix result (nr, nr, 1.0);
  183.       retval = tree_constant (result);
  184.     }
  185.       else
  186.     {
  187. // Too much copying?
  188. // XXX FIXME XXX -- we shouldn\'t do this if the exponent is large...
  189.       Matrix atmp;
  190.       if (btmp < 0)
  191.         {
  192.           btmp = -btmp;
  193.           atmp = a.inverse ();
  194.         }
  195.       else
  196.         atmp = a;
  197.  
  198.       Matrix result (atmp);
  199.       for (int i = 1; i < btmp; i++)
  200.         result = result * atmp;
  201.  
  202.       retval = tree_constant (result);
  203.     }
  204.     }
  205.   else
  206.     {
  207.       EIG a_eig (a);
  208.       ComplexColumnVector lambda (a_eig.eigenvalues ());
  209.       ComplexMatrix Q (a_eig.eigenvectors ());
  210.  
  211.       for (int i = 0; i < nr; i++)
  212.     lambda.elem (i) = pow (lambda.elem (i), b);
  213.  
  214.       ComplexDiagMatrix D (lambda);
  215.  
  216.       ComplexMatrix result = Q * D * Q.inverse ();
  217.       retval = tree_constant (result);
  218.     }
  219.  
  220.   return retval;
  221. }
  222.  
  223. // -*- 6 -*-
  224. tree_constant
  225. xpow (const Matrix& a, const Complex& b)
  226. {
  227.   int nr = a.rows ();
  228.   int nc = a.columns ();
  229.  
  230.   if (nr == 0 || nc == 0 || nr != nc)
  231.     {
  232.       error ("for A^b, A must be square");
  233.       return tree_constant ();
  234.     }
  235.  
  236.   EIG a_eig (a);
  237.   ComplexColumnVector lambda (a_eig.eigenvalues ());
  238.   ComplexMatrix Q (a_eig.eigenvectors ());
  239.  
  240.   for (int i = 0; i < nr; i++)
  241.     lambda.elem (i) = pow (lambda.elem (i), b);
  242.  
  243.   ComplexDiagMatrix D (lambda);
  244.  
  245.   ComplexMatrix result = Q * D * Q.inverse ();
  246.  
  247.   return tree_constant (result);
  248. }
  249.  
  250. // -*- 7 -*-
  251. tree_constant
  252. xpow (const Complex& a, double b)
  253. {
  254.   Complex result;
  255.   result = pow (a, b);
  256.   return tree_constant (result);
  257. }
  258.  
  259. // -*- 8 -*-
  260. tree_constant
  261. xpow (const Complex& a, const Matrix& b)
  262. {
  263.   tree_constant retval;
  264.  
  265.   int nr = b.rows ();
  266.   int nc = b.columns ();
  267.  
  268.   if (nr == 0 || nc == 0 || nr != nc)
  269.     {
  270.       error ("for x^A, A must be square");
  271.     }
  272.   else
  273.     {
  274.       EIG b_eig (b);
  275.       ComplexColumnVector lambda (b_eig.eigenvalues ());
  276.       ComplexMatrix Q (b_eig.eigenvectors ());
  277.  
  278.       for (int i = 0; i < nr; i++)
  279.     {
  280.       Complex elt = lambda.elem (i);
  281.       if (imag (elt) == 0.0)
  282.         lambda.elem (i) = pow (a, real (elt));
  283.       else
  284.         lambda.elem (i) = pow (a, elt);
  285.     }
  286.       ComplexDiagMatrix D (lambda);
  287.  
  288.       ComplexMatrix result = Q * D * Q.inverse ();
  289.       retval = tree_constant (result);
  290.     }
  291.  
  292.   return retval;
  293. }
  294.  
  295. // -*- 9 -*-
  296. tree_constant
  297. xpow (const Complex& a, const Complex& b)
  298. {
  299.   Complex result;
  300.   result = pow (a, b);
  301.   return tree_constant (result);
  302. }
  303.  
  304. // -*- 10 -*-
  305. tree_constant
  306. xpow (const Complex& a, const ComplexMatrix& b)
  307. {
  308.   tree_constant retval;
  309.  
  310.   int nr = b.rows ();
  311.   int nc = b.columns ();
  312.  
  313.   if (nr == 0 || nc == 0 || nr != nc)
  314.     error ("for x^A, A must be square");
  315.   else
  316.     {
  317.       EIG b_eig (b);
  318.       ComplexColumnVector lambda (b_eig.eigenvalues ());
  319.       ComplexMatrix Q (b_eig.eigenvectors ());
  320.  
  321.       for (int i = 0; i < nr; i++)
  322.     {
  323.       Complex elt = lambda.elem (i);
  324.       if (imag (elt) == 0.0)
  325.         lambda.elem (i) = pow (a, real (elt));
  326.       else
  327.         lambda.elem (i) = pow (a, elt);
  328.     }
  329.       ComplexDiagMatrix D (lambda);
  330.  
  331.       ComplexMatrix result = Q * D * Q.inverse ();
  332.       retval = tree_constant (result);
  333.     }
  334.  
  335.   return retval;
  336. }
  337.  
  338. // -*- 11 -*-
  339. tree_constant
  340. xpow (const ComplexMatrix& a, double b)
  341. {
  342.   tree_constant retval;
  343.  
  344.   int nr = a.rows ();
  345.   int nc = a.columns ();
  346.  
  347.   if (nr == 0 || nc == 0 || nr != nc)
  348.     {
  349.       error ("for A^b, A must be square");
  350.       return retval;
  351.     }
  352.  
  353.   if ((int) b == b)
  354.     {
  355.       int btmp = (int) b;
  356.       if (btmp == 0)
  357.     {
  358.       DiagMatrix result (nr, nr, 1.0);
  359.       retval = tree_constant (result);
  360.     }
  361.       else
  362.     {
  363. // Too much copying?
  364. // XXX FIXME XXX -- we shouldn\'t do this if the exponent is large...
  365.       ComplexMatrix atmp;
  366.       if (btmp < 0)
  367.         {
  368.           btmp = -btmp;
  369.           atmp = a.inverse ();
  370.         }
  371.       else
  372.         atmp = a;
  373.  
  374.       ComplexMatrix result (atmp);
  375.       for (int i = 1; i < btmp; i++)
  376.         result = result * atmp;
  377.  
  378.       retval = tree_constant (result);
  379.     }
  380.     }
  381.   else
  382.     {
  383.       EIG a_eig (a);
  384.       ComplexColumnVector lambda (a_eig.eigenvalues ());
  385.       ComplexMatrix Q (a_eig.eigenvectors ());
  386.  
  387.       for (int i = 0; i < nr; i++)
  388.     lambda.elem (i) = pow (lambda.elem (i), b);
  389.  
  390.       ComplexDiagMatrix D (lambda);
  391.  
  392.       ComplexMatrix result = Q * D * Q.inverse ();
  393.       retval = tree_constant (result);
  394.     }
  395.  
  396.   return retval;
  397. }
  398.  
  399. // -*- 12 -*-
  400. tree_constant
  401. xpow (const ComplexMatrix& a, const Complex& b)
  402. {
  403.   int nr = a.rows ();
  404.   int nc = a.columns ();
  405.  
  406.   if (nr == 0 || nc == 0 || nr != nc)
  407.     {
  408.       error ("for A^b, A must be square");
  409.       return tree_constant ();
  410.     }
  411.  
  412.   EIG a_eig (a);
  413.   ComplexColumnVector lambda (a_eig.eigenvalues ());
  414.   ComplexMatrix Q (a_eig.eigenvectors ());
  415.  
  416.   for (int i = 0; i < nr; i++)
  417.     lambda.elem (i) = pow (lambda.elem (i), b);
  418.  
  419.   ComplexDiagMatrix D (lambda);
  420.  
  421.   ComplexMatrix result = Q * D * Q.inverse ();
  422.  
  423.   return tree_constant (result);
  424. }
  425.  
  426. // Safer pow functions that work elementwise for matrices.
  427. //
  428. //       op2 \ op1:   s   m   cs   cm
  429. //            +--   +---+---+----+----+
  430. //   scalar   |     | * | 3 |  * |  9 |
  431. //                  +---+---+----+----+
  432. //   matrix         | 1 | 4 |  7 | 10 |
  433. //                  +---+---+----+----+
  434. //   complex_scalar | * | 5 |  * | 11 |
  435. //                  +---+---+----+----+
  436. //   complex_matrix | 2 | 6 |  8 | 12 |
  437. //                  +---+---+----+----+
  438. //
  439. //   * -> not needed.
  440.  
  441. // -*- 1 -*-
  442. tree_constant
  443. elem_xpow (double a, const Matrix& b)
  444. {
  445.   tree_constant retval;
  446.  
  447.   int nr = b.rows ();
  448.   int nc = b.columns ();
  449.  
  450. // For now, assume the worst.
  451.   if (a < 0.0)
  452.     {
  453.       Complex atmp (a);
  454.       ComplexMatrix result (nr, nc);
  455.       for (int j = 0; j < nc; j++)
  456.     for (int i = 0; i < nr; i++)
  457.       result.elem (i, j) = pow (atmp, b.elem (i, j));
  458.  
  459.       retval = tree_constant (result);
  460.     }
  461.   else
  462.     {
  463.       Matrix result (nr, nc);
  464.       for (int j = 0; j < nc; j++)
  465.     for (int i = 0; i < nr; i++)
  466.       result.elem (i, j) = pow (a, b.elem (i, j)); 
  467.  
  468.       retval = tree_constant (result);
  469.     }
  470.  
  471.   return retval;
  472. }
  473.  
  474. // -*- 2 -*-
  475. tree_constant
  476. elem_xpow (double a, const ComplexMatrix& b)
  477. {
  478.   int nr = b.rows ();
  479.   int nc = b.columns ();
  480.  
  481.   ComplexMatrix result (nr, nc);
  482.   for (int j = 0; j < nc; j++)
  483.     for (int i = 0; i < nr; i++)
  484.       result.elem (i, j) = pow (a, b.elem (i, j));
  485.  
  486.   return tree_constant (result);
  487. }
  488.  
  489. // -*- 3 -*-
  490. tree_constant
  491. elem_xpow (const Matrix& a, double b)
  492. {
  493.   tree_constant retval;
  494.  
  495.   int nr = a.rows ();
  496.   int nc = a.columns ();
  497.  
  498.   if ((int) b != b && any_element_is_negative (a))
  499.     {
  500.       ComplexMatrix result (nr, nc);
  501.       for (int j = 0; j < nc; j++)
  502.     for (int i = 0; i < nr; i++)
  503.       {
  504.         Complex atmp (a.elem (i, j));
  505.         result.elem (i, j) = pow (atmp, b);
  506.       }
  507.  
  508.       retval = tree_constant (result);
  509.     }
  510.   else
  511.     {
  512.       Matrix result (nr, nc);
  513.       for (int j = 0; j < nc; j++)
  514.     for (int i = 0; i < nr; i++)
  515.       result.elem (i, j) = pow (a.elem (i, j), b);
  516.  
  517.       retval = tree_constant (result);
  518.     }
  519.  
  520.   return retval;
  521. }
  522.  
  523. // -*- 4 -*-
  524. tree_constant
  525. elem_xpow (const Matrix& a, const Matrix& b)
  526. {
  527.   int nr = a.rows ();
  528.   int nc = a.columns ();
  529.  
  530.   assert (nr == b.rows () && nc == b.columns ());
  531.  
  532.   int convert_to_complex = 0;
  533.   int i;
  534.   for (int j = 0; j < nc; j++)
  535.     for (i = 0; i < nr; i++)
  536.       {
  537.     double atmp = a.elem (i, j);
  538.     double btmp = b.elem (i, j);
  539.     if (atmp < 0.0 && (int) btmp != btmp)
  540.       {
  541.         convert_to_complex = 1;
  542.         goto done;
  543.       }
  544.       }
  545.  
  546.  done:
  547.  
  548.   if (convert_to_complex)
  549.     {
  550.       ComplexMatrix complex_result (nr, nc);
  551.  
  552.       for (j = 0; j < nc; j++)
  553.     for (i = 0; i < nr; i++)
  554.       {
  555.         Complex atmp (a.elem (i, j));
  556.         Complex btmp (b.elem (i, j));
  557.         complex_result.elem (i, j) = pow (atmp, btmp);
  558.       }
  559.       return tree_constant (complex_result);
  560.     }
  561.   else
  562.     {
  563.       Matrix result (nr, nc);
  564.  
  565.       for (j = 0; j < nc; j++)
  566.     for (i = 0; i < nr; i++)
  567.       result.elem (i, j) = pow (a.elem (i, j), b.elem (i, j));
  568.  
  569.       return tree_constant (result);
  570.     }
  571. }
  572.  
  573. // -*- 5 -*-
  574. tree_constant
  575. elem_xpow (const Matrix& a, const Complex& b)
  576. {
  577.   int nr = a.rows ();
  578.   int nc = a.columns ();
  579.  
  580.   ComplexMatrix result (nr, nc);
  581.   for (int j = 0; j < nc; j++)
  582.     for (int i = 0; i < nr; i++)
  583.       result.elem (i, j) = pow (a.elem (i, j), b);
  584.  
  585.   return tree_constant (result);
  586. }
  587.  
  588. // -*- 6 -*-
  589. tree_constant
  590. elem_xpow (const Matrix& a, const ComplexMatrix& b)
  591. {
  592.   int nr = a.rows ();
  593.   int nc = a.columns ();
  594.  
  595.   assert (nr == b.rows () && nc == b.columns ());
  596.  
  597.   ComplexMatrix result (nr, nc);
  598.   for (int j = 0; j < nc; j++)
  599.     for (int i = 0; i < nr; i++)
  600.       result.elem (i, j) = pow (a.elem (i, j), b.elem (i, j));
  601.  
  602.   return tree_constant (result);
  603. }
  604.  
  605. // -*- 7 -*-
  606. tree_constant
  607. elem_xpow (const Complex& a, const Matrix& b)
  608. {
  609.   int nr = b.rows ();
  610.   int nc = b.columns ();
  611.  
  612.   ComplexMatrix result (nr, nc);
  613.   for (int j = 0; j < nc; j++)
  614.     for (int i = 0; i < nr; i++)
  615.       result.elem (i, j) = pow (a, b.elem (i, j));
  616.  
  617.   return tree_constant (result);
  618. }
  619.  
  620. // -*- 8 -*-
  621. tree_constant
  622. elem_xpow (const Complex& a, const ComplexMatrix& b)
  623. {
  624.   int nr = b.rows ();
  625.   int nc = b.columns ();
  626.  
  627.   ComplexMatrix result (nr, nc);
  628.   for (int j = 0; j < nc; j++)
  629.     for (int i = 0; i < nr; i++)
  630.       result.elem (i, j) = pow (a, b.elem (i, j));
  631.  
  632.   return tree_constant (result);
  633. }
  634.  
  635. // -*- 9 -*-
  636. tree_constant
  637. elem_xpow (const ComplexMatrix& a, double b)
  638. {
  639.   int nr = a.rows ();
  640.   int nc = a.columns ();
  641.  
  642.   ComplexMatrix result (nr, nc);
  643.   for (int j = 0; j < nc; j++)
  644.     for (int i = 0; i < nr; i++)
  645.       result.elem (i, j) = pow (a.elem (i, j), b);
  646.  
  647.   return tree_constant (result);
  648. }
  649.  
  650. // -*- 10 -*-
  651. tree_constant
  652. elem_xpow (const ComplexMatrix& a, const Matrix& b)
  653. {
  654.   int nr = a.rows ();
  655.   int nc = a.columns ();
  656.  
  657.   assert (nr == b.rows () && nc == b.columns ());
  658.  
  659.   ComplexMatrix result (nr, nc);
  660.   for (int j = 0; j < nc; j++)
  661.     for (int i = 0; i < nr; i++)
  662.       result.elem (i, j) = pow (a.elem (i, j), b.elem (i, j));
  663.  
  664.   return tree_constant (result);
  665. }
  666.  
  667. // -*- 11 -*-
  668. tree_constant
  669. elem_xpow (const ComplexMatrix& a, const Complex& b)
  670. {
  671.   int nr = a.rows ();
  672.   int nc = a.columns ();
  673.  
  674.   ComplexMatrix result (nr, nc);
  675.   for (int j = 0; j < nc; j++)
  676.     for (int i = 0; i < nr; i++)
  677.       result.elem (i, j) = pow (a.elem (i, j), b);
  678.  
  679.   return tree_constant (result);
  680. }
  681.  
  682. // -*- 12 -*-
  683. tree_constant
  684. elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b)
  685. {
  686.   int nr = a.rows ();
  687.   int nc = a.columns ();
  688.  
  689.   ComplexMatrix result (nr, nc);
  690.  
  691.   for (int j = 0; j < nc; j++)
  692.     for (int i = 0; i < nr; i++)
  693.       result.elem (i, j) = pow (a.elem (i, j), b.elem (i, j));
  694.  
  695.   return tree_constant (result);
  696. }
  697.  
  698. /*
  699. ;;; Local Variables: ***
  700. ;;; mode: C++ ***
  701. ;;; page-delimiter: "^/\\*" ***
  702. ;;; End: ***
  703. */
  704.