← Back to branch summary

~vivek-cs-iitr/helenos/gcc

« back to all changes in this revision

Viewing changes to uspace/app/gcc/dep/libmpc/src/acos.c

  • Committer: Vivek Prakash
  • Date: 2012-06-17 20:36:16 UTC
  • Revision ID: vivekprakash21@gmail.com-20120617203616-f4wl5llksht97f34
add libgmp, libmpfr and libmpc original source

Show diffs side-by-side

added added

removed removed

Lines of Context:
uspace/app/gcc/dep/libmpc/src/acos.c
 
1
/* mpc_acos -- arccosine of a complex number.
 
2
 
 
3
Copyright (C) INRIA, 2009, 2010
 
4
 
 
5
This file is part of the MPC Library.
 
6
 
 
7
The MPC Library is free software; you can redistribute it and/or modify
 
8
it under the terms of the GNU Lesser General Public License as published by
 
9
the Free Software Foundation; either version 2.1 of the License, or (at your
 
10
option) any later version.
 
11
 
 
12
The MPC Library is distributed in the hope that it will be useful, but
 
13
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
14
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
15
License for more details.
 
16
 
 
17
You should have received a copy of the GNU Lesser General Public License
 
18
along with the MPC Library; see the file COPYING.LIB.  If not, write to
 
19
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
20
MA 02111-1307, USA. */
 
21
 
 
22
#include <stdio.h>    /* for MPC_ASSERT */
 
23
#include "mpc-impl.h"
 
24
 
 
25
int
 
26
mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
 
27
{
 
28
  int inex_re, inex_im, inex;
 
29
  mpfr_prec_t p_re, p_im, p;
 
30
  mpc_t z1;
 
31
  mpfr_t pi_over_2;
 
32
  mpfr_exp_t e1, e2;
 
33
  mpfr_rnd_t rnd_im;
 
34
  mpc_rnd_t rnd1;
 
35
 
 
36
  inex_re = 0;
 
37
  inex_im = 0;
 
38
 
 
39
  /* special values */
 
40
  if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op)))
 
41
    {
 
42
      if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op)))
 
43
        {
 
44
          mpfr_set_inf (MPC_IM (rop), mpfr_signbit (MPC_IM (op)) ? +1 : -1);
 
45
          mpfr_set_nan (MPC_RE (rop));
 
46
        }
 
47
      else if (mpfr_zero_p (MPC_RE (op)))
 
48
        {
 
49
          inex_re = set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd));
 
50
          mpfr_set_nan (MPC_IM (rop));
 
51
        }
 
52
      else
 
53
        {
 
54
          mpfr_set_nan (MPC_RE (rop));
 
55
          mpfr_set_nan (MPC_IM (rop));
 
56
        }
 
57
 
 
58
      return MPC_INEX (inex_re, 0);
 
59
    }
 
60
 
 
61
  if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op)))
 
62
    {
 
63
      if (mpfr_inf_p (MPC_RE (op)))
 
64
        {
 
65
          if (mpfr_inf_p (MPC_IM (op)))
 
66
            {
 
67
              if (mpfr_sgn (MPC_RE (op)) > 0)
 
68
                {
 
69
                  inex_re =
 
70
                    set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd));
 
71
                  mpfr_div_2ui (MPC_RE (rop), MPC_RE (rop), 1, GMP_RNDN);
 
72
                }
 
73
              else
 
74
                {
 
75
 
 
76
                  /* the real part of the result is 3*pi/4
 
77
                     a = o(pi)  error(a) < 1 ulp(a)
 
78
                     b = o(3*a) error(b) < 2 ulp(b)
 
79
                     c = b/4    exact
 
80
                     thus 1 bit is lost */
 
81
                  mpfr_t x;
 
82
                  mpfr_prec_t prec;
 
83
                  int ok;
 
84
                  mpfr_init (x);
 
85
                  prec = mpfr_get_prec (MPC_RE (rop));
 
86
                  p = prec;
 
87
 
 
88
                  do
 
89
                    {
 
90
                      p += mpc_ceil_log2 (p);
 
91
                      mpfr_set_prec (x, p);
 
92
                      mpfr_const_pi (x, GMP_RNDD);
 
93
                      mpfr_mul_ui (x, x, 3, GMP_RNDD);
 
94
                      ok =
 
95
                        mpfr_can_round (x, p - 1, GMP_RNDD, MPC_RND_RE (rnd),
 
96
                                        prec+(MPC_RND_RE (rnd) == GMP_RNDN));
 
97
 
 
98
                    } while (ok == 0);
 
99
                  inex_re =
 
100
                    mpfr_div_2ui (MPC_RE (rop), x, 2, MPC_RND_RE (rnd));
 
101
                  mpfr_clear (x);
 
102
                }
 
103
            }
 
