hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

some more work on the fmpcb module; test code for multiplication

Browse source
This commit is contained in:
Fredrik Johansson 2012-11-07 11:48:22 +01:00
parent ea5f128274
commit 615ea353a5
6 changed files with 388 additions and 61 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
doc/source/fmpz_holonomic.rst
View file

@ -325,7 +325,7 @@ Sequence evaluation
.. function:: void fmpz_holonomic_get_nth_fmpq(fmpq_t res, const fmpz_holonomic_t op, const fmpq * initial, long n0, long n)
.. function:: mp_limb_t fmpz_holonomic_get_nth_nmod(const fmpz_holonomic_t op, mp_srcptr initial, long n0, long n, nmod_t mod)
.. function:: mp_limb_t fmpz_holonomic_get_nth_nmod(const fmpz_holonomic_t op, mp_srcptr initial, ulong n0, ulong n, nmod_t mod)
Computes element `c(n)` in the sequence annihilated by the
difference operator *op*, given the

21
doc/source/history.rst
View file

@ -4,6 +4,27 @@ History and changes
For more details, view the detailed commit log
in the git repository https://github.com/fredrik-johansson/arb
* version 0.3-git
* converted documentation to sphinx
* initial support for complex numbers and polynomials
* root isolation for complex polynomials
* new module fmpz_holonomic for functions/sequences
defined by linear differential/difference equations
with polynomial coefficients
* functions for creating various special sequences and functions
* some closure properties for sequences
* Taylor series expansion for differential equations
* computing the nth entry of a sequence using binary splitting
* computing the nth entry mod p using fast multipoint evaluation
* generic code for evaluating hypergeometric series
* matrix powering
* various other helper functions
* 2012-09-29 - version 0.2
* code for computing the gamma function (Karatsuba, Stirling's series)

80
fmpcb.h
View file

@ -125,6 +125,13 @@ fmpcb_swap(fmpcb_t z, fmpcb_t x)
fmprb_swap(fmpcb_imagref(z), fmpcb_imagref(x));
}
static __inline__ int
fmpcb_equal(const fmpcb_t x, const fmpcb_t y)
{
return fmprb_equal(fmpcb_realref(x), fmpcb_realref(y)) &&
fmprb_equal(fmpcb_imagref(x), fmpcb_imagref(y));
}
static __inline__ int
fmpcb_overlaps(const fmpcb_t x, const fmpcb_t y)
{
@ -199,66 +206,7 @@ fmpcb_mul_onei(fmpcb_t z, const fmpcb_t x)
}
}
static __inline__ void
fmpcb_mul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec)
{
#define a fmpcb_realref(x)
#define b fmpcb_imagref(x)
#define c fmpcb_realref(y)
#define d fmpcb_imagref(y)
/*
TODO: use these optimizations, but fix aliasing issues
if (fmprb_is_zero(b))
{
fmpcb_mul_fmprb(z, y, a, prec);
}
else if (fmprb_is_zero(d))
{
fmpcb_mul_fmprb(z, x, c, prec);
}
else if (fmprb_is_zero(a))
{
fmpcb_mul_fmprb(z, y, b, prec);
fmpcb_mul_i(z, z);
}
else if (fmprb_is_zero(c))
{
fmpcb_mul_fmprb(z, x, d, prec);
fmpcb_mul_i(z, z);
}
else
*/
{
fmprb_t t, u, v;
fmprb_init(t);
fmprb_init(u);
fmprb_init(v);
fmprb_add(t, a, b, prec);
fmprb_add(u, c, d, prec);
fmprb_mul(v, t, u, prec);
fmprb_mul(t, a, c, prec);
fmprb_mul(u, b, d, prec);
fmprb_sub(fmpcb_realref(z), t, u, prec);
fmprb_sub(fmpcb_imagref(z), v, t, prec);
fmprb_sub(fmpcb_imagref(z), fmpcb_imagref(z), u, prec);
fmprb_clear(t);
fmprb_clear(u);
fmprb_clear(v);
}
#undef a
#undef b
#undef c
#undef d
}
void fmpcb_mul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec);
static __inline__ void
fmpcb_inv(fmpcb_t z, const fmpcb_t x, long prec)
@ -300,6 +248,18 @@ _fmpcb_vec_set(fmpcb_struct * res, const fmpcb_struct * vec, long len)
fmpcb_set(res + i, vec + i);
}
static __inline__ void
fmpcb_print(const fmpcb_t x)
{
printf("(");
fmprb_print(fmpcb_realref(x));
printf(", ");
fmprb_print(fmpcb_imagref(x));
printf(")");
}
void fmpcb_printd(const fmpcb_t z, long digits);
void fmpcb_randtest(fmpcb_t z, flint_rand_t state, long prec, long mag_bits);
#endif

91
fmpcb/mul.c Normal file
View file

