arb/acb/rising_ui_get_mag.c

/*=============================================================================

    This file is part of ARB.

    ARB is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    ARB is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with ARB; if not, write to the Free Software
    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

=============================================================================*/
/******************************************************************************

    Copyright (C) 2014 Fredrik Johansson

******************************************************************************/

#include "acb.h"

void mag_add_ui_2exp_si(mag_t, const mag_t, ulong, long);

static void
acb_rising_get_mag2_right(mag_t bound, const arb_t a, const arb_t b, ulong n)
{
    mag_t t, u;
    ulong k;

    mag_init(t);
    mag_init(u);

    arb_get_mag(t, a);
    arb_get_mag(u, b);

    mag_mul(bound, t, t);
    mag_addmul(bound, u, u);
    mag_set(u, bound);
    mag_mul_2exp_si(t, t, 1);

    for (k = 1; k < n; k++)
    {
        mag_add_ui_2exp_si(u, u, 2 * k - 1, 0);
        mag_add(u, u, t);
        mag_mul(bound, bound, u);
    }

    mag_clear(t);
    mag_clear(u);
}

void
acb_rising_ui_get_mag(mag_t bound, const acb_t s, ulong n)
{
    if (n == 0)
    {
        mag_one(bound);
        return;
    }

    if (n == 1)
    {
        acb_get_mag(bound, s);
        return;
    }

    if (!acb_is_finite(s))
    {
        mag_inf(bound);
        return;
    }

    if (arf_sgn(arb_midref(acb_realref(s))) >= 0)
    {
        acb_rising_get_mag2_right(bound, acb_realref(s), acb_imagref(s), n);
    }
    else
    {
        arb_t a;
        long k;
        mag_t bound2, t, u;

        arb_init(a);
        mag_init(bound2);
        mag_init(t);
        mag_init(u);

        arb_get_mag(u, acb_imagref(s));
        mag_mul(u, u, u);
        mag_one(bound);

        for (k = 0; k < n; k++)
        {
            arb_add_ui(a, acb_realref(s), k, MAG_BITS);

            if (arf_sgn(arb_midref(a)) >= 0)
            {
                acb_rising_get_mag2_right(bound2, a, acb_imagref(s), n - k);
                mag_mul(bound, bound, bound2);
                break;
            }
            else
            {
                arb_get_mag(t, a);
                mag_mul(t, t, t);
                mag_add(t, t, u);
                mag_mul(bound, bound, t);
            }
        }

        arb_clear(a);
        mag_clear(bound2);
        mag_clear(u);
    }

    mag_sqrt(bound, bound);
}