104
          else
 
105
            {
 
106
              if (mpfr_sgn (MPC_RE (op)) > 0)
 
107
                mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN);
 
108
              else
 
109
                inex_re = mpfr_const_pi (MPC_RE (rop), MPC_RND_RE (rnd));
 
110
            }
 
111
        }
 
112
      else
 
113
        inex_re = set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd));
 
114
 
 
115
      mpfr_set_inf (MPC_IM (rop), mpfr_signbit (MPC_IM (op)) ? +1 : -1);
 
116
 
 
117
      return MPC_INEX (inex_re, 0);
 
118
    }
 
119
 
 
120
  /* pure real argument */
 
121
  if (mpfr_zero_p (MPC_IM (op)))
 
122
    {
 
123
      int s_im;
 
124
      s_im = mpfr_signbit (MPC_IM (op));
 
125
 
 
126
      if (mpfr_cmp_ui (MPC_RE (op), 1) > 0)
 
127
        {
 
128
          if (s_im)
 
129
            inex_im = mpfr_acosh (MPC_IM (rop), MPC_RE (op),
 
130
                                  MPC_RND_IM (rnd));
 
131
          else
 
132
            inex_im = -mpfr_acosh (MPC_IM (rop), MPC_RE (op),
 
133
                                   INV_RND (MPC_RND_IM (rnd)));
 
134
 
 
135
          mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN);
 
136
        }
 
137
      else if (mpfr_cmp_si (MPC_RE (op), -1) < 0)
 
138
        {
 
139
          mpfr_t minus_op_re;
 
140
          minus_op_re[0] = MPC_RE (op)[0];
 
141
          MPFR_CHANGE_SIGN (minus_op_re);
 
142
 
 
143
          if (s_im)
 
144
            inex_im = mpfr_acosh (MPC_IM (rop), minus_op_re,
 
145
                                  MPC_RND_IM (rnd));
 
146
          else
 
147
            inex_im = -mpfr_acosh (MPC_IM (rop), minus_op_re,
 
148
                                   INV_RND (MPC_RND_IM (rnd)));
 
149
          inex_re = mpfr_const_pi (MPC_RE (rop), MPC_RND_RE (rnd));
 
150
        }
 
151
      else
 
152
        {
 
153
          inex_re = mpfr_acos (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd));
 
154
          mpfr_set_ui (MPC_IM (rop), 0, MPC_RND_IM (rnd));
 
155
        }
 
156
 
 
157
      if (!s_im)
 
158
        mpc_conj (rop, rop, MPC_RNDNN);
 
159
 
 
160
      return MPC_INEX (inex_re, inex_im);
 
161
    }
 
162
 
 
163
  /* pure imaginary argument */
 
164
  if (mpfr_zero_p (MPC_RE (op)))
 
165
    {
 
166
      inex_re = set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd));
 
167
      inex_im = -mpfr_asinh (MPC_IM (rop), MPC_IM (op),
 
168
                             INV_RND (MPC_RND_IM (rnd)));
 
