.file "erfc.s"// Copyright (c) 2001 - 2005, Intel Corporation// All rights reserved.//// Contributed 2001 by the Intel Numerics Group, Intel Corporation//// Redistribution and use in source and binary forms, with or without// modification, are permitted provided that the following conditions are// met://// * Redistributions of source code must retain the above copyright// notice, this list of conditions and the following disclaimer.//// * Redistributions in binary form must reproduce the above copyright// notice, this list of conditions and the following disclaimer in the// documentation and/or other materials provided with the distribution.//// * The name of Intel Corporation may not be used to endorse or promote// products derived from this software without specific prior written// permission.// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // // Intel Corporation is the author of this code, and requests that all// problem reports or change requests be submitted to it directly at // http://www.intel.com/software/products/opensource/libraries/num.htm.//// History//==============================================================// 11/12/01  Initial version// 05/20/02  Cleaned up namespace and sf0 syntax// 02/06/03  Reordered header: .section, .global, .proc, .align// 03/31/05  Reformatted delimiters between data tables//// API//==============================================================// double erfc(double)//// Overview of operation//==============================================================// 1. 0 <= x <= 28.0//    //    erfc(x)  = P14(z) * exp( -x^2 ), z = x - x(i).////    Comment:////    Let x(i) = -1.0 + 2^(i/4),i=0,...19. So we have 20 unequal//    argument intervals [x(i),x(i+1)] with length ratio q = 2^(1/4).//    Values x(i) we have in the table erfc_xb_table.// //    Let x(i)<= x < x(i+1).//    We can find i as exponent of number (x + 1)^4.// //    Let P14(z) - polynomial approximation of degree 14 for function//    erfc(z+x(i)) * exp( (z+x(i))^2) and 0 <= z <= x(i+1)-x(i).//    Polynomial coeffitients we have in the table erfc_p_table.////    So we can find result for erfc(x) as above.//    Algorithm description for exp function see below.   //                     // 2. -6 <= x < 0////      erfc(x)  = 2.0 - erfc(-x)//// 3. x > 28.0//    erfc(x)  ~=~ 0.0                      //// 4. x < -6.0            //    erfc(x)  ~=~ 2.0                      // Special values //==============================================================// erfc(+0)    = 1.0// erfc(-0)    = 1.0// erfc(+qnan) = +qnan // erfc(-qnan) = -qnan // erfc(+snan) = +qnan // erfc(-snan) = -qnan // erfc(-inf)  = 2.0 // erfc(+inf)  = +0//==============================================================// Take double exp(double) from libm_64.//// Overview of operation//==============================================================// Take the input x. w is "how many log2/128 in x?"//  w = x * 128/log2//  n = int(w)//  x = n log2/128 + r + delta//  n = 128M + index_1 + 2^4 index_2//  x = M log2 + (log2/128) index_1 + (log2/8) index_2 + r + delta//  exp(x) = 2^M  2^(index_1/128)  2^(index_2/8) exp(r) exp(delta)//       Construct 2^M//       Get 2^(index_1/128) from table_1;//       Get 2^(index_2/8)   from table_2;//       Calculate exp(r) by series//          r = x - n (log2/128)_high//          delta = - n (log2/128)_low//       Calculate exp(delta) as 1 + delta//==============================================================//  Comment for exp for erfc:////  We use quad precision for calculate input argument -x^2 and add//  result low bits to value delta in exp.  // Registers used//==============================================================// Floating Point registers used: // f8, input// f9 -> f15,  f32 -> f93// General registers used: // r32 -> r68 // Predicate registers used:// p6 -> p15// Assembly macros//==============================================================exp_GR_rshf            = r33EXP_AD_TB1             = r34EXP_AD_TB2             = r35EXP_AD_P               = r36exp_GR_N               = r37exp_GR_index_1         = r38exp_GR_index_2_16      = r39exp_GR_biased_M        = r40EXP_AD_T1              = r41EXP_AD_T2              = r42exp_GR_sig_inv_ln2     = r43exp_GR_17ones          = r44exp_TB1_size           = r45exp_TB2_size           = r46exp_GR_rshf_2to56      = r47exp_GR_exp_2tom56      = r48// GR for erfc(x)//==============================================================GR_POS_ARG_ASYMP       = r49GR_NEG_ARG_ASYMP       = r50GR_ARG_ASYMP           = r51GR_ERFC_XB_TB          = r52GR_ERFC_P_TB           = r53GR_IndxPlusBias        = r54GR_BIAS                = r55GR_P_A12               = r56GR_P_A13               = r57GR_AbsArg              = r58GR_ShftXBi             = r59GR_ShftPi              = r60GR_mBIAS               = r61GR_ShftPi_bias         = r62GR_ShftXBi_bias        = r63GR_ShftA12             = r64GR_ShftA13             = r65GR_EpsNorm             = r66GR_0x1                 = r67GR_ShftPi_8            = r68// GR for __libm_support call//==============================================================GR_SAVE_B0             = r61GR_SAVE_PFS            = r62GR_SAVE_GP             = r63GR_SAVE_SP             = r64GR_Parameter_X         = r65GR_Parameter_Y         = r66GR_Parameter_RESULT    = r67GR_Parameter_TAG       = r68// FR for exp(-x^2)//==============================================================FR_X                   = f10FR_Y                   = f1FR_RESULT              = f8EXP_2TOM56             = f6EXP_INV_LN2_2TO63      = f7EXP_W_2TO56_RSH        = f9EXP_RSHF_2TO56         = f10exp_P4                 = f11 exp_P3                 = f12 exp_P2                 = f13 exp_P1                 = f14 exp_ln2_by_128_hi      = f15 exp_ln2_by_128_lo      = f32 EXP_RSHF               = f33EXP_Nfloat             = f34 exp_r                  = f35exp_f                  = f36exp_rsq                = f37exp_rcube              = f38EXP_2M                 = f39exp_S1                 = f40exp_T1                 = f41exp_rP4pP3             = f42exp_P_lo               = f43exp_P_hi               = f44exp_P                  = f45exp_S                  = f46EXP_NORM_f8            = f47   exp_S2                 = f48exp_T2                 = f49// FR for erfc(x)//==============================================================FR_AbsArg              = f50FR_Tmp                 = f51FR_Xb                  = f52FR_A0                  = f53FR_A1                  = f54FR_A2                  = f55FR_A3                  = f56FR_A4                  = f57FR_A5                  = f58FR_A6                  = f59FR_A7                  = f60FR_A8                  = f61FR_A9                  = f62FR_A10                 = f63FR_A11                 = f64FR_A12                 = f65FR_A13                 = f66FR_A14                 = f67FR_P14_0_1             = f68FR_P14_0_2             = f69FR_P14_1_1             = f70FR_P14_1_2             = f71FR_P14_2_1             = f72FR_P14_2_2             = f73FR_P14_3_1             = f74FR_P14_3_2             = f75FR_P14_6_1             = f76FR_P14_7_1             = f77FR_P14_7_2             = f78FR_P14_8_1             = f79FR_P14_8_2             = f80FR_P14_12_1            = f81FR_P14_13_1            = f82FR_P14_13_2            = f83FR_Pol                 = f84FR_Exp                 = f85FR_2                   = f86f8_sq_lo               = f87FR_LocArg              = f88FR_Tmpf                = f89FR_Tmp1                = f90FR_EpsNorm             = f91FR_UnfBound            = f92FR_NormX               = f93// Data tables//==============================================================RODATA.align 16// ************* DO NOT CHANGE ORDER OF THESE TABLES ********************LOCAL_OBJECT_START(exp_table_1)data8 0x403a8b12fc6e4892 , 0          // underflow boundarydata8 0xb17217f7d1cf79ab , 0x00003ff7 // ln2/128 hidata8 0xc9e3b39803f2f6af , 0x00003fb7 // ln2/128 lo//// Table 1 is 2^(index_1/128) where// index_1 goes from 0 to 15//data8 0x8000000000000000 , 0x00003FFFdata8 0x80B1ED4FD999AB6C , 0x00003FFFdata8 0x8164D1F3BC030773 , 0x00003FFFdata8 0x8218AF4373FC25EC , 0x00003FFFdata8 0x82CD8698AC2BA1D7 , 0x00003FFFdata8 0x8383594EEFB6EE37 , 0x00003FFF