|
| 1 | +//========================================================== |
| 2 | + |
| 3 | +// DATE : 2023-01-06 |
| 4 | +//========================================================== |
| 5 | + |
| 6 | +#include <math.h> |
| 7 | +#include <endian.h> |
| 8 | + |
| 9 | +namespace ModuleBase |
| 10 | +{ |
| 11 | +namespace libm |
| 12 | +{ |
| 13 | + |
| 14 | +typedef int int4; |
| 15 | +typedef union { unsigned int u[2]; int4 i[2]; double x; double d; } mynumber; |
| 16 | + |
| 17 | +#define max(x, y) (((y) > (x)) ? (y) : (x)) |
| 18 | +#define min(x, y) (((y) < (x)) ? (y) : (x)) |
| 19 | + |
| 20 | +#if (__BYTE_ORDER == __BIG_ENDIAN) |
| 21 | + |
| 22 | +#define HIGH_HALF 0 |
| 23 | +#define LOW_HALF 1 |
| 24 | + |
| 25 | +static const mynumber |
| 26 | + |
| 27 | +/**/ t576 = {{0x63f00000, 0x00000000}}, /* 2 ^ 576 */ |
| 28 | +/**/ tm600 = {{0x1a700000, 0x00000000}}, /* 2 ^- 600 */ |
| 29 | +/**/ tm24 = {{0x3e700000, 0x00000000}}, /* 2 ^- 24 */ |
| 30 | +/**/ big = {{0x43380000, 0x00000000}}, /* 6755399441055744 */ |
| 31 | +/**/ big1 = {{0x43580000, 0x00000000}}, /* 27021597764222976 */ |
| 32 | +/**/ hp0 = {{0x3FF921FB, 0x54442D18}} ,/* 1.5707963267948966 */ |
| 33 | +/**/ hp1 = {{0x3C91A626, 0x33145C07}} ,/* 6.123233995736766e-17 */ |
| 34 | +/**/ mp1 = {{0x3FF921FB, 0x58000000}}, /* 1.5707963407039642 */ |
| 35 | +/**/ mp2 = {{0xBE4DDE97, 0x40000000}}; /*-1.3909067675399456e-08 */ |
| 36 | + |
| 37 | +#endif |
| 38 | + |
| 39 | +#if (__BYTE_ORDER == __LITTLE_ENDIAN) |
| 40 | + |
| 41 | +#define HIGH_HALF 1 |
| 42 | +#define LOW_HALF 0 |
| 43 | + |
| 44 | +static const mynumber |
| 45 | + |
| 46 | +/**/ t576 = {{0x00000000, 0x63f00000}}, /* 2 ^ 576 */ |
| 47 | +/**/ tm600 = {{0x00000000, 0x1a700000}}, /* 2 ^- 600 */ |
| 48 | +/**/ tm24 = {{0x00000000, 0x3e700000}}, /* 2 ^- 24 */ |
| 49 | +/**/ big = {{0x00000000, 0x43380000}}, /* 6755399441055744 */ |
| 50 | +/**/ big1 = {{0x00000000, 0x43580000}}, /* 27021597764222976 */ |
| 51 | +/**/ hp0 = {{0x54442D18, 0x3FF921FB}}, /* 1.5707963267948966 */ |
| 52 | +/**/ hp1 = {{0x33145C07, 0x3C91A626}}, /* 6.123233995736766e-17 */ |
| 53 | +/**/ mp1 = {{0x58000000, 0x3FF921FB}}, /* 1.5707963407039642 */ |
| 54 | +/**/ mp2 = {{0x40000000, 0xBE4DDE97}}; /*-1.3909067675399456e-08 */ |
| 55 | + |
| 56 | +#endif |
| 57 | + |
| 58 | +static const double toverp[75] = { /* 2/ PI base 24*/ |
| 59 | + 10680707.0, 7228996.0, 1387004.0, 2578385.0, 16069853.0, |
| 60 | + 12639074.0, 9804092.0, 4427841.0, 16666979.0, 11263675.0, |
| 61 | + 12935607.0, 2387514.0, 4345298.0, 14681673.0, 3074569.0, |
| 62 | + 13734428.0, 16653803.0, 1880361.0, 10960616.0, 8533493.0, |
| 63 | + 3062596.0, 8710556.0, 7349940.0, 6258241.0, 3772886.0, |
| 64 | + 3769171.0, 3798172.0, 8675211.0, 12450088.0, 3874808.0, |
| 65 | + 9961438.0, 366607.0, 15675153.0, 9132554.0, 7151469.0, |
| 66 | + 3571407.0, 2607881.0, 12013382.0, 4155038.0, 6285869.0, |
| 67 | + 7677882.0, 13102053.0, 15825725.0, 473591.0, 9065106.0, |
| 68 | + 15363067.0, 6271263.0, 9264392.0, 5636912.0, 4652155.0, |
| 69 | + 7056368.0, 13614112.0, 10155062.0, 1944035.0, 9527646.0, |
| 70 | + 15080200.0, 6658437.0, 6231200.0, 6832269.0, 16767104.0, |
| 71 | + 5075751.0, 3212806.0, 1398474.0, 7579849.0, 6349435.0, |
| 72 | + 12618859.0, 4703257.0, 12806093.0, 14477321.0, 2786137.0, |
| 73 | + 12875403.0, 9837734.0, 14528324.0, 13719321.0, 343717.0 }; |
| 74 | + |
| 75 | +/* CN = 1+2**27 = '41a0000002000000' IEEE double format. Use it to split a |
| 76 | + double for better accuracy. */ |
| 77 | +#define CN 134217729.0 |
| 78 | +static const double split = CN; /* 2^27 + 1 */ |
| 79 | + |
| 80 | +/*******************************************************************/ |
| 81 | +/* Routine branred() performs range reduction of a double number */ |
| 82 | +/* x into Double length number a+aa,such that */ |
| 83 | +/* x=n*pi/2+(a+aa), abs(a+aa)<pi/4, n=0,+-1,+-2,.... */ |
| 84 | +/* Routine return integer (n mod 4) */ |
| 85 | +/*******************************************************************/ |
| 86 | +int |
| 87 | +__branred(double x, double *a, double *aa) |
| 88 | +{ |
| 89 | + int i,k; |
| 90 | + mynumber u,gor; |
| 91 | + double r[6],s,t,sum,b,bb,sum1,sum2,b1,bb1,b2,bb2,x1,x2,t1,t2; |
| 92 | + |
| 93 | + x*=tm600.x; |
| 94 | + t=x*split; /* split x to two numbers */ |
| 95 | + x1=t-(t-x); |
| 96 | + x2=x-x1; |
| 97 | + sum=0; |
| 98 | + u.x = x1; |
| 99 | + k = (u.i[HIGH_HALF]>>20)&2047; |
| 100 | + k = (k-450)/24; |
| 101 | + if (k<0) |
| 102 | + k=0; |
| 103 | + gor.x = t576.x; |
| 104 | + gor.i[HIGH_HALF] -= ((k*24)<<20); |
| 105 | + for (i=0;i<6;i++) |
| 106 | + { r[i] = x1*toverp[k+i]*gor.x; gor.x *= tm24.x; } |
| 107 | + for (i=0;i<3;i++) { |
| 108 | + s=(r[i]+big.x)-big.x; |
| 109 | + sum+=s; |
| 110 | + r[i]-=s; |
| 111 | + } |
| 112 | + t=0; |
| 113 | + for (i=0;i<6;i++) |
| 114 | + t+=r[5-i]; |
| 115 | + bb=(((((r[0]-t)+r[1])+r[2])+r[3])+r[4])+r[5]; |
| 116 | + s=(t+big.x)-big.x; |
| 117 | + sum+=s; |
| 118 | + t-=s; |
| 119 | + b=t+bb; |
| 120 | + bb=(t-b)+bb; |
| 121 | + s=(sum+big1.x)-big1.x; |
| 122 | + sum-=s; |
| 123 | + b1=b; |
| 124 | + bb1=bb; |
| 125 | + sum1=sum; |
| 126 | + sum=0; |
| 127 | + |
| 128 | + u.x = x2; |
| 129 | + k = (u.i[HIGH_HALF]>>20)&2047; |
| 130 | + k = (k-450)/24; |
| 131 | + if (k<0) |
| 132 | + k=0; |
| 133 | + gor.x = t576.x; |
| 134 | + gor.i[HIGH_HALF] -= ((k*24)<<20); |
| 135 | + for (i=0;i<6;i++) |
| 136 | + { r[i] = x2*toverp[k+i]*gor.x; gor.x *= tm24.x; } |
| 137 | + for (i=0;i<3;i++) { |
| 138 | + s=(r[i]+big.x)-big.x; |
| 139 | + sum+=s; |
| 140 | + r[i]-=s; |
| 141 | + } |
| 142 | + t=0; |
| 143 | + for (i=0;i<6;i++) |
| 144 | + t+=r[5-i]; |
| 145 | + bb=(((((r[0]-t)+r[1])+r[2])+r[3])+r[4])+r[5]; |
| 146 | + s=(t+big.x)-big.x; |
| 147 | + sum+=s; |
| 148 | + t-=s; |
| 149 | + b=t+bb; |
| 150 | + bb=(t-b)+bb; |
| 151 | + s=(sum+big1.x)-big1.x; |
| 152 | + sum-=s; |
| 153 | + |
| 154 | + b2=b; |
| 155 | + bb2=bb; |
| 156 | + sum2=sum; |
| 157 | + |
| 158 | + sum=sum1+sum2; |
| 159 | + b=b1+b2; |
| 160 | + bb = (fabs(b1)>fabs(b2))? (b1-b)+b2 : (b2-b)+b1; |
| 161 | + if (b > 0.5) |
| 162 | + {b-=1.0; sum+=1.0;} |
| 163 | + else if (b < -0.5) |
| 164 | + {b+=1.0; sum-=1.0;} |
| 165 | + s=b+(bb+bb1+bb2); |
| 166 | + t=((b-s)+bb)+(bb1+bb2); |
| 167 | + b=s*split; |
| 168 | + t1=b-(b-s); |
| 169 | + t2=s-t1; |
| 170 | + b=s*hp0.x; |
| 171 | + bb=(((t1*mp1.x-b)+t1*mp2.x)+t2*mp1.x)+(t2*mp2.x+s*hp1.x+t*hp0.x); |
| 172 | + s=b+bb; |
| 173 | + t=(b-s)+bb; |
| 174 | + *a=s; |
| 175 | + *aa=t; |
| 176 | + return ((int) sum)&3; /* return quater of unit circle */ |
| 177 | +} |
| 178 | + |
| 179 | +}; |
| 180 | +}; |
0 commit comments