@ -0,0 +1,91 @@
/*=============================================================================
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 "fmpcb.h"
void
fmpcb_mul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec)
{
#define a fmpcb_realref(x)
#define b fmpcb_imagref(x)
#define c fmpcb_realref(y)
#define d fmpcb_imagref(y)
#define e fmpcb_realref(z)
#define f fmpcb_imagref(z)
if (fmprb_is_zero(b))
{
fmprb_mul(f, d, a, prec);
fmprb_mul(e, c, a, prec);
}
else if (fmprb_is_zero(d))
{
fmprb_mul(f, b, c, prec);
fmprb_mul(e, a, c, prec);
}
else if (fmprb_is_zero(a))
{
fmprb_mul(e, c, b, prec);
fmprb_mul(f, d, b, prec);
fmpcb_mul_onei(z, z);
}
else if (fmprb_is_zero(c))
{
fmprb_mul(e, a, d, prec);
fmprb_mul(f, b, d, prec);
fmpcb_mul_onei(z, z);
}
else
{
fmprb_t t, u, v;
fmprb_init(t);
fmprb_init(u);
fmprb_init(v);
fmprb_add(t, a, b, prec);
fmprb_add(u, c, d, prec);
fmprb_mul(v, t, u, prec);
fmprb_mul(t, a, c, prec);
fmprb_mul(u, b, d, prec);
fmprb_sub(e, t, u, prec);
fmprb_sub(f, v, t, prec);
fmprb_sub(f, f, u, prec);
fmprb_clear(t);
fmprb_clear(u);
fmprb_clear(v);
}
#undef a
#undef b
#undef c
#undef d
#undef e
#undef f
}

34
fmpcb/randtest.c Normal file
View file

@ -0,0 +1,34 @@
/*=============================================================================
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 "fmpcb.h"
void
fmpcb_randtest(fmpcb_t z, flint_rand_t state, long prec, long mag_bits)
{
fmprb_randtest(fmpcb_realref(z), state, prec, mag_bits);
fmprb_randtest(fmpcb_imagref(z), state, prec, mag_bits);
}

221
fmpcb/test/t-mul.c Normal file
View file

@ -0,0 +1,221 @@
/*=============================================================================
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 "fmpcb.h"
int main()
{
long iter;
flint_rand_t state;
printf("mul....");
fflush(stdout);
flint_randinit(state);
/* test aliasing of c and a */
for (iter = 0; iter < 100000; iter++)
{
fmpcb_t a, b, c;
long prec;
fmpcb_init(a);
fmpcb_init(b);
fmpcb_init(c);
fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10);
prec = 2 + n_randint(state, 200);
fmpcb_mul(c, a, b, prec);
fmpcb_mul(a, a, b, prec);
if (!fmpcb_equal(a, c))
{
printf("FAIL: aliasing c, a\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
abort();
}
fmpcb_clear(a);
fmpcb_clear(b);
fmpcb_clear(c);
}
/* test aliasing of c and b */
for (iter = 0; iter < 100000; iter++)
{
fmpcb_t a, b, c;
long prec;
fmpcb_init(a);
fmpcb_init(b);
fmpcb_init(c);
fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10);
prec = 2 + n_randint(state, 200);
fmpcb_mul(c, a, b, prec);
fmpcb_mul(b, a, b, prec);
if (!fmpcb_equal(b, c))
{
printf("FAIL: aliasing b, a\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
abort();
}
fmpcb_clear(a);
fmpcb_clear(b);
fmpcb_clear(c);
}
/* test aliasing a, a */
for (iter = 0; iter < 100000; iter++)
{
fmpcb_t a, b, c, d;
long prec;
fmpcb_init(a);
fmpcb_init(b);
fmpcb_init(c);
fmpcb_init(d);
fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10);
prec = 2 + n_randint(state, 200);
fmpcb_set(b, a);
fmpcb_mul(c, a, a, prec);
fmpcb_mul(d, a, b, prec);
if (!fmpcb_equal(c, d))
{
printf("FAIL: aliasing a, a\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
printf("d = "); fmpcb_print(d); printf("\n\n");
abort();
}
fmpcb_clear(a);
fmpcb_clear(b);
fmpcb_clear(c);
fmpcb_clear(d);
}
/* test aliasing a, a, a */
for (iter = 0; iter < 100000; iter++)
{
fmpcb_t a, b, c;
long prec;
fmpcb_init(a);
fmpcb_init(b);
fmpcb_init(c);
fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10);
prec = 2 + n_randint(state, 200);
fmpcb_set(b, a);
fmpcb_mul(c, a, b, prec);
fmpcb_mul(a, a, a, prec);
if (!fmpcb_equal(a, c))
{
printf("FAIL: aliasing a, a, a\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
abort();
}
fmpcb_clear(a);
fmpcb_clear(b);
fmpcb_clear(c);
}
/* test a*(b+c) = a*b + a*c */
for (iter = 0; iter < 100000; iter++)
{
fmpcb_t a, b, c, d, e, f;
fmpcb_init(a);
fmpcb_init(b);
fmpcb_init(c);
fmpcb_init(d);
fmpcb_init(e);
fmpcb_init(f);
fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10);
fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10);
fmpcb_add(d, b, c, 2 + n_randint(state, 200));
fmpcb_mul(e, a, d, 2 + n_randint(state, 200));
fmpcb_mul(d, a, b, 2 + n_randint(state, 200));
fmpcb_mul(f, a, c, 2 + n_randint(state, 200));
fmpcb_add(f, d, f, 2 + n_randint(state, 200));
if (!fmpcb_overlaps(e, f))
{
printf("FAIL: a*(b+c) = a*b + a*c\n\n");
printf("a = "); fmpcb_print(a); printf("\n\n");
printf("b = "); fmpcb_print(b); printf("\n\n");
printf("c = "); fmpcb_print(c); printf("\n\n");
printf("e = "); fmpcb_print(e); printf("\n\n");
printf("f = "); fmpcb_print(f); printf("\n\n");
abort();
}
fmpcb_clear(a);
fmpcb_clear(b);
fmpcb_clear(c);
fmpcb_clear(d);
fmpcb_clear(e);
fmpcb_clear(f);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}