← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/gmp/mpfr/asin.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2006-05-17 02:46:26 UTC
  • Revision ID: james.westby@ubuntu.com-20060517024626-lljr08ftv9g9vefl
Tags: upstream-0.9h-20060510
Import upstream version 0.9h-20060510

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/mpfr/asin.c
 
1
/* mpfr_asin -- arc-sinus of a floating-point number
 
2
 
 
3
Copyright 2001 Free Software Foundation.
 
4
 
 
5
This file is part of the MPFR Library, and was contributed by Mathieu Dutour.
 
6
 
 
7
The MPFR 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 MPFR 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 MPFR 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>
 
23
#include <stdlib.h>
 
24
#include "gmp.h"
 
25
#include "gmp-impl.h"
 
26
#include "mpfr.h"
 
27
#include "mpfr-impl.h"
 
28
 
 
29
int
 
30
mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode)
 
31
{
 
32
  mpfr_t xp;
 
33
  mpfr_t arcs;
 
34
 
 
35
  int signe, suplement;
 
36
 
 
37
  mpfr_t tmp;
 
38
  int Prec;
 
39
  int prec_asin;
 
40
  int good = 0;
 
41
  int realprec;
 
42
  int estimated_delta;
 
43
  int compared; 
 
44
 
 
45
  /* Trivial cases */
 
46
  if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
 
47
    {
 
48
      MPFR_SET_NAN(asin);
 
49
      MPFR_RET_NAN;
 
50
    }
 
51
 
 
52
  /* Set x_p=|x| */
 
53
  signe = MPFR_SIGN(x);
 
54
  mpfr_init2 (xp, MPFR_PREC(x));
 
55
  mpfr_set (xp, x, rnd_mode);
 
56
  if (signe == -1)
 
57
    MPFR_CHANGE_SIGN(xp);
 
58
 
 
59
  compared = mpfr_cmp_ui (xp, 1);
 
60
 
 
61
  if (compared > 0) /* asin(x) = NaN for |x| > 1 */
 
62
    {
 
63
      MPFR_SET_NAN(asin);
 
64
      mpfr_clear (xp);
 
65
      MPFR_RET_NAN;
 
66
    }
 
67
 
 
68
  if (compared == 0) /* x = 1 or x = -1 */
 
69
    {
 
70
      if (signe > 0) /* asin(+1) = Pi/2 */
 
71
        mpfr_const_pi (asin, rnd_mode);
 
72
      else /* asin(-1) = -Pi/2 */
 
73
        {
 
74
          if (rnd_mode == GMP_RNDU)
 
75
            rnd_mode = GMP_RNDD;
 
76
          else if (rnd_mode == GMP_RNDD)
 
77
            rnd_mode = GMP_RNDU;
 
78
          mpfr_const_pi (asin, rnd_mode);
 
79
          mpfr_neg (asin, asin, rnd_mode);
 
80
        }
 
81
      MPFR_EXP(asin)--;
 
82
      mpfr_clear (xp);
 
83
      return 1; /* inexact */
 
84
    }
 
85
 
 
86
  if (MPFR_IS_ZERO(x)) /* x = 0 */
 
87
    {
 
88
      mpfr_set_ui (asin, 0, GMP_RNDN);
 
89
      mpfr_clear(xp);
 
90
      return 0; /* exact result */
 
91
    }
 
92
 
 
93
  prec_asin = MPFR_PREC(asin);
 
94
  mpfr_ui_sub (xp, 1, xp, GMP_RNDD);
 
95
  
 
96
  suplement = 2 - MPFR_EXP(xp);
 
97
#ifdef DEBUG
 
98
  printf("suplement=%d\n", suplement);
 
99
#endif
 
100
  realprec = prec_asin + 10;
 
101
 
 
102
  while (!good)
 
103
    {
 
104
      estimated_delta = 1 + suplement;
 
105
      Prec = realprec+estimated_delta;
 
106
 
 
107
      /* Initialisation    */
 
108
      mpfr_init2 (tmp, Prec);
 
109
      mpfr_init2 (arcs, Prec);
 
110
 
 
111
#ifdef DEBUG
 
112
      printf("Prec=%d\n", Prec);
 
113
      printf("              x=");
 
114
      mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);
 
115
      printf ("\n");
 
116
#endif
 
117
      mpfr_mul (tmp, x, x, GMP_RNDN);
 
118
#ifdef DEBUG
 
119
      printf("            x^2=");
 
120
      mpfr_out_str (stdout, 2, 0, tmp, GMP_RNDN);
 
121
      printf ("\n");
 
122
#endif
 
123
      mpfr_ui_sub (tmp, 1, tmp, GMP_RNDN);
 
124
#ifdef DEBUG
 
125
      printf("          1-x^2=");
 
126
      mpfr_out_str (stdout, 2, 0, tmp, GMP_RNDN);
 
127
      printf ("\n");
 
128
      printf("10:          1-x^2=");
 
129
      mpfr_out_str (stdout, 10, 0, tmp, GMP_RNDN);
 
130
      printf ("\n");
 
131
#endif
 
132
      mpfr_sqrt (tmp, tmp, GMP_RNDN);
 
133
#ifdef DEBUG
 
134
      printf("  sqrt(1-x^2)=");
 
135
      mpfr_out_str (stdout, 2, 0, tmp, GMP_RNDN);
 
136
      printf ("\n");
 
137
      printf("10:  sqrt(1-x^2)=");
 
138
      mpfr_out_str (stdout, 10, 0, tmp, GMP_RNDN);
 
139
      printf ("\n");
 
140
#endif
 
141
      mpfr_div (tmp, x, tmp, GMP_RNDN);
 
142
#ifdef DEBUG
 
143
      printf("x/sqrt(1-x^2)=");
 
144
      mpfr_out_str (stdout, 2, 0, tmp, GMP_RNDN);
 
145
      printf ("\n");
 
146
#endif
 
147
      mpfr_atan (arcs, tmp, GMP_RNDN);
 
148
#ifdef DEBUG
 
149
      printf("atan(x/..x^2)=");
 
150
      mpfr_out_str (stdout, 2, 0, arcs, GMP_RNDN);
 
151
      printf ("\n");
 
152
#endif
 
153
      if (mpfr_can_round (arcs, realprec, GMP_RNDN, rnd_mode, MPFR_PREC(asin)))
 
154
        {
 
155
          mpfr_set (asin, arcs, rnd_mode);
 
156
#ifdef DEBUG
 
157
          printf("asin         =");
 
158
          mpfr_out_str (stdout, 2, prec_asin, asin, GMP_RNDN);
 
159
          printf ("\n");
 
160
#endif
 
161
          good = 1;
 
162
        }
 
163
      else
 
164
        {
 
165
          realprec += _mpfr_ceil_log2 ((double) realprec);
 
166
#ifdef DEBUG
 
167
          printf("RETRY\n");
 
168
#endif
 
169
        }
 
170
      mpfr_clear (tmp);
 
171
      mpfr_clear (arcs);
 
172
  }
 
173
 
 
174
  mpfr_clear (xp);
 
175
 
 
176
  return 1; /* inexact result */
 
177
}
 
178
 
 
179
 
 
180
 
 
181
 
 
182