diff options
author | Wincent Balin <wincent@rockbox.org> | 2009-08-04 15:17:34 +0000 |
---|---|---|
committer | Wincent Balin <wincent@rockbox.org> | 2009-08-04 15:17:34 +0000 |
commit | 1f9675227ac7a47bf231f5095ce93b48f7ba5fbe (patch) | |
tree | b59d82eef6f58f176163a9c4c9c690dba339964c /apps/plugins/pdbox/pdbox-func.c | |
parent | 010bb8e6bad91c2c03a37f2a5f51820910500bdc (diff) | |
download | rockbox-1f9675227ac7a47bf231f5095ce93b48f7ba5fbe.tar.gz rockbox-1f9675227ac7a47bf231f5095ce93b48f7ba5fbe.zip |
PDBox: Undoing incorrect changes of math functions.
git-svn-id: svn://svn.rockbox.org/rockbox/trunk@22159 a1c6a512-1295-4272-9138-f99709370657
Diffstat (limited to 'apps/plugins/pdbox/pdbox-func.c')
-rw-r--r-- | apps/plugins/pdbox/pdbox-func.c | 856 |
1 files changed, 815 insertions, 41 deletions
diff --git a/apps/plugins/pdbox/pdbox-func.c b/apps/plugins/pdbox/pdbox-func.c index 42427a69c3..b019959790 100644 --- a/apps/plugins/pdbox/pdbox-func.c +++ b/apps/plugins/pdbox/pdbox-func.c | |||
@@ -582,55 +582,284 @@ float rb_log10(float x) | |||
582 | } | 582 | } |
583 | 583 | ||
584 | 584 | ||
585 | /* Power function, taken from glibc-2.8 and dietlibc-0.32 */ | 585 | /* Power function, |
586 | Taken from glibc-2.8 */ | ||
587 | |||
588 | int rb_isinf(float x) | ||
589 | { | ||
590 | int32_t ix, t; | ||
591 | GET_FLOAT_WORD(ix,x); | ||
592 | t = ix & 0x7fffffff; | ||
593 | t ^= 0x7f800000; | ||
594 | t |= -t; | ||
595 | return ~(t >> 31) & (ix >> 30); | ||
596 | } | ||
597 | |||
598 | float rb_copysignf(float x, float y) | ||
599 | { | ||
600 | uint32_t ix, iy; | ||
601 | GET_FLOAT_WORD(ix,x); | ||
602 | GET_FLOAT_WORD(iy,y); | ||
603 | SET_FLOAT_WORD(x,(ix&0x7fffffff)|(iy&0x80000000)); | ||
604 | return x; | ||
605 | } | ||
606 | |||
586 | static const float | 607 | static const float |
587 | huge = 1.0e+30, | 608 | huge = 1.0e+30, |
588 | tiny = 1.0e-30, | 609 | tiny = 1.0e-30, |
589 | one = 1.0f; | 610 | twom25 = 2.9802322388e-08; /* 0x33000000 */ |
590 | 611 | ||
591 | float rb_pow(float x, float y) | 612 | float rb_scalbnf(float x, int n) |
592 | { | 613 | { |
593 | unsigned int e; | 614 | int32_t k, ix; |
594 | float result; | 615 | GET_FLOAT_WORD(ix,x); |
595 | 616 | k = (ix&0x7f800000)>>23; /* extract exponent */ | |
596 | /* Special cases 0^x */ | 617 | if (k==0) { /* 0 or subnormal x */ |
597 | if(x == 0.0f) | 618 | if ((ix&0x7fffffff)==0) return x; /* +-0 */ |
598 | { | 619 | x *= two25; |
599 | if(y > 0.0f) | 620 | GET_FLOAT_WORD(ix,x); |
600 | return 0.0f; | 621 | k = ((ix&0x7f800000)>>23) - 25; |
601 | else if(y == 0.0f) | ||
602 | return 1.0f; | ||
603 | else | ||
604 | return 1.0f / x; | ||
605 | } | 622 | } |
623 | if (k==0xff) return x+x; /* NaN or Inf */ | ||
624 | k = k+n; | ||
625 | if (n> 50000 || k > 0xfe) | ||
626 | return huge*rb_copysignf(huge,x); /* overflow */ | ||
627 | if (n< -50000) | ||
628 | return tiny*rb_copysignf(tiny,x); /*underflow*/ | ||
629 | if (k > 0) /* normal result */ | ||
630 | {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;} | ||
631 | if (k <= -25) | ||
632 | return tiny*rb_copysignf(tiny,x); /*underflow*/ | ||
633 | k += 25; /* subnormal result */ | ||
634 | SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); | ||
635 | return x*twom25; | ||
636 | } | ||
606 | 637 | ||
607 | /* Special case x^n where n is integer */ | ||
608 | if(y == (int) (e = (int) y)) | ||
609 | { | ||
610 | if((int) e < 0) | ||
611 | { | ||
612 | e = -e; | ||
613 | x = 1.0f / x; | ||
614 | } | ||
615 | 638 | ||
616 | result = 1.0f; | 639 | static const float |
640 | bp[] = {1.0, 1.5,}, | ||
641 | dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */ | ||
642 | dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */ | ||
643 | one = 1.0, | ||
644 | two = 2.0, | ||
645 | two24 = 16777216.0, /* 0x4b800000 */ | ||
646 | /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ | ||
647 | L1 = 6.0000002384e-01, /* 0x3f19999a */ | ||
648 | L2 = 4.2857143283e-01, /* 0x3edb6db7 */ | ||
649 | L3 = 3.3333334327e-01, /* 0x3eaaaaab */ | ||
650 | L4 = 2.7272811532e-01, /* 0x3e8ba305 */ | ||
651 | L5 = 2.3066075146e-01, /* 0x3e6c3255 */ | ||
652 | L6 = 2.0697501302e-01, /* 0x3e53f142 */ | ||
653 | P1 = 1.6666667163e-01, /* 0x3e2aaaab */ | ||
654 | P2 = -2.7777778450e-03, /* 0xbb360b61 */ | ||
655 | P3 = 6.6137559770e-05, /* 0x388ab355 */ | ||
656 | P4 = -1.6533901999e-06, /* 0xb5ddea0e */ | ||
657 | P5 = 4.1381369442e-08; /* 0x3331bb4c */ | ||
617 | 658 | ||
618 | while(1) | 659 | static const float |
619 | { | 660 | lg2 = 6.9314718246e-01, /* 0x3f317218 */ |
620 | if(e & 1) | 661 | lg2_h = 6.93145752e-01, /* 0x3f317200 */ |
621 | result *= x; | 662 | lg2_l = 1.42860654e-06, /* 0x35bfbe8c */ |
663 | ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */ | ||
664 | cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */ | ||
665 | cp_h = 9.6179199219e-01, /* 0x3f763800 =head of cp */ | ||
666 | cp_l = 4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */ | ||
667 | ivln2 = 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */ | ||
668 | ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/ | ||
669 | ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ | ||
622 | 670 | ||
623 | if((e >>= 1) == 0) | 671 | float rb_pow(float x, float y) |
624 | break; | 672 | { |
673 | float z, ax, z_h, z_l, p_h, p_l; | ||
674 | float y1, t1, t2, r, s, t, u, v, w; | ||
675 | int32_t i, j, k, yisint, n; | ||
676 | int32_t hx, hy, ix, iy, is; | ||
677 | |||
678 | GET_FLOAT_WORD(hx,x); | ||
679 | GET_FLOAT_WORD(hy,y); | ||
680 | ix = hx&0x7fffffff; | ||
681 | iy = hy&0x7fffffff; | ||
682 | |||
683 | /* y==zero: x**0 = 1 */ | ||
684 | if(iy==0) return one; | ||
685 | |||
686 | /* x==+-1 */ | ||
687 | if(x == 1.0) return one; | ||
688 | if(x == -1.0 && rb_isinf(y)) return one; | ||
689 | |||
690 | /* +-NaN return x+y */ | ||
691 | if(ix > 0x7f800000 || iy > 0x7f800000) | ||
692 | return x+y; | ||
693 | |||
694 | /* determine if y is an odd int when x < 0 | ||
695 | * yisint = 0 ... y is not an integer | ||
696 | * yisint = 1 ... y is an odd int | ||
697 | * yisint = 2 ... y is an even int | ||
698 | */ | ||
699 | yisint = 0; | ||
700 | if(hx<0) { | ||
701 | if(iy>=0x4b800000) yisint = 2; /* even integer y */ | ||
702 | else if(iy>=0x3f800000) { | ||
703 | k = (iy>>23)-0x7f; /* exponent */ | ||
704 | j = iy>>(23-k); | ||
705 | if((j<<(23-k))==iy) yisint = 2-(j&1); | ||
706 | } | ||
707 | } | ||
625 | 708 | ||
626 | x *= x; | 709 | /* special value of y */ |
710 | if (iy==0x7f800000) { /* y is +-inf */ | ||
711 | if (ix==0x3f800000) | ||
712 | return y - y; /* inf**+-1 is NaN */ | ||
713 | else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */ | ||
714 | return (hy>=0)? y: zero; | ||
715 | else /* (|x|<1)**-,+inf = inf,0 */ | ||
716 | return (hy<0)?-y: zero; | ||
717 | } | ||
718 | if(iy==0x3f800000) { /* y is +-1 */ | ||
719 | if(hy<0) return one/x; else return x; | ||
720 | } | ||
721 | if(hy==0x40000000) return x*x; /* y is 2 */ | ||
722 | if(hy==0x3f000000) { /* y is 0.5 */ | ||
723 | if(hx>=0) /* x >= +0 */ | ||
724 | return rb_sqrt(x); | ||
725 | } | ||
726 | |||
727 | ax = rb_fabs(x); | ||
728 | /* special value of x */ | ||
729 | if(ix==0x7f800000||ix==0||ix==0x3f800000){ | ||
730 | z = ax; /*x is +-0,+-inf,+-1*/ | ||
731 | if(hy<0) z = one/z; /* z = (1/|x|) */ | ||
732 | if(hx<0) { | ||
733 | if(((ix-0x3f800000)|yisint)==0) { | ||
734 | z = (z-z)/(z-z); /* (-1)**non-int is NaN */ | ||
735 | } else if(yisint==1) | ||
736 | z = -z; /* (x<0)**odd = -(|x|**odd) */ | ||
627 | } | 737 | } |
738 | return z; | ||
739 | } | ||
628 | 740 | ||
629 | return result; | 741 | /* (x<0)**(non-int) is NaN */ |
742 | if(((((uint32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x); | ||
743 | |||
744 | /* |y| is huge */ | ||
745 | if(iy>0x4d000000) { /* if |y| > 2**27 */ | ||
746 | /* over/underflow if x is not close to one */ | ||
747 | if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny; | ||
748 | if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny; | ||
749 | /* now |1-x| is tiny <= 2**-20, suffice to compute | ||
750 | log(x) by x-x^2/2+x^3/3-x^4/4 */ | ||
751 | t = x-1; /* t has 20 trailing zeros */ | ||
752 | w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25)); | ||
753 | u = ivln2_h*t; /* ivln2_h has 16 sig. bits */ | ||
754 | v = t*ivln2_l-w*ivln2; | ||
755 | t1 = u+v; | ||
756 | GET_FLOAT_WORD(is,t1); | ||
757 | SET_FLOAT_WORD(t1,is&0xfffff000); | ||
758 | t2 = v-(t1-u); | ||
759 | } else { | ||
760 | float s2, s_h, s_l, t_h, t_l; | ||
761 | n = 0; | ||
762 | /* take care subnormal number */ | ||
763 | if(ix<0x00800000) | ||
764 | {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); } | ||
765 | n += ((ix)>>23)-0x7f; | ||
766 | j = ix&0x007fffff; | ||
767 | /* determine interval */ | ||
768 | ix = j|0x3f800000; /* normalize ix */ | ||
769 | if(j<=0x1cc471) k=0; /* |x|<sqrt(3/2) */ | ||
770 | else if(j<0x5db3d7) k=1; /* |x|<sqrt(3) */ | ||
771 | else {k=0;n+=1;ix -= 0x00800000;} | ||
772 | SET_FLOAT_WORD(ax,ix); | ||
773 | |||
774 | /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ | ||
775 | u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ | ||
776 | v = one/(ax+bp[k]); | ||
777 | s = u*v; | ||
778 | s_h = s; | ||
779 | GET_FLOAT_WORD(is,s_h); | ||
780 | SET_FLOAT_WORD(s_h,is&0xfffff000); | ||
781 | /* t_h=ax+bp[k] High */ | ||
782 | SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21)); | ||
783 | t_l = ax - (t_h-bp[k]); | ||
784 | s_l = v*((u-s_h*t_h)-s_h*t_l); | ||
785 | /* compute log(ax) */ | ||
786 | s2 = s*s; | ||
787 | r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); | ||
788 | r += s_l*(s_h+s); | ||
789 | s2 = s_h*s_h; | ||
790 | t_h = (float)3.0+s2+r; | ||
791 | GET_FLOAT_WORD(is,t_h); | ||
792 | SET_FLOAT_WORD(t_h,is&0xfffff000); | ||
793 | t_l = r-((t_h-(float)3.0)-s2); | ||
794 | /* u+v = s*(1+...) */ | ||
795 | u = s_h*t_h; | ||
796 | v = s_l*t_h+t_l*s; | ||
797 | /* 2/(3log2)*(s+...) */ | ||
798 | p_h = u+v; | ||
799 | GET_FLOAT_WORD(is,p_h); | ||
800 | SET_FLOAT_WORD(p_h,is&0xfffff000); | ||
801 | p_l = v-(p_h-u); | ||
802 | z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ | ||
803 | z_l = cp_l*p_h+p_l*cp+dp_l[k]; | ||
804 | /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ | ||
805 | t = (float)n; | ||
806 | t1 = (((z_h+z_l)+dp_h[k])+t); | ||
807 | GET_FLOAT_WORD(is,t1); | ||
808 | SET_FLOAT_WORD(t1,is&0xfffff000); | ||
809 | t2 = z_l-(((t1-t)-dp_h[k])-z_h); | ||
630 | } | 810 | } |
631 | 811 | ||
632 | /* Normal case */ | 812 | s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ |
633 | return rb_exp(rb_log(x) * y); | 813 | if(((((uint32_t)hx>>31)-1)|(yisint-1))==0) |
814 | s = -one; /* (-ve)**(odd int) */ | ||
815 | |||
816 | /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ | ||
817 | GET_FLOAT_WORD(is,y); | ||
818 | SET_FLOAT_WORD(y1,is&0xfffff000); | ||
819 | p_l = (y-y1)*t1+y*t2; | ||
820 | p_h = y1*t1; | ||
821 | z = p_l+p_h; | ||
822 | GET_FLOAT_WORD(j,z); | ||
823 | if (j>0x43000000) /* if z > 128 */ | ||
824 | return s*huge*huge; /* overflow */ | ||
825 | else if (j==0x43000000) { /* if z == 128 */ | ||
826 | if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ | ||
827 | } | ||
828 | else if ((j&0x7fffffff)>0x43160000) /* z <= -150 */ | ||
829 | return s*tiny*tiny; /* underflow */ | ||
830 | else if ((uint32_t) j==0xc3160000){ /* z == -150 */ | ||
831 | if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ | ||
832 | } | ||
833 | /* | ||
834 | * compute 2**(p_h+p_l) | ||
835 | */ | ||
836 | i = j&0x7fffffff; | ||
837 | k = (i>>23)-0x7f; | ||
838 | n = 0; | ||
839 | if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */ | ||
840 | n = j+(0x00800000>>(k+1)); | ||
841 | k = ((n&0x7fffffff)>>23)-0x7f; /* new k for n */ | ||
842 | SET_FLOAT_WORD(t,n&~(0x007fffff>>k)); | ||
843 | n = ((n&0x007fffff)|0x00800000)>>(23-k); | ||
844 | if(j<0) n = -n; | ||
845 | p_h -= t; | ||
846 | } | ||
847 | t = p_l+p_h; | ||
848 | GET_FLOAT_WORD(is,t); | ||
849 | SET_FLOAT_WORD(t,is&0xfffff000); | ||
850 | u = t*lg2_h; | ||
851 | v = (p_l-(t-p_h))*lg2+t*lg2_l; | ||
852 | z = u+v; | ||
853 | w = v-(z-u); | ||
854 | t = z*z; | ||
855 | t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); | ||
856 | r = (z*t1)/(t1-two)-(w+z*w); | ||
857 | z = one-(r-z); | ||
858 | GET_FLOAT_WORD(j,z); | ||
859 | j += (n<<23); | ||
860 | if((j>>23)<=0) z = rb_scalbnf(z,n); /* subnormal output */ | ||
861 | else SET_FLOAT_WORD(z,j); | ||
862 | return s*z; | ||
634 | } | 863 | } |
635 | 864 | ||
636 | 865 | ||
@@ -701,10 +930,14 @@ float rb_sqrt(float x) | |||
701 | return z; | 930 | return z; |
702 | } | 931 | } |
703 | 932 | ||
704 | /* Absolute value, simple calculus */ | 933 | /* Absolute value, |
934 | taken from glibc-2.8 */ | ||
705 | float rb_fabs(float x) | 935 | float rb_fabs(float x) |
706 | { | 936 | { |
707 | return (x < 0.0f) ? -x : x; | 937 | uint32_t ix; |
938 | GET_FLOAT_WORD(ix,x); | ||
939 | SET_FLOAT_WORD(x,ix&0x7fffffff); | ||
940 | return x; | ||
708 | } | 941 | } |
709 | 942 | ||
710 | /* Arc tangent, | 943 | /* Arc tangent, |
@@ -860,18 +1093,559 @@ float rb_atan2(float x, float y) | |||
860 | } | 1093 | } |
861 | 1094 | ||
862 | 1095 | ||
863 | /* Sine hyperbolic, taken from dietlibc-0.32 */ | 1096 | /* Sine hyperbolic, |
1097 | taken from glibc-2.8 */ | ||
1098 | |||
1099 | static const float | ||
1100 | o_threshold = 8.8721679688e+01,/* 0x42b17180 */ | ||
1101 | invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */ | ||
1102 | /* scaled coefficients related to expm1 */ | ||
1103 | Q1 = -3.3333335072e-02, /* 0xbd088889 */ | ||
1104 | Q2 = 1.5873016091e-03, /* 0x3ad00d01 */ | ||
1105 | Q3 = -7.9365076090e-05, /* 0xb8a670cd */ | ||
1106 | Q4 = 4.0082177293e-06, /* 0x36867e54 */ | ||
1107 | Q5 = -2.0109921195e-07; /* 0xb457edbb */ | ||
1108 | |||
1109 | float rb_expm1(float x) | ||
1110 | { | ||
1111 | float y,hi,lo,c=0,t,e,hxs,hfx,r1; | ||
1112 | int32_t k,xsb; | ||
1113 | uint32_t hx; | ||
1114 | |||
1115 | GET_FLOAT_WORD(hx,x); | ||
1116 | xsb = hx&0x80000000; /* sign bit of x */ | ||
1117 | if(xsb==0) y=x; else y= -x; /* y = |x| */ | ||
1118 | hx &= 0x7fffffff; /* high word of |x| */ | ||
1119 | |||
1120 | /* filter out huge and non-finite argument */ | ||
1121 | if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */ | ||
1122 | if(hx >= 0x42b17218) { /* if |x|>=88.721... */ | ||
1123 | if(hx>0x7f800000) | ||
1124 | return x+x; /* NaN */ | ||
1125 | if(hx==0x7f800000) | ||
1126 | return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ | ||
1127 | if(x > o_threshold) return huge*huge; /* overflow */ | ||
1128 | } | ||
1129 | if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */ | ||
1130 | if(x+tiny<(float)0.0) /* raise inexact */ | ||
1131 | return tiny-one; /* return -1 */ | ||
1132 | } | ||
1133 | } | ||
1134 | |||
1135 | /* argument reduction */ | ||
1136 | if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ | ||
1137 | if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ | ||
1138 | if(xsb==0) | ||
1139 | {hi = x - ln2_hi; lo = ln2_lo; k = 1;} | ||
1140 | else | ||
1141 | {hi = x + ln2_hi; lo = -ln2_lo; k = -1;} | ||
1142 | } else { | ||
1143 | k = invln2*x+((xsb==0)?(float)0.5:(float)-0.5); | ||
1144 | t = k; | ||
1145 | hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ | ||
1146 | lo = t*ln2_lo; | ||
1147 | } | ||
1148 | x = hi - lo; | ||
1149 | c = (hi-x)-lo; | ||
1150 | } | ||
1151 | else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ | ||
1152 | t = huge+x; /* return x with inexact flags when x!=0 */ | ||
1153 | return x - (t-(huge+x)); | ||
1154 | } | ||
1155 | else k = 0; | ||
1156 | |||
1157 | /* x is now in primary range */ | ||
1158 | hfx = (float)0.5*x; | ||
1159 | hxs = x*hfx; | ||
1160 | r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))); | ||
1161 | t = (float)3.0-r1*hfx; | ||
1162 | e = hxs*((r1-t)/((float)6.0 - x*t)); | ||
1163 | if(k==0) return x - (x*e-hxs); /* c is 0 */ | ||
1164 | else { | ||
1165 | e = (x*(e-c)-c); | ||
1166 | e -= hxs; | ||
1167 | if(k== -1) return (float)0.5*(x-e)-(float)0.5; | ||
1168 | if(k==1) { | ||
1169 | if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5)); | ||
1170 | else return one+(float)2.0*(x-e); | ||
1171 | } | ||
1172 | if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ | ||
1173 | int32_t i; | ||
1174 | y = one-(e-x); | ||
1175 | GET_FLOAT_WORD(i,y); | ||
1176 | SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ | ||
1177 | return y-one; | ||
1178 | } | ||
1179 | t = one; | ||
1180 | if(k<23) { | ||
1181 | int32_t i; | ||
1182 | SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ | ||
1183 | y = t-(e-x); | ||
1184 | GET_FLOAT_WORD(i,y); | ||
1185 | SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ | ||
1186 | } else { | ||
1187 | int32_t i; | ||
1188 | SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */ | ||
1189 | y = x-(e+t); | ||
1190 | y += one; | ||
1191 | GET_FLOAT_WORD(i,y); | ||
1192 | SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ | ||
1193 | } | ||
1194 | } | ||
1195 | return y; | ||
1196 | } | ||
1197 | |||
1198 | static const float shuge = 1.0e37; | ||
1199 | |||
864 | float rb_sinh(float x) | 1200 | float rb_sinh(float x) |
865 | { | 1201 | { |
866 | float y = rb_exp(x); | 1202 | float t,w,h; |
867 | return (y - 1.0f / y) * 0.5f; | 1203 | int32_t ix,jx; |
1204 | |||
1205 | GET_FLOAT_WORD(jx,x); | ||
1206 | ix = jx&0x7fffffff; | ||
1207 | |||
1208 | /* x is INF or NaN */ | ||
1209 | if(ix>=0x7f800000) return x+x; | ||
1210 | |||
1211 | h = 0.5; | ||
1212 | if (jx<0) h = -h; | ||
1213 | /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ | ||
1214 | if (ix < 0x41b00000) { /* |x|<22 */ | ||
1215 | if (ix<0x31800000) /* |x|<2**-28 */ | ||
1216 | if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ | ||
1217 | t = rb_expm1(rb_fabs(x)); | ||
1218 | if(ix<0x3f800000) return h*((float)2.0*t-t*t/(t+one)); | ||
1219 | return h*(t+t/(t+one)); | ||
1220 | } | ||
1221 | |||
1222 | /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ | ||
1223 | if (ix < 0x42b17180) return h*rb_exp(rb_fabs(x)); | ||
1224 | |||
1225 | /* |x| in [log(maxdouble), overflowthresold] */ | ||
1226 | if (ix<=0x42b2d4fc) { | ||
1227 | w = rb_exp((float)0.5*rb_fabs(x)); | ||
1228 | t = h*w; | ||
1229 | return t*w; | ||
1230 | } | ||
1231 | |||
1232 | /* |x| > overflowthresold, sinh(x) overflow */ | ||
1233 | return x*shuge; | ||
1234 | } | ||
1235 | |||
1236 | |||
1237 | /* Tangent, | ||
1238 | taken from glibc-2.8 */ | ||
1239 | |||
1240 | static const float | ||
1241 | pio4 = 7.8539812565e-01, /* 0x3f490fda */ | ||
1242 | pio4lo= 3.7748947079e-08, /* 0x33222168 */ | ||
1243 | T[] = { | ||
1244 | 3.3333334327e-01, /* 0x3eaaaaab */ | ||
1245 | 1.3333334029e-01, /* 0x3e088889 */ | ||
1246 | 5.3968254477e-02, /* 0x3d5d0dd1 */ | ||
1247 | 2.1869488060e-02, /* 0x3cb327a4 */ | ||
1248 | 8.8632395491e-03, /* 0x3c11371f */ | ||
1249 | 3.5920790397e-03, /* 0x3b6b6916 */ | ||
1250 | 1.4562094584e-03, /* 0x3abede48 */ | ||
1251 | 5.8804126456e-04, /* 0x3a1a26c8 */ | ||
1252 | 2.4646313977e-04, /* 0x398137b9 */ | ||
1253 | 7.8179444245e-05, /* 0x38a3f445 */ | ||
1254 | 7.1407252108e-05, /* 0x3895c07a */ | ||
1255 | -1.8558637748e-05, /* 0xb79bae5f */ | ||
1256 | 2.5907305826e-05, /* 0x37d95384 */ | ||
1257 | }; | ||
1258 | |||
1259 | float kernel_tan(float x, float y, int iy) | ||
1260 | { | ||
1261 | float z,r,v,w,s; | ||
1262 | int32_t ix,hx; | ||
1263 | GET_FLOAT_WORD(hx,x); | ||
1264 | ix = hx&0x7fffffff; /* high word of |x| */ | ||
1265 | if(ix<0x31800000) /* x < 2**-28 */ | ||
1266 | {if((int)x==0) { /* generate inexact */ | ||
1267 | if((ix|(iy+1))==0) return one/rb_fabs(x); | ||
1268 | else return (iy==1)? x: -one/x; | ||
1269 | } | ||
1270 | } | ||
1271 | if(ix>=0x3f2ca140) { /* |x|>=0.6744 */ | ||
1272 | if(hx<0) {x = -x; y = -y;} | ||
1273 | z = pio4-x; | ||
1274 | w = pio4lo-y; | ||
1275 | x = z+w; y = 0.0; | ||
1276 | } | ||
1277 | z = x*x; | ||
1278 | w = z*z; | ||
1279 | /* Break x^5*(T[1]+x^2*T[2]+...) into | ||
1280 | * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) + | ||
1281 | * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12])) | ||
1282 | */ | ||
1283 | r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11])))); | ||
1284 | v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12]))))); | ||
1285 | s = z*x; | ||
1286 | r = y + z*(s*(r+v)+y); | ||
1287 | r += T[0]*s; | ||
1288 | w = x+r; | ||
1289 | if(ix>=0x3f2ca140) { | ||
1290 | v = (float)iy; | ||
1291 | return (float)(1-((hx>>30)&2))*(v-(float)2.0*(x-(w*w/(w+v)-r))); | ||
1292 | } | ||
1293 | if(iy==1) return w; | ||
1294 | else { /* if allow error up to 2 ulp, | ||
1295 | simply return -1.0/(x+r) here */ | ||
1296 | /* compute -1.0/(x+r) accurately */ | ||
1297 | float a,t; | ||
1298 | int32_t i; | ||
1299 | z = w; | ||
1300 | GET_FLOAT_WORD(i,z); | ||
1301 | SET_FLOAT_WORD(z,i&0xfffff000); | ||
1302 | v = r-(z - x); /* z+v = r+x */ | ||
1303 | t = a = -(float)1.0/w; /* a = -1.0/w */ | ||
1304 | GET_FLOAT_WORD(i,t); | ||
1305 | SET_FLOAT_WORD(t,i&0xfffff000); | ||
1306 | s = (float)1.0+t*z; | ||
1307 | return t+a*(s+t*v); | ||
1308 | } | ||
868 | } | 1309 | } |
869 | 1310 | ||
870 | 1311 | ||
871 | /* Tangent, simple calculus solution. */ | 1312 | |
1313 | static const int init_jk[] = {4,7,9}; /* initial value for jk */ | ||
1314 | |||
1315 | static const float PIo2[] = { | ||
1316 | 1.5703125000e+00, /* 0x3fc90000 */ | ||
1317 | 4.5776367188e-04, /* 0x39f00000 */ | ||
1318 | 2.5987625122e-05, /* 0x37da0000 */ | ||
1319 | 7.5437128544e-08, /* 0x33a20000 */ | ||
1320 | 6.0026650317e-11, /* 0x2e840000 */ | ||
1321 | 7.3896444519e-13, /* 0x2b500000 */ | ||
1322 | 5.3845816694e-15, /* 0x27c20000 */ | ||
1323 | 5.6378512969e-18, /* 0x22d00000 */ | ||
1324 | 8.3009228831e-20, /* 0x1fc40000 */ | ||
1325 | 3.2756352257e-22, /* 0x1bc60000 */ | ||
1326 | 6.3331015649e-25, /* 0x17440000 */ | ||
1327 | }; | ||
1328 | |||
1329 | static const float | ||
1330 | two8 = 2.5600000000e+02, /* 0x43800000 */ | ||
1331 | twon8 = 3.9062500000e-03; /* 0x3b800000 */ | ||
1332 | |||
1333 | int kernel_rem_pio2(float *x, float *y, int e0, int nx, int prec, const int32_t *ipio2) | ||
1334 | { | ||
1335 | int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; | ||
1336 | float z,fw,f[20],fq[20],q[20]; | ||
1337 | |||
1338 | /* initialize jk*/ | ||
1339 | jk = init_jk[prec]; | ||
1340 | jp = jk; | ||
1341 | |||
1342 | /* determine jx,jv,q0, note that 3>q0 */ | ||
1343 | jx = nx-1; | ||
1344 | jv = (e0-3)/8; if(jv<0) jv=0; | ||
1345 | q0 = e0-8*(jv+1); | ||
1346 | |||
1347 | /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ | ||
1348 | j = jv-jx; m = jx+jk; | ||
1349 | for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (float) ipio2[j]; | ||
1350 | |||
1351 | /* compute q[0],q[1],...q[jk] */ | ||
1352 | for (i=0;i<=jk;i++) { | ||
1353 | for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw; | ||
1354 | } | ||
1355 | |||
1356 | jz = jk; | ||
1357 | recompute: | ||
1358 | /* distill q[] into iq[] reversingly */ | ||
1359 | for(i=0,j=jz,z=q[jz];j>0;i++,j--) { | ||
1360 | fw = (float)((int32_t)(twon8* z)); | ||
1361 | iq[i] = (int32_t)(z-two8*fw); | ||
1362 | z = q[j-1]+fw; | ||
1363 | } | ||
1364 | |||
1365 | /* compute n */ | ||
1366 | z = rb_scalbnf(z,q0); /* actual value of z */ | ||
1367 | z -= (float)8.0*rb_floor(z*(float)0.125); /* trim off integer >= 8 */ | ||
1368 | n = (int32_t) z; | ||
1369 | z -= (float)n; | ||
1370 | ih = 0; | ||
1371 | if(q0>0) { /* need iq[jz-1] to determine n */ | ||
1372 | i = (iq[jz-1]>>(8-q0)); n += i; | ||
1373 | iq[jz-1] -= i<<(8-q0); | ||
1374 | ih = iq[jz-1]>>(7-q0); | ||
1375 | } | ||
1376 | else if(q0==0) ih = iq[jz-1]>>8; | ||
1377 | else if(z>=(float)0.5) ih=2; | ||
1378 | |||
1379 | if(ih>0) { /* q > 0.5 */ | ||
1380 | n += 1; carry = 0; | ||
1381 | for(i=0;i<jz ;i++) { /* compute 1-q */ | ||
1382 | j = iq[i]; | ||
1383 | if(carry==0) { | ||
1384 | if(j!=0) { | ||
1385 | carry = 1; iq[i] = 0x100- j; | ||
1386 | } | ||
1387 | } else iq[i] = 0xff - j; | ||
1388 | } | ||
1389 | if(q0>0) { /* rare case: chance is 1 in 12 */ | ||
1390 | switch(q0) { | ||
1391 | case 1: | ||
1392 | iq[jz-1] &= 0x7f; break; | ||
1393 | case 2: | ||
1394 | iq[jz-1] &= 0x3f; break; | ||
1395 | } | ||
1396 | } | ||
1397 | if(ih==2) { | ||
1398 | z = one - z; | ||
1399 | if(carry!=0) z -= rb_scalbnf(one,q0); | ||
1400 | } | ||
1401 | } | ||
1402 | |||
1403 | /* check if recomputation is needed */ | ||
1404 | if(z==zero) { | ||
1405 | j = 0; | ||
1406 | for (i=jz-1;i>=jk;i--) j |= iq[i]; | ||
1407 | if(j==0) { /* need recomputation */ | ||
1408 | for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */ | ||
1409 | |||
1410 | for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */ | ||
1411 | f[jx+i] = (float) ipio2[jv+i]; | ||
1412 | for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; | ||
1413 | q[i] = fw; | ||
1414 | } | ||
1415 | jz += k; | ||
1416 | goto recompute; | ||
1417 | } | ||
1418 | } | ||
1419 | |||
1420 | /* chop off zero terms */ | ||
1421 | if(z==(float)0.0) { | ||
1422 | jz -= 1; q0 -= 8; | ||
1423 | while(iq[jz]==0) { jz--; q0-=8;} | ||
1424 | } else { /* break z into 8-bit if necessary */ | ||
1425 | z = rb_scalbnf(z,-q0); | ||
1426 | if(z>=two8) { | ||
1427 | fw = (float)((int32_t)(twon8*z)); | ||
1428 | iq[jz] = (int32_t)(z-two8*fw); | ||
1429 | jz += 1; q0 += 8; | ||
1430 | iq[jz] = (int32_t) fw; | ||
1431 | } else iq[jz] = (int32_t) z ; | ||
1432 | } | ||
1433 | |||
1434 | /* convert integer "bit" chunk to floating-point value */ | ||
1435 | fw = rb_scalbnf(one,q0); | ||
1436 | for(i=jz;i>=0;i--) { | ||
1437 | q[i] = fw*(float)iq[i]; fw*=twon8; | ||
1438 | } | ||
1439 | |||
1440 | /* compute PIo2[0,...,jp]*q[jz,...,0] */ | ||
1441 | for(i=jz;i>=0;i--) { | ||
1442 | for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k]; | ||
1443 | fq[jz-i] = fw; | ||
1444 | } | ||
1445 | |||
1446 | /* compress fq[] into y[] */ | ||
1447 | switch(prec) { | ||
1448 | case 0: | ||
1449 | fw = 0.0; | ||
1450 | for (i=jz;i>=0;i--) fw += fq[i]; | ||
1451 | y[0] = (ih==0)? fw: -fw; | ||
1452 | break; | ||
1453 | case 1: | ||
1454 | case 2: | ||
1455 | fw = 0.0; | ||
1456 | for (i=jz;i>=0;i--) fw += fq[i]; | ||
1457 | y[0] = (ih==0)? fw: -fw; | ||
1458 | fw = fq[0]-fw; | ||
1459 | for (i=1;i<=jz;i++) fw += fq[i]; | ||
1460 | y[1] = (ih==0)? fw: -fw; | ||
1461 | break; | ||
1462 | case 3: /* painful */ | ||
1463 | for (i=jz;i>0;i--) { | ||
1464 | fw = fq[i-1]+fq[i]; | ||
1465 | fq[i] += fq[i-1]-fw; | ||
1466 | fq[i-1] = fw; | ||
1467 | } | ||
1468 | for (i=jz;i>1;i--) { | ||
1469 | fw = fq[i-1]+fq[i]; | ||
1470 | fq[i] += fq[i-1]-fw; | ||
1471 | fq[i-1] = fw; | ||
1472 | } | ||
1473 | for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; | ||
1474 | if(ih==0) { | ||
1475 | y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; | ||
1476 | } else { | ||
1477 | y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; | ||
1478 | } | ||
1479 | } | ||
1480 | return n&7; | ||
1481 | } | ||
1482 | |||
1483 | |||
1484 | static const int32_t two_over_pi[] = { | ||
1485 | 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC, | ||
1486 | 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62, | ||
1487 | 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63, | ||
1488 | 0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A, | ||
1489 | 0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09, | ||
1490 | 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29, | ||
1491 | 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44, | ||
1492 | 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41, | ||
1493 | 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C, | ||
1494 | 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8, | ||
1495 | 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11, | ||
1496 | 0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF, | ||
1497 | 0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E, | ||
1498 | 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5, | ||
1499 | 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92, | ||
1500 | 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08, | ||
1501 | 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0, | ||
1502 | 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3, | ||
1503 | 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85, | ||
1504 | 0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80, | ||
1505 | 0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA, | ||
1506 | 0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B, | ||
1507 | }; | ||
1508 | |||
1509 | static const int32_t npio2_hw[] = { | ||
1510 | 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00, | ||
1511 | 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00, | ||
1512 | 0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100, | ||
1513 | 0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00, | ||
1514 | 0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00, | ||
1515 | 0x4242c700, 0x42490f00 | ||
1516 | }; | ||
1517 | |||
1518 | /* | ||
1519 | * invpio2: 24 bits of 2/pi | ||
1520 | * pio2_1: first 17 bit of pi/2 | ||
1521 | * pio2_1t: pi/2 - pio2_1 | ||
1522 | * pio2_2: second 17 bit of pi/2 | ||
1523 | * pio2_2t: pi/2 - (pio2_1+pio2_2) | ||
1524 | * pio2_3: third 17 bit of pi/2 | ||
1525 | * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) | ||
1526 | */ | ||
1527 | |||
1528 | static const float | ||
1529 | half = 5.0000000000e-01, /* 0x3f000000 */ | ||
1530 | invpio2 = 6.3661980629e-01, /* 0x3f22f984 */ | ||
1531 | pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */ | ||
1532 | pio2_1t = 1.0804334124e-05, /* 0x37354443 */ | ||
1533 | pio2_2 = 1.0804273188e-05, /* 0x37354400 */ | ||
1534 | pio2_2t = 6.0770999344e-11, /* 0x2e85a308 */ | ||
1535 | pio2_3 = 6.0770943833e-11, /* 0x2e85a300 */ | ||
1536 | pio2_3t = 6.1232342629e-17; /* 0x248d3132 */ | ||
1537 | |||
1538 | int32_t rem_pio2(float x, float *y) | ||
1539 | { | ||
1540 | float z,w,t,r,fn; | ||
1541 | float tx[3]; | ||
1542 | int32_t e0,i,j,nx,n,ix,hx; | ||
1543 | |||
1544 | GET_FLOAT_WORD(hx,x); | ||
1545 | ix = hx&0x7fffffff; | ||
1546 | if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */ | ||
1547 | {y[0] = x; y[1] = 0; return 0;} | ||
1548 | if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */ | ||
1549 | if(hx>0) { | ||
1550 | z = x - pio2_1; | ||
1551 | if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */ | ||
1552 | y[0] = z - pio2_1t; | ||
1553 | y[1] = (z-y[0])-pio2_1t; | ||
1554 | } else { /* near pi/2, use 24+24+24 bit pi */ | ||
1555 | z -= pio2_2; | ||
1556 | y[0] = z - pio2_2t; | ||
1557 | y[1] = (z-y[0])-pio2_2t; | ||
1558 | } | ||
1559 | return 1; | ||
1560 | } else { /* negative x */ | ||
1561 | z = x + pio2_1; | ||
1562 | if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */ | ||
1563 | y[0] = z + pio2_1t; | ||
1564 | y[1] = (z-y[0])+pio2_1t; | ||
1565 | } else { /* near pi/2, use 24+24+24 bit pi */ | ||
1566 | z += pio2_2; | ||
1567 | y[0] = z + pio2_2t; | ||
1568 | y[1] = (z-y[0])+pio2_2t; | ||
1569 | } | ||
1570 | return -1; | ||
1571 | } | ||
1572 | } | ||
1573 | if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */ | ||
1574 | t = rb_fabs(x); | ||
1575 | n = (int32_t) (t*invpio2+half); | ||
1576 | fn = (float)n; | ||
1577 | r = t-fn*pio2_1; | ||
1578 | w = fn*pio2_1t; /* 1st round good to 40 bit */ | ||
1579 | if(n<32&&(int32_t)(ix&0xffffff00)!=npio2_hw[n-1]) { | ||
1580 | y[0] = r-w; /* quick check no cancellation */ | ||
1581 | } else { | ||
1582 | uint32_t high; | ||
1583 | j = ix>>23; | ||
1584 | y[0] = r-w; | ||
1585 | GET_FLOAT_WORD(high,y[0]); | ||
1586 | i = j-((high>>23)&0xff); | ||
1587 | if(i>8) { /* 2nd iteration needed, good to 57 */ | ||
1588 | t = r; | ||
1589 | w = fn*pio2_2; | ||
1590 | r = t-w; | ||
1591 | w = fn*pio2_2t-((t-r)-w); | ||
1592 | y[0] = r-w; | ||
1593 | GET_FLOAT_WORD(high,y[0]); | ||
1594 | i = j-((high>>23)&0xff); | ||
1595 | if(i>25) { /* 3rd iteration need, 74 bits acc */ | ||
1596 | t = r; /* will cover all possible cases */ | ||
1597 | w = fn*pio2_3; | ||
1598 | r = t-w; | ||
1599 | w = fn*pio2_3t-((t-r)-w); | ||
1600 | y[0] = r-w; | ||
1601 | } | ||
1602 | } | ||
1603 | } | ||
1604 | y[1] = (r-y[0])-w; | ||
1605 | if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} | ||
1606 | else return n; | ||
1607 | } | ||
1608 | /* | ||
1609 | * all other (large) arguments | ||
1610 | */ | ||
1611 | if(ix>=0x7f800000) { /* x is inf or NaN */ | ||
1612 | y[0]=y[1]=x-x; return 0; | ||
1613 | } | ||
1614 | /* set z = scalbn(|x|,ilogb(x)-7) */ | ||
1615 | e0 = (ix>>23)-134; /* e0 = ilogb(z)-7; */ | ||
1616 | SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23))); | ||
1617 | for(i=0;i<2;i++) { | ||
1618 | tx[i] = (float)((int32_t)(z)); | ||
1619 | z = (z-tx[i])*two8; | ||
1620 | } | ||
1621 | tx[2] = z; | ||
1622 | nx = 3; | ||
1623 | while(tx[nx-1]==zero) nx--; /* skip zero term */ | ||
1624 | n = kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); | ||
1625 | if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} | ||
1626 | return n; | ||
1627 | } | ||
1628 | |||
872 | float rb_tan(float x) | 1629 | float rb_tan(float x) |
873 | { | 1630 | { |
874 | return rb_sin(x) / rb_cos(x); | 1631 | float y[2],z=0.0; |
1632 | int32_t n, ix; | ||
1633 | |||
1634 | GET_FLOAT_WORD(ix,x); | ||
1635 | |||
1636 | /* |x| ~< pi/4 */ | ||
1637 | ix &= 0x7fffffff; | ||
1638 | if(ix <= 0x3f490fda) return kernel_tan(x,z,1); | ||
1639 | |||
1640 | /* tan(Inf or NaN) is NaN */ | ||
1641 | else if (ix>=0x7f800000) return x-x; /* NaN */ | ||
1642 | |||
1643 | /* argument reduction needed */ | ||
1644 | else { | ||
1645 | n = rem_pio2(x,y); | ||
1646 | return kernel_tan(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even | ||
1647 | -1 -- n odd */ | ||
1648 | } | ||
875 | } | 1649 | } |
876 | 1650 | ||
877 | 1651 | ||