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

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/complex/k_atan2.c
          +++ new/usr/src/lib/libm/common/complex/k_atan2.c
↓ open down ↓ 10 lines elided ↑ open up ↑
  11   11   * and limitations under the License.
  12   12   *
  13   13   * When distributing Covered Code, include this CDDL HEADER in each
  14   14   * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15   15   * If applicable, add the following below this CDDL HEADER, with the
  16   16   * fields enclosed by brackets "[]" replaced with your own identifying
  17   17   * information: Portions Copyright [yyyy] [name of copyright owner]
  18   18   *
  19   19   * CDDL HEADER END
  20   20   */
       21 +
  21   22  /*
  22   23   * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23   24   */
       25 +
  24   26  /*
  25   27   * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26   28   * Use is subject to license terms.
  27   29   */
  28   30  
  29      -#include "libm.h"               /* __k_atan2 */
       31 +#include "libm.h"                       /* __k_atan2 */
  30   32  #include "complex_wrapper.h"
  31   33  
  32   34  /*
  33   35   * double __k_atan2(double y, double x, double *e)
  34   36   *
  35   37   * Compute atan2 with error terms.
  36   38   *
  37   39   * Important formula:
  38   40   *                  3       5
  39   41   *                 x       x
↓ open down ↓ 22 lines elided ↑ open up ↑
  62   64   * return
  63   65   *      e = a  + z*(a  + z*(a  + ... z*(a  + e)...))
  64   66   *           0       2       4           2n
  65   67   * Note:
  66   68   * 1.   e and coefficient ai are represented by two double numbers.
  67   69   *      For e, the first one contain the leading 24 bits rounded, and the
  68   70   *      second one contain the remaining 53 bits (total 77 bits accuracy).
  69   71   *      For ai, the first one contian the leading 53 bits rounded, and the
  70   72   *      second is the remaining 53 bits (total 106 bits accuracy).
  71   73   * 2.   z is an array of three doubles.
  72      - *      z[0] :  the rounded value of Z (the intended value of z)
  73      - *      z[1] :  the leading 24 bits of Z rounded
  74      - *      z[2] :  the remaining 53 bits of Z
       74 + *      z[0] :  the rounded value of Z (the intended value of z)
       75 + *      z[1] :  the leading 24 bits of Z rounded
       76 + *      z[2] :  the remaining 53 bits of Z
  75   77   *              Note that z[0] = z[1]+z[2] rounded.
  76   78   *
  77   79   */
  78      -
  79   80  static void
  80      -mx_poly(const double *z, const double *a, double *e, int n) {
       81 +mx_poly(const double *z, const double *a, double *e, int n)
       82 +{
  81   83          double r, s, t, p_h, p_l, z_h, z_l, p;
  82   84          int i;
  83   85  
  84   86          n = n + n;
  85   87          p = e[0] + a[n];
  86   88          p_l = a[n + 1];
  87      -        p_h = (double) ((float) p);
  88      -        p   = a[n - 2] + z[0] * p;
  89      -        z_h = z[1]; z_l = z[2];
       89 +        p_h = (double)((float)p);
       90 +        p = a[n - 2] + z[0] * p;
       91 +        z_h = z[1];
       92 +        z_l = z[2];
  90   93          p_l += e[0] - (p_h - a[n]);
  91   94  
  92   95          for (i = n - 2; i >= 2; i -= 2) {
  93      -        /* compute p = ai + z * p */
  94      -                t   = z_h * p_h;
  95      -                s   = z[0] * p_l + p_h * z_l;
  96      -                p_h = (double) ((float) p);
  97      -                s  += a[i + 1];
  98      -                r   = t - (p_h - a[i]);
  99      -                p   = a[i - 2] + z[0] * p;
       96 +                /* compute p = ai + z * p */
       97 +                t = z_h * p_h;
       98 +                s = z[0] * p_l + p_h * z_l;
       99 +                p_h = (double)((float)p);
      100 +                s += a[i + 1];
      101 +                r = t - (p_h - a[i]);
      102 +                p = a[i - 2] + z[0] * p;
 100  103                  p_l = r + s;
 101  104          }
 102      -        e[0] = (double)((float) p);
 103      -        t   = z_h * p_h;
 104      -        s   = z[0] * p_l + p_h * z_l;
 105      -        r   = t - (e[0] - a[0]);
      105 +
      106 +        e[0] = (double)((float)p);
      107 +        t = z_h * p_h;
      108 +        s = z[0] * p_l + p_h * z_l;
      109 +        r = t - (e[0] - a[0]);
 106  110          e[1] = r + s;
 107  111  }
 108  112  
 109  113  /*
 110  114   * Table of constants for atan from 0.125 to 8
 111  115   *      0.125 -- 0x3fc00000  --- (increment at bit 16)
 112  116   *               0x3fc10000
 113  117   *               0x3fc20000
 114  118   *      ...     ...
 115  119   *               0x401f0000
 116  120   *      8.000 -- 0x40200000      (total: 97)
 117  121   * By K.C. Ng, March 9, 1989
 118  122   */
 119  123  
 120  124  static const double TBL_atan_hi[] = {
 121      -1.243549945467614382e-01, 1.320397616146387620e-01, 1.397088742891636204e-01,
 122      -1.473614810886516302e-01, 1.549967419239409727e-01, 1.626138285979485676e-01,
 123      -1.702119252854744080e-01, 1.777902289926760471e-01, 1.853479499956947607e-01,
 124      -1.928843122579746439e-01, 2.003985538258785115e-01, 2.078899272022629863e-01,
 125      -2.153576996977380476e-01, 2.228011537593945213e-01, 2.302195872768437179e-01,
 126      -2.376123138654712419e-01, 2.449786631268641435e-01, 2.596296294082575118e-01,
 127      -2.741674511196587893e-01, 2.885873618940774099e-01, 3.028848683749714166e-01,
 128      -3.170557532091470287e-01, 3.310960767041321029e-01, 3.450021772071051318e-01,
 129      -3.587706702705721895e-01, 3.723984466767542023e-01, 3.858826693980737521e-01,
 130      -3.992207695752525431e-01, 4.124104415973872673e-01, 4.254496373700422662e-01,
 131      -4.383365598579578304e-01, 4.510696559885234436e-01, 4.636476090008060935e-01,
 132      -4.883339510564055352e-01, 5.123894603107377321e-01, 5.358112379604637043e-01,
 133      -5.585993153435624414e-01, 5.807563535676704136e-01, 6.022873461349641522e-01,
 134      -6.231993299340659043e-01, 6.435011087932843710e-01, 6.632029927060932861e-01,
 135      -6.823165548747480713e-01, 7.008544078844501923e-01, 7.188299996216245269e-01,
 136      -7.362574289814280970e-01, 7.531512809621944138e-01, 7.695264804056582975e-01,
 137      -7.853981633974482790e-01, 8.156919233162234217e-01, 8.441539861131710509e-01,
 138      -8.709034570756529758e-01, 8.960553845713439269e-01, 9.197196053504168578e-01,
 139      -9.420000403794636101e-01, 9.629943306809362058e-01, 9.827937232473290541e-01,
 140      -1.001483135694234639e+00, 1.019141344266349725e+00, 1.035841253008800145e+00,
 141      -1.051650212548373764e+00, 1.066630365315743623e+00, 1.080839000541168327e+00,
 142      -1.094328907321189925e+00, 1.107148717794090409e+00, 1.130953743979160375e+00,
 143      -1.152571997215667610e+00, 1.172273881128476303e+00, 1.190289949682531656e+00,
 144      -1.206817370285252489e+00, 1.222025323210989667e+00, 1.236059489478081863e+00,
 145      -1.249045772398254428e+00, 1.261093382252440387e+00, 1.272297395208717319e+00,
 146      -1.282740879744270757e+00, 1.292496667789785336e+00, 1.301628834009196156e+00,
 147      -1.310193935047555547e+00, 1.318242051016837113e+00, 1.325817663668032553e+00,
 148      -1.339705659598999565e+00, 1.352127380920954636e+00, 1.363300100359693845e+00,
 149      -1.373400766945015894e+00, 1.382574821490125894e+00, 1.390942827002418447e+00,
 150      -1.398605512271957618e+00, 1.405647649380269870e+00, 1.412141064608495311e+00,
 151      -1.418146998399631542e+00, 1.423717971406494032e+00, 1.428899272190732761e+00,
 152      -1.433730152484709031e+00, 1.438244794498222623e+00, 1.442473099109101931e+00,
 153      -1.446441332248135092e+00,
      125 +        1.243549945467614382e-01, 1.320397616146387620e-01,
      126 +        1.397088742891636204e-01, 1.473614810886516302e-01,
      127 +        1.549967419239409727e-01, 1.626138285979485676e-01,
      128 +        1.702119252854744080e-01, 1.777902289926760471e-01,
      129 +        1.853479499956947607e-01, 1.928843122579746439e-01,
      130 +        2.003985538258785115e-01, 2.078899272022629863e-01,
      131 +        2.153576996977380476e-01, 2.228011537593945213e-01,
      132 +        2.302195872768437179e-01, 2.376123138654712419e-01,
      133 +        2.449786631268641435e-01, 2.596296294082575118e-01,
      134 +        2.741674511196587893e-01, 2.885873618940774099e-01,
      135 +        3.028848683749714166e-01, 3.170557532091470287e-01,
      136 +        3.310960767041321029e-01, 3.450021772071051318e-01,
      137 +        3.587706702705721895e-01, 3.723984466767542023e-01,
      138 +        3.858826693980737521e-01, 3.992207695752525431e-01,
      139 +        4.124104415973872673e-01, 4.254496373700422662e-01,
      140 +        4.383365598579578304e-01, 4.510696559885234436e-01,
      141 +        4.636476090008060935e-01, 4.883339510564055352e-01,
      142 +        5.123894603107377321e-01, 5.358112379604637043e-01,
      143 +        5.585993153435624414e-01, 5.807563535676704136e-01,
      144 +        6.022873461349641522e-01, 6.231993299340659043e-01,
      145 +        6.435011087932843710e-01, 6.632029927060932861e-01,
      146 +        6.823165548747480713e-01, 7.008544078844501923e-01,
      147 +        7.188299996216245269e-01, 7.362574289814280970e-01,
      148 +        7.531512809621944138e-01, 7.695264804056582975e-01,
      149 +        7.853981633974482790e-01, 8.156919233162234217e-01,
      150 +        8.441539861131710509e-01, 8.709034570756529758e-01,
      151 +        8.960553845713439269e-01, 9.197196053504168578e-01,
      152 +        9.420000403794636101e-01, 9.629943306809362058e-01,
      153 +        9.827937232473290541e-01, 1.001483135694234639e+00,
      154 +        1.019141344266349725e+00, 1.035841253008800145e+00,
      155 +        1.051650212548373764e+00, 1.066630365315743623e+00,
      156 +        1.080839000541168327e+00, 1.094328907321189925e+00,
      157 +        1.107148717794090409e+00, 1.130953743979160375e+00,
      158 +        1.152571997215667610e+00, 1.172273881128476303e+00,
      159 +        1.190289949682531656e+00, 1.206817370285252489e+00,
      160 +        1.222025323210989667e+00, 1.236059489478081863e+00,
      161 +        1.249045772398254428e+00, 1.261093382252440387e+00,
      162 +        1.272297395208717319e+00, 1.282740879744270757e+00,
      163 +        1.292496667789785336e+00, 1.301628834009196156e+00,
      164 +        1.310193935047555547e+00, 1.318242051016837113e+00,
      165 +        1.325817663668032553e+00, 1.339705659598999565e+00,
      166 +        1.352127380920954636e+00, 1.363300100359693845e+00,
      167 +        1.373400766945015894e+00, 1.382574821490125894e+00,
      168 +        1.390942827002418447e+00, 1.398605512271957618e+00,
      169 +        1.405647649380269870e+00, 1.412141064608495311e+00,
      170 +        1.418146998399631542e+00, 1.423717971406494032e+00,
      171 +        1.428899272190732761e+00, 1.433730152484709031e+00,
      172 +        1.438244794498222623e+00, 1.442473099109101931e+00,
      173 +        1.446441332248135092e+00,
 154  174  };
 155  175  
 156  176  static const double TBL_atan_lo[] = {
 157      --3.125324142453938311e-18, -1.276925400709959526e-17, 2.479758919089733066e-17,
 158      -5.409599147666297957e-18, 9.585415594114323829e-18, 7.784470643106252464e-18,
 159      --3.541164079802125137e-18, 2.372599351477449041e-17, 4.180692268843078977e-18,
 160      -2.034098543938166622e-17, 3.139954287184449286e-18, 7.333160666520898500e-18,
 161      -4.738160130078732886e-19, -5.498822172446843173e-18, 1.231340452914270316e-17,
 162      -1.058231431371112987e-17, 1.069875561873445139e-17, 1.923875492461530410e-17,
 163      -8.261353575163771936e-18, -1.428369957377257085e-17, -1.101082790300136900e-17,
 164      --1.893928924292642146e-17, -7.952610375793798701e-18, -2.293880475557830393e-17,
 165      -3.088733564861919217e-17, 1.961231150484565340e-17, 2.378822732491940868e-17,
 166      -2.246598105617042065e-17, 3.963462895355093301e-17, 2.331553074189288466e-17,
 167      --2.494277030626540909e-17, 3.280735600183735558e-17, 2.269877745296168709e-17,
 168      --1.137323618932958456e-17, -2.546278147285580353e-17, -4.063795683482557497e-18,
 169      --5.455630548591626394e-18, -1.441464378193066908e-17, 2.950430737228402307e-17,
 170      -2.672403885140095079e-17, 1.583478505144428617e-17, -3.076054864429649001e-17,
 171      -6.943223671560007740e-18, -1.987626234335816123e-17, -2.147838844445698302e-17,
 172      -3.473937648299456719e-17, -2.425693465918206812e-17, -3.704991905602721293e-17,
 173      -3.061616997868383018e-17, -1.071456562778743077e-17, -4.841337011934916763e-17,
 174      --2.269823590747287052e-17, 2.923876285774304890e-17, -4.057439412852767923e-17,
 175      -5.460837485846687627e-17, -3.986660595210752445e-18, 1.390331103123099845e-17,
 176      -9.438308023545392000e-17, 1.000401886936679889e-17, 3.194313981784503706e-17,
 177      --9.650564731467513515e-17, -5.956589637160374564e-17, -1.567632251135907253e-17,
 178      --5.490676155022364226e-18, 9.404471373566379412e-17, 7.123833804538446299e-17,
 179      --9.159738508900378819e-17, 8.385188614028674371e-17, 7.683333629842068806e-17,
 180      -4.172467638861439118e-17, -2.979162864892849274e-17, 7.879752739459421280e-17,
 181      --2.196203799612310905e-18, 3.242139621534960503e-17, 2.245875015034507026e-17,
 182      --9.283188754266129476e-18, -6.830804768926660334e-17, -1.236918499824626670e-17,
 183      -8.745413734780278834e-17, -6.319394031144676258e-17, -8.824429373951136321e-17,
 184      --2.599011860304134377e-17, 2.147674250751150961e-17, 1.093246171526936217e-16,
 185      --3.307710355769516504e-17, -3.561490438648230100e-17, -9.843712133488842595e-17,
 186      --2.324061182591627982e-17, -8.922630138234492386e-17, -9.573807110557223276e-17,
 187      --8.263883782511013632e-17, 8.721870922223967507e-17, -6.457134743238754385e-17,
 188      --4.396204466767636187e-17, -2.493019910264565554e-17, -1.105119435430315713e-16,
 189      -9.211323971545051565e-17,
      177 +        -3.125324142453938311e-18, -1.276925400709959526e-17,
      178 +        2.479758919089733066e-17, 5.409599147666297957e-18,
      179 +        9.585415594114323829e-18, 7.784470643106252464e-18,
      180 +        -3.541164079802125137e-18, 2.372599351477449041e-17,
      181 +        4.180692268843078977e-18, 2.034098543938166622e-17,
      182 +        3.139954287184449286e-18, 7.333160666520898500e-18,
      183 +        4.738160130078732886e-19, -5.498822172446843173e-18,
      184 +        1.231340452914270316e-17, 1.058231431371112987e-17,
      185 +        1.069875561873445139e-17, 1.923875492461530410e-17,
      186 +        8.261353575163771936e-18, -1.428369957377257085e-17,
      187 +        -1.101082790300136900e-17, -1.893928924292642146e-17,
      188 +        -7.952610375793798701e-18, -2.293880475557830393e-17,
      189 +        3.088733564861919217e-17, 1.961231150484565340e-17,
      190 +        2.378822732491940868e-17, 2.246598105617042065e-17,
      191 +        3.963462895355093301e-17, 2.331553074189288466e-17,
      192 +        -2.494277030626540909e-17, 3.280735600183735558e-17,
      193 +        2.269877745296168709e-17, -1.137323618932958456e-17,
      194 +        -2.546278147285580353e-17, -4.063795683482557497e-18,
      195 +        -5.455630548591626394e-18, -1.441464378193066908e-17,
      196 +        2.950430737228402307e-17, 2.672403885140095079e-17,
      197 +        1.583478505144428617e-17, -3.076054864429649001e-17,
      198 +        6.943223671560007740e-18, -1.987626234335816123e-17,
      199 +        -2.147838844445698302e-17, 3.473937648299456719e-17,
      200 +        -2.425693465918206812e-17, -3.704991905602721293e-17,
      201 +        3.061616997868383018e-17, -1.071456562778743077e-17,
      202 +        -4.841337011934916763e-17, -2.269823590747287052e-17,
      203 +        2.923876285774304890e-17, -4.057439412852767923e-17,
      204 +        5.460837485846687627e-17, -3.986660595210752445e-18,
      205 +        1.390331103123099845e-17, 9.438308023545392000e-17,
      206 +        1.000401886936679889e-17, 3.194313981784503706e-17,
      207 +        -9.650564731467513515e-17, -5.956589637160374564e-17,
      208 +        -1.567632251135907253e-17, -5.490676155022364226e-18,
      209 +        9.404471373566379412e-17, 7.123833804538446299e-17,
      210 +        -9.159738508900378819e-17, 8.385188614028674371e-17,
      211 +        7.683333629842068806e-17, 4.172467638861439118e-17,
      212 +        -2.979162864892849274e-17, 7.879752739459421280e-17,
      213 +        -2.196203799612310905e-18, 3.242139621534960503e-17,
      214 +        2.245875015034507026e-17, -9.283188754266129476e-18,
      215 +        -6.830804768926660334e-17, -1.236918499824626670e-17,
      216 +        8.745413734780278834e-17, -6.319394031144676258e-17,
      217 +        -8.824429373951136321e-17, -2.599011860304134377e-17,
      218 +        2.147674250751150961e-17, 1.093246171526936217e-16,
      219 +        -3.307710355769516504e-17, -3.561490438648230100e-17,
      220 +        -9.843712133488842595e-17, -2.324061182591627982e-17,
      221 +        -8.922630138234492386e-17, -9.573807110557223276e-17,
      222 +        -8.263883782511013632e-17, 8.721870922223967507e-17,
      223 +        -6.457134743238754385e-17, -4.396204466767636187e-17,
      224 +        -2.493019910264565554e-17, -1.105119435430315713e-16,
      225 +        9.211323971545051565e-17,
 190  226  };
 191  227  
 192  228  /*
 193  229   * mx_atan(x,err)
 194  230   * Table look-up algorithm
 195  231   * By K.C. Ng, March 9, 1989
 196  232   *
 197  233   * Algorithm.
 198  234   *
 199  235   * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
 200  236   * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
 201  237   * error (relative)
 202      - *      |(atan(x)-poly1(x))/x|<= 2^-83.41
      238 + *      |(atan(x)-poly1(x))/x|<= 2^-83.41
 203  239   *
 204  240   * and use poly2(x) to approximate atan(x) for x in [0,1/65] with
 205  241   * error
 206  242   *      |atan(x)-poly2(x)|<= 2^-86.8
 207  243   *
 208  244   * Here poly1 and poly2 are odd polynomial with the following form:
 209  245   *              x + x^3*(a1+x^2*(a2+...))
 210  246   *
 211  247   * (0). Purge off Inf and NaN and 0
 212  248   * (1). Reduce x to positive by atan(x) = -atan(-x).
↓ open down ↓ 15 lines elided ↑ open up ↑
 228  264   *              quad   : j = (iy - 0x3ffc0000) >> 12
 229  265   *
 230  266   *      Let s = (x-y)/(1+x*y). Then
 231  267   *      atan(x) = atan(y) + poly1(s)
 232  268   *              = _TBL_atan_hi[j] + (_TBL_atan_lo[j] + poly2(s) )
 233  269   *
 234  270   *      Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
 235  271   *
 236  272   */
 237  273  
 238      -#define P1 p[2]
 239      -#define P4 p[8]
 240      -#define P5 p[9]
 241      -#define P6 p[10]
 242      -#define P7 p[11]
 243      -#define P8 p[12]
 244      -#define P9 p[13]
      274 +#define P1              p[2]
      275 +#define P4              p[8]
      276 +#define P5              p[9]
      277 +#define P6              p[10]
      278 +#define P7              p[11]
      279 +#define P8              p[12]
      280 +#define P9              p[13]
      281 +
 245  282  static const double p[] = {
 246  283          1.0,
 247  284          0.0,
 248  285          -3.33333333333333314830e-01,    /* p1   = BFD55555 55555555 */
 249  286          -1.85030852238476921863e-17,    /* p1_l = BC755525 9783A49C */
 250  287          2.00000000000000011102e-01,     /* p2   = 3FC99999 9999999A */
 251  288          -1.27263196576150347368e-17,    /* p2_l = BC6D584B 0D874007 */
 252  289          -1.42857142857141405923e-01,    /* p3   = BFC24924 9249245E */
 253  290          -1.34258204847170493327e-17,    /* p3_l = BC6EF534 A112500D */
 254  291          1.11111111110486909803e-01,     /* p4   = 3FBC71C7 1C71176A */
 255  292          -9.09090907557387889470e-02,    /* p5   = BFB745D1 73B47A7D */
 256  293          7.69230541541713053189e-02,     /* p6   = 3FB3B13A B1E68DE6 */
 257  294          -6.66645815401964159097e-02,    /* p7   = BFB110EE 1584446A */
 258  295          5.87081768778560317279e-02,     /* p8   = 3FAE0EFF 87657733 */
 259  296          -4.90818147456113240690e-02,    /* p9   = BFA92140 6A524B5C */
 260  297  };
 261      -#define Q1 q[2]
 262      -#define Q3 q[6]
 263      -#define Q4 q[7]
 264      -#define Q5 q[8]
      298 +
      299 +#define Q1              q[2]
      300 +#define Q3              q[6]
      301 +#define Q4              q[7]
      302 +#define Q5              q[8]
      303 +
 265  304  static const double q[] = {
 266  305          1.0,
 267  306          0.0,
 268  307          -3.33333333333333314830e-01,    /* q1   = BFD55555 55555555 */
 269      -        -1.85022941571278638733e-17,    /* q1_l = BC7554E9 D20EFA66 */
 270      -        1.99999999999999927836e-01,     /* q2   = 3FC99999 99999997 */
 271      -        -1.28782564407438833398e-17,    /* q2_l = BC6DB1FB 17217417 */
 272      -        -1.42857142855492280642e-01,    /* q3   = BFC24924 92483C46 */
 273      -        1.11111097130183356096e-01,     /* q4   = 3FBC71C6 E06595CC */
 274      -        -9.08553303569109294013e-02,    /* q5   = BFB7424B 808CDA76 */
      308 +        -1.85022941571278638733e-17,            /* q1_l = BC7554E9 D20EFA66 */
      309 +        1.99999999999999927836e-01,             /* q2   = 3FC99999 99999997 */
      310 +        -1.28782564407438833398e-17,            /* q2_l = BC6DB1FB 17217417 */
      311 +        -1.42857142855492280642e-01,            /* q3   = BFC24924 92483C46 */
      312 +        1.11111097130183356096e-01,             /* q4   = 3FBC71C6 E06595CC */
      313 +        -9.08553303569109294013e-02,            /* q5   = BFB7424B 808CDA76 */
 275  314  };
 276      -static const double
 277      -one = 1.0,
 278      -pio2hi = 1.570796326794896558e+00,
 279      -pio2lo = 6.123233995736765886e-17;
      315 +
      316 +static const double one = 1.0,
      317 +        pio2hi = 1.570796326794896558e+00,
      318 +        pio2lo = 6.123233995736765886e-17;
 280  319  
 281  320  static double
 282      -mx_atan(double x, double *err) {
 283      -        double y, z, r, s, t, w, s_h, s_l, x_h, x_l, zz[3], ee[2], z_h,
 284      -                z_l, r_h, r_l, u, v;
      321 +mx_atan(double x, double *err)
      322 +{
      323 +        double y, z, r, s, t, w, s_h, s_l, x_h, x_l, zz[3], ee[2], z_h, z_l,
      324 +            r_h, r_l, u, v;
 285  325          int ix, iy, sign, j;
 286  326  
 287      -        ix = ((int *) &x)[HIWORD];
      327 +        ix = ((int *)&x)[HIWORD];
 288  328          sign = ix & 0x80000000;
 289  329          ix ^= sign;
 290  330  
 291  331          /* for |x| < 1/8 */
 292  332          if (ix < 0x3fc00000) {
 293      -                if (ix < 0x3f300000) {  /* when |x| < 2**-12 */
      333 +                if (ix < 0x3f300000) {          /* when |x| < 2**-12 */
 294  334                          if (ix < 0x3d800000) {  /* if |x| < 2**-39 */
 295      -                                *err = (double) ((int) x);
      335 +                                *err = (double)((int)x);
 296  336                                  return (x);
 297  337                          }
      338 +
 298  339                          z = x * x;
 299  340                          t = x * z * (q[2] + z * (q[4] + z * q[6]));
 300  341                          r = x + t;
 301  342                          *err = t - (r - x);
 302  343                          return (r);
 303  344                  }
      345 +
 304  346                  z = x * x;
 305  347  
 306  348                  /* use double precision at p4 and on */
 307      -                ee[0] = z *
 308      -                        (P4 + z *
 309      -                        (P5 + z * (P6 + z * (P7 + z * (P8 + z * P9)))));
      349 +                ee[0] = z * (P4 + z * (P5 + z * (P6 + z * (P7 + z * (P8 + z *
      350 +                    P9)))));
 310  351  
 311      -                x_h = (double) ((float) x);
 312      -                z_h = (double) ((float) z);
      352 +                x_h = (double)((float)x);
      353 +                z_h = (double)((float)z);
 313  354                  x_l = x - x_h;
 314  355                  z_l = (x_h * x_h - z_h);
 315  356                  zz[0] = z;
 316  357                  zz[1] = z_h;
 317  358                  zz[2] = z_l + x_l * (x + x_h);
 318  359  
 319  360                  /*
 320  361                   * compute (1+z*(p1+z*(p2+z*(p3+e)))) by call
 321  362                   * mx_poly
 322  363                   */
 323  364  
 324  365                  mx_poly(zz, p, ee, 3);
 325  366  
 326  367                  /* finally x*(1+z*(p1+...)) */
 327  368                  r = x_h * ee[0];
 328  369                  t = x * ee[1] + x_l * ee[0];
 329  370                  s = t + r;
 330  371                  *err = t - (s - r);
 331  372                  return (s);
 332  373          }
      374 +
 333  375          /* for |x| >= 8.0 */
 334      -        if (ix >= 0x40200000) { /* x >=  8 */
      376 +        if (ix >= 0x40200000) {                 /* x >=  8 */
 335  377                  x = fabs(x);
 336      -                if (ix >= 0x42600000) { /* x >=  2**39 */
 337      -                        if (ix >= 0x44c00000) { /* x >=  2**77 */
      378 +
      379 +                if (ix >= 0x42600000) {         /* x >=  2**39 */
      380 +                        if (ix >= 0x44c00000)   /* x >=  2**77 */
 338  381                                  y = -pio2lo;
 339      -                        } else
      382 +                        else
 340  383                                  y = one / x - pio2lo;
      384 +
 341  385                          if (sign == 0) {
 342  386                                  t = pio2hi - y;
 343  387                                  *err = -(y - (pio2hi - t));
 344  388                          } else {
 345  389                                  t = y - pio2hi;
 346  390                                  *err = y - (pio2hi + t);
 347  391                          }
      392 +
 348  393                          return (t);
 349  394                  } else {
 350  395                          /* compute r = 1/x */
 351  396                          r = one / x;
 352  397                          z = r * r;
 353      -                        if (ix < 0x40504000) {  /* 8 <  x <  65 */
 354  398  
      399 +                        if (ix < 0x40504000) {  /* 8 <  x <  65 */
 355  400                                  /* use double precision at p4 and on */
 356      -                                ee[0] = z *
 357      -                                        (P4 + z *
 358      -                                        (P5 + z *
 359      -                                        (P6 + z * (P7 + z * (P8 + z * P9)))));
 360      -                                x_h = (double) ((float) x);
 361      -                                r_h = (double) ((float) r);
 362      -                                z_h = (double) ((float) z);
      401 +                                ee[0] = z * (P4 + z * (P5 + z * (P6 + z * (P7 +
      402 +                                    z * (P8 + z * P9)))));
      403 +                                x_h = (double)((float)x);
      404 +                                r_h = (double)((float)r);
      405 +                                z_h = (double)((float)z);
 363  406                                  r_l = r * ((x_h - x) * r_h - (x_h * r_h - one));
 364  407                                  z_l = (r_h * r_h - z_h);
 365  408                                  zz[0] = z;
 366  409                                  zz[1] = z_h;
 367  410                                  zz[2] = z_l + r_l * (r + r_h);
      411 +
 368  412                                  /*
 369  413                                   * compute (1+z*(p1+z*(p2+z*(p3+e)))) by call
 370  414                                   * mx_poly
 371  415                                   */
 372  416                                  mx_poly(zz, p, ee, 3);
 373      -                        } else { /* x < 65 < 2**39 */
      417 +                        } else {        /* x < 65 < 2**39 */
 374  418                                  /* use double precision at q3 and on */
 375  419                                  ee[0] = z * (Q3 + z * (Q4 + z * Q5));
 376      -                                x_h = (double) ((float) x);
 377      -                                r_h = (double) ((float) r);
 378      -                                z_h = (double) ((float) z);
      420 +                                x_h = (double)((float)x);
      421 +                                r_h = (double)((float)r);
      422 +                                z_h = (double)((float)z);
 379  423                                  r_l = r * ((x_h - x) * r_h - (x_h * r_h - one));
 380  424                                  z_l = (r_h * r_h - z_h);
 381  425                                  zz[0] = z;
 382  426                                  zz[1] = z_h;
 383  427                                  zz[2] = z_l + r_l * (r + r_h);
      428 +
 384  429                                  /*
 385  430                                   * compute (1+z*(q1+z*(q2+e))) by call
 386  431                                   * mx_poly
 387  432                                   */
 388  433                                  mx_poly(zz, q, ee, 2);
 389  434                          }
      435 +
 390  436                          /* pio2 - r*(1+...) */
 391  437                          v = r_h * ee[0];
 392  438                          t = pio2lo - (r * ee[1] + r_l * ee[0]);
      439 +
 393  440                          if (sign == 0) {
 394  441                                  s = pio2hi - v;
 395  442                                  t -= (v - (pio2hi - s));
 396  443                          } else {
 397  444                                  s = v - pio2hi;
 398  445                                  t = -(t - (v - (s + pio2hi)));
 399  446                          }
      447 +
 400  448                          w = s + t;
 401  449                          *err = t - (w - s);
 402  450                          return (w);
 403  451                  }
 404  452          }
      453 +
 405  454          /* now x is between 1/8 and 8 */
 406      -        ((int *) &x)[HIWORD] = ix;
      455 +        ((int *)&x)[HIWORD] = ix;
 407  456          iy = (ix + 0x00008000) & 0x7fff0000;
 408      -        ((int *) &y)[HIWORD] = iy;
 409      -        ((int *) &y)[LOWORD] = 0;
      457 +        ((int *)&y)[HIWORD] = iy;
      458 +        ((int *)&y)[LOWORD] = 0;
 410  459          j = (iy - 0x3fc00000) >> 16;
 411  460  
 412  461          w = (x - y);
 413  462          v = 1 / (one + x * y);
 414  463          s = w * v;
 415  464          z = s * s;
 416  465          /* use double precision at q3 and on */
 417  466          ee[0] = z * (Q3 + z * (Q4 + z * Q5));
 418      -        s_h = (double) ((float) s);
 419      -        z_h = (double) ((float) z);
 420      -        x_h = (double) ((float) x);
 421      -        t = (double) ((float) (one + x * y));
      467 +        s_h = (double)((float)s);
      468 +        z_h = (double)((float)z);
      469 +        x_h = (double)((float)x);
      470 +        t = (double)((float)(one + x * y));
 422  471          r = -((x_h - x) * y - (x_h * y - (t - one)));
 423  472          s_l = -v * (s_h * r - (w - s_h * t));
 424  473          z_l = (s_h * s_h - z_h);
 425  474          zz[0] = z;
 426  475          zz[1] = z_h;
 427  476          zz[2] = z_l + s_l * (s + s_h);
 428  477          /* compute (1+z*(q1+z*(q2+e))) by call mx_poly */
 429  478          mx_poly(zz, q, ee, 2);
 430  479          v = s_h * ee[0];
 431  480          t = TBL_atan_lo[j] + (s * ee[1] + s_l * ee[0]);
 432  481          u = TBL_atan_hi[j];
 433  482          s = u + v;
 434  483          t += (v - (s - u));
 435  484          w = s + t;
 436  485          *err = t - (w - s);
      486 +
 437  487          if (sign != 0) {
 438  488                  w = -w;
 439  489                  *err = -*err;
 440  490          }
      491 +
 441  492          return (w);
 442  493  }
 443  494  
 444      -static const double
 445      -        twom768 = 6.441148769597133308e-232,    /* 2^-768 */
 446      -        two768  = 1.552518092300708935e+231,    /* 2^768 */
      495 +static const double twom768 = 6.441148769597133308e-232,        /* 2^-768 */
      496 +        two768 = 1.552518092300708935e+231,                     /* 2^768 */
 447  497          pi = 3.1415926535897931159979634685,
 448  498          pi_lo = 1.224646799147353177e-16,
 449  499          pio2 = 1.570796326794896558e+00,
 450  500          pio2_lo = 6.123233995736765886e-17,
 451  501          pio4 = 0.78539816339744827899949,
 452  502          pio4_lo = 3.061616997868382943e-17,
 453  503          pi3o4 = 2.356194490192344836998,
 454  504          pi3o4_lo = 9.184850993605148829195e-17;
 455  505  
 456  506  double
 457      -__k_atan2(double y, double x, double *w) {
      507 +__k_atan2(double y, double x, double *w)
      508 +{
 458  509          double t, xh, th, t1, t2, w1, w2;
 459  510          int ix, iy, hx, hy, lx, ly;
 460  511  
 461      -        hy = ((int *) &y)[HIWORD];
 462      -        ly = ((int *) &y)[LOWORD];
      512 +        hy = ((int *)&y)[HIWORD];
      513 +        ly = ((int *)&y)[LOWORD];
 463  514          iy = hy & ~0x80000000;
 464  515  
 465      -        hx = ((int *) &x)[HIWORD];
 466      -        lx = ((int *) &x)[LOWORD];
      516 +        hx = ((int *)&x)[HIWORD];
      517 +        lx = ((int *)&x)[LOWORD];
 467  518          ix = hx & ~0x80000000;
 468  519  
 469  520          *w = 0.0;
      521 +
 470  522          if (ix >= 0x7ff00000 || iy >= 0x7ff00000) {     /* ignore inexact */
 471      -                if (isnan(x) || isnan(y))
      523 +                if (isnan(x) || isnan(y)) {
 472  524                          return (x * y);
 473      -                else if (iy < 0x7ff00000) {
      525 +                } else if (iy < 0x7ff00000) {
 474  526                          if (hx >= 0) {  /* ATAN2(+-finite, +inf) is +-0 */
 475  527                                  *w *= y;
 476  528                                  return (*w);
 477  529                          } else {        /* ATAN2(+-finite, -inf) is +-pi */
 478  530                                  *w = copysign(pi_lo, y);
 479  531                                  return (copysign(pi, y));
 480  532                          }
 481  533                  } else if (ix < 0x7ff00000) {
 482      -                                        /* ATAN2(+-inf, finite) is +-pi/2 */
 483      -                        *w = (hy >= 0)? pio2_lo : -pio2_lo;
 484      -                        return ((hy >= 0)? pio2 : -pio2);
      534 +                        /* ATAN2(+-inf, finite) is +-pi/2 */
      535 +                        *w = (hy >= 0) ? pio2_lo : -pio2_lo;
      536 +                        return ((hy >= 0) ? pio2 : -pio2);
 485  537                  } else if (hx > 0) {    /* ATAN2(+-INF,+INF) = +-pi/4 */
 486      -                        *w = (hy >= 0)? pio4_lo : -pio4_lo;
 487      -                        return ((hy >= 0)? pio4 : -pio4);
      538 +                        *w = (hy >= 0) ? pio4_lo : -pio4_lo;
      539 +                        return ((hy >= 0) ? pio4 : -pio4);
 488  540                  } else {                /* ATAN2(+-INF,-INF) = +-3pi/4 */
 489      -                        *w = (hy >= 0)? pi3o4_lo : -pi3o4_lo;
 490      -                        return ((hy >= 0)? pi3o4 : -pi3o4);
      541 +                        *w = (hy >= 0) ? pi3o4_lo : -pi3o4_lo;
      542 +                        return ((hy >= 0) ? pi3o4 : -pi3o4);
 491  543                  }
 492  544          } else if ((ix | lx) == 0 || (iy | ly) == 0) {
 493  545                  if ((iy | ly) == 0) {
 494      -                        if (hx >= 0) /* ATAN2(+-0, +(0 <= x <= inf)) is +-0 */
      546 +                        if (hx >= 0) { /* ATAN2(+-0, +(0 <= x <= inf)) is +-0 */
 495  547                                  return (y);
 496      -                        else {  /* ATAN2(+-0, -(0 <= x <= inf)) is +-pi */
 497      -                                *w = (hy >= 0)? pi_lo : -pi_lo;
 498      -                                return ((hy >= 0)? pi : -pi);
      548 +                        } else { /* ATAN2(+-0, -(0 <= x <= inf)) is +-pi */
      549 +                                *w = (hy >= 0) ? pi_lo : -pi_lo;
      550 +                                return ((hy >= 0) ? pi : -pi);
 499  551                          }
 500      -                } else { /* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2 */
 501      -                        *w = (hy >= 0)? pio2_lo : -pio2_lo;
 502      -                        return ((hy >= 0)? pio2 : -pio2);
      552 +                } else {  /* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2 */
      553 +                        *w = (hy >= 0) ? pio2_lo : -pio2_lo;
      554 +                        return ((hy >= 0) ? pio2 : -pio2);
 503  555                  }
 504      -        } else if (iy - ix > 0x06400000) { /* |x/y| < 2 ** -100 */
 505      -                *w = (hy >= 0)? pio2_lo : -pio2_lo;
 506      -                return ((hy >= 0)? pio2 : -pio2);
 507      -        } else if (ix - iy > 0x06400000) { /* |y/x| < 2 ** -100 */
      556 +        } else if (iy - ix > 0x06400000) {      /* |x/y| < 2 ** -100 */
      557 +                *w = (hy >= 0) ? pio2_lo : -pio2_lo;
      558 +                return ((hy >= 0) ? pio2 : -pio2);
      559 +        } else if (ix - iy > 0x06400000) {      /* |y/x| < 2 ** -100 */
 508  560                  if (hx < 0) {
 509      -                        *w = (hy >= 0)? pi_lo : -pi_lo;
 510      -                        return ((hy >= 0)? pi : -pi);
      561 +                        *w = (hy >= 0) ? pi_lo : -pi_lo;
      562 +                        return ((hy >= 0) ? pi : -pi);
 511  563                  } else {
 512  564                          t = y / x;
 513  565                          th = t;
 514      -                        ((int *) &th)[LOWORD] &= 0xf8000000;
      566 +                        ((int *)&th)[LOWORD] &= 0xf8000000;
 515  567                          xh = x;
 516      -                        ((int *) &xh)[LOWORD] &= 0xf8000000;
      568 +                        ((int *)&xh)[LOWORD] &= 0xf8000000;
 517  569                          t1 = (x - xh) * t + xh * (t - th);
 518  570                          t2 = y - xh * th;
 519  571                          *w = (t2 - t1) / x;
 520  572                          return (t);
 521  573                  }
 522  574          } else {
 523  575                  if (ix >= 0x5f300000) {
 524  576                          x *= twom768;
 525  577                          y *= twom768;
 526  578                  } else if (ix < 0x23d00000) {
 527  579                          x *= two768;
 528  580                          y *= two768;
 529  581                  }
      582 +
 530  583                  y = fabs(y);
 531  584                  x = fabs(x);
 532  585                  t = y / x;
 533  586                  th = t;
 534      -                ((int *) &th)[LOWORD] &= 0xf8000000;
      587 +                ((int *)&th)[LOWORD] &= 0xf8000000;
 535  588                  xh = x;
 536      -                ((int *) &xh)[LOWORD] &= 0xf8000000;
      589 +                ((int *)&xh)[LOWORD] &= 0xf8000000;
 537  590                  t1 = (x - xh) * t + xh * (t - th);
 538  591                  t2 = y - xh * th;
 539  592                  w1 = mx_atan(t, &w2);
 540  593                  w2 += (t2 - t1) / (x + y * t);
      594 +
 541  595                  if (hx < 0) {
 542  596                          t1 = pi - w1;
 543  597                          t2 = pi - t1;
 544  598                          w2 = (pi_lo - w2) - (w1 - t2);
 545  599                          w1 = t1;
 546  600                  }
 547      -                *w = (hy >= 0)? w2 : -w2;
 548      -                return ((hy >= 0)? w1 : -w1);
      601 +
      602 +                *w = (hy >= 0) ? w2 : -w2;
      603 +                return ((hy >= 0) ? w1 : -w1);
 549  604          }
 550  605  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX