25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * double __k_atan2(double y, double x, double *e)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Compute atan2 with error terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Important formula:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * atan(x) = x - ----- + ----- - ... (for x <= 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * = --- - --- + --- - ... (for x > 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Arg(x + y i) = sign(y) * atan2(|y|, x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * = sign(y) * atan(|y|/x) (for x > 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * sign(y) * (PI - atan(|y|/|x|)) (for x < 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Thus if x >> y (IEEE double: EXP(x) - EXP(y) >= 60):
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1. (x > 0): atan2(y,x) ~ y/x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2. (x < 0): atan2(y,x) ~ sign(y) (PI - |y/x|))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Otherwise if x << y:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * atan2(y,x) ~ sign(y)*PI/2 - x/y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * __k_atan2 call static functions mx_poly, mx_atan
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (void) mx_poly (double *z, double *a, double *e, int n)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * e = a + z*(a + z*(a + ... z*(a + e)...))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1. e and coefficient ai are represented by two double numbers.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * For e, the first one contain the leading 24 bits rounded, and the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * second one contain the remaining 53 bits (total 77 bits accuracy).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * For ai, the first one contian the leading 53 bits rounded, and the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * second is the remaining 53 bits (total 106 bits accuracy).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2. z is an array of three doubles.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * z[0] : the rounded value of Z (the intended value of z)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * z[1] : the leading 24 bits of Z rounded
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * z[2] : the remaining 53 bits of Z
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Note that z[0] = z[1]+z[2] rounded.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtismx_poly(const double *z, const double *a, double *e, int n) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis p = e[0] + a[n];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis p_h = (double) ((float) p);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis p = a[n - 2] + z[0] * p;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* compute p = ai + z * p */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis p_h = (double) ((float) p);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis s += a[i + 1];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis r = t - (p_h - a[i]);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis p = a[i - 2] + z[0] * p;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis e[0] = (double)((float) p);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis r = t - (e[0] - a[0]);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Table of constants for atan from 0.125 to 8
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 0.125 -- 0x3fc00000 --- (increment at bit 16)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 8.000 -- 0x40200000 (total: 97)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * By K.C. Ng, March 9, 1989
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.243549945467614382e-01, 1.320397616146387620e-01, 1.397088742891636204e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.473614810886516302e-01, 1.549967419239409727e-01, 1.626138285979485676e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.702119252854744080e-01, 1.777902289926760471e-01, 1.853479499956947607e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.928843122579746439e-01, 2.003985538258785115e-01, 2.078899272022629863e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis2.153576996977380476e-01, 2.228011537593945213e-01, 2.302195872768437179e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis2.376123138654712419e-01, 2.449786631268641435e-01, 2.596296294082575118e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis2.741674511196587893e-01, 2.885873618940774099e-01, 3.028848683749714166e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis3.170557532091470287e-01, 3.310960767041321029e-01, 3.450021772071051318e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis3.587706702705721895e-01, 3.723984466767542023e-01, 3.858826693980737521e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis3.992207695752525431e-01, 4.124104415973872673e-01, 4.254496373700422662e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis4.383365598579578304e-01, 4.510696559885234436e-01, 4.636476090008060935e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis4.883339510564055352e-01, 5.123894603107377321e-01, 5.358112379604637043e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis5.585993153435624414e-01, 5.807563535676704136e-01, 6.022873461349641522e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis6.231993299340659043e-01, 6.435011087932843710e-01, 6.632029927060932861e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis6.823165548747480713e-01, 7.008544078844501923e-01, 7.188299996216245269e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis7.362574289814280970e-01, 7.531512809621944138e-01, 7.695264804056582975e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis7.853981633974482790e-01, 8.156919233162234217e-01, 8.441539861131710509e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis8.709034570756529758e-01, 8.960553845713439269e-01, 9.197196053504168578e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis9.420000403794636101e-01, 9.629943306809362058e-01, 9.827937232473290541e-01,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.001483135694234639e+00, 1.019141344266349725e+00, 1.035841253008800145e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.051650212548373764e+00, 1.066630365315743623e+00, 1.080839000541168327e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.094328907321189925e+00, 1.107148717794090409e+00, 1.130953743979160375e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.152571997215667610e+00, 1.172273881128476303e+00, 1.190289949682531656e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.206817370285252489e+00, 1.222025323210989667e+00, 1.236059489478081863e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.249045772398254428e+00, 1.261093382252440387e+00, 1.272297395208717319e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.282740879744270757e+00, 1.292496667789785336e+00, 1.301628834009196156e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.310193935047555547e+00, 1.318242051016837113e+00, 1.325817663668032553e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.339705659598999565e+00, 1.352127380920954636e+00, 1.363300100359693845e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.373400766945015894e+00, 1.382574821490125894e+00, 1.390942827002418447e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.398605512271957618e+00, 1.405647649380269870e+00, 1.412141064608495311e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.418146998399631542e+00, 1.423717971406494032e+00, 1.428899272190732761e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.433730152484709031e+00, 1.438244794498222623e+00, 1.442473099109101931e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.446441332248135092e+00,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-3.125324142453938311e-18, -1.276925400709959526e-17, 2.479758919089733066e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis5.409599147666297957e-18, 9.585415594114323829e-18, 7.784470643106252464e-18,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-3.541164079802125137e-18, 2.372599351477449041e-17, 4.180692268843078977e-18,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis2.034098543938166622e-17, 3.139954287184449286e-18, 7.333160666520898500e-18,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis4.738160130078732886e-19, -5.498822172446843173e-18, 1.231340452914270316e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis1.058231431371112987e-17, 1.069875561873445139e-17, 1.923875492461530410e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis8.261353575163771936e-18, -1.428369957377257085e-17, -1.101082790300136900e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-1.893928924292642146e-17, -7.952610375793798701e-18, -2.293880475557830393e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis3.088733564861919217e-17, 1.961231150484565340e-17, 2.378822732491940868e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis2.246598105617042065e-17, 3.963462895355093301e-17, 2.331553074189288466e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-2.494277030626540909e-17, 3.280735600183735558e-17, 2.269877745296168709e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-1.137323618932958456e-17, -2.546278147285580353e-17, -4.063795683482557497e-18,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-5.455630548591626394e-18, -1.441464378193066908e-17, 2.950430737228402307e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis2.672403885140095079e-17, 1.583478505144428617e-17, -3.076054864429649001e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis6.943223671560007740e-18, -1.987626234335816123e-17, -2.147838844445698302e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis3.473937648299456719e-17, -2.425693465918206812e-17, -3.704991905602721293e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis3.061616997868383018e-17, -1.071456562778743077e-17, -4.841337011934916763e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-2.269823590747287052e-17, 2.923876285774304890e-17, -4.057439412852767923e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis5.460837485846687627e-17, -3.986660595210752445e-18, 1.390331103123099845e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis9.438308023545392000e-17, 1.000401886936679889e-17, 3.194313981784503706e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-9.650564731467513515e-17, -5.956589637160374564e-17, -1.567632251135907253e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-5.490676155022364226e-18, 9.404471373566379412e-17, 7.123833804538446299e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-9.159738508900378819e-17, 8.385188614028674371e-17, 7.683333629842068806e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis4.172467638861439118e-17, -2.979162864892849274e-17, 7.879752739459421280e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-2.196203799612310905e-18, 3.242139621534960503e-17, 2.245875015034507026e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-9.283188754266129476e-18, -6.830804768926660334e-17, -1.236918499824626670e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis8.745413734780278834e-17, -6.319394031144676258e-17, -8.824429373951136321e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-2.599011860304134377e-17, 2.147674250751150961e-17, 1.093246171526936217e-16,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-3.307710355769516504e-17, -3.561490438648230100e-17, -9.843712133488842595e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-2.324061182591627982e-17, -8.922630138234492386e-17, -9.573807110557223276e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-8.263883782511013632e-17, 8.721870922223967507e-17, -6.457134743238754385e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis-4.396204466767636187e-17, -2.493019910264565554e-17, -1.105119435430315713e-16,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis9.211323971545051565e-17,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * mx_atan(x,err)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Table look-up algorithm
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * By K.C. Ng, March 9, 1989
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * error (relative)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * |(atan(x)-poly1(x))/x|<= 2^-83.41
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and use poly2(x) to approximate atan(x) for x in [0,1/65] with
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * |atan(x)-poly2(x)|<= 2^-86.8
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Here poly1 and poly2 are odd polynomial with the following form:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * x + x^3*(a1+x^2*(a2+...))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (0). Purge off Inf and NaN and 0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (1). Reduce x to positive by atan(x) = -atan(-x).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (2). For x <= 1/8, use
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (2.1) if x < 2^(-prec/2), atan(x) = x with inexact flag raised
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (2.2) Otherwise
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * atan(x) = poly1(x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (3). For x >= 8 then (prec = 78)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (3.1) if x >= 2^prec, atan(x) = atan(inf) - pio2lo
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (3.2) if x >= 2^(prec/3), atan(x) = atan(inf) - 1/x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (3.3) if x > 65, atan(x) = atan(inf) - poly2(1/x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (3.4) Otherwise, atan(x) = atan(inf) - poly1(1/x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (4). Now x is in (0.125, 8)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Find y that match x to 4.5 bit after binary (easy).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If iy is the high word of y, then
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * single : j = (iy - 0x3e000000) >> 19
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * double : j = (iy - 0x3fc00000) >> 16
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * quad : j = (iy - 0x3ffc0000) >> 12
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Let s = (x-y)/(1+x*y). Then
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * atan(x) = atan(y) + poly1(s)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * = _TBL_atan_hi[j] + (_TBL_atan_lo[j] + poly2(s) )
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double p[] = {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -3.33333333333333314830e-01, /* p1 = BFD55555 55555555 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -1.85030852238476921863e-17, /* p1_l = BC755525 9783A49C */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 2.00000000000000011102e-01, /* p2 = 3FC99999 9999999A */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -1.27263196576150347368e-17, /* p2_l = BC6D584B 0D874007 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -1.42857142857141405923e-01, /* p3 = BFC24924 9249245E */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -1.34258204847170493327e-17, /* p3_l = BC6EF534 A112500D */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 1.11111111110486909803e-01, /* p4 = 3FBC71C7 1C71176A */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -9.09090907557387889470e-02, /* p5 = BFB745D1 73B47A7D */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 7.69230541541713053189e-02, /* p6 = 3FB3B13A B1E68DE6 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -6.66645815401964159097e-02, /* p7 = BFB110EE 1584446A */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 5.87081768778560317279e-02, /* p8 = 3FAE0EFF 87657733 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -4.90818147456113240690e-02, /* p9 = BFA92140 6A524B5C */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double q[] = {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -3.33333333333333314830e-01, /* q1 = BFD55555 55555555 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -1.85022941571278638733e-17, /* q1_l = BC7554E9 D20EFA66 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 1.99999999999999927836e-01, /* q2 = 3FC99999 99999997 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -1.28782564407438833398e-17, /* q2_l = BC6DB1FB 17217417 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -1.42857142855492280642e-01, /* q3 = BFC24924 92483C46 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 1.11111097130183356096e-01, /* q4 = 3FBC71C6 E06595CC */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis -9.08553303569109294013e-02, /* q5 = BFB7424B 808CDA76 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis double y, z, r, s, t, w, s_h, s_l, x_h, x_l, zz[3], ee[2], z_h,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* for |x| < 1/8 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ix < 0x3f300000) { /* when |x| < 2**-12 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ix < 0x3d800000) { /* if |x| < 2**-39 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *err = (double) ((int) x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = x * z * (q[2] + z * (q[4] + z * q[6]));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *err = t - (r - x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* use double precision at p4 and on */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis (P5 + z * (P6 + z * (P7 + z * (P8 + z * P9)))));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x_h = (double) ((float) x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z_h = (double) ((float) z);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * compute (1+z*(p1+z*(p2+z*(p3+e)))) by call
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* finally x*(1+z*(p1+...)) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *err = t - (s - r);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* for |x| >= 8.0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* compute r = 1/x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* use double precision at p4 and on */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x_h = (double) ((float) x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis r_h = (double) ((float) r);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z_h = (double) ((float) z);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis r_l = r * ((x_h - x) * r_h - (x_h * r_h - one));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * compute (1+z*(p1+z*(p2+z*(p3+e)))) by call
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else { /* x < 65 < 2**39 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* use double precision at q3 and on */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x_h = (double) ((float) x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis r_h = (double) ((float) r);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z_h = (double) ((float) z);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis r_l = r * ((x_h - x) * r_h - (x_h * r_h - one));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * compute (1+z*(q1+z*(q2+e))) by call
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* pio2 - r*(1+...) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t -= (v - (pio2hi - s));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = -(t - (v - (s + pio2hi)));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *err = t - (w - s);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* now x is between 1/8 and 8 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ((int *) &y)[LOWORD] = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* use double precision at q3 and on */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis s_h = (double) ((float) s);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z_h = (double) ((float) z);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x_h = (double) ((float) x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = (double) ((float) (one + x * y));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis r = -((x_h - x) * y - (x_h * y - (t - one)));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* compute (1+z*(q1+z*(q2+e))) by call mx_poly */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = TBL_atan_lo[j] + (s * ee[1] + s_l * ee[0]);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t += (v - (s - u));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *err = t - (w - s);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis twom768 = 6.441148769597133308e-232, /* 2^-768 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis two768 = 1.552518092300708935e+231, /* 2^768 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ix >= 0x7ff00000 || iy >= 0x7ff00000) { /* ignore inexact */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (x * y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (hx >= 0) { /* ATAN2(+-finite, +inf) is +-0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else { /* ATAN2(+-finite, -inf) is +-pi */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* ATAN2(+-inf, finite) is +-pi/2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (hx > 0) { /* ATAN2(+-INF,+INF) = +-pi/4 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else { /* ATAN2(+-INF,-INF) = +-3pi/4 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if ((ix | lx) == 0 || (iy | ly) == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (hx >= 0) /* ATAN2(+-0, +(0 <= x <= inf)) is +-0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else { /* ATAN2(+-0, -(0 <= x <= inf)) is +-pi */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else { /* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (iy - ix > 0x06400000) { /* |x/y| < 2 ** -100 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (ix - iy > 0x06400000) { /* |y/x| < 2 ** -100 */