sysdeps/powerpc/powerpc64/fpu/s_nearbyintl.S

   1 /* nearbyint long double.
   2    IBM extended format long double version.
   3    Copyright (C) 2004, 2006 Free Software Foundation, Inc.
   4    This file is part of the GNU C Library.
   5
   6    The GNU C Library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Lesser General Public
   8    License as published by the Free Software Foundation; either
   9    version 2.1 of the License, or (at your option) any later version.
  10
  11    The GNU C Library is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14    Lesser General Public License for more details.
  15
  16    You should have received a copy of the GNU Lesser General Public
  17    License along with the GNU C Library; if not, see
  18    <http://www.gnu.org/licenses/>.  */
  19
  20 #include <sysdep.h>
  21 #include <math_ldbl_opt.h>
  22
  23         .section        ".toc","aw"
  24 .LC0:   /* 2**52 */
  25         .tc FD_43300000_0[TC],0x4330000000000000
  26         .section        ".text"
  27
  28 /* long double [fp1,fp2] nearbyintl (long double x [fp1,fp2])
  29    IEEE 1003.1 nearbyintl function.  nearbyintl is simular to the rintl
  30    but does raise the "inexact" exception.  This implementation is
  31    based on rintl but explicitly maskes the inexact exception on entry
  32    and clears any pending inexact before restoring the exception mask
  33    on exit.
  34
  35    PowerPC64 long double uses the IBM extended format which is
  36    represented two 64-floating point double values. The values are
  37    non-overlapping giving an effective precision of 106 bits. The first
  38    double contains the high order bits of mantisa and is always rounded
  39    to represent a normal rounding of long double to double. Since the
  40    long double value is sum of the high and low values, the low double
  41    normally has the opposite sign to compensate for the this rounding.
  42
  43    For long double there are two cases:
  44    1) |x| < 2**52, all the integer bits are in the high double.
  45       floor the high double and set the low double to -0.0.
  46    2) |x| >= 2**52, Rounding involves both doubles.
  47       See the comment before lable .L2 for details.
  48    */
  49 ENTRY (__nearbyintl)
  50         mffs    fp11            /* Save current FPSCR.  */
  51         lfd     fp13,.LC0@toc(2)
  52         fabs    fp0,fp1
  53         mtfsb0  28              /* Disable "inexact" exceptions.  */
  54         fsub    fp12,fp13,fp13  /* generate 0.0  */
  55         fabs    fp9,fp2
  56         fcmpu   cr7,fp0,fp13    /* if (fabs(x) > TWO52)  */
  57         fcmpu   cr6,fp1,fp12    /* if (x > 0.0)  */
  58         bnl-    cr7,.L2
  59         fmr     fp2,fp12
  60         bng-    cr6,.L4
  61         fadd    fp1,fp1,fp13    /* x+= TWO52;  */
  62         fsub    fp1,fp1,fp13    /* x-= TWO52;  */
  63         b       .L9
  64 .L4:
  65         bnl-    cr6,.L9         /* if (x < 0.0)  */
  66         fsub    fp1,fp13,fp1    /* x = TWO52 - x;  */
  67         fsub    fp0,fp1,fp13    /* x = - (x - TWO52);  */
  68         fneg    fp1,fp0
  69 .L9:
  70         mtfsb0  6               /* Clear any pending "inexact" exceptions.  */
  71         mtfsf   0x01,fp11       /* restore exception mask.  */
  72         blr
  73
  74 /* The high double is > TWO52 so we need to round the low double and
  75    perhaps the high double.  This gets a bit tricky so we use the
  76    following algorithm:
  77
  78    tau = floor(x_high/TWO52);
  79    x0 = x_high - tau;
  80    x1 = x_low + tau;
  81    r1 = nearbyint(x1);
  82    y_high = x0 + r1;
  83    y_low = r1 - tau;
  84    return y;  */
  85 .L2:
  86         fcmpu   cr7,fp9,fp13    /* if (|x_low| > TWO52)  */
  87         fcmpu   cr0,fp9,fp12    /* || (|x_low| == 0.0)  */
  88         bge-    cr7,.L9         /*   return x;  */
  89         beq-  cr0,.L9
  90         fdiv    fp8,fp1,fp13    /* x_high/TWO52  */
  91         fctidz  fp0,fp8
  92         fcfid   fp8,fp0         /* tau = floor(x_high/TWO52);  */
  93         fsub    fp3,fp1,fp8     /* x0 = x_high - tau;  */
  94         fadd    fp4,fp2,fp8     /* x1 = x_low + tau;  */
  95
  96         fcmpu   cr6,fp4,fp12    /* if (x1 > 0.0)  */
  97         bng-    cr6,.L8
  98         fadd    fp5,fp4,fp13    /* r1 = x1 + TWO52;  */
  99         fsub    fp5,fp5,fp13    /* r1 = r1 - TWO52;  */
 100         b       .L6
 101 .L8:
 102         fmr     fp5,fp4
 103         bge-    cr6,.L6         /* if (x1 < 0.0)  */
 104         fsub    fp5,fp13,fp4    /* r1 = TWO52 - x1;  */
 105         fsub    fp0,fp5,fp13    /* r1 = - (r1 - TWO52);  */
 106         fneg    fp5,fp0
 107 .L6:
 108         fadd    fp1,fp3,fp5     /* y_high = x0 + r1;  */
 109         fsub    fp2,fp5,fp8     /* y_low = r1 - tau;  */
 110         b       .L9
 111 END (__nearbyintl)
 112
 113 long_double_symbol (libm, __nearbyintl, nearbyintl)