arb/fmprb_poly/log_series.c
Fredrik Johansson 2bec8fb919 add log_series
2012-09-08 17:57:41 +02:00

100 lines
2.9 KiB
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) 2012 Fredrik Johansson
******************************************************************************/
#include "fmprb_poly.h"
void
_fmprb_poly_log_series(fmprb_struct * res, fmprb_struct * f, long n, long prec)
{
/* log(c+f) = log(c*(1+f/c)) = log(c) + log(1+f/c) */
if (!fmprb_is_one(f))
{
fmprb_struct * h;
long i;
h = _fmprb_vec_init(n);
fmprb_one(h);
for (i = 1; i < n; i++)
fmprb_div(h + i, f + i, f, prec);
_fmprb_poly_log_series(res, h, n, prec);
fmprb_log(res, f, prec);
_fmprb_vec_clear(h, n);
}
else
{
fmprb_struct * f_diff;
fmprb_struct * f_inv;
f_diff = _fmprb_vec_init(n);
f_inv = _fmprb_vec_init(n);
_fmprb_poly_derivative(f_diff, f, n, prec);
fmprb_zero(f_diff + n - 1);
_fmprb_poly_inv_series(f_inv, f, n, prec);
_fmprb_poly_mullow(res, f_diff, n - 1, f_inv, n - 1, n - 1, prec);
_fmprb_poly_integral(res, res, n, prec);
_fmprb_vec_clear(f_diff, n);
_fmprb_vec_clear(f_inv, n);
}
}
void
fmprb_poly_log_series(fmprb_poly_t res, const fmprb_poly_t f, long n, long prec)
{
if (f->length < 1 || n < 1)
{
printf("Exception: log_series: input must be nonzero\n");
abort();
}
if (f->length < n || res == f)
{
long i;
fmprb_poly_t t;
fmprb_poly_init(t);
fmprb_poly_fit_length(t, n);
for (i = 0; i < FLINT_MIN(n, f->length); i++)
fmprb_set(t->coeffs + i, f->coeffs + i);
for (i = FLINT_MIN(n, f->length); i < n; i++)
fmprb_zero(t->coeffs + i);
_fmprb_poly_set_length(t, n);
fmprb_poly_log_series(res, t, n, prec);
fmprb_poly_clear(t);
return;
}
fmprb_poly_fit_length(res, n);
_fmprb_poly_log_series(res->coeffs, f->coeffs, n, prec);
_fmprb_poly_set_length(res, n);
_fmprb_poly_normalise(res);
}