Print this page
11210 libm should be cstyle(1ONBLD) clean


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */

  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak __casin = casin
  31 
  32 /* INDENT OFF */
  33 /*
  34  * dcomplex casin(dcomplex z);
  35  *
  36  * Alogrithm
  37  * (based on T.E.Hull, Thomas F. Fairgrieve and Ping Tak Peter Tang's
  38  * paper "Implementing the Complex Arcsine and Arccosine Functins Using
  39  * Exception Handling", ACM TOMS, Vol 23, pp 299-335)
  40  *
  41  * The principal value of complex inverse sine function casin(z),
  42  * where z = x+iy, can be defined by
  43  *
  44  *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
  45  *
  46  * where the log function is the natural log, and
  47  *             ____________           ____________
  48  *       1    /     2    2      1    /     2    2
  49  *  A = ---  / (x+1)  + y   +  ---  / (x-1)  + y
  50  *       2 \/                   2 \/
  51  *             ____________           ____________
  52  *       1    /     2    2      1    /     2    2


 174  *         A ~ sqrt(x*x+y*y)
 175  *         B ~ x/sqrt(x*x+y*y).
 176  *      Thus
 177  *         real part = asin(B) = atan(x/y),
 178  *         imag part = log(A+sqrt(A*A-1)) ~ log(2A)
 179  *                   = log(2) + 0.5*log(x*x+y*y)
 180  *                   = log(2) + log(y) + 0.5*log(1+(x/y)^2)
 181  *
 182  *  case 6. x < 4 sqrt(u). In this case, we have
 183  *          A ~ sqrt(1+y*y), B = x/sqrt(1+y*y).
 184  *      Since B is tiny, we have
 185  *          real part = asin(B) ~ B = x/sqrt(1+y*y)
 186  *          imag part = log(A+sqrt(A*A-1)) = log (A+sqrt(y*y))
 187  *                    = log(y+sqrt(1+y*y))
 188  *                    = 0.5*log(y^2+2ysqrt(1+y^2)+1+y^2)
 189  *                    = 0.5*log(1+2y(y+sqrt(1+y^2)));
 190  *                    = 0.5*log1p(2y(y+A));
 191  *
 192  *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
 193  */
 194 /* INDENT ON */
 195 
 196 #include "libm.h"               /* asin/atan/fabs/log/log1p/sqrt */
 197 #include "complex_wrapper.h"
 198 
 199 /* INDENT OFF */
 200 static const double
 201         zero = 0.0,
 202         one = 1.0,
 203         E = 1.11022302462515654042e-16,                 /* 2**-53 */
 204         ln2 = 6.93147180559945286227e-01,
 205         pi_2 = 1.570796326794896558e+00,
 206         pi_2_l = 6.123233995736765886e-17,
 207         pi_4 = 7.85398163397448278999e-01,
 208         Foursqrtu = 5.96667258496016539463e-154,        /* 2**(-509) */
 209         Acrossover = 1.5,
 210         Bcrossover = 0.6417,
 211         half = 0.5;
 212 /* INDENT ON */
 213 
 214 dcomplex
 215 casin(dcomplex z) {

 216         double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx;
 217         int ix, iy, hx, hy;
 218         unsigned lx, ly;
 219         dcomplex ans;
 220 
 221         x = D_RE(z);
 222         y = D_IM(z);
 223         hx = HI_WORD(x);
 224         lx = LO_WORD(x);
 225         hy = HI_WORD(y);
 226         ly = LO_WORD(y);
 227         ix = hx & 0x7fffffff;
 228         iy = hy & 0x7fffffff;
 229         x = fabs(x);
 230         y = fabs(y);
 231 
 232         /* special cases */
 233 
 234         /* x is inf or NaN */
 235         if (ix >= 0x7ff00000) {      /* x is inf or NaN */
 236                 if (ISINF(ix, lx)) {    /* x is INF */
 237                         D_IM(ans) = x;

 238                         if (iy >= 0x7ff00000) {
 239                                 if (ISINF(iy, ly))
 240                                         /* casin(inf + i inf) = pi/4 + i inf */
 241                                         D_RE(ans) = pi_4;
 242                                 else    /* casin(inf + i NaN) = NaN  + i inf  */
 243                                         D_RE(ans) = y + y;
 244                         } else  /* casin(inf + iy) = pi/2 + i inf */
 245                                 D_RE(ans) = pi_2;

 246                 } else {                /* x is NaN */
 247                         if (iy >= 0x7ff00000) {
 248                                 /* INDENT OFF */
 249                                 /*
 250                                  * casin(NaN + i inf) = NaN + i inf
 251                                  * casin(NaN + i NaN) = NaN + i NaN
 252                                  */
 253                                 /* INDENT ON */
 254                                 D_IM(ans) = y + y;
 255                                 D_RE(ans) = x + x;
 256                         } else {
 257                                 /* casin(NaN + i y ) = NaN  + i NaN */
 258                                 D_IM(ans) = D_RE(ans) = x + y;
 259                         }
 260                 }

 261                 if (hx < 0)
 262                         D_RE(ans) = -D_RE(ans);

 263                 if (hy < 0)
 264                         D_IM(ans) = -D_IM(ans);

 265                 return (ans);
 266         }
 267 
 268         /* casin(+0 + i 0  ) =  0   + i 0. */
 269         if ((ix | lx | iy | ly) == 0)
 270                 return (z);
 271 
 272         if (iy >= 0x7ff00000) {      /* y is inf or NaN */
 273                 if (ISINF(iy, ly)) {    /* casin(x + i inf) =  0   + i inf */
 274                         D_IM(ans) = y;
 275                         D_RE(ans) = zero;
 276                 } else {                /* casin(x + i NaN) = NaN  + i NaN */
 277                         D_IM(ans) = x + y;

 278                         if ((ix | lx) == 0)
 279                                 D_RE(ans) = x;
 280                         else
 281                                 D_RE(ans) = y;
 282                 }

 283                 if (hx < 0)
 284                         D_RE(ans) = -D_RE(ans);

 285                 if (hy < 0)
 286                         D_IM(ans) = -D_IM(ans);

 287                 return (ans);
 288         }
 289 
 290         if ((iy | ly) == 0) {   /* region 1: y=0 */
 291                 if (ix < 0x3ff00000) {       /* |x| < 1 */
 292                         D_RE(ans) = asin(x);
 293                         D_IM(ans) = zero;
 294                 } else {
 295                         D_RE(ans) = pi_2;
 296                         if (ix >= 0x43500000)        /* |x| >= 2**54 */

 297                                 D_IM(ans) = ln2 + log(x);
 298                         else if (ix >= 0x3ff80000)   /* x > Acrossover */
 299                                 D_IM(ans) = log(x + sqrt((x - one) * (x +
 300                                         one)));
 301                         else {
 302                                 xm1 = x - one;
 303                                 D_IM(ans) = log1p(xm1 + sqrt(xm1 * (x + one)));
 304                         }
 305                 }
 306         } else if (y <= E * fabs(x - one)) { /* region 2: y < tiny*|x-1| */
 307                 if (ix < 0x3ff00000) {       /* x < 1 */
 308                         D_RE(ans) = asin(x);
 309                         D_IM(ans) = y / sqrt((one + x) * (one - x));
 310                 } else {
 311                         D_RE(ans) = pi_2;
 312                         if (ix >= 0x43500000) {      /* |x| >= 2**54 */

 313                                 D_IM(ans) = ln2 + log(x);
 314                         } else if (ix >= 0x3ff80000) /* x > Acrossover */
 315                                 D_IM(ans) = log(x + sqrt((x - one) * (x +
 316                                         one)));
 317                         else
 318                                 D_IM(ans) = log1p((x - one) + sqrt((x - one) *
 319                                         (x + one)));
 320                 }
 321         } else if (y < Foursqrtu) {  /* region 3 */
 322                 t = sqrt(y);
 323                 D_RE(ans) = pi_2 - (t - pi_2_l);
 324                 D_IM(ans) = t;
 325         } else if (E * y - one >= x) {       /* region 4 */
 326                 D_RE(ans) = x / y;      /* need to fix underflow cases */
 327                 D_IM(ans) = ln2 + log(y);
 328         } else if (ix >= 0x5fc00000 || iy >= 0x5fc00000) {        /* x,y>2**509 */
 329                 /* region 5: x+1 or y is very large (>= sqrt(max)/8) */
 330                 t = x / y;
 331                 D_RE(ans) = atan(t);
 332                 D_IM(ans) = ln2 + log(y) + half * log1p(t * t);
 333         } else if (x < Foursqrtu) {
 334                 /* region 6: x is very small, < 4sqrt(min) */
 335                 A = sqrt(one + y * y);
 336                 D_RE(ans) = x / A;      /* may underflow */

 337                 if (iy >= 0x3ff80000)        /* if y > Acrossover */
 338                         D_IM(ans) = log(y + A);
 339                 else
 340                         D_IM(ans) = half * log1p((y + y) * (y + A));
 341         } else {        /* safe region */
 342                 y2 = y * y;
 343                 xp1 = x + one;
 344                 xm1 = x - one;
 345                 R = sqrt(xp1 * xp1 + y2);
 346                 S = sqrt(xm1 * xm1 + y2);
 347                 A = half * (R + S);
 348                 B = x / A;
 349 
 350                 if (B <= Bcrossover)
 351                         D_RE(ans) = asin(B);
 352                 else {          /* use atan and an accurate approx to a-x */
 353                         Apx = A + x;

 354                         if (x <= one)
 355                                 D_RE(ans) = atan(x / sqrt(half * Apx * (y2 /
 356                                         (R + xp1) + (S - xm1))));
 357                         else
 358                                 D_RE(ans) = atan(x / (y * sqrt(half * (Apx /
 359                                         (R + xp1) + Apx / (S + xm1)))));
 360                 }

 361                 if (A <= Acrossover) {
 362                         /* use log1p and an accurate approx to A-1 */
 363                         if (x < one)
 364                                 Am1 = half * (y2 / (R + xp1) + y2 / (S - xm1));
 365                         else
 366                                 Am1 = half * (y2 / (R + xp1) + (S + xm1));

 367                         D_IM(ans) = log1p(Am1 + sqrt(Am1 * (A + one)));
 368                 } else {
 369                         D_IM(ans) = log(A + sqrt(A * A - one));
 370                 }
 371         }
 372 
 373         if (hx < 0)
 374                 D_RE(ans) = -D_RE(ans);

 375         if (hy < 0)
 376                 D_IM(ans) = -D_IM(ans);
 377 
 378         return (ans);
 379 }


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __casin = casin
  32 
  33 
  34 /*
  35  * dcomplex casin(dcomplex z);
  36  *
  37  * Alogrithm
  38  * (based on T.E.Hull, Thomas F. Fairgrieve and Ping Tak Peter Tang's
  39  * paper "Implementing the Complex Arcsine and Arccosine Functins Using
  40  * Exception Handling", ACM TOMS, Vol 23, pp 299-335)
  41  *
  42  * The principal value of complex inverse sine function casin(z),
  43  * where z = x+iy, can be defined by
  44  *
  45  *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
  46  *
  47  * where the log function is the natural log, and
  48  *             ____________           ____________
  49  *       1    /     2    2      1    /     2    2
  50  *  A = ---  / (x+1)  + y   +  ---  / (x-1)  + y
  51  *       2 \/                   2 \/
  52  *             ____________           ____________
  53  *       1    /     2    2      1    /     2    2


 175  *         A ~ sqrt(x*x+y*y)
 176  *         B ~ x/sqrt(x*x+y*y).
 177  *      Thus
 178  *         real part = asin(B) = atan(x/y),
 179  *         imag part = log(A+sqrt(A*A-1)) ~ log(2A)
 180  *                   = log(2) + 0.5*log(x*x+y*y)
 181  *                   = log(2) + log(y) + 0.5*log(1+(x/y)^2)
 182  *
 183  *  case 6. x < 4 sqrt(u). In this case, we have
 184  *          A ~ sqrt(1+y*y), B = x/sqrt(1+y*y).
 185  *      Since B is tiny, we have
 186  *          real part = asin(B) ~ B = x/sqrt(1+y*y)
 187  *          imag part = log(A+sqrt(A*A-1)) = log (A+sqrt(y*y))
 188  *                    = log(y+sqrt(1+y*y))
 189  *                    = 0.5*log(y^2+2ysqrt(1+y^2)+1+y^2)
 190  *                    = 0.5*log(1+2y(y+sqrt(1+y^2)));
 191  *                    = 0.5*log1p(2y(y+A));
 192  *
 193  *      casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)),
 194  */

 195 
 196 #include "libm.h"                       /* asin/atan/fabs/log/log1p/sqrt */
 197 #include "complex_wrapper.h"
 198 
 199 static const double zero = 0.0,


 200         one = 1.0,
 201         E = 1.11022302462515654042e-16, /* 2**-53 */
 202         ln2 = 6.93147180559945286227e-01,
 203         pi_2 = 1.570796326794896558e+00,
 204         pi_2_l = 6.123233995736765886e-17,
 205         pi_4 = 7.85398163397448278999e-01,
 206         Foursqrtu = 5.96667258496016539463e-154, /* 2**(-509) */
 207         Acrossover = 1.5,
 208         Bcrossover = 0.6417,
 209         half = 0.5;
 210 
 211 
 212 dcomplex
 213 casin(dcomplex z)
 214 {
 215         double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx;
 216         int ix, iy, hx, hy;
 217         unsigned lx, ly;
 218         dcomplex ans;
 219 
 220         x = D_RE(z);
 221         y = D_IM(z);
 222         hx = HI_WORD(x);
 223         lx = LO_WORD(x);
 224         hy = HI_WORD(y);
 225         ly = LO_WORD(y);
 226         ix = hx & 0x7fffffff;
 227         iy = hy & 0x7fffffff;
 228         x = fabs(x);
 229         y = fabs(y);
 230 
 231         /* special cases */
 232 
 233         /* x is inf or NaN */
 234         if (ix >= 0x7ff00000) {              /* x is inf or NaN */
 235                 if (ISINF(ix, lx)) {    /* x is INF */
 236                         D_IM(ans) = x;
 237 
 238                         if (iy >= 0x7ff00000) {
 239                                 if (ISINF(iy, ly))
 240                                         /* casin(inf + i inf) = pi/4 + i inf */
 241                                         D_RE(ans) = pi_4;
 242                                 else    /* casin(inf + i NaN) = NaN  + i inf  */
 243                                         D_RE(ans) = y + y;
 244                         } else { /* casin(inf + iy) = pi/2 + i inf */
 245                                 D_RE(ans) = pi_2;
 246                         }
 247                 } else {                /* x is NaN */
 248                         if (iy >= 0x7ff00000) {
 249 
 250                                 /*
 251                                  * casin(NaN + i inf) = NaN + i inf
 252                                  * casin(NaN + i NaN) = NaN + i NaN
 253                                  */

 254                                 D_IM(ans) = y + y;
 255                                 D_RE(ans) = x + x;
 256                         } else {
 257                                 /* casin(NaN + i y ) = NaN  + i NaN */
 258                                 D_IM(ans) = D_RE(ans) = x + y;
 259                         }
 260                 }
 261 
 262                 if (hx < 0)
 263                         D_RE(ans) = -D_RE(ans);
 264 
 265                 if (hy < 0)
 266                         D_IM(ans) = -D_IM(ans);
 267 
 268                 return (ans);
 269         }
 270 
 271         /* casin(+0 + i 0  ) =  0   + i 0. */
 272         if ((ix | lx | iy | ly) == 0)
 273                 return (z);
 274 
 275         if (iy >= 0x7ff00000) {              /* y is inf or NaN */
 276                 if (ISINF(iy, ly)) {    /* casin(x + i inf) =  0   + i inf */
 277                         D_IM(ans) = y;
 278                         D_RE(ans) = zero;
 279                 } else {                /* casin(x + i NaN) = NaN  + i NaN */
 280                         D_IM(ans) = x + y;
 281 
 282                         if ((ix | lx) == 0)
 283                                 D_RE(ans) = x;
 284                         else
 285                                 D_RE(ans) = y;
 286                 }
 287 
 288                 if (hx < 0)
 289                         D_RE(ans) = -D_RE(ans);
 290 
 291                 if (hy < 0)
 292                         D_IM(ans) = -D_IM(ans);
 293 
 294                 return (ans);
 295         }
 296 
 297         if ((iy | ly) == 0) {           /* region 1: y=0 */
 298                 if (ix < 0x3ff00000) {       /* |x| < 1 */
 299                         D_RE(ans) = asin(x);
 300                         D_IM(ans) = zero;
 301                 } else {
 302                         D_RE(ans) = pi_2;
 303 
 304                         if (ix >= 0x43500000) {              /* |x| >= 2**54 */
 305                                 D_IM(ans) = ln2 + log(x);
 306                         } else if (ix >= 0x3ff80000) {       /* x > Acrossover */
 307                                 D_IM(ans) = log(x + sqrt((x - one) * (x +
 308                                     one)));
 309                         } else {
 310                                 xm1 = x - one;
 311                                 D_IM(ans) = log1p(xm1 + sqrt(xm1 * (x + one)));
 312                         }
 313                 }
 314         } else if (y <= E * fabs(x - one)) { /* region 2: y < tiny*|x-1| */
 315                 if (ix < 0x3ff00000) {               /* x < 1 */
 316                         D_RE(ans) = asin(x);
 317                         D_IM(ans) = y / sqrt((one + x) * (one - x));
 318                 } else {
 319                         D_RE(ans) = pi_2;
 320 
 321                         if (ix >= 0x43500000)                /* |x| >= 2**54 */
 322                                 D_IM(ans) = ln2 + log(x);
 323                         else if (ix >= 0x3ff80000)   /* x > Acrossover */
 324                                 D_IM(ans) = log(x + sqrt((x - one) * (x +
 325                                     one)));
 326                         else
 327                                 D_IM(ans) = log1p((x - one) + sqrt((x - one) *
 328                                     (x + one)));
 329                 }
 330         } else if (y < Foursqrtu) {  /* region 3 */
 331                 t = sqrt(y);
 332                 D_RE(ans) = pi_2 - (t - pi_2_l);
 333                 D_IM(ans) = t;
 334         } else if (E * y - one >= x) { /* region 4 */
 335                 D_RE(ans) = x / y; /* need to fix underflow cases */
 336                 D_IM(ans) = ln2 + log(y);
 337         } else if (ix >= 0x5fc00000 || iy >= 0x5fc00000) {        /* x,y>2**509 */
 338                 /* region 5: x+1 or y is very large (>= sqrt(max)/8) */
 339                 t = x / y;
 340                 D_RE(ans) = atan(t);
 341                 D_IM(ans) = ln2 + log(y) + half * log1p(t * t);
 342         } else if (x < Foursqrtu) {
 343                 /* region 6: x is very small, < 4sqrt(min) */
 344                 A = sqrt(one + y * y);
 345                 D_RE(ans) = x / A;      /* may underflow */
 346 
 347                 if (iy >= 0x3ff80000)        /* if y > Acrossover */
 348                         D_IM(ans) = log(y + A);
 349                 else
 350                         D_IM(ans) = half * log1p((y + y) * (y + A));
 351         } else {                        /* safe region */
 352                 y2 = y * y;
 353                 xp1 = x + one;
 354                 xm1 = x - one;
 355                 R = sqrt(xp1 * xp1 + y2);
 356                 S = sqrt(xm1 * xm1 + y2);
 357                 A = half * (R + S);
 358                 B = x / A;
 359 
 360                 if (B <= Bcrossover) {
 361                         D_RE(ans) = asin(B);
 362                 } else {        /* use atan and an accurate approx to a-x */
 363                         Apx = A + x;
 364 
 365                         if (x <= one)
 366                                 D_RE(ans) = atan(x / sqrt(half * Apx * (y2 /
 367                                     (R + xp1) + (S - xm1))));
 368                         else
 369                                 D_RE(ans) = atan(x / (y * sqrt(half * (Apx /
 370                                     (R + xp1) + Apx / (S + xm1)))));
 371                 }
 372 
 373                 if (A <= Acrossover) {
 374                         /* use log1p and an accurate approx to A-1 */
 375                         if (x < one)
 376                                 Am1 = half * (y2 / (R + xp1) + y2 / (S - xm1));
 377                         else
 378                                 Am1 = half * (y2 / (R + xp1) + (S + xm1));
 379 
 380                         D_IM(ans) = log1p(Am1 + sqrt(Am1 * (A + one)));
 381                 } else {
 382                         D_IM(ans) = log(A + sqrt(A * A - one));
 383                 }
 384         }
 385 
 386         if (hx < 0)
 387                 D_RE(ans) = -D_RE(ans);
 388 
 389         if (hy < 0)
 390                 D_IM(ans) = -D_IM(ans);
 391 
 392         return (ans);
 393 }