usr/src/lib/libm/common/complex/k_atan2l.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   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  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23  */
  24 /*
  25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28
  29 #include "libm.h"               /* __k_atan2l */
  30 #include "complex_wrapper.h"
  31
  32 #if defined(__sparc)
  33 #define HALF(x)  ((int *) &x)[3] = 0; ((int *) &x)[2] &= 0xfe000000
  34 #elif defined(__x86)
  35 #define HALF(x)  ((int *) &x)[0] = 0
  36 #endif
  37
  38 /*
  39  * long double __k_atan2l(long double y, long double x, long double *e)
  40  *
  41  * Compute atan2l with error terms.
  42  *
  43  * Important formula:
  44  *                  3       5
  45  *                 x       x
  46  * atan(x) = x - ----- + ----- - ...    (for x <= 1)
  47  *                 3       5
  48  *
  49  *           pi     1     1
  50  *         = --- - --- + --- - ...      (for x > 1)
  51  *                         3
  52  *            2     x    3x
  53  *
  54  * Arg(x + y i) = sign(y) * atan2(|y|, x)
  55  *              = sign(y) * atan(|y|/x)           (for x > 0)
  56  *                sign(y) * (PI - atan(|y|/|x|))  (for x < 0)
  57  * Thus if x >> y (IEEE double: EXP(x) - EXP(y) >= 60):
  58  *      1. (x > 0): atan2(y,x) ~ y/x
  59  *      2. (x < 0): atan2(y,x) ~ sign(y) (PI - |y/x|))
  60  * Otherwise if x << y:
  61  *      atan2(y,x) ~ sign(y)*PI/2 - x/y
  62  *
  63  * __k_atan2l call static functions mx_polyl, mx_atanl
  64  */
  65
  66
  67 /*
  68  * (void) mx_polyl (long double *z, long double *a, long double *e, int n)
  69  * return
  70  *      e = a  + z*(a  + z*(a  + ... z*(a  + e)...))
  71  *           0       2       4           2n
  72  * Note:
  73  * 1.   e and coefficient ai are represented by two long double numbers.
  74  *      For e, the first one contain the leading 53 bits (30 for x86 exteneded)
  75  *      and the second one contain the remaining 113 bits (64 for x86 extended).
  76  *      For ai, the first one contian the leading 53 bits (or 30 for x86)
  77  *      rounded, and the second is the remaining 113 bits (or 64 for x86).
  78  * 2.   z is an array of three doubles.
  79  *      z[0] :  the rounded value of Z (the intended value of z)
  80  *      z[1] :  the leading 32 (or 56) bits of Z rounded
  81  *      z[2] :  the remaining 113 (or 64) bits of Z
  82  *              Note that z[0] = z[1]+z[2] rounded.
  83  *
  84  */
  85
  86 static void
  87 mx_polyl(const long double *z, const long double *a, long double *e, int n) {
  88         long double r, s, t, p_h, p_l, z_h, z_l, p, w;
  89         int i;
  90         n = n + n;
  91         p = e[0] + a[n];
  92         p_l = a[n + 1];
  93         w = p; HALF(w);
  94         p_h = w;
  95         p   = a[n - 2] + z[0] * p;
  96         z_h = z[1]; z_l = z[2];
  97         p_l += e[0] - (p_h - a[n]);
  98
  99         for (i = n - 2; i >= 2; i -= 2) {
 100
 101         /* compute p = ai + z * p */
 102                 t = z_h * p_h;
 103                 s = z[0] * p_l + p_h * z_l;
 104                 w = p;  HALF(w);
 105                 p_h = w;
 106                 s  += a[i + 1];
 107                 r   = t - (p_h - a[i]);
 108                 p   = a[i - 2] + z[0] * p;
 109                 p_l = r + s;
 110         }
 111         w = p; HALF(w);
 112         e[0] = w;
 113         t   = z_h * p_h;
 114         s   = z[0] * p_l + p_h * z_l;
 115         r   = t - (e[0] - a[0]);
 116         e[1] = r + s;
 117 }
 118
 119 /*
 120  * Table of constants for atan from 0.125 to 8
 121  *      0.125 -- 0x3ffc0000  --- (increment at bit 12)
 122  *               0x3ffc1000
 123  *               0x3ffc2000
 124  *      ...     ...
 125  *               0x4001f000
 126  *      8.000 -- 0x40020000      (total: 97)
 127  */
 128
 129 static const long double TBL_atan_hil[] = {
 130 #if defined(__sparc)
 131 1.2435499454676143503135484916387102416568e-01L,
 132 1.3203976161463874927468440652656953226250e-01L,
 133 1.3970887428916364518336777673909505681607e-01L,
 134 1.4736148108865163560980276039684551821066e-01L,
 135 1.5499674192394098230371437493349219133371e-01L,
 136 1.6261382859794857537364156376155780062019e-01L,
 137 1.7021192528547440449049660709976171369543e-01L,
 138 1.7779022899267607079662479921582468899456e-01L,
 139 1.8534794999569476488602596122854464667261e-01L,
 140 1.9288431225797466419705871069022730349878e-01L,
 141 2.0039855382587851465394578503437838446153e-01L,
 142 2.0788992720226299360533498310299432475629e-01L,
 143 2.1535769969773804802445962716648964165745e-01L,
 144 2.2280115375939451577103212214043255525024e-01L,
 145 2.3021958727684373024017095967980299065551e-01L,
 146 2.3761231386547125247388363432563777919892e-01L,
 147 2.4497866312686415417208248121127580641959e-01L,
 148 2.5962962940825753102994644318397190560106e-01L,
 149 2.7416745111965879759937189834217578592444e-01L,
 150 2.8858736189407739562361141995821834504332e-01L,
 151 3.0288486837497140556055609450555821812277e-01L,
 152 3.1705575320914700980901557667446732975852e-01L,
 153 3.3109607670413209494433878775694455421259e-01L,
 154 3.4500217720710510886768128690005168408290e-01L,
 155 3.5877067027057222039592006392646052215363e-01L,
 156 3.7239844667675422192365503828370182641413e-01L,
 157 3.8588266939807377589769548460723139638186e-01L,
 158 3.9922076957525256561471669615886476491104e-01L,
 159 4.1241044159738730689979128966712694260920e-01L,
 160 4.2544963737004228954226360518079233013817e-01L,
 161 4.3833655985795780544561604921477130895882e-01L,
 162 4.5106965598852347637563925728219344073798e-01L,
 163 4.6364760900080611621425623146121439713344e-01L,
 164 4.8833395105640552386716496074706484459644e-01L,
 165 5.1238946031073770666660102058425923805558e-01L,
 166 5.3581123796046370026908506870769144698471e-01L,
 167 5.5859931534356243597150821640166122875873e-01L,
 168 5.8075635356767039920327447500150082375122e-01L,
 169 6.0228734613496418168212269420423291922459e-01L,
 170 6.2319932993406593099247534906037459367793e-01L,
 171 6.4350110879328438680280922871732260447265e-01L,
 172 6.6320299270609325536325431023827583417226e-01L,
 173 6.8231655487474807825642998171115298784729e-01L,
 174 7.0085440788445017245795128178675127318623e-01L,
 175 7.1882999962162450541701415152590469891043e-01L,
 176 7.3625742898142813174283527108914662479274e-01L,
 177 7.5315128096219438952473937026902888600575e-01L,
 178 7.6952648040565826040682003598565401726598e-01L,
 179 7.8539816339744830961566084581987569936977e-01L,
 180 8.1569192331622341102146083874564582672284e-01L,
 181 8.4415398611317100251784414827164746738632e-01L,
 182 8.7090345707565295314017311259781407291650e-01L,
 183 8.9605538457134395617480071802993779546602e-01L,
 184 9.1971960535041681722860345482108940969311e-01L,
 185 9.4200004037946366473793717053459362115891e-01L,
 186 9.6299433068093620181519583599709989677298e-01L,
 187 9.8279372324732906798571061101466603762572e-01L,
 188 1.0014831356942347329183295953014374896343e+00L,
 189 1.0191413442663497346383429170230636212354e+00L,
 190 1.0358412530088001765846944703254440735476e+00L,
 191 1.0516502125483736674598673120862999026920e+00L,
 192 1.0666303653157435630791763474202799086015e+00L,
 193 1.0808390005411683108871567292171997859003e+00L,
 194 1.0943289073211899198927883146102352763033e+00L,
 195 1.1071487177940905030170654601785370497543e+00L,
 196 1.1309537439791604464709335155363277560026e+00L,
 197 1.1525719972156675180401498626127514672834e+00L,
 198 1.1722738811284763866005949441337046006865e+00L,
 199 1.1902899496825317329277337748293182803384e+00L,
 200 1.2068173702852525303955115800565576625682e+00L,
 201 1.2220253232109896370417417439225704120294e+00L,
 202 1.2360594894780819419094519711090786146210e+00L,
 203 1.2490457723982544258299170772810900483550e+00L,
 204 1.2610933822524404193139408812473357640124e+00L,
 205 1.2722973952087173412961937498224805746463e+00L,
 206 1.2827408797442707473628852511364955164072e+00L,
 207 1.2924966677897852679030914214070816723528e+00L,
 208 1.3016288340091961438047858503666855024453e+00L,
 209 1.3101939350475556342564376891719053437537e+00L,
 210 1.3182420510168370498593302023271363040427e+00L,
 211 1.3258176636680324650592392104284756886164e+00L,
 212 1.3397056595989995393283037525895557850243e+00L,
 213 1.3521273809209546571891479413898127598774e+00L,
 214 1.3633001003596939542892985278250991560269e+00L,
 215 1.3734007669450158608612719264449610604836e+00L,
 216 1.3825748214901258580599674177685685163955e+00L,
 217 1.3909428270024183486427686943836432395486e+00L,
 218 1.3986055122719575950126700816114282727858e+00L,
 219 1.4056476493802697809521934019958080664406e+00L,
 220 1.4121410646084952153676136718584890852820e+00L,
 221 1.4181469983996314594038603039700988632607e+00L,
 222 1.4237179714064941189018190466107297108905e+00L,
 223 1.4288992721907326964184700745371984001389e+00L,
 224 1.4337301524847089866404719096698873880264e+00L,
 225 1.4382447944982225979614042479354816039669e+00L,
 226 1.4424730991091018200252920599377291810352e+00L,
 227 1.4464413322481351841999668424758803866109e+00L,
 228 #elif defined(__x86)
 229 1.243549945356789976358413696289e-01L,  1.320397615781985223293304443359e-01L,
 230 1.397088742814958095550537109375e-01L,  1.473614810383878648281097412109e-01L,
 231 1.549967419123277068138122558594e-01L,  1.626138285500928759574890136719e-01L,
 232 1.702119252295233309268951416016e-01L,  1.777902289759367704391479492188e-01L,
 233 1.853479499695822596549987792969e-01L,  1.928843122441321611404418945312e-01L,
 234 2.003985538030974566936492919922e-01L,  2.078899272019043564796447753906e-01L,
 235 2.153576996643096208572387695312e-01L,  2.228011537226848304271697998047e-01L,
 236 2.302195872762240469455718994141e-01L,  2.376123138237744569778442382812e-01L,
 237 2.449786631041206419467926025391e-01L,  2.596296293195337057113647460938e-01L,
 238 2.741674510762095451354980468750e-01L,  2.885873618070036172866821289062e-01L,
 239 3.028848683461546897888183593750e-01L,  3.170557531993836164474487304688e-01L,
 240 3.310960766393691301345825195312e-01L,  3.450021771714091300964355468750e-01L,
 241 3.587706702528521418571472167969e-01L,  3.723984466632828116416931152344e-01L,
 242 3.858826693613082170486450195312e-01L,  3.992207695264369249343872070312e-01L,
 243 4.124104415532201528549194335938e-01L,  4.254496373469009995460510253906e-01L,
 244 4.383365598041564226150512695312e-01L,  4.510696559445932507514953613281e-01L,
 245 4.636476089945062994956970214844e-01L,  4.883339509833604097366333007812e-01L,
 246 5.123894601128995418548583984375e-01L,  5.358112377580255270004272460938e-01L,
 247 5.585993151180446147918701171875e-01L,  5.807563534472137689590454101562e-01L,
 248 6.022873460315167903900146484375e-01L,  6.231993297114968299865722656250e-01L,
 249 6.435011087451130151748657226562e-01L,  6.632029926404356956481933593750e-01L,
 250 6.823165547102689743041992187500e-01L,  7.008544078562408685684204101562e-01L,
 251 7.188299994450062513351440429688e-01L,  7.362574287690222263336181640625e-01L,
 252 7.531512808054685592651367187500e-01L,  7.695264802314341068267822265625e-01L,
 253 7.853981633670628070831298828125e-01L,  8.156919232569634914398193359375e-01L,
 254 8.441539860796183347702026367188e-01L,  8.709034570492804050445556640625e-01L,
 255 8.960553845390677452087402343750e-01L,  9.197196052409708499908447265625e-01L,
 256 9.420000403188169002532958984375e-01L,  9.629943305626511573791503906250e-01L,
 257 9.827937232330441474914550781250e-01L,  1.001483135391026735305786132812e+00L,
 258 1.019141343887895345687866210938e+00L,  1.035841252654790878295898437500e+00L,
 259 1.051650212146341800689697265625e+00L,  1.066630364861339330673217773438e+00L,
 260 1.080839000176638364791870117188e+00L,  1.094328907318413257598876953125e+00L,
 261 1.107148717623203992843627929688e+00L,  1.130953743588179349899291992188e+00L,
 262 1.152571997139602899551391601562e+00L,  1.172273880802094936370849609375e+00L,
 263 1.190289949532598257064819335938e+00L,  1.206817369908094406127929687500e+00L,
 264 1.222025323193520307540893554688e+00L,  1.236059489194303750991821289062e+00L,
 265 1.249045772012323141098022460938e+00L,  1.261093381792306900024414062500e+00L,
 266 1.272297394927591085433959960938e+00L,  1.282740879338234663009643554688e+00L,
 267 1.292496667709201574325561523438e+00L,  1.301628833636641502380371093750e+00L,
 268 1.310193934943526983261108398438e+00L,  1.318242050707340240478515625000e+00L,
 269 1.325817663222551345825195312500e+00L,  1.339705659542232751846313476562e+00L,
 270 1.352127380669116973876953125000e+00L,  1.363300099968910217285156250000e+00L,
 271 1.373400766868144273757934570312e+00L,  1.382574821356683969497680664062e+00L,
 272 1.390942826867103576660156250000e+00L,  1.398605511989444494247436523438e+00L,
 273 1.405647648964077234268188476562e+00L,  1.412141064181923866271972656250e+00L,
 274 1.418146998155862092971801757812e+00L,  1.423717970959842205047607421875e+00L,
 275 1.428899271879345178604125976562e+00L,  1.433730152435600757598876953125e+00L,
 276 1.438244794495403766632080078125e+00L,  1.442473099101334810256958007812e+00L,
 277 1.446441331878304481506347656250e+00L,
 278 #endif
 279 };
 280 static const long double TBL_atan_lol[] = {
 281 #if defined(__sparc)
 282 1.4074869197628063802317202820414310039556e-36L,
 283 -4.9596961594739925555730439437999675295505e-36L,
 284 8.9527745625194648873931213446361849472788e-36L,
 285 1.1880437423207895718180765843544965589427e-35L,
 286 -2.7810278112045145378425375128234365381448e-37L,
 287 1.4797220377023800327295536234315147262387e-36L,
 288 -4.2169561400548198732870384801849639863829e-36L,
 289 7.2431229666913484649930323656316023494680e-36L,
 290 -2.1573430089839170299895679353790663182462e-36L,
 291 -9.9515745405126723554452367298128605186305e-36L,
 292 -3.9065558992324838181617569730397882363067e-36L,
 293 5.5260292271793726813211980664661124518807e-36L,
 294 8.8415722215914321807682254318036452043689e-36L,
 295 -8.1767728791586179254193323628285599800711e-36L,
 296 -1.3344123034656142243797113823028330070762e-36L,
 297 -4.4927331207813382908930733924681325892188e-36L,
 298 4.4945511471812490393201824336762495687730e-36L,
 299 -1.6688081504279223555776724459648440567274e-35L,
 300 1.5629757586107955769461086568937329684113e-35L,
 301 -2.2389835563308078552507970385331510848109e-35L,
 302 -4.8312321745547311551870450671182151367050e-36L,
 303 -1.4336172352905832876958926610980698844309e-35L,
 304 -8.7440181998899932802989174170960593316080e-36L,
 305 5.9284636008529837445780360785464550143016e-36L,
 306 -2.2376651248436241276061055295043514993630e-35L,
 307 6.0745837599336105414280310756677442136480e-36L,
 308 1.5372187110451949677792344762029967023093e-35L,
 309 2.0976068056751156241657121582478790247159e-35L,
 310 -5.5623956405495438060726862202622807523700e-36L,
 311 1.9697366707832471841858411934897351901523e-35L,
 312 2.1070311964479488509034733639424887543697e-35L,
 313 -2.3027356362982001602256518510854229844561e-35L,
 314 4.8950964225733349266861843522029764772843e-36L,
 315 -7.2380143477794458213872723050820253166391e-36L,
 316 1.6365648865703614031637443396049568858105e-35L,
 317 -3.9885811958234530793729129919803234197399e-35L,
 318 4.1587722120912613510417783923227421336929e-35L,
 319 3.8347421454556472153684687377337135027394e-35L,
 320 -9.2251178933638721723515896465489002497864e-36L,
 321 1.4094619690455989526175736741854656192178e-36L,
 322 3.3568857805472235270612851425810803679451e-35L,
 323 3.9090991055522552395018106803232118803401e-35L,
 324 5.2956416979654208140521862707297033857956e-36L,
 325 -5.0960846819945514367847063923662507136721e-36L,
 326 -4.4959014425277615858329680393918315204998e-35L,
 327 3.8039226544551634266566857615962609653834e-35L,
 328 -4.4056522872895512108308642196611689657618e-36L,
 329 1.6025024192482161076223807753425619076948e-36L,
 330 2.1679525325309452561992610065108380635264e-35L,
 331 1.9844038013515422125715362925736754104066e-35L,
 332 3.9139619471799746834505227353568432457241e-35L,
 333 2.1113443807975453505518453436799561854730e-35L,
 334 3.1558557277444692755039816944392770185432e-35L,
 335 1.6295044520355461408265585619500238335614e-35L,
 336 -3.5087245209270305856151230356171213582305e-35L,
 337 2.9041041864282855679591055270946117300088e-35L,
 338 -2.3128843453818356590931995209806627233282e-35L,
 339 -7.7124923181471578439967973820714857839953e-35L,
 340 2.7539027829886922429092063590445808781462e-35L,
 341 -9.4500899453181308951084545990839335972452e-35L,
 342 -7.3061755302032092337594946001641651543473e-35L,
 343 -4.1736144813953752193952770157406952602798e-35L,
 344 3.4369948356256407045344855262863733571105e-35L,
 345 -6.3790243492298090907302084924276831116460e-35L,
 346 -9.6842943816353261291004127866079538980649e-36L,
 347 4.8746757539138870909275958326700072821615e-35L,
 348 -8.7533886477084190884511601368582548254655e-35L,
 349 1.4284743992327918892692551138086727754845e-35L,
 350 5.7262776211073389542565625693479173445042e-35L,
 351 -3.2254883148780411245594822270747948565684e-35L,
 352 7.8853548190609877325965525252380833808405e-35L,
 353 8.4081736739037194097515038365370730251333e-35L,
 354 7.4722870357563683815078242981933587273670e-35L,
 355 7.9977202825793435289434813600890494256112e-36L,
 356 -8.0577840773362139054848492346292673645405e-35L,
 357 1.4217746753670583065490040209048757624336e-35L,
 358 1.2232486914221205004109743560319090913328e-35L,
 359 8.9696055070830036447361957217943988339065e-35L,
 360 -3.1480394435081884410686066739846269858951e-35L,
 361 -5.0927146040715345013240642517608928352977e-35L,
 362 -5.7431997715924136568133859432702789493569e-35L,
 363 -4.3920451405083770279099766080476485439987e-35L,
 364 9.1106753984907715563018666776308759323326e-35L,
 365 -3.7032569014272841009512400773061537538358e-35L,
 366 8.8167419429746714276909825405131416764489e-35L,
 367 -3.8389341696028352503752312861740895209678e-36L,
 368 -3.3462959341960891546340895508017603408404e-35L,
 369 -3.9212626776786074383916188498955828634947e-35L,
 370 -7.8340397396377867255864494568594088378648e-35L,
 371 7.4681018632456986520600640340627309824469e-35L,
 372 8.9110918618956918451135594876165314884113e-35L,
 373 3.9418160632271890530431797145664308529115e-35L,
 374 -4.1048114088580104820193435638327617443913e-35L,
 375 -2.3165419451582153326383944756220900454330e-35L,
 376 -1.8428312581525319409399330203703211113843e-35L,
 377 7.1477316546709482345411712017906842769961e-35L,
 378 2.9914501578435874662153637707016094237004e-35L,
 379 #elif defined(__x86)
 380 1.108243739551347953496477557317e-11L,  3.644022694535396219063202730280e-11L,
 381 7.667835628314065801595065768845e-12L,  5.026377078169301918590803009109e-11L,
 382 1.161327548990211907411719105561e-11L,  4.785569941615255008968280209991e-11L,
 383 5.595107356360146549819920947848e-11L,  1.673930035747684999707469623769e-11L,
 384 2.611250523102718193166964451527e-11L,  1.384250305661681615897729354721e-11L,
 385 2.278105796029649304219088055497e-11L,  3.586371256902077123693302823191e-13L,
 386 3.342842716722085763523965049902e-11L,  3.670968534386232233574504707347e-11L,
 387 6.196832945990602657404893210974e-13L,  4.169679549603939604438777470618e-11L,
 388 2.274351222528987867221331091414e-11L,  8.872382531858169709022188891298e-11L,
 389 4.344925246387385146717580155420e-11L,  8.707377833692929105196832265348e-11L,
 390 2.881671577173773513055821329154e-11L,  9.763393361566846205717315422347e-12L,
 391 6.476296480975626822569454546857e-11L,  3.569597877124574002505169001136e-11L,
 392 1.772007853877284712958549977698e-11L,  1.347141028196192304932683248872e-11L,
 393 3.676555884905046507598141175404e-11L,  4.881564068032948912761478588710e-11L,
 394 4.416715404487185607337693704681e-11L,  2.314128999621257979016734983553e-11L,
 395 5.380138283056477968352133002913e-11L,  4.393022562414389595406841771063e-11L,
 396 6.299816718559209976839402028537e-12L,  7.304511413053165996581483735843e-11L,
 397 1.978381648117426221467592544212e-10L,  2.024381732686578226139414070989e-10L,
 398 2.255178211796380992141612703464e-10L,  1.204566302442290648452508620986e-10L,
 399 1.034473912921080457667329099995e-10L,  2.225691010059030834353745950874e-10L,
 400 4.817137162794350606107263804151e-11L,  6.565755971506095086327587326326e-11L,
 401 1.644791039522307629611529931429e-10L,  2.820930388953087163050126809014e-11L,
 402 1.766182540818701085571546539514e-10L,  2.124059054092171070266466628320e-10L,
 403 1.567258302596026515190288816001e-10L,  1.742241535800378094231540188685e-10L,
 404 3.038550253253096300737572104929e-11L,  5.925991958164150280814584656688e-11L,
 405 3.355266774764151155289750652594e-11L,  2.637254809561744853531409402995e-11L,
 406 3.227621096606048365493782702458e-11L,  1.094459672377587282585894259882e-10L,
 407 6.064676448464127209709358607166e-11L,  1.182850444360454453720999258140e-10L,
 408 1.428492049425553288966601449688e-11L,  3.032079976125434624889374125094e-10L,
 409 3.784543889504767060855636487744e-10L,  3.540092982887960328254439790467e-10L,
 410 4.020318667701700464612998296302e-10L,  4.544042324059585739827798668654e-10L,
 411 3.645299460952866120296998202703e-10L,  2.776662293911361485235212513020e-12L,
 412 1.708865101734375304910370400700e-10L,  3.909810965716415233488278047493e-10L,
 413 7.606461848875826105025137974947e-11L,  3.263814502297453347587046149712e-10L,
 414 1.499334758629144388918183376012e-10L,  3.771581242675818925565576303133e-10L,
 415 1.746932950084818923507049088298e-11L,  2.837781909176306820465786987027e-10L,
 416 3.859312847318946163435901230778e-10L,  4.601335192895268187473357720101e-10L,
 417 2.811262558622337888849804940684e-10L,  4.060360843532416964489955306249e-10L,
 418 8.058369357752989796958168458531e-11L,  3.725546414244147566166855921414e-10L,
 419 1.040286509953292907344053122733e-10L,  3.094968093808145773271362531155e-10L,
 420 4.454811192340438979284756311844e-10L,  5.676678748199027602705574110388e-11L,
 421 2.518376833121948163898128509842e-10L,  3.907837370041422778250991189943e-10L,
 422 7.687158710333735613246114865100e-11L,  1.334418885622867537060685125566e-10L,
 423 1.353147719826124443836432060856e-10L,  2.825131007652335581739282335732e-10L,
 424 4.161925466840049254333079881002e-10L,  4.265713490956410156084891599630e-10L,
 425 2.437693664320585461575989523716e-10L,  4.466519138542116247357297503086e-10L,
 426 3.113875178143440979746983590908e-10L,  4.910822904159495654488736486097e-11L,
 427 2.818831329324169810481585538618e-12L,  7.767009768334052125229252512543e-12L,
 428 3.698307026936191862258804165254e-10L,
 429 #endif
 430 };
 431
 432 /*
 433  * mx_atanl(x, err)
 434  * Table look-up algorithm
 435  * By K.C. Ng, March 9, 1989
 436  *
 437  * Algorithm.
 438  *
 439  * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
 440  * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
 441  * error (relative)
 442  *      |(atan(x)-poly1(x))/x|<= 2^-140
 443  *
 444  * and use poly2(x) to approximate atan(x) for x in [0,1/65] with
 445  * error
 446  *      |atan(x)-poly2(x)|<= 2^-143.7
 447  *
 448  * Here poly1 and poly2 are odd polynomial with the following form:
 449  *              x + x^3*(a1+x^2*(a2+...))
 450  *
 451  * (0). Purge off Inf and NaN and 0
 452  * (1). Reduce x to positive by atan(x) = -atan(-x).
 453  * (2). For x <= 1/8, use
 454  *      (2.1) if x < 2^(-prec/2), atan(x) =  x  with inexact flag raised
 455  *      (2.2) Otherwise
 456  *              atan(x) = poly1(x)
 457  * (3). For x >= 8 then (prec = 78)
 458  *      (3.1) if x >= 2^prec,     atan(x) = atan(inf) - pio2_lo
 459  *      (3.2) if x >= 2^(prec/3), atan(x) = atan(inf) - 1/x
 460  *      (3.3) if x >  65,         atan(x) = atan(inf) - poly2(1/x)
 461  *      (3.4) Otherwise,          atan(x) = atan(inf) - poly1(1/x)
 462  *
 463  * (4). Now x is in (0.125, 8)
 464  *      Find y that match x to 4.5 bit after binary (easy).
 465  *      If iy is the high word of y, then
 466  *              single : j = (iy - 0x3e000000) >> 19
 467  *              double : j = (iy - 0x3fc00000) >> 16
 468  *              quad   : j = (iy - 0x3ffc0000) >> 12
 469  *
 470  *      Let s = (x-y)/(1+x*y). Then
 471  *      atan(x) = atan(y) + poly1(s)
 472  *              = _TBL_atan_hi[j] + (_TBL_atan_lo[j] + poly2(s) )
 473  *
 474  *      Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
 475  *
 476  */
 477
 478 /*
 479  * p[0] - p[16] for atan(x) =
 480  *              x + x^3*(p1+x^2*(p2+...))
 481  */
 482 static const long double pe[] = {
 483         1.0L,
 484         0.0L,
 485 #if defined(__sparc)
 486         -0.33333333333333332870740406406184774823L,
 487         -4.62592926927148558508441072595508240609e-18L,
 488         0.19999999999999999722444243843710864894L,
 489         2.77555756156289124602047010782090464486e-18L,
 490         -0.14285714285714285615158658515611023176L,
 491         -9.91270557700756738621231719241800559409e-19L,
 492 #elif defined(__x86)
 493         -0.33333333325572311878204345703125L,
 494         -7.76102145512898763020833333192787755766644373e-11L,
 495         0.19999999995343387126922607421875L,
 496         4.65661287307739257812498949613909375938538636e-11L,
 497         -0.142857142840512096881866455078125L,
 498         -1.66307602609906877787419703858463013035681375e-11L,
 499 #endif
 500 };
 501
 502 static const long double p[] = { /* p[0] - p[16] */
 503         1.0L,
 504         -3.33333333333333333333333333333333333319278775586e-0001L,
 505         1.99999999999999999999999999999999894961390937601e-0001L,
 506         -1.42857142857142857142857142856866970385846301312e-0001L,
 507         1.11111111111111111111111110742899094415954427738e-0001L,
 508         -9.09090909090909090909087972707015549231951421806e-0002L,
 509         7.69230769230769230767699003016385628597359717046e-0002L,
 510         -6.66666666666666666113842763495291228025226575259e-0002L,
 511         5.88235294117646915706902204947653640091126695962e-0002L,
 512         -5.26315789473657016886225044679594035524579379810e-0002L,
 513         4.76190476186633969331771169790375592681525481267e-0002L,
 514         -4.34782608290146274616081389793141896576997370161e-0002L,
 515         3.99999968161267722260103962788865225205057218988e-0002L,
 516         -3.70368536844778256320786172745225703228683638328e-0002L,
 517         3.44752320396524479494062858284036892703898522150e-0002L,
 518         -3.20491216046653214683721787776813360591233428081e-0002L,
 519         2.67632651033434456758550618122802167256870856514e-0002L,
 520 };
 521
 522 /* q[0] - q[9] */
 523 static const long double qe[] = {
 524         1.0L,
 525         0.0L,
 526 #if defined(__sparc)
 527         -0.33333333333333332870740406406184774823486804962158203125L,
 528         -4.625929269271485585069345465471207312531868714634217630e-18L,
 529         0.19999999999999999722444243843710864894092082977294921875L,
 530         2.7755575615628864268260553912956813621977220359134667560e-18L,
 531 #elif defined(__x86)
 532         -0.33333333325572311878204345703125L,
 533         -7.76102145512898763020833333042135150927893e-11L,
 534         0.19999999995343387126922607421875L,
 535         4.656612873077392578124507576697622106863058e-11L,
 536 #endif
 537 };
 538
 539 static const long double q[] = {        /* q[0] - q[9] */
 540         -3.33333333333333333333333333333333333304213515094e-0001L,
 541         1.99999999999999999999999999999995075766976221077e-0001L,
 542         -1.42857142857142857142857142570379604317921113079e-0001L,
 543         1.11111111111111111111102923861900979127978214077e-0001L,
 544         -9.09090909090909089586854075816999506863320031460e-0002L,
 545         7.69230769230756334929213246003824644696974730368e-0002L,
 546         -6.66666666589192433974402013508912138168133579856e-0002L,
 547         5.88235013696778007696800252045588307023299350858e-0002L,
 548         -5.25754959898164576495303840687699583228444695685e-0002L,
 549 };
 550
 551 static const long double
 552 two8700  = 9.140338438955067659002088492701e+2618L, /* 2^8700 */
 553 twom8700 = 1.094051392821643668051436593760e-2619L, /* 2^-8700 */
 554 one   = 1.0L,
 555 zero  = 0.0L,
 556 pi    = 3.1415926535897932384626433832795028841971693993751L,
 557 pio2  = 1.57079632679489661923132169163975144209858469968755L,
 558 pio4  = 0.785398163397448309615660845819875721049292349843776L,
 559 pi3o4 = 2.356194490192344928846982537459627163147877049531329L,
 560 #if defined(__sparc)
 561 pi_lo    = 8.67181013012378102479704402604335196876232e-35L,
 562 pio2_lo  = 4.33590506506189051239852201302167598438116e-35L,
 563 pio4_lo  = 2.16795253253094525619926100651083799219058e-35L,
 564 pi3o4_lo = 6.50385759759283576859778301953251397657174e-35L;
 565 #elif defined(__x86)
 566 pi_lo    = -5.01655761266833202355732708e-20L,
 567 pio2_lo  = -2.50827880633416601177866354e-20L,
 568 pio4_lo  = -1.25413940316708300588933177e-20L,
 569 pi3o4_lo = -9.18342907192877118770525931e-20L;
 570 #endif
 571
 572 static long double
 573 mx_atanl(long double x, long double *err) {
 574         long double y, z, r, s, t, w, s_h, s_l, x_h, x_l, zz[3], ee[2], z_h,
 575                 z_l, r_h, r_l, u, v;
 576         int ix, iy, hx, i, j;
 577         float fx;
 578
 579         hx = HI_XWORD(x);
 580         ix = hx & (~0x80000000);
 581
 582         /* for |x| < 1/8 */
 583         if (ix < 0x3ffc0000) {
 584                 if (ix < 0x3ff30000) {  /* when |x| < 2**-12 */
 585                         if (ix < 0x3fc60000) {  /* if |x| < 2**-prec/2 */
 586                                 *err = (long double) ((int) x);
 587                                 return (x);
 588                         }
 589                         z = x * x;
 590                         t = q[8];
 591                         for (i = 7; i >= 0; i--) t = q[i] + z * t;
 592                         t *= x * z;
 593                         r = x + t;
 594                         *err = t - (r - x);
 595                         return (r);
 596                 }
 597                 z = x * x;
 598
 599                 /* use long double precision at p4 and on */
 600                 t = p[16];
 601                 for (i = 15; i >= 4; i--) t = p[i] + z * t;
 602                 ee[0] = z * t;
 603
 604                 x_h = x;  HALF(x_h);
 605                 z_h = z;  HALF(z_h);
 606                 x_l = x - x_h;
 607                 z_l = (x_h * x_h - z_h);
 608                 zz[0] = z;
 609                 zz[1] = z_h;
 610                 zz[2] = z_l + x_l * (x + x_h);
 611
 612                 /* compute (1+z*(p1+z*(p2+z*(p3+e)))) */
 613
 614                 mx_polyl(zz, pe, ee, 3);
 615
 616                 /* finally x*(1+z*(p1+...)) */
 617                 r = x_h * ee[0];
 618                 t = x * ee[1] + x_l * ee[0];
 619                 s = t + r;
 620                 *err = t - (s - r);
 621                 return (s);
 622         }
 623         /* for |x| >= 8.0 */
 624         if (ix >= 0x40020000) { /* x >=  8 */
 625                 x = fabsl(x);
 626                 if (ix >= 0x402e0000) { /* x >=  2**47 */
 627                         if (ix >= 0x408b0000) { /* x >=  2**140 */
 628                                 y = -pio2_lo;
 629                         } else
 630                                 y = one / x - pio2_lo;
 631                         if (hx >= 0) {
 632                                 t = pio2 - y;
 633                                 *err = -(y - (pio2 - t));
 634                         } else {
 635                                 t = y - pio2;
 636                                 *err = y - (pio2 + t);
 637                         }
 638                         return (t);
 639                 } else {
 640                         /* compute r = 1/x */
 641                         r = one / x;
 642                         z = r * r;
 643                         x_h = x; HALF(x_h);
 644                         r_h = r; HALF(r_h);
 645                         z_h = z; HALF(z_h);
 646                         r_l = r * ((x_h - x) * r_h - (x_h * r_h - one));
 647                         z_l = (r_h * r_h - z_h);
 648                         zz[0] = z;
 649                         zz[1] = z_h;
 650                         zz[2] = z_l + r_l * (r + r_h);
 651                         if (ix < 0x40050400) {  /* 8 <  x <  65 */
 652                                 /* use double precision at p4 and on */
 653                                 t = p[16];
 654                                 for (i = 15; i >= 4; i--) t = p[i] + z * t;
 655                                 ee[0] = z * t;
 656                                 /* compute (1+z*(p1+z*(p2+z*(p3+e)))) */
 657                                 mx_polyl(zz, pe, ee, 3);
 658                         } else { /* x < 65 < 2**47 */
 659                                 /* use long double at q3 and on */
 660                                 t = q[8];
 661                                 for (i = 7; i >= 2; i--) t = q[i] + z * t;
 662                                 ee[0] = z * t;
 663                                 /* compute (1+z*(q1+z*(q2+e))) */
 664                                 mx_polyl(zz, qe, ee, 2);
 665                         }
 666                         /* pio2 - r*(1+...) */
 667                         v = r_h * ee[0];
 668                         t = pio2_lo - (r * ee[1] + r_l * ee[0]);
 669                         if (hx >= 0) {
 670                                 s = pio2 - v;
 671                                 t -= (v - (pio2 - s));
 672                         } else {
 673                                 s = v - pio2;
 674                                 t = -(t - (v - (s + pio2)));
 675                         }
 676                         w = s + t;
 677                         *err = t - (w - s);
 678                         return (w);
 679                 }
 680         }
 681         /* now x is between 1/8 and 8 */
 682         iy = (ix + 0x00000800) & 0x7ffff000;
 683         j = (iy - 0x3ffc0000) >> 12;
 684         ((int *) &fx)[0] = 0x3e000000 + (j << 19);
 685         y = (long double) fx;
 686         x = fabsl(x);
 687
 688         w = (x - y);
 689         v = 1.0L / (one + x * y);
 690         s = w * v;
 691         z = s * s;
 692         /* use long double precision at q3 and on */
 693         t = q[8];
 694         for (i = 7; i >= 2; i--) t = q[i] + z * t;
 695         ee[0] = z * t;
 696         s_h = s; HALF(s_h);
 697         z_h = z; HALF(z_h);
 698         x_h = x; HALF(x_h);
 699         t = one + x * y; HALF(t);
 700         r = -((x_h - x) * y - (x_h * y - (t - one)));
 701         s_l = -v * (s_h * r - (w - s_h * t));
 702         z_l = (s_h * s_h - z_h);
 703         zz[0] = z;
 704         zz[1] = z_h;
 705         zz[2] = z_l + s_l * (s + s_h);
 706         /* compute (1+z*(q1+z*(q2+e))) by call mx_poly */
 707         mx_polyl(zz, qe, ee, 2);
 708         v = s_h * ee[0];
 709         t = TBL_atan_lol[j] + (s * ee[1] + s_l * ee[0]);
 710         u = TBL_atan_hil[j];
 711         s = u + v;
 712         t += (v - (s - u));
 713         w = s + t;
 714         *err = t - (w - s);
 715         if (hx < 0) {
 716                 w = -w;
 717                 *err = -*err;
 718         }
 719         return (w);
 720 }
 721
 722 long double
 723 __k_atan2l(long double y, long double x, long double *w) {
 724         long double t, xh, th, t1, t2, w1, w2;
 725         int ix, iy, hx, hy;
 726
 727         hy = HI_XWORD(y);
 728         hx = HI_XWORD(x);
 729         iy = hy & ~0x80000000;
 730         ix = hx & ~0x80000000;
 731
 732         *w = 0.0;
 733         if (ix >= 0x7fff0000 || iy >= 0x7fff0000) {     /* ignore inexact */
 734                 if (isnanl(x) || isnanl(y))
 735                         return (x * y);
 736                 else if (iy < 0x7fff0000) {
 737                         if (hx >= 0) {  /* ATAN2(+-finite, +inf) is +-0 */
 738                                 *w *= y;
 739                                 return (*w);
 740                         } else {        /* ATAN2(+-finite, -inf) is +-pi */
 741                                 *w = copysignl(pi_lo, y);
 742                                 return (copysignl(pi, y));
 743                         }
 744                 } else if (ix < 0x7fff0000) {
 745                                         /* ATAN2(+-inf, finite) is +-pi/2 */
 746                         *w = (hy >= 0)? pio2_lo : -pio2_lo;
 747                         return ((hy >= 0)? pio2 : -pio2);
 748                 } else if (hx > 0) {    /* ATAN2(+-INF,+INF) = +-pi/4 */
 749                         *w = (hy >= 0)? pio4_lo : -pio4_lo;
 750                         return ((hy >= 0)? pio4 : -pio4);
 751                 } else {                /* ATAN2(+-INF,-INF) = +-3pi/4 */
 752                         *w = (hy >= 0)? pi3o4_lo : -pi3o4_lo;
 753                         return ((hy >= 0)? pi3o4 : -pi3o4);
 754                 }
 755         } else if (x == zero || y == zero) {
 756                 if (y == zero) {
 757                         if (hx >= 0) /* ATAN2(+-0, +(0 <= x <= inf)) is +-0 */
 758                                 return (y);
 759                         else {  /* ATAN2(+-0, -(0 <= x <= inf)) is +-pi */
 760                                 *w = (hy >= 0)? pi_lo : -pi_lo;
 761                                 return ((hy >= 0)? pi : -pi);
 762                         }
 763                 } else { /* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2 */
 764                         *w = (hy >= 0)? pio2_lo : -pio2_lo;
 765                         return ((hy >= 0)? pio2 : -pio2);
 766                 }
 767         } else if (iy - ix > 0x00640000) { /* |x/y| < 2 ** -100 */
 768                 *w = (hy >= 0)? pio2_lo : -pio2_lo;
 769                 return ((hy >= 0)? pio2 : -pio2);
 770         } else if (ix - iy > 0x00640000) { /* |y/x| < 2 ** -100 */
 771                 if (hx < 0) {
 772                         *w = (hy >= 0)? pi_lo : -pi_lo;
 773                         return ((hy >= 0)? pi : -pi);
 774                 } else {
 775                         t = y / x;
 776                         th = t; HALF(th);
 777                         xh = x; HALF(xh);
 778                         t1 = (x - xh) * t + xh * (t - th);
 779                         t2 = y - xh * th;
 780                         *w = (t2 - t1) / x;
 781                         return (t);
 782                 }
 783         } else {
 784                 if (ix >= 0x5fff3000) {
 785                         x *= twom8700;
 786                         y *= twom8700;
 787                 } else if (ix < 0x203d0000) {
 788                         x *= two8700;
 789                         y *= two8700;
 790                 }
 791                 y = fabsl(y);
 792                 x = fabsl(x);
 793                 t = y / x;
 794                 th = t; HALF(th);
 795                 xh = x; HALF(xh);
 796                 t1 = (x - xh) * t + xh * (t - th);
 797                 t2 = y - xh * th;
 798                 w1 = mx_atanl(t, &w2);
 799                 w2 += (t2 - t1) / (x + y * t);
 800                 if (hx < 0) {
 801                         t1 = pi - w1;
 802                         t2 = pi - t1;
 803                         w2 = (pi_lo - w2) - (w1 - t2);
 804                         w1 = t1;
 805                 }
 806                 *w = (hy >= 0)? w2 : -w2;
 807                 return ((hy >= 0)? w1 : -w1);
 808         }
 809 }