home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATMATH.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  13.8 KB  |  540 lines
  1. /**************************************************/
  2. /*  matmath.c source for transcendental functions */
  3. /**************************************************/
  4.  
  5.  
  6. /**************************************************/
  7. /*            MatClass Source File                */
  8. /*       Copyright of C. R. Birchenhall           */
  9. /*       University of Manchester, UK.            */
  10. /*   MatClass is freeware. This file should be    */
  11. /* made freely available to users of any software */
  12. /* whose creation is wholly or partly dependent   */
  13. /*                on this file.                   */
  14. /**************************************************/
  15.  
  16. #include "matmath.hpp"
  17.  
  18. #include <math.h>
  19.  
  20. /***********************************************/
  21. /*            Transcendental Functions         */
  22. /***********************************************/
  23.  
  24. enum transcendental {
  25.    LOG, LOG10, EXP, SIN, COS, TAN, ASIN, ACOS, ATAN,
  26.    SINH, COSH, TANH
  27. } ; // transcendental
  28.  
  29. static char* transcNames[] = {
  30.    "Log", "Log10", "Exp",
  31.    "Sin", "Cos", "Tan",
  32.    "ArcSin", "ArcCos", "ArcTan",
  33.    "Sinh", "Cosh", "Tanh"
  34. } ;
  35.  
  36. void transcend( matrix& z, const matrix& x, int fn )
  37. {
  38.    static char *mName = "transcend" ;
  39.    matFunc func( mName ) ; x.debugInfo( func ) ;
  40.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  41.    z.reset( nr, nc ) ;
  42.    int error ;
  43.    REAL eps = matrixEps(), maxExp = -log( eps ) ;
  44.    // do tests ---- these need more work !!!
  45.    switch ( fn ) {
  46.       case LOG :
  47.       case LOG10 :
  48.          error = x.compare( eps ) != GREATER ;
  49.          break ;
  50.       case EXP :
  51.       case SINH :
  52.       case COSH :
  53.       case TANH :
  54.          error = x.compare( maxExp ) != LESS || 
  55.                  x.compare( -maxExp ) != GREATER ;
  56.          break ;
  57.       case ASIN :
  58.       case ACOS :
  59.          error = x.compare( 1.0 ) & GREATER  ||
  60.                  x.compare( -1.0 ) & LESS ;
  61.          break ;
  62.       case TAN :     // ????
  63.       default :
  64.          error = FALSE ;
  65.          break ;
  66.    } // switch
  67.    if ( error )
  68.       x.errorExit( transcNames[ fn ], TRANS ) ;
  69.    /*****************************/
  70.    /* Loop can be optimised if  */
  71.    /* willing to exploit col    */
  72.    /* major form of matrices    */
  73.    /*****************************/
  74.    for ( j = 1 ; j <= nc ; j++ ) {
  75.       switch ( fn ) {
  76.          case LOG :
  77.             for ( i = 1 ; i <= nr ; i++ )
  78.                z(i,j) = (REAL) log( x(i,j) ) ;
  79.             break ;
  80.          case LOG10 :
  81.             for ( i = 1 ; i <= nr ; i++ )
  82.                z(i,j) = (REAL) log10( x(i,j) ) ;
  83.             break ;
  84.          case EXP :
  85.             for ( i = 1 ; i <= nr ; i++ )
  86.                z(i,j) = (REAL) exp( x(i,j) ) ;
  87.             break ;
  88.          case SIN :
  89.             for ( i = 1 ; i <= nr ; i++ )
  90.                z(i,j) = (REAL) sin( x(i,j) ) ;
  91.             break ;
  92.          case COS :
  93.             for ( i = 1 ; i <= nr ; i++ )
  94.                z(i,j) = (REAL) cos( x(i,j) ) ;
  95.             break ;
  96.          case TAN :
  97.             for ( i = 1 ; i <= nr ; i++ )
  98.                z(i,j) = (REAL) tan( x(i,j) ) ;
  99.             break ;
  100.          case ASIN :
  101.             for ( i = 1 ; i <= nr ; i++ )
  102.                z(i,j) = (REAL) asin( x(i,j) ) ;
  103.             break ;
  104.          case ACOS :
  105.             for ( i = 1 ; i <= nr ; i++ )
  106.                z(i,j) = (REAL) acos( x(i,j) ) ;
  107.             break ;
  108.          case ATAN :
  109.             for ( i = 1 ; i <= nr ; i++ )
  110.                z(i,j) = (REAL) atan( x(i,j) ) ;
  111.             break ;
  112.          case SINH :
  113.             for ( i = 1 ; i <= nr ; i++ )
  114.                z(i,j) = (REAL) sinh( x(i,j) ) ;
  115.             break ;
  116.          case COSH :
  117.             for ( i = 1 ; i <= nr ; i++ )
  118.                z(i,j) = (REAL) cosh( x(i,j) ) ;
  119.             break ;
  120.          case TANH :
  121.             for ( i = 1 ; i <= nr ; i++ )
  122.                z(i,j) = (REAL) tanh( x(i,j) ) ;
  123.             break ;
  124.       } // switch
  125.    } // for i
  126.    return ;
  127. } // matrix transcend
  128.  
  129. #ifdef __cplusplus
  130.  
  131. matrix abs( matrix& x )
  132. {
  133.    static char *mName = "abs" ;
  134.    matFunc func( mName ) ; x.debugInfo( func ) ;
  135.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  136.    matrix z( nr, nc ) ;
  137.    for ( i = 1 ; i <= nr ; i++ )
  138.       for ( j = 1 ; j <= nc ; j++ )
  139.          z(i,j) = (REAL) fabs( x(i,j) ) ;
  140.    return z ;
  141. } // matrix abs
  142.  
  143.  
  144. matrix log( const matrix& x )
  145. { 
  146.    static char *mName = "log" ;
  147.    matFunc func( mName ) ; x.debugInfo( func ) ;
  148.    matrix z( x.nRows(), x.nCols() ) ;
  149.    transcend( z, x, LOG ) ;
  150.    return z ;
  151. } // log
  152.  
  153. matrix exp( const matrix& x )
  154. { 
  155.    static char *mName = "exp" ;
  156.    matFunc func( mName ) ; x.debugInfo( func ) ;
  157.    matrix z( x.nRows(), x.nCols() ) ;
  158.    transcend( z, x, EXP ) ;
  159.    return z ;
  160. } // exp
  161.  
  162. matrix log10( const matrix& x )
  163. { 
  164.    static char *mName = "log10" ;
  165.    matFunc func( mName ) ; x.debugInfo( func ) ;
  166.    matrix z( x.nRows(), x.nCols() ) ;
  167.    transcend( z, x, LOG10 ) ;
  168.    return z ;
  169. } // log10m
  170.  
  171. matrix ln( const matrix& x )
  172. { 
  173.    static char *mName = "ln" ;
  174.    matFunc func( mName ) ; x.debugInfo( func ) ;
  175.    matrix z( x.nRows(), x.nCols() ) ;
  176.    transcend( z, x, LOG ) ;
  177.    return z ;
  178. } // ln
  179.  
  180. matrix sin( const matrix& x )
  181. { 
  182.    static char *mName = "sin" ;
  183.    matFunc func( mName ) ; x.debugInfo( func ) ;
  184.    matrix z( x.nRows(), x.nCols() ) ;
  185.    transcend( z, x, SIN ) ;
  186.    return z ;
  187. } //sin
  188.  
  189. matrix cos( const matrix& x )
  190. { 
  191.    static char *mName = "cos" ;
  192.    matFunc func( mName ) ; x.debugInfo( func ) ;
  193.    matrix z( x.nRows(), x.nCols() ) ;
  194.    transcend( z, x, COS ) ;
  195.    return z ;
  196. } // cos
  197.  
  198. matrix tan( const matrix& x )
  199. { 
  200.    static char *mName = "tan" ;
  201.    matFunc func( mName ) ; x.debugInfo( func ) ;
  202.    matrix z( x.nRows(), x.nCols() ) ;
  203.    transcend( z, x, TAN ) ;
  204.    return z ;
  205. } // tan
  206.  
  207. matrix asin( const matrix& x )
  208. { 
  209.    static char *mName = "asin" ;
  210.    matFunc func( mName ) ; x.debugInfo( func ) ;
  211.    matrix z( x.nRows(), x.nCols() ) ;
  212.    transcend( z, x, ASIN ) ;
  213.    return z ;
  214. } // asinm
  215.  
  216. matrix acos( const matrix& x )
  217. { 
  218.    static char *mName = "acos" ;
  219.    matFunc func( mName ) ; x.debugInfo( func ) ;
  220.    matrix z( x.nRows(), x.nCols() ) ;
  221.    transcend( z, x, ACOS ) ;
  222.    return z ;
  223. } // acos
  224.  
  225. matrix atan( const matrix& x )
  226. { 
  227.    static char *mName = "atan" ;
  228.    matFunc func( mName ) ; x.debugInfo( func ) ;
  229.    matrix z( x.nRows(), x.nCols() ) ;
  230.    transcend( z, x, ATAN ) ;
  231.    return z ;
  232. } // atan
  233.  
  234. matrix sinh( const matrix& x )
  235. { 
  236.    static char *mName = "sinh" ;
  237.    matFunc func( mName ) ; x.debugInfo( func ) ;
  238.    matrix z( x.nRows(), x.nCols() ) ;
  239.    transcend( z, x, SINH ) ;
  240.    return z ;
  241. } // sinh
  242.  
  243. matrix cosh( const matrix& x )
  244. { 
  245.    static char *mName = "cosh" ;
  246.    matFunc func( mName ) ; x.debugInfo( func ) ;
  247.    matrix z( x.nRows(), x.nCols() ) ;
  248.    transcend( z, x, COSH ) ;
  249.    return z ;
  250. } // cosh
  251.  
  252. matrix tanh( const matrix& x )
  253. { 
  254.    static char *mName = "tanh" ;
  255.    matFunc func( mName ) ; x.debugInfo( func ) ;
  256.    matrix z( x.nRows(), x.nCols() ) ;
  257.    transcend( z, x, TANH ) ;
  258.    return z ;
  259. } // tanh
  260.  
  261. #else
  262.  
  263. matrix absm( matrix& x )
  264. {
  265.    static char *mName = "absm" ;
  266.    matFunc func( mName ) ; x.debugInfo( func ) ;
  267.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  268.    matrix z( nr, nc ) ;
  269.    for ( i = 1 ; i <= nr ; i++ )
  270.       for ( j = 1 ; j <= nc ; j++ )
  271.          z(i,j) = (REAL) fabs( x(i,j) ) ;
  272.    return z ;
  273. } // matrix absm
  274.  
  275.  
  276. matrix logm( const matrix& x )
  277. { 
  278.    static char *mName = "log" ;
  279.    matFunc func( mName ) ; x.debugInfo( func ) ;
  280.    matrix z( x.nRows(), x.nCols() ) ;
  281.    transcend( z, x, LOG ) ;
  282.    return z ;
  283. } // log
  284.  
  285. matrix expm( const matrix& x )
  286. { 
  287.    static char *mName = "exp" ;
  288.    matFunc func( mName ) ; x.debugInfo( func ) ;
  289.    matrix z( x.nRows(), x.nCols() ) ;
  290.    transcend( z, x, EXP ) ;
  291.    return z ;
  292. } // exp
  293.  
  294. matrix log10m( const matrix& x )
  295. { 
  296.    static char *mName = "log10" ;
  297.    matFunc func( mName ) ; x.debugInfo( func ) ;
  298.    matrix z( x.nRows(), x.nCols() ) ;
  299.    transcend( z, x, LOG10 ) ;
  300.    return z ;
  301. } // log10m
  302.  
  303. matrix lnm( const matrix& x )
  304. { 
  305.    static char *mName = "ln" ;
  306.    matFunc func( mName ) ; x.debugInfo( func ) ;
  307.    matrix z( x.nRows(), x.nCols() ) ;
  308.    transcend( z, x, LOG ) ;
  309.    return z ;
  310. } // ln
  311.  
  312. matrix sinm( const matrix& x )
  313. { 
  314.    static char *mName = "sin" ;
  315.    matFunc func( mName ) ; x.debugInfo( func ) ;
  316.    matrix z( x.nRows(), x.nCols() ) ;
  317.    transcend( z, x, SIN ) ;
  318.    return z ;
  319. } //sin
  320.  
  321. matrix cosm( const matrix& x )
  322. { 
  323.    static char *mName = "cos" ;
  324.    matFunc func( mName ) ; x.debugInfo( func ) ;
  325.    matrix z( x.nRows(), x.nCols() ) ;
  326.    transcend( z, x, COS ) ;
  327.    return z ;
  328. } // cos
  329.  
  330. matrix tanm( const matrix& x )
  331. { 
  332.    static char *mName = "tan" ;
  333.    matFunc func( mName ) ; x.debugInfo( func ) ;
  334.    matrix z( x.nRows(), x.nCols() ) ;
  335.    transcend( z, x, TAN ) ;
  336.    return z ;
  337. } // tan
  338.  
  339. matrix asinm( const matrix& x )
  340. { 
  341.    static char *mName = "asin" ;
  342.    matFunc func( mName ) ; x.debugInfo( func ) ;
  343.    matrix z( x.nRows(), x.nCols() ) ;
  344.    transcend( z, x, ASIN ) ;
  345.    return z ;
  346. } // asinm
  347.  
  348. matrix acosm( const matrix& x )
  349. { 
  350.    static char *mName = "acos" ;
  351.    matFunc func( mName ) ; x.debugInfo( func ) ;
  352.    matrix z( x.nRows(), x.nCols() ) ;
  353.    transcend( z, x, ACOS ) ;
  354.    return z ;
  355. } // acos
  356.  
  357. matrix atanm( const matrix& x )
  358. { 
  359.    static char *mName = "atan" ;
  360.    matFunc func( mName ) ; x.debugInfo( func ) ;
  361.    matrix z( x.nRows(), x.nCols() ) ;
  362.    transcend( z, x, ATAN ) ;
  363.    return z ;
  364. } // atan
  365.  
  366. matrix sinhm( const matrix& x )
  367. { 
  368.    static char *mName = "sinh" ;
  369.    matFunc func( mName ) ; x.debugInfo( func ) ;
  370.    matrix z( x.nRows(), x.nCols() ) ;
  371.    transcend( z, x, SINH ) ;
  372.    return z ;
  373. } // sinh
  374.  
  375. matrix coshm( const matrix& x )
  376. { 
  377.    static char *mName = "cosh" ;
  378.    matFunc func( mName ) ; x.debugInfo( func ) ;
  379.    matrix z( x.nRows(), x.nCols() ) ;
  380.    transcend( z, x, COSH ) ;
  381.    return z ;
  382. } // cosh
  383.  
  384. matrix tanhm( const matrix& x )
  385. { 
  386.    static char *mName = "tanh" ;
  387.    matFunc func( mName ) ; x.debugInfo( func ) ;
  388.    matrix z( x.nRows(), x.nCols() ) ;
  389.    transcend( z, x, TANH ) ;
  390.    return z ;
  391. } // tanh
  392.  
  393. #endif
  394.  
  395.  
  396. matrix& matrix::roundOf( const matrix& x, truncation option )
  397. {
  398.    static char *mName = "roundOf" ;
  399.    matFunc func( mName ) ; debugInfo( func ) ;
  400.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  401.    if ( x.isNull() )
  402.       x.errorExit( mName, NULRF ) ;
  403.    reset( nr, nc ) ;
  404.    DOUBLE d, tmp ;
  405.    for ( i = 1 ; i <= nr ; i++ ) {
  406.       for ( j = 1 ; j <= nc ; j++ ) {
  407.          d = (DOUBLE) x(i,j) ;
  408.          switch ( option ) {
  409.             case ROUND :
  410.                tmp = floor( d ) ;
  411.                d = ( d - tmp > 0.5 ) ? tmp + 1.0 : tmp ;
  412.                break ;
  413.             case FIX :
  414.                d = ( d < 0 ) ? ceil( d ) : floor( d ) ;
  415.                break ;
  416.             case FLOOR :
  417.                d = floor( d ) ;
  418.                break ;
  419.             case CEIL :
  420.                d = ceil( d ) ;
  421.                break ;
  422.          } // switch
  423.          mat(i,j) = (REAL) d ;
  424.       } // for j
  425.    } // for i
  426.    return *this ;
  427. } // matrix roundOf
  428.  
  429. #ifdef __cplusplus
  430. matrix sign( const matrix& x )
  431. #else
  432. matrix signm( const matrix& x )
  433. #endif
  434. {
  435.    static char *mName = "signm(matrix)" ;
  436.    matFunc func( mName ) ; x.debugInfo( func ) ;
  437.    matrix z( x.nRows(), x.nCols() ) ;
  438.    z.signOf( x ) ;
  439.    return z ;
  440. } // sign
  441.  
  442. #ifdef __cplusplus
  443. matrix round( const matrix& x )
  444. #else
  445. matrix roundm( const matrix& x )
  446. #endif
  447. {
  448.    static char *mName = "roundm(matrix)" ;
  449.    matFunc func( mName ) ; x.debugInfo( func ) ;
  450.    matrix z( x.nRows(), x.nCols() ) ;
  451.    z.roundOf( x, ROUND ) ;
  452.    return z ;
  453. } // round
  454.  
  455. #ifdef __cplusplus
  456. matrix fix( const matrix& x )
  457. #else
  458. matrix fixm( const matrix& x )
  459. #endif
  460. {
  461.    static char *mName = "fixm(matrix)" ;
  462.    matFunc func( mName ) ; x.debugInfo( func ) ;
  463.    matrix z( x.nRows(), x.nCols() ) ;
  464.    z.roundOf( x, FIX ) ;
  465.    return z ;
  466. } //  fix
  467.  
  468. #ifdef __cplusplus
  469. matrix floor( const matrix& x )
  470. #else
  471. matrix floorm( const matrix& x )
  472. #endif
  473. {
  474.    static char *mName = "floorm(matrix)" ;
  475.    matFunc func( mName ) ; x.debugInfo( func ) ;
  476.    matrix z( x.nRows(), x.nCols() ) ;
  477.    z.roundOf( x, FLOOR ) ;
  478.    return z ;
  479. } // floor
  480.  
  481. #ifdef __cplusplus
  482. matrix ceil( const matrix& x )
  483. #else
  484. matrix ceilm( const matrix& x )
  485. #endif
  486. {
  487.    static char *mName = "ceilm(matrix)" ;
  488.    matFunc func( mName ) ; x.debugInfo( func ) ;
  489.    matrix z( x.nRows(), x.nCols() ) ;
  490.    z.roundOf( x, CEIL ) ;
  491.    return z ;
  492. } // ceil
  493.  
  494. matrix& matrix::remOf( const matrix& x, const matrix& y )
  495. {
  496.    char *mName = "remOf" ;
  497.    matFunc func( mName ) ; debugInfo( func ) ;
  498.    x.checkDims( y, mName ) ;
  499.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  500.    reset( nr, nc ) ;
  501.    DOUBLE d ;
  502.    for ( i = 1 ; i <= nr ; i++ ) {
  503.       for ( j = 1 ; j <= nc ; j++ ) {
  504.          if ( y(i,j) == 0 )
  505.             errorExit( mName, ZEROD ) ;
  506.          d = x(i,j) / y(i,j) ;
  507.          d = ( d < 0.0 ) ? ceil( d ) : ( ( d > 0.0 ) ? floor( d ) : d ) ;
  508.          mat(i,j) = x(i,j) - y(i,j) * d ;
  509.       } // for j
  510.    } // for j
  511.    return *this ;
  512. } // matrix remOf
  513.  
  514. matrix matrix::rem( const matrix& y ) M_CONST
  515. {
  516.    char *mName = "rem" ;
  517.    matFunc func( mName ) ; debugInfo( func ) ;
  518.    matrix z( nRows(), nCols() ) ;
  519.    z.remOf( *this, y ) ;
  520.    return z ;
  521. } // matrix::rem
  522.  
  523. matrix& matrix::signOf( const matrix& x )
  524. {
  525.    static char *mName = "signOf" ;
  526.    matFunc func( mName ) ; debugInfo( func ) ;
  527.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  528.    if ( x.isNull() )
  529.       x.errorExit( mName, NULRF ) ;
  530.    reset( nr, nc ) ;
  531.    REAL d ;
  532.    for ( i = 1 ; i <= nr ; i++ ) {
  533.       for ( j = 1 ; j <= nc ; j++ ) {
  534.          d = x(i,j) ;
  535.          mat(i,j) = ( d < 0.0 ) ? -1.0 : ( ( d > 0.0 ) ? 1.0 : 0.0 ) ;
  536.       } // for j
  537.    } // for j
  538.    return *this ;
  539. } // matrix signOf
  540.