Codebase list flint-arb / HEAD acb_poly / mullow_transpose.c
HEAD

Tree @HEAD (Download .tar.gz)

mullow_transpose.c @HEAD — raw · history · blame 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
/*
    Copyright (C) 2012 Fredrik Johansson

    This file is part of Arb.

    Arb is free software: you can redistribute it and/or modify it under
    the terms of the GNU Lesser General Public License (LGPL) as published
    by the Free Software Foundation; either version 2.1 of the License, or
    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
*/

#include "acb_poly.h"

void
_acb_poly_mullow_transpose(acb_ptr res,
    acb_srcptr poly1, slong len1,
    acb_srcptr poly2, slong len2, slong n, slong prec)
{
    arb_ptr a, b, c, d, e, f, w;
    arb_ptr t;
    slong i;

    len1 = FLINT_MIN(len1, n);
    len2 = FLINT_MIN(len2, n);

    w = flint_malloc(sizeof(arb_struct) * (2 * (len1 + len2 + n)));
    a = w;
    b = a + len1;
    c = b + len1;
    d = c + len2;
    e = d + len2;
    f = e + n;

    /* (e+fi) = (a+bi)(c+di) = (ac - bd) + (ad + bc)i */
    t = _arb_vec_init(n);

    for (i = 0; i < len1; i++)
    {
        a[i] = *acb_realref(poly1 + i);
        b[i] = *acb_imagref(poly1 + i);
    }

    for (i = 0; i < len2; i++)
    {
        c[i] = *acb_realref(poly2 + i);
        d[i] = *acb_imagref(poly2 + i);
    }

    for (i = 0; i < n; i++)
    {
        e[i] = *acb_realref(res + i);
        f[i] = *acb_imagref(res + i);
    }

    _arb_poly_mullow(e, a, len1, c, len2, n, prec);
    _arb_poly_mullow(t, b, len1, d, len2, n, prec);
    _arb_vec_sub(e, e, t, n, prec);

    _arb_poly_mullow(f, a, len1, d, len2, n, prec);
    /* squaring */
    if (poly1 == poly2 && len1 == len2)
    {
        _arb_vec_scalar_mul_2exp_si(f, f, n, 1);
    }
    else
    {
        _arb_poly_mullow(t, b, len1, c, len2, n, prec);
        _arb_vec_add(f, f, t, n, prec);
    }

    for (i = 0; i < n; i++)
    {
        *acb_realref(res + i) = e[i];
        *acb_imagref(res + i) = f[i];
    }

    _arb_vec_clear(t, n);
    flint_free(w);
}

void
acb_poly_mullow_transpose(acb_poly_t res, const acb_poly_t poly1,
                                            const acb_poly_t poly2,
                                                slong n, slong prec)
{
    slong len1, len2;

    len1 = poly1->length;
    len2 = poly2->length;

    if (len1 == 0 || len2 == 0 || n == 0)
    {
        acb_poly_zero(res);
        return;
    }

    n = FLINT_MIN((len1 + len2 - 1), n);
    len1 = FLINT_MIN(len1, n);
    len2 = FLINT_MIN(len2, n);

    if (res == poly1 || res == poly2)
    {
        acb_poly_t t;
        acb_poly_init2(t, n);
        _acb_poly_mullow_transpose(t->coeffs, poly1->coeffs, len1,
                                poly2->coeffs, len2, n, prec);
        acb_poly_swap(res, t);
        acb_poly_clear(t);
    }
    else
    {
        acb_poly_fit_length(res, n);
        _acb_poly_mullow_transpose(res->coeffs, poly1->coeffs, len1,
                                poly2->coeffs, len2, n, prec);
    }

    _acb_poly_set_length(res, n);
    _acb_poly_normalise(res);
}