477 return (-1.0L / x);
478 }
479 } else if (ny == 0) { /* y subnormal, scale it */
480 y *= two240;
481 iy = HI_XWORD(y);
482 ny = (iy >> 16) - 240;
483 }
484 n = nx - ny;
485 /*
486 * When y is zero or when x >> y, i.e., n > 62, 78, 122 for DBLE,
487 * EXTENDED, QUAD respectively,
488 * log(x) = log(sqrt(x * x + y * y)) to 27 extra bits.
489 */
490
491 #if defined(__x86)
492 if (n > 78 || y == 0.0L) {
493 #else
494 if (n > 122 || y == 0.0L) {
495 #endif
496
497 XFSCALE(x, 0x3fff - (ix >> 16));
498 i = ((ix & 0xffff) + 0x100) >> 9; /* 7.5 bits of x */
499 zk = 1.0L + ((long double) i) * 0.0078125L;
500 z = x - zk;
501 dk = (double)z;
502
503 #if defined(__x86)
504 ((unsigned *)&dk)[LOWORD] &= 0xfffe0000;
505 #endif
506
507 zh = (long double)dk;
508 k = i & 0x7f; /* index of zk */
509 n = nx - 0x3fff;
510 *er = z - zh;
511 if (i == 0x80) { /* if zk = 2.0, adjust scaling */
512 n += 1;
513 zh *= 0.5L; *er *= 0.5L;
514 }
515 w = k_log_NKzl(n, k, zh, er);
516 } else {
517 /*
518 * compute z = x*x + y*y
519 */
520 XFSCALE(x, 0x3fff - (ix >> 16));
521 XFSCALE(y, 0x3fff - n - (iy >> 16));
522 ix = (ix & 0xffff) | 0x3fff0000;
523 iy = (iy & 0xffff) | (0x3fff0000 - (n << 16));
524 nx -= 0x3fff;
525 t1 = x * x; t2 = y * y;
526 wh = x;
527
528 /* split x into correctly rounded half */
529 #if defined(__x86)
530 ((unsigned *)&wh)[0] = 0; /* 32 bits chopped */
531 #else
532 lx = ((unsigned *)&wh)[2]; /* 56 rounded */
533 j = ((lx >> 24) + 1) >> 1;
534 ((unsigned *)&wh)[2] = (j << 25);
535 lx = ((unsigned *)&wh)[1];
536 ly = lx + (j >> 7);
537 ((unsigned *)&wh)[1] = ly;
538 ((unsigned *)&wh)[0] += (ly == 0 && lx != 0);
539 ((unsigned *)&wh)[3] = 0;
540 #endif
541
|
477 return (-1.0L / x);
478 }
479 } else if (ny == 0) { /* y subnormal, scale it */
480 y *= two240;
481 iy = HI_XWORD(y);
482 ny = (iy >> 16) - 240;
483 }
484 n = nx - ny;
485 /*
486 * When y is zero or when x >> y, i.e., n > 62, 78, 122 for DBLE,
487 * EXTENDED, QUAD respectively,
488 * log(x) = log(sqrt(x * x + y * y)) to 27 extra bits.
489 */
490
491 #if defined(__x86)
492 if (n > 78 || y == 0.0L) {
493 #else
494 if (n > 122 || y == 0.0L) {
495 #endif
496
497 XFSCALE(x, (0x3fff - (ix >> 16)));
498 i = ((ix & 0xffff) + 0x100) >> 9; /* 7.5 bits of x */
499 zk = 1.0L + ((long double) i) * 0.0078125L;
500 z = x - zk;
501 dk = (double)z;
502
503 #if defined(__x86)
504 ((unsigned *)&dk)[LOWORD] &= 0xfffe0000;
505 #endif
506
507 zh = (long double)dk;
508 k = i & 0x7f; /* index of zk */
509 n = nx - 0x3fff;
510 *er = z - zh;
511 if (i == 0x80) { /* if zk = 2.0, adjust scaling */
512 n += 1;
513 zh *= 0.5L; *er *= 0.5L;
514 }
515 w = k_log_NKzl(n, k, zh, er);
516 } else {
517 /*
518 * compute z = x*x + y*y
519 */
520 XFSCALE(x, (0x3fff - (ix >> 16)));
521 XFSCALE(y, (0x3fff - n - (iy >> 16)));
522 ix = (ix & 0xffff) | 0x3fff0000;
523 iy = (iy & 0xffff) | (0x3fff0000 - (n << 16));
524 nx -= 0x3fff;
525 t1 = x * x; t2 = y * y;
526 wh = x;
527
528 /* split x into correctly rounded half */
529 #if defined(__x86)
530 ((unsigned *)&wh)[0] = 0; /* 32 bits chopped */
531 #else
532 lx = ((unsigned *)&wh)[2]; /* 56 rounded */
533 j = ((lx >> 24) + 1) >> 1;
534 ((unsigned *)&wh)[2] = (j << 25);
535 lx = ((unsigned *)&wh)[1];
536 ly = lx + (j >> 7);
537 ((unsigned *)&wh)[1] = ly;
538 ((unsigned *)&wh)[0] += (ly == 0 && lx != 0);
539 ((unsigned *)&wh)[3] = 0;
540 #endif
541
|