169
      mpc_conj (rop,rop, MPC_RNDNN);
 
170
 
 
171
      return MPC_INEX (inex_re, inex_im);
 
172
    }
 
173
 
 
174
  /* regular complex argument: acos(z) = Pi/2 - asin(z) */
 
175
  p_re = mpfr_get_prec (MPC_RE(rop));
 
176
  p_im = mpfr_get_prec (MPC_IM(rop));
 
177
  p = p_re;
 
178
  mpc_init3 (z1, p, p_im); /* we round directly the imaginary part to p_im,
 
179
                              with rounding mode opposite to rnd_im */
 
180
  rnd_im = MPC_RND_IM(rnd);
 
181
  /* the imaginary part of asin(z) has the same sign as Im(z), thus if
 
182
     Im(z) > 0 and rnd_im = RNDZ, we want to round the Im(asin(z)) to -Inf
 
183
     so that -Im(asin(z)) is rounded to zero */
 
184
  if (rnd_im == GMP_RNDZ)
 
185
    rnd_im = mpfr_sgn (MPC_IM(op)) > 0 ? GMP_RNDD : GMP_RNDU;
 
186
  else
 
187
    rnd_im = rnd_im == GMP_RNDU ? GMP_RNDD
 
188
      : rnd_im == GMP_RNDD ? GMP_RNDU
 
189
      : rnd_im; /* both RNDZ and RNDA map to themselves for -asin(z) */
 
190
  rnd1 = RNDC(GMP_RNDN, rnd_im);
 
191
  mpfr_init2 (pi_over_2, p);
 
192
  for (;;)
 
193
    {
 
194
      p += mpc_ceil_log2 (p) + 3;
 
195
 
 
196
      mpfr_set_prec (MPC_RE(z1), p);
 
197
      mpfr_set_prec (pi_over_2, p);
 
198
 
 
199
      mpfr_const_pi (pi_over_2, GMP_RNDN);
 
200
      mpfr_div_2exp (pi_over_2, pi_over_2, 1, GMP_RNDN); /* Pi/2 */
 
201
      e1 = 1; /* Exp(pi_over_2) */
 
202
      inex = mpc_asin (z1, op, rnd1); /* asin(z) */
 
203
      MPC_ASSERT (mpfr_sgn (MPC_IM(z1)) * mpfr_sgn (MPC_IM(op)) > 0);
 
204
      inex_im = MPC_INEX_IM(inex); /* inex_im is in {-1, 0, 1} */
 
205
      e2 = mpfr_get_exp (MPC_RE(z1));
 
206
      mpfr_sub (MPC_RE(z1), pi_over_2, MPC_RE(z1), GMP_RNDN);
 
207
      if (!mpfr_zero_p (MPC_RE(z1)))
 
208
        {
 
209
          /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) +
 
210
             2^(e2-p-1) */
 
211
          e1 = e1 >= e2 ? e1 + 1 : e2 + 1;
 
212
          /* the error on x is bounded by 1/2 ulp(x) + 2^(e1-p-1) */
 
213
          e1 -= mpfr_get_exp (MPC_RE(z1));
 
214
          /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */
 
215
          e1 = e1 <= 0 ? 0 : e1;
 
216
          /* the error on x is bounded by 2^e1 * ulp(x) */
 
217
          mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN); /* exact */
 
218
          inex_im = -inex_im;
 
219
          if (mpfr_can_round (MPC_RE(z1), p - e1, GMP_RNDN, GMP_RNDZ,
 
220
                              p_re + (MPC_RND_RE(rnd) == GMP_RNDN)))
 
221
            break;
 
222
        }
 
223
    }
 
224
  inex = mpc_set (rop, z1, rnd);
 
225
  inex_re = MPC_INEX_RE(inex);
 
226
  mpc_clear (z1);
 
227
  mpfr_clear (pi_over_2);
 
228
 
 
229
  return MPC_INEX(inex_re, inex_im);
 
230
}