-mpf_t pi, half_pi, two_pi, e;
-
-/* The gfc_(integer|real)_kinds[] structures have everything the front
- end needs to know about integers and real numbers on the target.
- Other entries of the structure are calculated from these values.
- The first entry is the default kind, the second entry of the real
- structure is the default double kind. */
-
-#define MPZ_NULL {{0,0,0}}
-#define MPF_NULL {{0,0,0,0}}
-
-#define DEF_GFC_INTEGER_KIND(KIND,RADIX,DIGITS,BIT_SIZE) \
- {KIND, RADIX, DIGITS, BIT_SIZE, 0, MPZ_NULL, MPZ_NULL, MPZ_NULL}
-
-#define DEF_GFC_LOGICAL_KIND(KIND,BIT_SIZE) \
- {KIND, BIT_SIZE}
-
-#define DEF_GFC_REAL_KIND(KIND,RADIX,DIGITS,MIN_EXP, MAX_EXP) \
- {KIND, RADIX, DIGITS, MIN_EXP, MAX_EXP, \
- 0, 0, MPF_NULL, MPF_NULL, MPF_NULL}
-
-gfc_integer_info gfc_integer_kinds[] = {
- DEF_GFC_INTEGER_KIND (4, 2, 31, 32),
- DEF_GFC_INTEGER_KIND (8, 2, 63, 64),
- DEF_GFC_INTEGER_KIND (2, 2, 15, 16),
- DEF_GFC_INTEGER_KIND (1, 2, 7, 8),
- DEF_GFC_INTEGER_KIND (0, 0, 0, 0)
-};
-
-gfc_logical_info gfc_logical_kinds[] = {
- DEF_GFC_LOGICAL_KIND (4, 32),
- DEF_GFC_LOGICAL_KIND (8, 64),
- DEF_GFC_LOGICAL_KIND (2, 16),
- DEF_GFC_LOGICAL_KIND (1, 8),
- DEF_GFC_LOGICAL_KIND (0, 0)
-};
-
-gfc_real_info gfc_real_kinds[] = {
- DEF_GFC_REAL_KIND (4, 2, 24, -125, 128),
- DEF_GFC_REAL_KIND (8, 2, 53, -1021, 1024),
- DEF_GFC_REAL_KIND (0, 0, 0, 0, 0)
-};
-
-
-/* The integer kind to use for array indices. This will be set to the
- proper value based on target information from the backend. */
-
-int gfc_index_integer_kind;
-
-
-/* Compute the natural log of arg.
-
- We first get the argument into the range 0.5 to 1.5 by successive
- multiplications or divisions by e. Then we use the series:
-
- ln(x) = (x-1) - (x-1)^2/2 + (x-1)^3/3 - (x-1)^4/4 + ...
-
- Because we are expanding in powers of (x-1), and 0.5 < x < 1.5, we
- have -0.5 < (x-1) < 0.5. Ignoring the harmonic term, this means
- that each term is at most 1/(2^i), meaning one bit is gained per
- iteration.
-
- Not very efficient, but it doesn't have to be. */
-
-void
-natural_logarithm (mpf_t * arg, mpf_t * result)
-{
- mpf_t x, xp, t, log;
- int i, p;
-
- mpf_init_set (x, *arg);
- mpf_init (t);
-
- p = 0;
-
- mpf_set_str (t, "0.5", 10);
- while (mpf_cmp (x, t) < 0)
- {
- mpf_mul (x, x, e);
- p--;
- }
-
- mpf_set_str (t, "1.5", 10);
- while (mpf_cmp (x, t) > 0)
- {
- mpf_div (x, x, e);
- p++;
- }
-
- mpf_sub_ui (x, x, 1);
- mpf_init_set_ui (log, 0);
- mpf_init_set_ui (xp, 1);
-
- for (i = 1; i < GFC_REAL_BITS; i++)
- {
- mpf_mul (xp, xp, x);
- mpf_div_ui (t, xp, i);
-
- if (i % 2 == 0)
- mpf_sub (log, log, t);
- else
- mpf_add (log, log, t);
- }
-
- /* Add in the log (e^p) = p */
-
- if (p < 0)
- mpf_sub_ui (log, log, -p);
- else
- mpf_add_ui (log, log, p);
-
- mpf_clear (x);
- mpf_clear (xp);
- mpf_clear (t);
-
- mpf_set (*result, log);
- mpf_clear (log);
-}
-
-
-/* Calculate the common logarithm of arg. We use the natural
- logaritm of arg and of 10:
-
- log10(arg) = log(arg)/log(10) */
-
-void
-common_logarithm (mpf_t * arg, mpf_t * result)
-{
- mpf_t i10, log10;
-
- natural_logarithm (arg, result);
-
- mpf_init_set_ui (i10, 10);
- mpf_init (log10);
- natural_logarithm (&i10, &log10);
-
- mpf_div (*result, *result, log10);
- mpf_clear (i10);
- mpf_clear (log10);
-}
-
-/* Calculate exp(arg).
-
- We use a reduction of the form
-
- x = Nln2 + r
-
- Then we obtain exp(r) from the Maclaurin series.
- exp(x) is then recovered from the identity
-
- exp(x) = 2^N*exp(r). */
-
-void
-exponential (mpf_t * arg, mpf_t * result)
-{
- mpf_t two, ln2, power, q, r, num, denom, term, x, xp;
- int i;
- long n;
- unsigned long p, mp;
-
-
- mpf_init_set (x, *arg);
-
- if (mpf_cmp_ui (x, 0) == 0)
- {
- mpf_set_ui (*result, 1);
- }
- else if (mpf_cmp_ui (x, 1) == 0)
- {
- mpf_set (*result, e);
- }
- else
- {
- mpf_init_set_ui (two, 2);
- mpf_init (ln2);
- mpf_init (q);
- mpf_init (r);
- mpf_init (power);
- mpf_init (term);
-
- natural_logarithm (&two, &ln2);
-
- mpf_div (q, x, ln2);
- mpf_floor (power, q);
- mpf_mul (q, power, ln2);
- mpf_sub (r, x, q);
-
- mpf_init_set_ui (xp, 1);
- mpf_init_set_ui (num, 1);
- mpf_init_set_ui (denom, 1);
-
- for (i = 1; i <= GFC_REAL_BITS + 10; i++)
- {
- mpf_mul (num, num, r);
- mpf_mul_ui (denom, denom, i);
- mpf_div (term, num, denom);
- mpf_add (xp, xp, term);
- }
-
- /* Reconstruction step */
- n = (long) mpf_get_d (power);
-
- if (n > 0)
- {
- p = (unsigned int) n;
- mpf_mul_2exp (*result, xp, p);
- }
- else
- {
- mp = (unsigned int) (-n);
- mpf_div_2exp (*result, xp, mp);
- }
-
- mpf_clear (two);
- mpf_clear (ln2);
- mpf_clear (q);
- mpf_clear (r);
- mpf_clear (power);
- mpf_clear (num);
- mpf_clear (denom);
- mpf_clear (term);
- mpf_clear (xp);
- }
-
- mpf_clear (x);
-}
-
-
-/* Calculate sin(arg).
-
- We use a reduction of the form
-
- x= N*2pi + r
-
- Then we obtain sin(r) from the Maclaurin series. */