home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / ELFUN.DI$ / SINH.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.9 KB  |  63 lines
  1. function y = sinh(xin)
  2. %SINH    Hyperbolic sine.
  3. %    SINH(X) is the hyperbolic sine of the elements of X.
  4.  
  5. %    J.M. Hammond, 6-17-92
  6. %    Copyright (c) 1984-93 by the MathWorks, Inc.
  7.  
  8. %       21 Oct 92 (jmh) -- fixed for matrix inputs
  9.  
  10. % Complex argument
  11. if any(any(imag(xin) ~= 0))        % 21 Oct 92 (jmh)
  12.     x = real(xin);
  13.     y = imag(xin);
  14.     y = sinh(x).*cos(y) + i*cosh(x).*sin(y);
  15.  
  16. % Real argument
  17. else
  18.     % Reference: W. J. Cody and W. Waite, "Software Manual
  19.     % for the Elementary Functions", 1980, chapter 12.
  20.    
  21.     lnv  = 0.6931610107421875;     % v is close to 2, log(v) is exact
  22.     xsmall = 1.49e-8;              % approximately sqrt(eps)
  23.     xbar = 7.083964185322641e+02;  % log(1/realmin)
  24.     vf   = 0.138302778796019e-04;  % v/2-1, where v is as above
  25.     
  26.     x = abs(xin);
  27.     y = zeros(size(x));
  28.     
  29.     % sinh(x) = x to working precision
  30.     k = (xsmall >= x) | (isnan(x));
  31.     if any(any(k))            % 21 Oct 92 (jmh)
  32.         y(k) = x(k);
  33.     end
  34.     
  35.     % Use rational approximation
  36.     k = (xsmall < x) & (x <= 1);
  37.     if any(any(k))            % 21 Oct 92 (jmh)
  38.         p = [-0.35181283430177118e6 -0.11563521196851768e5 ...
  39.              -0.16375798202630751e3 -0.78966127417357099e0];
  40.         q = [-0.21108770058106271e7  0.36162723109421836e5 ...
  41.              -0.27773523119650702e3  1.0];
  42.         xx = x(k).^2;
  43.         y(k) = x(k) + x(k).*xx.*((((p(4)*xx + p(3)).*xx + p(2)).*xx + p(1))./...
  44.                                 (((xx + q(3)).*xx + q(2)).*xx + q(1)));
  45.     end
  46.     
  47.     % Use exp
  48.     k = (1 < x) & (x < xbar);
  49.     if any(any(k))            % 21 Oct 92 (jmh)
  50.         y(k) = (exp(x(k)) - exp(-x(k)))/2;
  51.     end
  52.     
  53.     % Dodge overflow problem for answers close to realmax
  54.     k = (xbar < x);
  55.     if any(any(k))            % 21 Oct 92 (jmh)
  56.         x(k) = x(k) - lnv;
  57.         y(k) = exp(x(k)) + vf .* exp(x(k));
  58.     end
  59.  
  60.     y = y.*sign(xin);
  61.  
  62. end
  63.