Commit bdee44d3100163b767bc92a2c426c5bda2abfafa - flint-arb

+5

-3

CMakeLists.txt less more

135	135	target_link_libraries(arb m)
136	136	endif()
137	137
	138	include(GNUInstallDirs)
	139
138	140	install(TARGETS arb
139		RUNTIME DESTINATION bin
140		ARCHIVE DESTINATION lib
141		LIBRARY DESTINATION lib
	141	RUNTIME DESTINATION "${CMAKE_INSTALL_FULL_BINDIR}"
	142	ARCHIVE DESTINATION "${CMAKE_INSTALL_FULL_LIBDIR}"
	143	LIBRARY DESTINATION "${CMAKE_INSTALL_FULL_LIBDIR}"
142	144	)
143	145
144	146	foreach (FOLDER ${FOLDERS})

+146

-0

acb/dot_fmpz.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	void
	14	acb_dot_fmpz(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec)
	15	{
	16	arb_ptr t;
	17	slong i, ssize, size, tmp_size;
	18	mp_ptr ztmp;
	19	fmpz v;
	20	ulong av, al;
	21	unsigned int bc;
	22	TMP_INIT;
	23
	24	/* todo: fast fma and fmma (len=2) code */
	25	if (len <= 1)
	26	{
	27	if (initial == NULL)
	28	{
	29	if (len <= 0)
	30	acb_zero(res);
	31	else
	32	{
	33	acb_mul_fmpz(res, x, y, prec);
	34	if (subtract)
	35	acb_neg(res, res);
	36	}
	37	return;
	38	}
	39	else if (len <= 0)
	40	{
	41	acb_set_round(res, initial, prec);
	42	return;
	43	}
	44	}
	45
	46	TMP_START;
	47	t = TMP_ALLOC(sizeof(arb_struct) * len);
	48
	49	tmp_size = 0;
	50	for (i = 0; i < len; i++)
	51	{
	52	v = y[i * ystep];
	53
	54	MAG_EXP(arb_radref(t + i)) = 0;
	55	MAG_MAN(arb_radref(t + i)) = 0;
	56
	57	if (v == 0)
	58	{
	59	ARF_XSIZE(arb_midref(t + i)) = 0;
	60	ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO;
	61	}
	62	else if (!COEFF_IS_MPZ(v))
	63	{
	64	av = FLINT_ABS(v);
	65	count_leading_zeros(bc, av);
	66
	67	ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
	68	ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc;
	69	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, v < 0);
	70	}
	71	else
	72	{
	73	__mpz_struct * z = COEFF_TO_PTR(v);
	74
	75	ssize = z->_mp_size;
	76	size = FLINT_ABS(ssize);
	77
	78	av = z->_mp_d[size - 1];
	79	count_leading_zeros(bc, av);
	80
	81	if (size == 1)
	82	{
	83	ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
	84	ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc;
	85	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, ssize < 0);
	86	}
	87	else if (size == 2)
	88	{
	89	al = z->_mp_d[0];
	90
	91	ARF_EXP(arb_midref(t + i)) = 2 * FLINT_BITS - bc;
	92
	93	if (bc != 0)
	94	{
	95	av = (av << bc) \| (al >> (FLINT_BITS - bc));
	96	al = al << bc;
	97	}
	98
	99	ARF_NOPTR_D(arb_midref(t + i))[0] = al;
	100	ARF_NOPTR_D(arb_midref(t + i))[1] = av;
	101	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(2, ssize < 0);
	102	}
	103	else
	104	{
	105	if (bc != 0)
	106	{
	107	tmp_size += size;
	108	/* use to flag tmp where we need tmp storage */
	109	MAG_MAN(arb_radref(t + i)) = bc;
	110	}
	111
	112	ARF_EXP(arb_midref(t + i)) = size * FLINT_BITS - bc;
	113	ARF_PTR_D(arb_midref(t + i)) = z->_mp_d;
	114	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(size, ssize < 0);
	115	}
	116	}
	117	}
	118
	119	if (tmp_size != 0)
	120	{
	121	ztmp = TMP_ALLOC(sizeof(mp_limb_t) * tmp_size);
	122
	123	for (i = 0; i < len; i++)
	124	{
	125	bc = MAG_MAN(arb_radref(t + i));
	126
	127	if (bc != 0)
	128	{
	129	size = ARF_SIZE(arb_midref(t + i));
	130
	131	mpn_lshift(ztmp, ARF_PTR_D(arb_midref(t + i)), size, bc);
	132	ARF_PTR_D(arb_midref(t + i)) = ztmp;
	133	ztmp += size;
	134	}
	135
	136	MAG_MAN(arb_radref(t + i)) = 0;
	137	}
	138	}
	139
	140	arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec);
	141	arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec);
	142
	143	TMP_END;
	144	}
	145

+77

-0

acb/dot_si.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	void
	14	acb_dot_si(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec)
	15	{
	16	arb_ptr t;
	17	slong i;
	18	slong v;
	19	ulong av;
	20	unsigned int bc;
	21	TMP_INIT;
	22
	23	/* todo: fast fma and fmma (len=2) code */
	24	if (len <= 1)
	25	{
	26	if (initial == NULL)
	27	{
	28	if (len <= 0)
	29	acb_zero(res);
	30	else
	31	{
	32	acb_mul_si(res, x, y[0], prec);
	33	if (subtract)
	34	acb_neg(res, res);
	35	}
	36	return;
	37	}
	38	else if (len <= 0)
	39	{
	40	acb_set_round(res, initial, prec);
	41	return;
	42	}
	43	}
	44
	45	TMP_START;
	46	t = TMP_ALLOC(sizeof(arb_struct) * len);
	47
	48	for (i = 0; i < len; i++)
	49	{
	50	v = y[i * ystep];
	51
	52	if (v == 0)
	53	{
	54	ARF_XSIZE(arb_midref(t + i)) = 0;
	55	ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO;
	56	}
	57	else
	58	{
	59	av = FLINT_ABS(v);
	60	count_leading_zeros(bc, av);
	61
	62	ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
	63	ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc;
	64	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, v < 0);
	65	}
	66
	67	MAG_EXP(arb_radref(t + i)) = 0;
	68	MAG_MAN(arb_radref(t + i)) = 0;
	69	}
	70
	71	arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec);
	72	arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec);
	73
	74	TMP_END;
	75	}
	76

+115

-0

acb/dot_siui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	static void
	14	arf_shallow_set_siui(arf_t res, ulong vhi, ulong vlo)
	15	{
	16	int negative;
	17	unsigned int bc;
	18
	19	negative = ((slong) vhi) < 0;
	20
	21	if (negative)
	22	{
	23	vhi = -vhi - (vlo != 0);
	24	vlo = -vlo;
	25	}
	26
	27	if (vhi == 0)
	28	{
	29	if (vlo == 0)
	30	{
	31	ARF_XSIZE(res) = 0;
	32	ARF_EXP(res) = ARF_EXP_ZERO;
	33	}
	34	else
	35	{
	36	count_leading_zeros(bc, vlo);
	37	ARF_EXP(res) = FLINT_BITS - bc;
	38	ARF_NOPTR_D(res)[0] = vlo << bc;
	39	ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative);
	40	}
	41	}
	42	else if (vlo == 0)
	43	{
	44	count_leading_zeros(bc, vhi);
	45	ARF_EXP(res) = 2 * FLINT_BITS - bc;
	46	ARF_NOPTR_D(res)[0] = vhi << bc;
	47	ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative);
	48	}
	49	else
	50	{
	51	count_leading_zeros(bc, vhi);
	52	ARF_EXP(res) = 2 * FLINT_BITS - bc;
	53	ARF_NOPTR_D(res)[0] = vlo << bc;
	54	if (bc == 0)
	55	ARF_NOPTR_D(res)[1] = vhi;
	56	else
	57	ARF_NOPTR_D(res)[1] = (vhi << bc) \| (vlo >> (FLINT_BITS - bc));
	58	ARF_XSIZE(res) = ARF_MAKE_XSIZE(2, negative);
	59	}
	60	}
	61
	62	void
	63	acb_dot_siui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	64	{
	65	arb_ptr t;
	66	slong i;
	67	ulong vhi, vlo;
	68	TMP_INIT;
	69
	70	/* todo: fast fma and fmma (len=2) code */
	71	if (len <= 1)
	72	{
	73	if (initial == NULL)
	74	{
	75	if (len <= 0)
	76	acb_zero(res);
	77	else
	78	{
	79	arf_t t;
	80	arf_shallow_set_siui(t, y[1], y[0]);
	81	arb_mul_arf(acb_realref(res), acb_realref(x), t, prec);
	82	arb_mul_arf(acb_imagref(res), acb_imagref(x), t, prec);
	83	if (subtract)
	84	acb_neg(res, res);
	85	}
	86	return;
	87	}
	88	else if (len <= 0)
	89	{
	90	acb_set_round(res, initial, prec);
	91	return;
	92	}
	93	}
	94
	95	TMP_START;
	96	t = TMP_ALLOC(sizeof(arb_struct) * len);
	97
	98	for (i = 0; i < len; i++)
	99	{
	100	vlo = y[2 * i * ystep];
	101	vhi = y[2 * i * ystep + 1];
	102
	103	arf_shallow_set_siui(arb_midref(t + i), vhi, vlo);
	104
	105	MAG_EXP(arb_radref(t + i)) = 0;
	106	MAG_MAN(arb_radref(t + i)) = 0;
	107	}
	108
	109	arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec);
	110	arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec);
	111
	112	TMP_END;
	113	}
	114

+75

-0

acb/dot_ui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	void
	14	acb_dot_ui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	15	{
	16	arb_ptr t;
	17	slong i;
	18	ulong v;
	19	unsigned int bc;
	20	TMP_INIT;
	21
	22	/* todo: fast fma and fmma (len=2) code */
	23	if (len <= 1)
	24	{
	25	if (initial == NULL)
	26	{
	27	if (len <= 0)
	28	acb_zero(res);
	29	else
	30	{
	31	acb_mul_ui(res, x, y[0], prec);
	32	if (subtract)
	33	acb_neg(res, res);
	34	}
	35	return;
	36	}
	37	else if (len <= 0)
	38	{
	39	acb_set_round(res, initial, prec);
	40	return;
	41	}
	42	}
	43
	44	TMP_START;
	45	t = TMP_ALLOC(sizeof(arb_struct) * len);
	46
	47	for (i = 0; i < len; i++)
	48	{
	49	v = y[i * ystep];
	50
	51	if (v == 0)
	52	{
	53	ARF_XSIZE(arb_midref(t + i)) = 0;
	54	ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO;
	55	}
	56	else
	57	{
	58	count_leading_zeros(bc, v);
	59
	60	ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
	61	ARF_NOPTR_D(arb_midref(t + i))[0] = v << bc;
	62	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, 0);
	63	}
	64
	65	MAG_EXP(arb_radref(t + i)) = 0;
	66	MAG_MAN(arb_radref(t + i)) = 0;
	67	}
	68
	69	arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec);
	70	arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec);
	71
	72	TMP_END;
	73	}
	74

+106

-0

acb/dot_uiui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	static void
	14	arf_shallow_set_uiui(arf_t res, ulong vhi, ulong vlo)
	15	{
	16	unsigned int bc;
	17
	18	if (vhi == 0)
	19	{
	20	if (vlo == 0)
	21	{
	22	ARF_XSIZE(res) = 0;
	23	ARF_EXP(res) = ARF_EXP_ZERO;
	24	}
	25	else
	26	{
	27	count_leading_zeros(bc, vlo);
	28	ARF_EXP(res) = FLINT_BITS - bc;
	29	ARF_NOPTR_D(res)[0] = vlo << bc;
	30	ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, 0);
	31	}
	32	}
	33	else if (vlo == 0)
	34	{
	35	count_leading_zeros(bc, vhi);
	36	ARF_EXP(res) = 2 * FLINT_BITS - bc;
	37	ARF_NOPTR_D(res)[0] = vhi << bc;
	38	ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, 0);
	39	}
	40	else
	41	{
	42	count_leading_zeros(bc, vhi);
	43	ARF_EXP(res) = 2 * FLINT_BITS - bc;
	44	ARF_NOPTR_D(res)[0] = vlo << bc;
	45	if (bc == 0)
	46	ARF_NOPTR_D(res)[1] = vhi;
	47	else
	48	ARF_NOPTR_D(res)[1] = (vhi << bc) \| (vlo >> (FLINT_BITS - bc));
	49	ARF_XSIZE(res) = ARF_MAKE_XSIZE(2, 0);
	50	}
	51	}
	52
	53	void
	54	acb_dot_uiui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	55	{
	56	arb_ptr t;
	57	slong i;
	58	ulong vhi, vlo;
	59	TMP_INIT;
	60
	61	/* todo: fast fma and fmma (len=2) code */
	62	if (len <= 1)
	63	{
	64	if (initial == NULL)
	65	{
	66	if (len <= 0)
	67	acb_zero(res);
	68	else
	69	{
	70	arf_t t;
	71	arf_shallow_set_uiui(t, y[1], y[0]);
	72	arb_mul_arf(acb_realref(res), acb_realref(x), t, prec);
	73	arb_mul_arf(acb_imagref(res), acb_imagref(x), t, prec);
	74	if (subtract)
	75	acb_neg(res, res);
	76	}
	77	return;
	78	}
	79	else if (len <= 0)
	80	{
	81	acb_set_round(res, initial, prec);
	82	return;
	83	}
	84	}
	85
	86	TMP_START;
	87	t = TMP_ALLOC(sizeof(arb_struct) * len);
	88
	89	for (i = 0; i < len; i++)
	90	{
	91	vlo = y[2 * i * ystep];
	92	vhi = y[2 * i * ystep + 1];
	93
	94	arf_shallow_set_uiui(arb_midref(t + i), vhi, vlo);
	95
	96	MAG_EXP(arb_radref(t + i)) = 0;
	97	MAG_MAN(arb_radref(t + i)) = 0;
	98	}
	99
	100	arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec);
	101	arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec);
	102
	103	TMP_END;
	104	}
	105

+102

-0

acb/test/t-dot_fmpz.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_siui....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	acb_ptr x, y;
	26	fmpz * w;
	27	acb_t s1, s2, z;
	28	slong i, len, prec;
	29	int initial, subtract, revx, revy;
	30
	31	len = n_randint(state, 5);
	32	prec = 2 + n_randint(state, 200);
	33
	34	initial = n_randint(state, 2);
	35	subtract = n_randint(state, 2);
	36	revx = n_randint(state, 2);
	37	revy = n_randint(state, 2);
	38
	39	x = _acb_vec_init(len);
	40	y = _acb_vec_init(len);
	41	w = _fmpz_vec_init(len);
	42	acb_init(s1);
	43	acb_init(s2);
	44	acb_init(z);
	45
	46	for (i = 0; i < len; i++)
	47	{
	48	acb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	49	fmpz_randtest(w + i, state, 1 + n_randint(state, 200));
	50	acb_set_fmpz(y + i, w + i);
	51	}
	52
	53	acb_randtest(s1, state, 200, 10);
	54	acb_randtest(s2, state, 200, 10);
	55	acb_randtest(z, state, 200, 10);
	56
	57	acb_dot(s1, initial ? z : NULL, subtract,
	58	revx ? (x + len - 1) : x, revx ? -1 : 1,
	59	revy ? (y + len - 1) : y, revy ? -1 : 1,
	60	len, prec);
	61
	62	acb_dot_fmpz(s2, initial ? z : NULL, subtract,
	63	revx ? (x + len - 1) : x, revx ? -1 : 1,
	64	revy ? (w + len - 1) : w, revy ? -1 : 1,
	65	len, prec);
	66
	67	if (!acb_overlaps(s1, s2))
	68	{
	69	flint_printf("FAIL\n\n");
	70	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	71
	72	if (initial)
	73	{
	74	flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
	75	}
	76
	77	for (i = 0; i < len; i++)
	78	{
	79	flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
	80	flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
	81	}
	82	flint_printf("\n\n");
	83	flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	84	flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	85	flint_abort();
	86	}
	87
	88	acb_clear(s1);
	89	acb_clear(s2);
	90	acb_clear(z);
	91	_acb_vec_clear(x, len);
	92	_acb_vec_clear(y, len);
	93	_fmpz_vec_clear(w, len);
	94	}
	95
	96	flint_randclear(state);
	97	flint_cleanup();
	98	flint_printf("PASS\n");
	99	return EXIT_SUCCESS;
	100	}
	101

+102

-0

acb/test/t-dot_si.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_si....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	acb_ptr x, y;
	26	slong * w;
	27	acb_t s1, s2, z;
	28	slong i, len, prec;
	29	int initial, subtract, revx, revy;
	30
	31	len = n_randint(state, 5);
	32	prec = 2 + n_randint(state, 200);
	33
	34	initial = n_randint(state, 2);
	35	subtract = n_randint(state, 2);
	36	revx = n_randint(state, 2);
	37	revy = n_randint(state, 2);
	38
	39	x = _acb_vec_init(len);
	40	y = _acb_vec_init(len);
	41	w = flint_malloc(sizeof(ulong) * len);
	42	acb_init(s1);
	43	acb_init(s2);
	44	acb_init(z);
	45
	46	for (i = 0; i < len; i++)
	47	{
	48	acb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	49	w[i] = n_randtest(state);
	50	acb_set_si(y + i, w[i]);
	51	}
	52
	53	acb_randtest(s1, state, 200, 10);
	54	acb_randtest(s2, state, 200, 10);
	55	acb_randtest(z, state, 200, 10);
	56
	57	acb_dot(s1, initial ? z : NULL, subtract,
	58	revx ? (x + len - 1) : x, revx ? -1 : 1,
	59	revy ? (y + len - 1) : y, revy ? -1 : 1,
	60	len, prec);
	61
	62	acb_dot_si(s2, initial ? z : NULL, subtract,
	63	revx ? (x + len - 1) : x, revx ? -1 : 1,
	64	revy ? (w + len - 1) : w, revy ? -1 : 1,
	65	len, prec);
	66
	67	if (!acb_equal(s1, s2))
	68	{
	69	flint_printf("FAIL\n\n");
	70	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	71
	72	if (initial)
	73	{
	74	flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
	75	}
	76
	77	for (i = 0; i < len; i++)
	78	{
	79	flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
	80	flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
	81	}
	82	flint_printf("\n\n");
	83	flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	84	flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	85	flint_abort();
	86	}
	87
	88	acb_clear(s1);
	89	acb_clear(s2);
	90	acb_clear(z);
	91	_acb_vec_clear(x, len);
	92	_acb_vec_clear(y, len);
	93	flint_free(w);
	94	}
	95
	96	flint_randclear(state);
	97	flint_cleanup();
	98	flint_printf("PASS\n");
	99	return EXIT_SUCCESS;
	100	}
	101

+107

-0

acb/test/t-dot_siui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_siui....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	acb_ptr x, y;
	26	ulong * w;
	27	acb_t s1, s2, z;
	28	fmpz_t c;
	29	slong i, len, prec;
	30	int initial, subtract, revx, revy;
	31
	32	len = n_randint(state, 5);
	33	prec = 2 + n_randint(state, 200);
	34
	35	initial = n_randint(state, 2);
	36	subtract = n_randint(state, 2);
	37	revx = n_randint(state, 2);
	38	revy = n_randint(state, 2);
	39
	40	x = _acb_vec_init(len);
	41	y = _acb_vec_init(len);
	42	w = flint_malloc(2 * sizeof(ulong) * len);
	43	acb_init(s1);
	44	acb_init(s2);
	45	acb_init(z);
	46	fmpz_init(c);
	47
	48	for (i = 0; i < len; i++)
	49	{
	50	acb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	51	w[2 * i] = n_randtest(state);
	52	w[2 * i + 1] = n_randtest(state);
	53	fmpz_set_signed_uiui(c, w[2 * i + 1], w[2 * i]);
	54	acb_set_fmpz(y + i, c);
	55	}
	56
	57	acb_randtest(s1, state, 200, 10);
	58	acb_randtest(s2, state, 200, 10);
	59	acb_randtest(z, state, 200, 10);
	60
	61	acb_dot(s1, initial ? z : NULL, subtract,
	62	revx ? (x + len - 1) : x, revx ? -1 : 1,
	63	revy ? (y + len - 1) : y, revy ? -1 : 1,
	64	len, prec);
	65
	66	acb_dot_siui(s2, initial ? z : NULL, subtract,
	67	revx ? (x + len - 1) : x, revx ? -1 : 1,
	68	revy ? (w + 2 * len - 2) : w, revy ? -1 : 1,
	69	len, prec);
	70
	71	if (!acb_overlaps(s1, s2))
	72	{
	73	flint_printf("FAIL\n\n");
	74	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	75
	76	if (initial)
	77	{
	78	flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
	79	}
	80
	81	for (i = 0; i < len; i++)
	82	{
	83	flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
	84	flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
	85	}
	86	flint_printf("\n\n");
	87	flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	88	flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	89	flint_abort();
	90	}
	91
	92	acb_clear(s1);
	93	acb_clear(s2);
	94	acb_clear(z);
	95	_acb_vec_clear(x, len);
	96	_acb_vec_clear(y, len);
	97	flint_free(w);
	98	fmpz_clear(c);
	99	}
	100
	101	flint_randclear(state);
	102	flint_cleanup();
	103	flint_printf("PASS\n");
	104	return EXIT_SUCCESS;
	105	}
	106

+102

-0

acb/test/t-dot_ui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_ui....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	acb_ptr x, y;
	26	ulong * w;
	27	acb_t s1, s2, z;
	28	slong i, len, prec;
	29	int initial, subtract, revx, revy;
	30
	31	len = n_randint(state, 5);
	32	prec = 2 + n_randint(state, 200);
	33
	34	initial = n_randint(state, 2);
	35	subtract = n_randint(state, 2);
	36	revx = n_randint(state, 2);
	37	revy = n_randint(state, 2);
	38
	39	x = _acb_vec_init(len);
	40	y = _acb_vec_init(len);
	41	w = flint_malloc(sizeof(ulong) * len);
	42	acb_init(s1);
	43	acb_init(s2);
	44	acb_init(z);
	45
	46	for (i = 0; i < len; i++)
	47	{
	48	acb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	49	w[i] = n_randtest(state);
	50	acb_set_ui(y + i, w[i]);
	51	}
	52
	53	acb_randtest(s1, state, 200, 10);
	54	acb_randtest(s2, state, 200, 10);
	55	acb_randtest(z, state, 200, 10);
	56
	57	acb_dot(s1, initial ? z : NULL, subtract,
	58	revx ? (x + len - 1) : x, revx ? -1 : 1,
	59	revy ? (y + len - 1) : y, revy ? -1 : 1,
	60	len, prec);
	61
	62	acb_dot_ui(s2, initial ? z : NULL, subtract,
	63	revx ? (x + len - 1) : x, revx ? -1 : 1,
	64	revy ? (w + len - 1) : w, revy ? -1 : 1,
	65	len, prec);
	66
	67	if (!acb_equal(s1, s2))
	68	{
	69	flint_printf("FAIL\n\n");
	70	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	71
	72	if (initial)
	73	{
	74	flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
	75	}
	76
	77	for (i = 0; i < len; i++)
	78	{
	79	flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
	80	flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
	81	}
	82	flint_printf("\n\n");
	83	flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	84	flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	85	flint_abort();
	86	}
	87
	88	acb_clear(s1);
	89	acb_clear(s2);
	90	acb_clear(z);
	91	_acb_vec_clear(x, len);
	92	_acb_vec_clear(y, len);
	93	flint_free(w);
	94	}
	95
	96	flint_randclear(state);
	97	flint_cleanup();
	98	flint_printf("PASS\n");
	99	return EXIT_SUCCESS;
	100	}
	101

+107

-0

acb/test/t-dot_uiui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_uiui....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	acb_ptr x, y;
	26	ulong * w;
	27	acb_t s1, s2, z;
	28	fmpz_t c;
	29	slong i, len, prec;
	30	int initial, subtract, revx, revy;
	31
	32	len = n_randint(state, 5);
	33	prec = 2 + n_randint(state, 200);
	34
	35	initial = n_randint(state, 2);
	36	subtract = n_randint(state, 2);
	37	revx = n_randint(state, 2);
	38	revy = n_randint(state, 2);
	39
	40	x = _acb_vec_init(len);
	41	y = _acb_vec_init(len);
	42	w = flint_malloc(2 * sizeof(ulong) * len);
	43	acb_init(s1);
	44	acb_init(s2);
	45	acb_init(z);
	46	fmpz_init(c);
	47
	48	for (i = 0; i < len; i++)
	49	{
	50	acb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	51	w[2 * i] = n_randtest(state);
	52	w[2 * i + 1] = n_randtest(state);
	53	fmpz_set_uiui(c, w[2 * i + 1], w[2 * i]);
	54	acb_set_fmpz(y + i, c);
	55	}
	56
	57	acb_randtest(s1, state, 200, 10);
	58	acb_randtest(s2, state, 200, 10);
	59	acb_randtest(z, state, 200, 10);
	60
	61	acb_dot(s1, initial ? z : NULL, subtract,
	62	revx ? (x + len - 1) : x, revx ? -1 : 1,
	63	revy ? (y + len - 1) : y, revy ? -1 : 1,
	64	len, prec);
	65
	66	acb_dot_uiui(s2, initial ? z : NULL, subtract,
	67	revx ? (x + len - 1) : x, revx ? -1 : 1,
	68	revy ? (w + 2 * len - 2) : w, revy ? -1 : 1,
	69	len, prec);
	70
	71	if (!acb_overlaps(s1, s2))
	72	{
	73	flint_printf("FAIL\n\n");
	74	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	75
	76	if (initial)
	77	{
	78	flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z));
	79	}
	80
	81	for (i = 0; i < len; i++)
	82	{
	83	flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i));
	84	flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i));
	85	}
	86	flint_printf("\n\n");
	87	flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	88	flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	89	flint_abort();
	90	}
	91
	92	acb_clear(s1);
	93	acb_clear(s2);
	94	acb_clear(z);
	95	_acb_vec_clear(x, len);
	96	_acb_vec_clear(y, len);
	97	flint_free(w);
	98	fmpz_clear(c);
	99	}
	100
	101	flint_randclear(state);
	102	flint_cleanup();
	103	flint_printf("PASS\n");
	104	return EXIT_SUCCESS;
	105	}
	106

+11

-0

acb.h less more

601	601	void acb_approx_dot(acb_t res, const acb_t initial, int subtract,
602	602	acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec);
603	603
	604	void acb_dot_ui(acb_t res, const acb_t initial, int subtract,
	605	acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec);
	606	void acb_dot_si(acb_t res, const acb_t initial, int subtract,
	607	acb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec);
	608	void acb_dot_uiui(acb_t res, const acb_t initial, int subtract,
	609	acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec);
	610	void acb_dot_siui(acb_t res, const acb_t initial, int subtract,
	611	acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec);
	612	void acb_dot_fmpz(acb_t res, const acb_t initial, int subtract,
	613	acb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec);
	614
604	615	void acb_inv(acb_t z, const acb_t x, slong prec);
605	616
606	617	void acb_div(acb_t z, const acb_t x, const acb_t y, slong prec);

+58

-0

acb_elliptic/p_prime.c less more

	0	/*
	1	Copyright (C) 2021 Daniel Schultz
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_elliptic.h"
	12	#include "acb_modular.h"
	13
	14	void
	15	acb_elliptic_p_prime(acb_t r, const acb_t z, const acb_t tau, slong prec)
	16	{
	17	acb_struct tz[4];
	18	acb_t t1, t2, t3;
	19	int i, real;
	20
	21	real = acb_is_real(z) && arb_is_int_2exp_si(acb_realref(tau), -1) &&
	22	arb_is_positive(acb_imagref(tau));
	23
	24	acb_init(t1);
	25	acb_init(t2);
	26	acb_init(t3);
	27
	28	for (i = 0; i < 4; i++)
	29	acb_init(tz + i);
	30
	31	acb_modular_theta(tz + 0, tz + 1, tz + 2, tz + 3, z, tau, prec);
	32
	33	/* (-2pieta^2/theta1)^3theta2theta3theta4 /
	34	acb_const_pi(t2, prec);
	35	acb_mul_2exp_si(t2, t2, 1);
	36	acb_neg(t2, t2);
	37	acb_modular_eta(t3, tau, prec);
	38	acb_mul(t1, t3, t3, prec);
	39	acb_mul(t3, t1, t2, prec);
	40	acb_div(t1, t3, tz + 0, prec);
	41	acb_mul(t2, t1, t1, prec);
	42	acb_mul(t3, t1, t2, prec);
	43	acb_mul(t1, tz + 1, tz + 2, prec);
	44	acb_mul(t2, t1, tz + 3, prec);
	45	acb_mul(r, t3, t2, prec);
	46
	47	if (real)
	48	arb_zero(acb_imagref(r));
	49
	50	acb_clear(t1);
	51	acb_clear(t2);
	52	acb_clear(t3);
	53
	54	for (i = 0; i < 4; i++)
	55	acb_clear(tz + i);
	56	}
	57

+98

-0

acb_elliptic/test/t-p_p_prime.c less more

	0	/*
	1	Copyright (C) 2021 Daniel Schultz
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_elliptic.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("p_p_prime....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
	24	{
	25	acb_struct pj[2];
	26	acb_t tau, z, p, pp, g2, g3, t;
	27	slong prec;
	28
	29	acb_init(tau);
	30	acb_init(z);
	31	acb_init(p);
	32	acb_init(pp);
	33	acb_init(pj + 0);
	34	acb_init(pj + 1);
	35	acb_init(g2);
	36	acb_init(g3);
	37	acb_init(t);
	38
	39	prec = 2 + n_randint(state, 400);
	40
	41	acb_randtest(z, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
	42	acb_randtest(tau, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
	43	if (arf_sgn(arb_midref(acb_imagref(tau))) < 0)
	44	acb_neg(tau, tau);
	45
	46	acb_elliptic_p(p, z, tau, prec);
	47	acb_elliptic_p_prime(pp, z, tau, prec);
	48	acb_elliptic_p_jet(pj, z, tau, 2, prec);
	49
	50	if (!acb_overlaps(p, pj + 0) \|\| !acb_overlaps(pp, pj + 1))
	51	{
	52	flint_printf("FAIL (overlap)\n");
	53	flint_printf("tau = "); acb_printd(tau, 30); flint_printf("\n\n");
	54	flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
	55	flint_printf("p = "); acb_printd(p, 30); flint_printf("\n\n");
	56	flint_printf("pp = "); acb_printd(pp, 30); flint_printf("\n\n");
	57	flint_printf("pj0 = "); acb_printd(pj + 0, 30); flint_printf("\n\n");
	58	flint_printf("pj1 = "); acb_printd(pj + 1, 30); flint_printf("\n\n");
	59	flint_abort();
	60	}
	61
	62	acb_elliptic_invariants(g2, g3, tau, prec);
	63	acb_pow_ui(pj + 0, pp, 2, prec);
	64
	65	acb_mul(t, p, g2, prec);
	66	acb_add(t, t, g3, prec);
	67	acb_pow_ui(pj + 1, p, 3, prec);
	68	acb_mul_ui(pj + 1, pj + 1, 4, prec);
	69	acb_sub(pj + 1, pj + 1, t, prec);
	70
	71	if (!acb_overlaps(pj + 0, pj + 1))
	72	{
	73	flint_printf("FAIL (check pp^2 = 4p^3-g2*p-g3)\n");
	74	flint_printf("tau = "); acb_printd(tau, 30); flint_printf("\n\n");
	75	flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
	76	flint_printf("p = "); acb_printd(p, 30); flint_printf("\n\n");
	77	flint_printf("pp = "); acb_printd(pp, 30); flint_printf("\n\n");
	78	flint_abort();
	79	}
	80
	81	acb_clear(tau);
	82	acb_clear(z);
	83	acb_clear(p);
	84	acb_clear(pp);
	85	acb_clear(pj + 0);
	86	acb_clear(pj + 1);
	87	acb_clear(g2);
	88	acb_clear(g3);
	89	acb_clear(t);
	90	}
	91
	92	flint_randclear(state);
	93	flint_cleanup();
	94	flint_printf("PASS\n");
	95	return EXIT_SUCCESS;
	96	}
	97

+2

-0

acb_elliptic.h less more

49	49
50	50	void acb_elliptic_p(acb_t r, const acb_t z, const acb_t tau, slong prec);
51	51
	52	void acb_elliptic_p_prime(acb_t r, const acb_t z, const acb_t tau, slong prec);
	53
52	54	void acb_elliptic_p_jet(acb_ptr r, const acb_t z, const acb_t tau, slong len, slong prec);
53	55
54	56	void _acb_elliptic_p_series(acb_ptr res, acb_srcptr z, slong zlen, const acb_t tau, slong len, slong prec);

+114

-0

acb_hypgeom/rising_ui_rs.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12	#include "acb_hypgeom.h"
	13
	14	void
	15	acb_hypgeom_rising_ui_rs(acb_t res, const acb_t x, ulong n, ulong m, slong prec)
	16	{
	17	slong i, k, l, m0, climbs, climbs_max, wp;
	18	acb_ptr xpow;
	19	acb_t t, u;
	20	mp_ptr c;
	21	TMP_INIT;
	22
	23	if (n <= 1)
	24	{
	25	if (n == 0)
	26	acb_one(res);
	27	else
	28	acb_set_round(res, x, prec);
	29	return;
	30	}
	31
	32	TMP_START;
	33
	34	if (m == 0 \|\| m == -1)
	35	{
	36	if (n <= 6)
	37	m = 2;
	38	else if (n <= 16)
	39	m = 4;
	40	else if (n <= 40)
	41	m = 6;
	42	else
	43	{
	44	m0 = n_sqrt(n);
	45	m = 8 + 0.27 * pow(FLINT_MAX(0, prec - 1024), 0.4);
	46	m = FLINT_MIN(m, m0);
	47	m = FLINT_MIN(m, 64);
	48	}
	49	}
	50
	51	wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
	52
	53	climbs_max = FLINT_BIT_COUNT(n - 1) * m;
	54	c = TMP_ALLOC(sizeof(mp_limb_t) * climbs_max * m);
	55
	56	xpow = _acb_vec_init(m + 1);
	57	_acb_vec_set_powers(xpow, x, m + 1, wp);
	58	acb_init(t);
	59	acb_init(u);
	60
	61	for (k = 0; k < n; k += m)
	62	{
	63	l = FLINT_MIN(m, n - k);
	64	climbs = FLINT_BIT_COUNT(k + l - 1) * l;
	65	climbs = (climbs + FLINT_BITS - 1) / FLINT_BITS;
	66
	67	/* assumes l >= 2 */
	68	if (l == 1)
	69	{
	70	acb_add_ui(u, x, k, wp);
	71	}
	72	else
	73	{
	74	if (climbs == 1)
	75	{
	76	_arb_hypgeom_rising_coeffs_1(c, k, l);
	77	acb_dot_ui(u, xpow + l, 0, xpow, 1, c, 1, l, wp);
	78	}
	79	else if (climbs == 2)
	80	{
	81	_arb_hypgeom_rising_coeffs_2(c, k, l);
	82	acb_dot_uiui(u, xpow + l, 0, xpow, 1, c, 1, l, wp);
	83	}
	84	else
	85	{
	86	fmpz * f = (fmpz *) c;
	87
	88	for (i = 0; i < l; i++)
	89	fmpz_init(f + i);
	90
	91	_arb_hypgeom_rising_coeffs_fmpz(f, k, l);
	92
	93	acb_dot_fmpz(u, xpow + l, 0, xpow, 1, f, 1, l, wp);
	94
	95	for (i = 0; i < l; i++)
	96	fmpz_clear(f + i);
	97	}
	98	}
	99
	100	if (k == 0)
	101	acb_swap(t, u);
	102	else
	103	acb_mul(t, t, u, wp);
	104	}
	105
	106	acb_set_round(res, t, prec);
	107
	108	acb_clear(t);
	109	acb_clear(u);
	110	_acb_vec_clear(xpow, m + 1);
	111	TMP_END;
	112	}
	113

+76

-0

acb_hypgeom/test/t-rising_ui_rs.c less more

	0	/*
	1	Copyright (C) 2018 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "acb_hypgeom.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("rising_ui_rs....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
	24	{
	25	acb_t x, xk, y, ya, yb, yayb;
	26	ulong k, n, m1, m2, m3;
	27	slong prec;
	28
	29	prec = 2 + n_randint(state, 200);
	30	k = n_randint(state, 10);
	31	n = n_randint(state, 50);
	32	m1 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n + k, 1));
	33	m2 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(k, 1));
	34	m3 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n, 1));
	35
	36	acb_init(x);
	37	acb_init(xk);
	38	acb_init(y);
	39	acb_init(ya);
	40	acb_init(yb);
	41	acb_init(yayb);
	42
	43	acb_randtest(x, state, prec, 10);
	44	acb_add_ui(xk, x, k, prec);
	45
	46	acb_hypgeom_rising_ui_rs(y, x, n + k, m1, prec);
	47	acb_hypgeom_rising_ui_rs(ya, x, k, m2, prec);
	48	acb_hypgeom_rising_ui_rs(yb, xk, n, m3, prec);
	49	acb_mul(yayb, ya, yb, prec);
	50
	51	if (!acb_overlaps(y, yayb))
	52	{
	53	flint_printf("FAIL\n\n");
	54	flint_printf("k = %wu, n = %wu, m1 = %wu, m2 = %wu, m3 = %wu\n\n", k, n, m1, m2, m3);
	55	flint_printf("x = "); acb_printn(x, 100, 0); flint_printf("\n\n");
	56	flint_printf("y = "); acb_printn(y, 100, 0); flint_printf("\n\n");
	57	flint_printf("ya = "); acb_printn(ya, 100, 0); flint_printf("\n\n");
	58	flint_printf("yb = "); acb_printn(yb, 100, 0); flint_printf("\n\n");
	59	flint_printf("yayb = "); acb_printn(yayb, 100, 0); flint_printf("\n\n");
	60	flint_abort();
	61	}
	62
	63	acb_clear(x);
	64	acb_clear(xk);
	65	acb_clear(y);
	66	acb_clear(ya);
	67	acb_clear(yb);
	68	acb_clear(yayb);
	69	}
	70
	71	flint_randclear(state);
	72	flint_cleanup();
	73	flint_printf("PASS\n");
	74	return EXIT_SUCCESS;
	75	}

+2

-0

acb_hypgeom.h less more

17	17	#ifdef __cplusplus
18	18	extern "C" {
19	19	#endif
	20
	21	void acb_hypgeom_rising_ui_rs(acb_t res, const acb_t x, ulong n, ulong m, slong prec);
20	22
21	23	void acb_hypgeom_pfq_bound_factor(mag_t C,
22	24	acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, ulong n);

+1

-1

acb_mat/approx_mul.c less more

37	37	acb_mat_t T;
38	38	acb_mat_init(T, ar, bc);
39	39	acb_mat_approx_mul_classical(T, A, B, prec);
40		acb_mat_swap(T, C);
	40	acb_mat_swap_entrywise(T, C);
41	41	acb_mat_clear(T);
42	42	return;
43	43	}

+1

-1

acb_mat/mul_classical.c less more

44	44	acb_mat_t T;
45	45	acb_mat_init(T, ar, bc);
46	46	acb_mat_mul_classical(T, A, B, prec);
47		acb_mat_swap(T, C);
	47	acb_mat_swap_entrywise(T, C);
48	48	acb_mat_clear(T);
49	49	return;
50	50	}

+1

-1

acb_mat/mul_threaded.c less more

86	86	acb_mat_t T;
87	87	acb_mat_init(T, ar, bc);
88	88	acb_mat_mul_threaded(T, A, B, prec);
89		acb_mat_swap(T, C);
	89	acb_mat_swap_entrywise(T, C);
90	90	acb_mat_clear(T);
91	91	return;
92	92	}

+10

-0

acb_mat.h less more

65	65	*mat2 = t;
66	66	}
67	67
	68	ACB_MAT_INLINE void
	69	acb_mat_swap_entrywise(acb_mat_t mat1, acb_mat_t mat2)
	70	{
	71	slong i, j;
	72
	73	for (i = 0; i < acb_mat_nrows(mat1); i++)
	74	for (j = 0; j < acb_mat_ncols(mat1); j++)
	75	acb_swap(acb_mat_entry(mat2, i, j), acb_mat_entry(mat1, i, j));
	76	}
	77
68	78	/* Window matrices */
69	79
70	80	void acb_mat_window_init(acb_mat_t window, const acb_mat_t mat, slong r1, slong c1, slong r2, slong c2);

+4

-4

acb_modular/epsilon_arg.c less more

11	11	#include "acb_modular.h"
12	12
13	13	static int
14		fmpz_kronecker(const fmpz_t a, const fmpz_t b)
	14	fmpz_kronecker1(const fmpz_t a, const fmpz_t b)
15	15	{
16	16	if (fmpz_sgn(b) < 0)
17	17	{

19	19	fmpz_t t;
20	20	fmpz_init(t);
21	21	fmpz_neg(t, b);
22		r = fmpz_kronecker(a, t);
	22	r = fmpz_kronecker1(a, t);
23	23	fmpz_clear(t);
24	24	return r;
25	25	}

57	57
58	58	if (cc % 2 == 1)
59	59	{
60		u = fmpz_kronecker(a, c);
	60	u = fmpz_kronecker1(a, c);
61	61	aa = aabb + 2aacc - 3cc + ccdd(1-aa*aa);
62	62	}
63	63	else
64	64	{
65		u = fmpz_kronecker(c, a);
	65	u = fmpz_kronecker1(c, a);
66	66	aa = aabb - aacc + 3aa - 3 + ccdd(1-aaaa);
67	67	}
68	68

+145

-0

arb/dot_fmpz.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	void
	14	arb_dot_fmpz(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec)
	15	{
	16	arb_ptr t;
	17	slong i, ssize, size, tmp_size;
	18	mp_ptr ztmp;
	19	fmpz v;
	20	ulong av, al;
	21	unsigned int bc;
	22	TMP_INIT;
	23
	24	/* todo: fast fma and fmma (len=2) code */
	25	if (len <= 1)
	26	{
	27	if (initial == NULL)
	28	{
	29	if (len <= 0)
	30	arb_zero(res);
	31	else
	32	{
	33	arb_mul_fmpz(res, x, y, prec);
	34	if (subtract)
	35	arb_neg(res, res);
	36	}
	37	return;
	38	}
	39	else if (len <= 0)
	40	{
	41	arb_set_round(res, initial, prec);
	42	return;
	43	}
	44	}
	45
	46	TMP_START;
	47	t = TMP_ALLOC(sizeof(arb_struct) * len);
	48
	49	tmp_size = 0;
	50	for (i = 0; i < len; i++)
	51	{
	52	v = y[i * ystep];
	53
	54	MAG_EXP(arb_radref(t + i)) = 0;
	55	MAG_MAN(arb_radref(t + i)) = 0;
	56
	57	if (v == 0)
	58	{
	59	ARF_XSIZE(arb_midref(t + i)) = 0;
	60	ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO;
	61	}
	62	else if (!COEFF_IS_MPZ(v))
	63	{
	64	av = FLINT_ABS(v);
	65	count_leading_zeros(bc, av);
	66
	67	ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
	68	ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc;
	69	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, v < 0);
	70	}
	71	else
	72	{
	73	__mpz_struct * z = COEFF_TO_PTR(v);
	74
	75	ssize = z->_mp_size;
	76	size = FLINT_ABS(ssize);
	77
	78	av = z->_mp_d[size - 1];
	79	count_leading_zeros(bc, av);
	80
	81	if (size == 1)
	82	{
	83	ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
	84	ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc;
	85	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, ssize < 0);
	86	}
	87	else if (size == 2)
	88	{
	89	al = z->_mp_d[0];
	90
	91	ARF_EXP(arb_midref(t + i)) = 2 * FLINT_BITS - bc;
	92
	93	if (bc != 0)
	94	{
	95	av = (av << bc) \| (al >> (FLINT_BITS - bc));
	96	al = al << bc;
	97	}
	98
	99	ARF_NOPTR_D(arb_midref(t + i))[0] = al;
	100	ARF_NOPTR_D(arb_midref(t + i))[1] = av;
	101	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(2, ssize < 0);
	102	}
	103	else
	104	{
	105	if (bc != 0)
	106	{
	107	tmp_size += size;
	108	/* use to flag tmp where we need tmp storage */
	109	MAG_MAN(arb_radref(t + i)) = bc;
	110	}
	111
	112	ARF_EXP(arb_midref(t + i)) = size * FLINT_BITS - bc;
	113	ARF_PTR_D(arb_midref(t + i)) = z->_mp_d;
	114	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(size, ssize < 0);
	115	}
	116	}
	117	}
	118
	119	if (tmp_size != 0)
	120	{
	121	ztmp = TMP_ALLOC(sizeof(mp_limb_t) * tmp_size);
	122
	123	for (i = 0; i < len; i++)
	124	{
	125	bc = MAG_MAN(arb_radref(t + i));
	126
	127	if (bc != 0)
	128	{
	129	size = ARF_SIZE(arb_midref(t + i));
	130
	131	mpn_lshift(ztmp, ARF_PTR_D(arb_midref(t + i)), size, bc);
	132	ARF_PTR_D(arb_midref(t + i)) = ztmp;
	133	ztmp += size;
	134	}
	135
	136	MAG_MAN(arb_radref(t + i)) = 0;
	137	}
	138	}
	139
	140	arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec);
	141
	142	TMP_END;
	143	}
	144

+76

-0

arb/dot_si.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	void
	14	arb_dot_si(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec)
	15	{
	16	arb_ptr t;
	17	slong i;
	18	slong v;
	19	ulong av;
	20	unsigned int bc;
	21	TMP_INIT;
	22
	23	/* todo: fast fma and fmma (len=2) code */
	24	if (len <= 1)
	25	{
	26	if (initial == NULL)
	27	{
	28	if (len <= 0)
	29	arb_zero(res);
	30	else
	31	{
	32	arb_mul_si(res, x, y[0], prec);
	33	if (subtract)
	34	arb_neg(res, res);
	35	}
	36	return;
	37	}
	38	else if (len <= 0)
	39	{
	40	arb_set_round(res, initial, prec);
	41	return;
	42	}
	43	}
	44
	45	TMP_START;
	46	t = TMP_ALLOC(sizeof(arb_struct) * len);
	47
	48	for (i = 0; i < len; i++)
	49	{
	50	v = y[i * ystep];
	51
	52	if (v == 0)
	53	{
	54	ARF_XSIZE(arb_midref(t + i)) = 0;
	55	ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO;
	56	}
	57	else
	58	{
	59	av = FLINT_ABS(v);
	60	count_leading_zeros(bc, av);
	61
	62	ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
	63	ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc;
	64	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, v < 0);
	65	}
	66
	67	MAG_EXP(arb_radref(t + i)) = 0;
	68	MAG_MAN(arb_radref(t + i)) = 0;
	69	}
	70
	71	arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec);
	72
	73	TMP_END;
	74	}
	75

+113

-0

arb/dot_siui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	static void
	14	arf_shallow_set_siui(arf_t res, ulong vhi, ulong vlo)
	15	{
	16	int negative;
	17	unsigned int bc;
	18
	19	negative = ((slong) vhi) < 0;
	20
	21	if (negative)
	22	{
	23	vhi = -vhi - (vlo != 0);
	24	vlo = -vlo;
	25	}
	26
	27	if (vhi == 0)
	28	{
	29	if (vlo == 0)
	30	{
	31	ARF_XSIZE(res) = 0;
	32	ARF_EXP(res) = ARF_EXP_ZERO;
	33	}
	34	else
	35	{
	36	count_leading_zeros(bc, vlo);
	37	ARF_EXP(res) = FLINT_BITS - bc;
	38	ARF_NOPTR_D(res)[0] = vlo << bc;
	39	ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative);
	40	}
	41	}
	42	else if (vlo == 0)
	43	{
	44	count_leading_zeros(bc, vhi);
	45	ARF_EXP(res) = 2 * FLINT_BITS - bc;
	46	ARF_NOPTR_D(res)[0] = vhi << bc;
	47	ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative);
	48	}
	49	else
	50	{
	51	count_leading_zeros(bc, vhi);
	52	ARF_EXP(res) = 2 * FLINT_BITS - bc;
	53	ARF_NOPTR_D(res)[0] = vlo << bc;
	54	if (bc == 0)
	55	ARF_NOPTR_D(res)[1] = vhi;
	56	else
	57	ARF_NOPTR_D(res)[1] = (vhi << bc) \| (vlo >> (FLINT_BITS - bc));
	58	ARF_XSIZE(res) = ARF_MAKE_XSIZE(2, negative);
	59	}
	60	}
	61
	62	void
	63	arb_dot_siui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	64	{
	65	arb_ptr t;
	66	slong i;
	67	ulong vhi, vlo;
	68	TMP_INIT;
	69
	70	/* todo: fast fma and fmma (len=2) code */
	71	if (len <= 1)
	72	{
	73	if (initial == NULL)
	74	{
	75	if (len <= 0)
	76	arb_zero(res);
	77	else
	78	{
	79	arf_t t;
	80	arf_shallow_set_siui(t, y[1], y[0]);
	81	arb_mul_arf(res, x, t, prec);
	82	if (subtract)
	83	arb_neg(res, res);
	84	}
	85	return;
	86	}
	87	else if (len <= 0)
	88	{
	89	arb_set_round(res, initial, prec);
	90	return;
	91	}
	92	}
	93
	94	TMP_START;
	95	t = TMP_ALLOC(sizeof(arb_struct) * len);
	96
	97	for (i = 0; i < len; i++)
	98	{
	99	vlo = y[2 * i * ystep];
	100	vhi = y[2 * i * ystep + 1];
	101
	102	arf_shallow_set_siui(arb_midref(t + i), vhi, vlo);
	103
	104	MAG_EXP(arb_radref(t + i)) = 0;
	105	MAG_MAN(arb_radref(t + i)) = 0;
	106	}
	107
	108	arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec);
	109
	110	TMP_END;
	111	}
	112

+74

-0

arb/dot_ui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	void
	14	arb_dot_ui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	15	{
	16	arb_ptr t;
	17	slong i;
	18	ulong v;
	19	unsigned int bc;
	20	TMP_INIT;
	21
	22	/* todo: fast fma and fmma (len=2) code */
	23	if (len <= 1)
	24	{
	25	if (initial == NULL)
	26	{
	27	if (len <= 0)
	28	arb_zero(res);
	29	else
	30	{
	31	arb_mul_ui(res, x, y[0], prec);
	32	if (subtract)
	33	arb_neg(res, res);
	34	}
	35	return;
	36	}
	37	else if (len <= 0)
	38	{
	39	arb_set_round(res, initial, prec);
	40	return;
	41	}
	42	}
	43
	44	TMP_START;
	45	t = TMP_ALLOC(sizeof(arb_struct) * len);
	46
	47	for (i = 0; i < len; i++)
	48	{
	49	v = y[i * ystep];
	50
	51	if (v == 0)
	52	{
	53	ARF_XSIZE(arb_midref(t + i)) = 0;
	54	ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO;
	55	}
	56	else
	57	{
	58	count_leading_zeros(bc, v);
	59
	60	ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc;
	61	ARF_NOPTR_D(arb_midref(t + i))[0] = v << bc;
	62	ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, 0);
	63	}
	64
	65	MAG_EXP(arb_radref(t + i)) = 0;
	66	MAG_MAN(arb_radref(t + i)) = 0;
	67	}
	68
	69	arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec);
	70
	71	TMP_END;
	72	}
	73

+104

-0

arb/dot_uiui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	static void
	14	arf_shallow_set_uiui(arf_t res, ulong vhi, ulong vlo)
	15	{
	16	unsigned int bc;
	17
	18	if (vhi == 0)
	19	{
	20	if (vlo == 0)
	21	{
	22	ARF_XSIZE(res) = 0;
	23	ARF_EXP(res) = ARF_EXP_ZERO;
	24	}
	25	else
	26	{
	27	count_leading_zeros(bc, vlo);
	28	ARF_EXP(res) = FLINT_BITS - bc;
	29	ARF_NOPTR_D(res)[0] = vlo << bc;
	30	ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, 0);
	31	}
	32	}
	33	else if (vlo == 0)
	34	{
	35	count_leading_zeros(bc, vhi);
	36	ARF_EXP(res) = 2 * FLINT_BITS - bc;
	37	ARF_NOPTR_D(res)[0] = vhi << bc;
	38	ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, 0);
	39	}
	40	else
	41	{
	42	count_leading_zeros(bc, vhi);
	43	ARF_EXP(res) = 2 * FLINT_BITS - bc;
	44	ARF_NOPTR_D(res)[0] = vlo << bc;
	45	if (bc == 0)
	46	ARF_NOPTR_D(res)[1] = vhi;
	47	else
	48	ARF_NOPTR_D(res)[1] = (vhi << bc) \| (vlo >> (FLINT_BITS - bc));
	49	ARF_XSIZE(res) = ARF_MAKE_XSIZE(2, 0);
	50	}
	51	}
	52
	53	void
	54	arb_dot_uiui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	55	{
	56	arb_ptr t;
	57	slong i;
	58	ulong vhi, vlo;
	59	TMP_INIT;
	60
	61	/* todo: fast fma and fmma (len=2) code */
	62	if (len <= 1)
	63	{
	64	if (initial == NULL)
	65	{
	66	if (len <= 0)
	67	arb_zero(res);
	68	else
	69	{
	70	arf_t t;
	71	arf_shallow_set_uiui(t, y[1], y[0]);
	72	arb_mul_arf(res, x, t, prec);
	73	if (subtract)
	74	arb_neg(res, res);
	75	}
	76	return;
	77	}
	78	else if (len <= 0)
	79	{
	80	arb_set_round(res, initial, prec);
	81	return;
	82	}
	83	}
	84
	85	TMP_START;
	86	t = TMP_ALLOC(sizeof(arb_struct) * len);
	87
	88	for (i = 0; i < len; i++)
	89	{
	90	vlo = y[2 * i * ystep];
	91	vhi = y[2 * i * ystep + 1];
	92
	93	arf_shallow_set_uiui(arb_midref(t + i), vhi, vlo);
	94
	95	MAG_EXP(arb_radref(t + i)) = 0;
	96	MAG_MAN(arb_radref(t + i)) = 0;
	97	}
	98
	99	arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec);
	100
	101	TMP_END;
	102	}
	103

+283

-7

arb/gamma.c less more

639	639	}
640	640	}
641	641
	642	slong _arb_compute_bs_exponents(slong * tab, slong n);
	643	slong _arb_get_exp_pos(const slong * tab, slong step);
	644
	645	static void
	646	bsplit2(arb_t P, arb_t Q, const fmpz_t zp, const fmpz_t zq,
	647	const slong * xexp, arb_srcptr xpow,
	648	ulong N, slong a, slong b, int cont, slong prec)
	649	{
	650	if (b - a == 1)
	651	{
	652	fmpz_t t;
	653	fmpz_init(t);
	654	fmpz_set(t, zp);
	655	fmpz_addmul_ui(t, zq, a + 1);
	656	arb_set_fmpz(P, t);
	657	arb_set(Q, P);
	658	fmpz_clear(t);
	659	}
	660	else
	661	{
	662	arb_t Pb, Qb;
	663	slong step, i, m;
	664
	665	arb_init(Pb);
	666	arb_init(Qb);
	667
	668	step = (b - a) / 2;
	669	m = a + step;
	670
	671	bsplit2(P, Q, zp, zq, xexp, xpow, N, a, m, 1, prec);
	672	bsplit2(Pb, Qb, zp, zq, xexp, xpow, N, m, b, 1, prec);
	673
	674	arb_mul(P, P, Pb, prec);
	675	arb_mul(Q, Q, Pb, prec);
	676
	677	i = (step == 1) ? 0 : _arb_get_exp_pos(xexp, step);
	678	arb_addmul(Q, Qb, xpow + i, prec);
	679
	680	arb_clear(Pb);
	681	arb_clear(Qb);
	682	}
	683	}
	684
	685	static void
	686	bsplit3(arb_t P, arb_t Q, const fmpz_t zp, const fmpz_t zq,
	687	const slong * xexp, arb_srcptr xpow,
	688	ulong N, slong a, slong b, int cont, slong prec)
	689	{
	690	if (b - a == 1)
	691	{
	692	fmpz_t t;
	693	fmpz_init(t);
	694	arb_set(P, xpow + 0); /* N zq */
	695	fmpz_set(t, zp);
	696	fmpz_submul_ui(t, zq, a + 1); /* zp - (a + 1) zq */
	697	arb_set_fmpz(Q, t);
	698	fmpz_clear(t);
	699	}
	700	else
	701	{
	702	arb_t Pb, Qb;
	703	slong step, i, m;
	704
	705	arb_init(Pb);
	706	arb_init(Qb);
	707
	708	step = (b - a) / 2;
	709	m = a + step;
	710
	711	bsplit3(P, Q, zp, zq, xexp, xpow, N, a, m, 1, prec);
	712	bsplit3(Pb, Qb, zp, zq, xexp, xpow, N, m, b, 1, prec);
	713
	714	i = _arb_get_exp_pos(xexp, m - a);
	715
	716	arb_mul(P, P, xpow + i, prec);
	717	if (b - m != m - a)
	718	arb_mul(P, P, xpow + 0, prec);
	719
	720	arb_addmul(P, Q, Pb, prec);
	721
	722	if (cont)
	723	{
	724	arb_mul(Q, Q, Qb, prec);
	725	}
	726	else
	727	{
	728	i = _arb_get_exp_pos(xexp, m - a);
	729
	730	arb_mul(Q, xpow + i, xpow + i, prec);
	731	if (b - m != m - a)
	732	arb_mul(Q, Q, xpow + 0, prec);
	733	}
	734
	735	arb_clear(Pb);
	736	arb_clear(Qb);
	737	}
	738	}
	739
	740	double d_lambertw_branch1(double x);
	741
	742	static ulong
	743	more_trailing_zeros(ulong N)
	744	{
	745	ulong bc, N2;
	746
	747	bc = FLINT_BIT_COUNT(N);
	748
	749	if (bc >= 8)
	750	{
	751	N2 = (N >> (bc - 5)) << (bc - 5);
	752	N2 += ((N2 != N) << (bc - 5));
	753	return N2;
	754	}
	755	else
	756	{
	757	return N;
	758	}
	759	}
	760
	761	#define C_LOG2 0.69314718055994530942
	762	#define C_EXP1 2.7182818284590452354
	763
	764	static void
	765	build_bsplit_power_table(arb_ptr xpow, const slong * xexp, slong len, slong prec)
	766	{
	767	slong i;
	768
	769	for (i = 1; i < len; i++)
	770	{
	771	if (xexp[i] == 2 * xexp[i-1])
	772	{
	773	arb_mul(xpow + i, xpow + i - 1, xpow + i - 1, prec);
	774	}
	775	else if (xexp[i] == 2 * xexp[i-2]) /* prefer squaring if possible */
	776	{
	777	arb_mul(xpow + i, xpow + i - 2, xpow + i - 2, prec);
	778	}
	779	else if (xexp[i] == 2 * xexp[i-1] + 1)
	780	{
	781	arb_mul(xpow + i, xpow + i - 1, xpow + i - 1, prec);
	782	arb_mul(xpow + i, xpow + i, xpow, prec);
	783	}
	784	else if (xexp[i] == 2 * xexp[i-2] + 1)
	785	{
	786	arb_mul(xpow + i, xpow + i - 2, xpow + i - 2, prec);
	787	arb_mul(xpow + i, xpow + i, xpow, prec);
	788	}
	789	else
	790	{
	791	flint_printf("power table has the wrong structure!\n");
	792	flint_abort();
	793	}
	794	}
	795	}
	796
	797	/* assumes z in [1, 2] */
	798	static void
	799	arb_gamma_fmpq_general_off1(arb_t res, const fmpq_t z, slong prec)
	800	{
	801	slong wp, N, n, n2, length, length2, wp2;
	802	double alpha;
	803	arb_t P, Q;
	804	slong xexp, xexp2;
	805	arb_ptr xpow;
	806	mag_t err, err2;
	807
	808	wp = prec + 30;
	809
	810	alpha = 0.52; /* tuning parameter between 0.5 and 1 */
	811
	812	N = alpha * C_LOG2 * wp;
	813	N = more_trailing_zeros(N);
	814	alpha = N / (C_LOG2 * wp);
	815
	816	/* Terms in convergent series */
	817	n = (1 - alpha) / d_lambertw((1 - alpha) / (alpha * C_EXP1)) * C_LOG2 * wp;
	818
	819	/* Precision and terms in asymptotic series */
	820	wp2 = wp * (1 - alpha);
	821	wp2 = FLINT_MAX(wp2, 30);
	822	n2 = (alpha - 1) / d_lambertw_branch1((alpha - 1) / (alpha * C_EXP1)) * C_LOG2 * wp;
	823	n2 = FLINT_MAX(n2, 2); /* binary splitting correctness */
	824
	825	mag_init(err);
	826	mag_init(err2);
	827	arb_init(P);
	828	arb_init(Q);
	829
	830	/* compute the powers of x = Nzq that will appear (at least x^1) /
	831	xexp = flint_calloc(2 * FLINT_BITS, sizeof(slong));
	832	xexp2 = flint_calloc(2 * FLINT_BITS, sizeof(slong));
	833
	834	length = _arb_compute_bs_exponents(xexp, n);
	835	length2 = _arb_compute_bs_exponents(xexp2, n2);
	836
	837	xpow = _arb_vec_init(FLINT_MAX(length, length2));
	838
	839	arb_set_fmpz(xpow + 0, fmpq_denref(z));
	840	arb_mul_ui(xpow + 0, xpow + 0, N, wp);
	841
	842	build_bsplit_power_table(xpow, xexp, length, wp);
	843
	844	/* 1F1(1, 1+z, N) */
	845	bsplit2(P, Q, fmpq_numref(z), fmpq_denref(z), xexp, xpow, N, 0, n, 0, wp);
	846	arb_div(P, Q, P, wp);
	847
	848	/* Convergent series error bound: N^n / n! (1 + (N/n) + ...) */
	849	mag_set_ui(err, N);
	850	mag_pow_ui(err, err, n);
	851	mag_rfac_ui(err2, n);
	852	mag_mul(err, err, err2);
	853	mag_set_ui(err2, N);
	854	mag_div_ui(err2, err2, n);
	855	mag_geom_series(err2, err2, 0);
	856	mag_mul(err, err, err2);
	857	arb_add_error_mag(P, err);
	858
	859	/* divide 1F1 by z */
	860	arb_mul_fmpz(P, P, fmpq_denref(z), wp);
	861	arb_div_fmpz(P, P, fmpq_numref(z), wp);
	862	arb_swap(res, P);
	863
	864	build_bsplit_power_table(xpow, xexp2, length2, wp2);
	865
	866	bsplit3(P, Q, fmpq_numref(z), fmpq_denref(z), xexp2, xpow, N, 0, n2, 0, wp2);
	867	arb_div(P, P, Q, wp2);
	868
	869	/* 2F0 error bound (bounded by first omitted term) */
	870	mag_fac_ui(err, n2);
	871	mag_set_ui_lower(err2, N);
	872	mag_pow_ui_lower(err2, err2, n2);
	873	mag_div(err, err, err2);
	874	arb_add_error_mag(P, err);
	875
	876	/* N^z * exp(-N) * (1F1/z + 2F0/N) */
	877	arb_div_ui(P, P, N, wp2);
	878
	879	arb_add(res, res, P, wp);
	880	arb_set_ui(Q, N);
	881	arb_pow_fmpq(Q, Q, z, wp);
	882	arb_mul(res, res, Q, wp);
	883
	884	arb_set_si(Q, -N);
	885	arb_exp(Q, Q, wp);
	886	arb_mul(res, res, Q, wp);
	887
	888	_arb_vec_clear(xpow, FLINT_MAX(length, length2));
	889	flint_free(xexp);
	890	flint_free(xexp2);
	891
	892	arb_clear(P);
	893	arb_clear(Q);
	894	mag_clear(err);
	895	mag_clear(err2);
	896	}
	897
	898	/* assumes z in (0, 1] */
	899	void
	900	arb_gamma_fmpq_hyp(arb_t res, const fmpq_t z, slong prec)
	901	{
	902	fmpq_t t;
	903	fmpq_init(t);
	904	fmpq_add_ui(t, z, 1);
	905	arb_gamma_fmpq_general_off1(res, t, prec);
	906	arb_mul_fmpz(res, res, fmpq_denref(z), prec + 30);
	907	arb_div_fmpz(res, res, fmpq_numref(z), prec);
	908	fmpq_clear(t);
	909	}
	910
642	911	void
643	912	arb_gamma_fmpq_outward(arb_t y, const fmpq_t x, slong prec)
644	913	{

683	952	}
684	953	else
685	954	{
686		flint_printf("arb_gamma_fmpq: invalid fraction\n");
687		flint_abort();
	955	arb_gamma_fmpq_hyp(t, a, prec);
688	956	}
689	957
690	958	/* argument reduction */

743	1011	}
744	1012	}
745	1013
	1014	if (q != 1 && prec > 7000 + 300 * fmpz_bits(fmpq_denref(x)) &&
	1015	(slong) fmpz_bits(fmpq_numref(x)) - (slong) fmpz_bits(fmpq_denref(x)) < FLINT_BITS &&
	1016	fabs(fmpq_get_d(x)) < 0.03 * prec * sqrt(prec))
	1017	{
	1018	arb_gamma_fmpq_outward(y, x, prec);
	1019	return;
	1020	}
	1021
746	1022	arb_gamma_fmpq_stirling(y, x, prec);
747	1023	}
748	1024

783	1059	}
784	1060	else
785	1061	{
786		/* fast gamma(n), gamma(n/2) or gamma(n/4) */
	1062	/* fast gamma(n), gamma(n/2) or gamma(n/4), ... */
787	1063	if (arf_cmpabs_2exp_si(mid, prec) < 0 &&
788		arf_is_int_2exp_si(mid, -2))
	1064	(arf_is_int_2exp_si(mid, -2) \|\| (prec > 1000 && arf_is_int_2exp_si(mid, -1000))))
789	1065	{
790	1066	fmpq_t a;
791	1067	fmpq_init(a);

916	1192	return;
917	1193	}
918	1194
919		/* fast gamma(n), gamma(n/2) or gamma(n/4) */
920		if (arb_is_exact(x) && arf_is_int_2exp_si(arb_midref(x), -2) &&
921		arf_cmpabs_2exp_si(arb_midref(x), prec) < 0)
	1195	/* fast gamma(n), gamma(n/2) or gamma(n/4), ... */
	1196	if (arb_is_exact(x) && arf_cmpabs_2exp_si(arb_midref(x), prec) < 0 &&
	1197	(arf_is_int_2exp_si(arb_midref(x), -2) \|\| (prec > 10000 && arf_is_int_2exp_si(arb_midref(x), -1000))))
922	1198	{
923	1199	fmpq_t a;
924	1200	fmpq_init(a);

+55

-35

arb/get_str.c less more

451	451
452	452	arb_get_str_parts(&negative, &mid_digits, mid_exp, &rad_digits, rad_exp, x, n, more);
453	453
454		skip_mid = mid_digits[0] == '0';
455		skip_rad = (rad_digits[0] == '0') \|\| ((flags & ARB_STR_NO_RADIUS) && !skip_mid);
456
457		_arb_digits_as_float_str(&mid_digits, mid_exp, -4, FLINT_MAX(6, n - 1));
458		_arb_digits_as_float_str(&rad_digits, rad_exp, -2, 2);
459
460		if (skip_rad)
461		{
462		res = flint_malloc(strlen(mid_digits) + 2);
463
464		if (negative)
465		strcpy(res, "-");
466		else
467		strcpy(res, "");
468
469		strcat(res, mid_digits);
470		}
471		else if (skip_mid)
472		{
473		res = flint_malloc(strlen(rad_digits) + 7);
474
475		strcpy(res, "[+/- ");
476		strcat(res, rad_digits);
477		strcat(res, "]");
	454	if ((flags & ARB_STR_NO_RADIUS) && mid_digits[0] == '0')
	455	{
	456	fmpz_add_ui(rad_exp, rad_exp, strlen(rad_digits));
	457
	458	res = flint_malloc(fmpz_sizeinbase(rad_exp, 10) + 4);
	459	res[0] = '0';
	460	res[1] = 'e';
	461	if (fmpz_sgn(rad_exp) >= 0)
	462	{
	463	res[2] = '+';
	464	fmpz_get_str(res + 3, 10, rad_exp);
	465	}
	466	else
	467	{
	468	fmpz_get_str(res + 2, 10, rad_exp);
	469	}
478	470	}
479	471	else
480	472	{
481		res = flint_malloc(strlen(mid_digits) + strlen(rad_digits) + 9);
482
483		strcpy(res, "[");
484
485		if (negative)
486		strcat(res, "-");
487
488		strcat(res, mid_digits);
489		strcat(res, " +/- ");
490		strcat(res, rad_digits);
491		strcat(res, "]");
	473	skip_mid = mid_digits[0] == '0';
	474	skip_rad = (rad_digits[0] == '0') \|\| ((flags & ARB_STR_NO_RADIUS) && !skip_mid);
	475
	476	_arb_digits_as_float_str(&mid_digits, mid_exp, -4, FLINT_MAX(6, n - 1));
	477	_arb_digits_as_float_str(&rad_digits, rad_exp, -2, 2);
	478
	479	if (skip_rad)
	480	{
	481	res = flint_malloc(strlen(mid_digits) + 2);
	482
	483	if (negative)
	484	strcpy(res, "-");
	485	else
	486	strcpy(res, "");
	487
	488	strcat(res, mid_digits);
	489	}
	490	else if (skip_mid)
	491	{
	492	res = flint_malloc(strlen(rad_digits) + 7);
	493
	494	strcpy(res, "[+/- ");
	495	strcat(res, rad_digits);
	496	strcat(res, "]");
	497	}
	498	else
	499	{
	500	res = flint_malloc(strlen(mid_digits) + strlen(rad_digits) + 9);
	501
	502	strcpy(res, "[");
	503
	504	if (negative)
	505	strcat(res, "-");
	506
	507	strcat(res, mid_digits);
	508	strcat(res, " +/- ");
	509	strcat(res, rad_digits);
	510	strcat(res, "]");
	511	}
492	512	}
493	513
494	514	if (condense)

+1

-1

arb/lambertw.c less more

62	62	1.0,-3556.4306263369027831,1.4761527435056145298e6,
63	63	-9.8425904825010893103e7,7.0373606710750560344e8 };
64	64
65		static double
	65	double
66	66	d_lambertw_branch1(double x)
67	67	{
68	68	double w, u;

+102

-0

arb/test/t-dot_fmpz.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_siui....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	arb_ptr x, y;
	26	fmpz * w;
	27	arb_t s1, s2, z;
	28	slong i, len, prec;
	29	int initial, subtract, revx, revy;
	30
	31	len = n_randint(state, 5);
	32	prec = 2 + n_randint(state, 200);
	33
	34	initial = n_randint(state, 2);
	35	subtract = n_randint(state, 2);
	36	revx = n_randint(state, 2);
	37	revy = n_randint(state, 2);
	38
	39	x = _arb_vec_init(len);
	40	y = _arb_vec_init(len);
	41	w = _fmpz_vec_init(len);
	42	arb_init(s1);
	43	arb_init(s2);
	44	arb_init(z);
	45
	46	for (i = 0; i < len; i++)
	47	{
	48	arb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	49	fmpz_randtest(w + i, state, 1 + n_randint(state, 200));
	50	arb_set_fmpz(y + i, w + i);
	51	}
	52
	53	arb_randtest(s1, state, 200, 10);
	54	arb_randtest(s2, state, 200, 10);
	55	arb_randtest(z, state, 200, 10);
	56
	57	arb_dot(s1, initial ? z : NULL, subtract,
	58	revx ? (x + len - 1) : x, revx ? -1 : 1,
	59	revy ? (y + len - 1) : y, revy ? -1 : 1,
	60	len, prec);
	61
	62	arb_dot_fmpz(s2, initial ? z : NULL, subtract,
	63	revx ? (x + len - 1) : x, revx ? -1 : 1,
	64	revy ? (w + len - 1) : w, revy ? -1 : 1,
	65	len, prec);
	66
	67	if (!arb_overlaps(s1, s2))
	68	{
	69	flint_printf("FAIL\n\n");
	70	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	71
	72	if (initial)
	73	{
	74	flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z));
	75	}
	76
	77	for (i = 0; i < len; i++)
	78	{
	79	flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i));
	80	flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i));
	81	}
	82	flint_printf("\n\n");
	83	flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	84	flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	85	flint_abort();
	86	}
	87
	88	arb_clear(s1);
	89	arb_clear(s2);
	90	arb_clear(z);
	91	_arb_vec_clear(x, len);
	92	_arb_vec_clear(y, len);
	93	_fmpz_vec_clear(w, len);
	94	}
	95
	96	flint_randclear(state);
	97	flint_cleanup();
	98	flint_printf("PASS\n");
	99	return EXIT_SUCCESS;
	100	}
	101

+102

-0

arb/test/t-dot_si.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_si....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	arb_ptr x, y;
	26	slong * w;
	27	arb_t s1, s2, z;
	28	slong i, len, prec;
	29	int initial, subtract, revx, revy;
	30
	31	len = n_randint(state, 5);
	32	prec = 2 + n_randint(state, 200);
	33
	34	initial = n_randint(state, 2);
	35	subtract = n_randint(state, 2);
	36	revx = n_randint(state, 2);
	37	revy = n_randint(state, 2);
	38
	39	x = _arb_vec_init(len);
	40	y = _arb_vec_init(len);
	41	w = flint_malloc(sizeof(ulong) * len);
	42	arb_init(s1);
	43	arb_init(s2);
	44	arb_init(z);
	45
	46	for (i = 0; i < len; i++)
	47	{
	48	arb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	49	w[i] = n_randtest(state);
	50	arb_set_si(y + i, w[i]);
	51	}
	52
	53	arb_randtest(s1, state, 200, 10);
	54	arb_randtest(s2, state, 200, 10);
	55	arb_randtest(z, state, 200, 10);
	56
	57	arb_dot(s1, initial ? z : NULL, subtract,
	58	revx ? (x + len - 1) : x, revx ? -1 : 1,
	59	revy ? (y + len - 1) : y, revy ? -1 : 1,
	60	len, prec);
	61
	62	arb_dot_si(s2, initial ? z : NULL, subtract,
	63	revx ? (x + len - 1) : x, revx ? -1 : 1,
	64	revy ? (w + len - 1) : w, revy ? -1 : 1,
	65	len, prec);
	66
	67	if (!arb_equal(s1, s2))
	68	{
	69	flint_printf("FAIL\n\n");
	70	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	71
	72	if (initial)
	73	{
	74	flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z));
	75	}
	76
	77	for (i = 0; i < len; i++)
	78	{
	79	flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i));
	80	flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i));
	81	}
	82	flint_printf("\n\n");
	83	flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	84	flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	85	flint_abort();
	86	}
	87
	88	arb_clear(s1);
	89	arb_clear(s2);
	90	arb_clear(z);
	91	_arb_vec_clear(x, len);
	92	_arb_vec_clear(y, len);
	93	flint_free(w);
	94	}
	95
	96	flint_randclear(state);
	97	flint_cleanup();
	98	flint_printf("PASS\n");
	99	return EXIT_SUCCESS;
	100	}
	101

+107

-0

arb/test/t-dot_siui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_siui....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	arb_ptr x, y;
	26	ulong * w;
	27	arb_t s1, s2, z;
	28	fmpz_t c;
	29	slong i, len, prec;
	30	int initial, subtract, revx, revy;
	31
	32	len = n_randint(state, 5);
	33	prec = 2 + n_randint(state, 200);
	34
	35	initial = n_randint(state, 2);
	36	subtract = n_randint(state, 2);
	37	revx = n_randint(state, 2);
	38	revy = n_randint(state, 2);
	39
	40	x = _arb_vec_init(len);
	41	y = _arb_vec_init(len);
	42	w = flint_malloc(2 * sizeof(ulong) * len);
	43	arb_init(s1);
	44	arb_init(s2);
	45	arb_init(z);
	46	fmpz_init(c);
	47
	48	for (i = 0; i < len; i++)
	49	{
	50	arb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	51	w[2 * i] = n_randtest(state);
	52	w[2 * i + 1] = n_randtest(state);
	53	fmpz_set_signed_uiui(c, w[2 * i + 1], w[2 * i]);
	54	arb_set_fmpz(y + i, c);
	55	}
	56
	57	arb_randtest(s1, state, 200, 10);
	58	arb_randtest(s2, state, 200, 10);
	59	arb_randtest(z, state, 200, 10);
	60
	61	arb_dot(s1, initial ? z : NULL, subtract,
	62	revx ? (x + len - 1) : x, revx ? -1 : 1,
	63	revy ? (y + len - 1) : y, revy ? -1 : 1,
	64	len, prec);
	65
	66	arb_dot_siui(s2, initial ? z : NULL, subtract,
	67	revx ? (x + len - 1) : x, revx ? -1 : 1,
	68	revy ? (w + 2 * len - 2) : w, revy ? -1 : 1,
	69	len, prec);
	70
	71	if (!arb_overlaps(s1, s2))
	72	{
	73	flint_printf("FAIL\n\n");
	74	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	75
	76	if (initial)
	77	{
	78	flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z));
	79	}
	80
	81	for (i = 0; i < len; i++)
	82	{
	83	flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i));
	84	flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i));
	85	}
	86	flint_printf("\n\n");
	87	flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	88	flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	89	flint_abort();
	90	}
	91
	92	arb_clear(s1);
	93	arb_clear(s2);
	94	arb_clear(z);
	95	_arb_vec_clear(x, len);
	96	_arb_vec_clear(y, len);
	97	flint_free(w);
	98	fmpz_clear(c);
	99	}
	100
	101	flint_randclear(state);
	102	flint_cleanup();
	103	flint_printf("PASS\n");
	104	return EXIT_SUCCESS;
	105	}
	106

+102

-0

arb/test/t-dot_ui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_ui....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	arb_ptr x, y;
	26	ulong * w;
	27	arb_t s1, s2, z;
	28	slong i, len, prec;
	29	int initial, subtract, revx, revy;
	30
	31	len = n_randint(state, 5);
	32	prec = 2 + n_randint(state, 200);
	33
	34	initial = n_randint(state, 2);
	35	subtract = n_randint(state, 2);
	36	revx = n_randint(state, 2);
	37	revy = n_randint(state, 2);
	38
	39	x = _arb_vec_init(len);
	40	y = _arb_vec_init(len);
	41	w = flint_malloc(sizeof(ulong) * len);
	42	arb_init(s1);
	43	arb_init(s2);
	44	arb_init(z);
	45
	46	for (i = 0; i < len; i++)
	47	{
	48	arb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	49	w[i] = n_randtest(state);
	50	arb_set_ui(y + i, w[i]);
	51	}
	52
	53	arb_randtest(s1, state, 200, 10);
	54	arb_randtest(s2, state, 200, 10);
	55	arb_randtest(z, state, 200, 10);
	56
	57	arb_dot(s1, initial ? z : NULL, subtract,
	58	revx ? (x + len - 1) : x, revx ? -1 : 1,
	59	revy ? (y + len - 1) : y, revy ? -1 : 1,
	60	len, prec);
	61
	62	arb_dot_ui(s2, initial ? z : NULL, subtract,
	63	revx ? (x + len - 1) : x, revx ? -1 : 1,
	64	revy ? (w + len - 1) : w, revy ? -1 : 1,
	65	len, prec);
	66
	67	if (!arb_equal(s1, s2))
	68	{
	69	flint_printf("FAIL\n\n");
	70	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	71
	72	if (initial)
	73	{
	74	flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z));
	75	}
	76
	77	for (i = 0; i < len; i++)
	78	{
	79	flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i));
	80	flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i));
	81	}
	82	flint_printf("\n\n");
	83	flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	84	flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	85	flint_abort();
	86	}
	87
	88	arb_clear(s1);
	89	arb_clear(s2);
	90	arb_clear(z);
	91	_arb_vec_clear(x, len);
	92	_arb_vec_clear(y, len);
	93	flint_free(w);
	94	}
	95
	96	flint_randclear(state);
	97	flint_cleanup();
	98	flint_printf("PASS\n");
	99	return EXIT_SUCCESS;
	100	}
	101

+107

-0

arb/test/t-dot_uiui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	int main()
	14	{
	15	slong iter;
	16	flint_rand_t state;
	17
	18	flint_printf("dot_uiui....");
	19	fflush(stdout);
	20
	21	flint_randinit(state);
	22
	23	for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
	24	{
	25	arb_ptr x, y;
	26	ulong * w;
	27	arb_t s1, s2, z;
	28	fmpz_t c;
	29	slong i, len, prec;
	30	int initial, subtract, revx, revy;
	31
	32	len = n_randint(state, 5);
	33	prec = 2 + n_randint(state, 200);
	34
	35	initial = n_randint(state, 2);
	36	subtract = n_randint(state, 2);
	37	revx = n_randint(state, 2);
	38	revy = n_randint(state, 2);
	39
	40	x = _arb_vec_init(len);
	41	y = _arb_vec_init(len);
	42	w = flint_malloc(2 * sizeof(ulong) * len);
	43	arb_init(s1);
	44	arb_init(s2);
	45	arb_init(z);
	46	fmpz_init(c);
	47
	48	for (i = 0; i < len; i++)
	49	{
	50	arb_randtest(x + i, state, 2 + n_randint(state, 200), 10);
	51	w[2 * i] = n_randtest(state);
	52	w[2 * i + 1] = n_randtest(state);
	53	fmpz_set_uiui(c, w[2 * i + 1], w[2 * i]);
	54	arb_set_fmpz(y + i, c);
	55	}
	56
	57	arb_randtest(s1, state, 200, 10);
	58	arb_randtest(s2, state, 200, 10);
	59	arb_randtest(z, state, 200, 10);
	60
	61	arb_dot(s1, initial ? z : NULL, subtract,
	62	revx ? (x + len - 1) : x, revx ? -1 : 1,
	63	revy ? (y + len - 1) : y, revy ? -1 : 1,
	64	len, prec);
	65
	66	arb_dot_uiui(s2, initial ? z : NULL, subtract,
	67	revx ? (x + len - 1) : x, revx ? -1 : 1,
	68	revy ? (w + 2 * len - 2) : w, revy ? -1 : 1,
	69	len, prec);
	70
	71	if (!arb_overlaps(s1, s2))
	72	{
	73	flint_printf("FAIL\n\n");
	74	flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec);
	75
	76	if (initial)
	77	{
	78	flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z));
	79	}
	80
	81	for (i = 0; i < len; i++)
	82	{
	83	flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i));
	84	flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i));
	85	}
	86	flint_printf("\n\n");
	87	flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n");
	88	flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n");
	89	flint_abort();
	90	}
	91
	92	arb_clear(s1);
	93	arb_clear(s2);
	94	arb_clear(z);
	95	_arb_vec_clear(x, len);
	96	_arb_vec_clear(y, len);
	97	flint_free(w);
	98	fmpz_clear(c);
	99	}
	100
	101	flint_randclear(state);
	102	flint_cleanup();
	103	flint_printf("PASS\n");
	104	return EXIT_SUCCESS;
	105	}
	106

+52

-0

arb/test/t-gamma_fmpq.c less more

71	71	fmpq_clear(q);
72	72	}
73	73
	74	for (iter = 0; iter < 50 * arb_test_multiplier(); iter++)
	75	{
	76	arb_t r, s;
	77	fmpq_t q;
	78	slong accuracy, prec;
	79
	80	prec = 2 + n_randint(state, 25000);
	81
	82	arb_init(r);
	83	arb_init(s);
	84	fmpq_init(q);
	85
	86	fmpz_randtest(fmpq_numref(q), state, 3 + n_randlimb(state) % 30);
	87	fmpz_randtest_not_zero(fmpq_denref(q), state, 3 + n_randlimb(state) % 30);
	88	fmpq_canonicalise(q);
	89
	90	arb_gamma_fmpq(r, q, prec);
	91
	92	arb_set_fmpq(s, q, prec);
	93	arb_gamma(s, s, prec);
	94
	95	if (!arb_overlaps(r, s))
	96	{
	97	flint_printf("FAIL: containment\n\n");
	98	flint_printf("prec = %wd\n", prec);
	99	flint_printf("q = "); fmpq_print(q); flint_printf("\n\n");
	100	flint_printf("r = "); arb_printd(r, prec / 3.33); flint_printf("\n\n");
	101	flint_printf("s = "); arb_printd(s, prec / 3.33); flint_printf("\n\n");
	102	flint_abort();
	103	}
	104
	105	if (!(fmpz_is_one(fmpq_denref(q)) && fmpz_sgn(fmpq_numref(q)) <= 0) &&
	106	fabs(fmpq_get_d(q)) < 10.0)
	107	{
	108	accuracy = arb_rel_accuracy_bits(r);
	109
	110	if (accuracy < prec - 6)
	111	{
	112	flint_printf("FAIL: poor accuracy\n\n");
	113	flint_printf("prec = %wd\n", prec);
	114	flint_printf("q = "); fmpq_print(q); flint_printf("\n\n");
	115	flint_printf("r = "); arb_printn(r, prec / 3.33, 0); flint_printf("\n\n");
	116	flint_printf("s = "); arb_printn(s, prec / 3.33, 0); flint_printf("\n\n");
	117	flint_abort();
	118	}
	119	}
	120
	121	arb_clear(r);
	122	arb_clear(s);
	123	fmpq_clear(q);
	124	}
	125
74	126	flint_randclear(state);
75	127	flint_cleanup();
76	128	flint_printf("PASS\n");

+60

-0

arb/test/t-get_str.c less more

74	74	arb_clear(y);
75	75	}
76	76
	77	/* test ARB_STR_NO_RADIUS */
	78	{
	79	arb_t x;
	80	char * s;
	81
	82	arb_init(x);
	83
	84	arb_set_str(x, "3.1415926535897932", 53);
	85	s = arb_get_str(x, 10, ARB_STR_NO_RADIUS);
	86	if (strcmp(s, "3.141592654"))
	87	{
	88	flint_printf("FAIL (ARB_STR_NO_RADIUS)\n");
	89	flint_printf("%s\n", s);
	90	flint_abort();
	91	}
	92	flint_free(s);
	93
	94	arb_set_str(x, "+/- 3.45e-10", 53);
	95	s = arb_get_str(x, 10, ARB_STR_NO_RADIUS);
	96	if (strcmp(s, "0e-9"))
	97	{
	98	flint_printf("FAIL (ARB_STR_NO_RADIUS)\n");
	99	flint_printf("%s\n", s);
	100	flint_abort();
	101	}
	102	flint_free(s);
	103
	104	arb_set_str(x, "+/- 3.45e10", 53);
	105	s = arb_get_str(x, 10, ARB_STR_NO_RADIUS);
	106	if (strcmp(s, "0e+11"))
	107	{
	108	flint_printf("FAIL (ARB_STR_NO_RADIUS)\n");
	109	flint_printf("%s\n", s);
	110	flint_abort();
	111	}
	112	flint_free(s);
	113
	114	arb_set_str(x, "5e10 +/- 6e10", 53);
	115	s = arb_get_str(x, 10, ARB_STR_NO_RADIUS);
	116	if (strcmp(s, "0e+12"))
	117	{
	118	flint_printf("FAIL (ARB_STR_NO_RADIUS)\n");
	119	flint_printf("%s\n", s);
	120	flint_abort();
	121	}
	122	flint_free(s);
	123
	124	arb_set_str(x, "5e-100000000000000000002 +/- 5e-100000000000000000002", 53);
	125	s = arb_get_str(x, 10, ARB_STR_NO_RADIUS);
	126	if (strcmp(s, "0e-100000000000000000000"))
	127	{
	128	flint_printf("FAIL (ARB_STR_NO_RADIUS)\n");
	129	flint_printf("%s\n", s);
	130	flint_abort();
	131	}
	132	flint_free(s);
	133
	134	arb_clear(x);
	135	}
	136
77	137	flint_randclear(state);
78	138	flint_cleanup();
79	139	flint_printf("PASS\n");

+85

-0

arb/test/t-urandom.c less more

	0	/*
	1	Copyright (C) 2021 Albin Ahlbäck
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	#define N 10000
	14
	15	int main()
	16	{
	17	slong iter;
	18	slong prec;
	19	flint_rand_t state;
	20	arb_ptr rand;
	21	arb_t m; /* mean */
	22	arb_t s; /* variance */
	23	arb_t mp;
	24	arb_t sp;
	25	arb_t tmp;
	26
	27	flint_printf("urandom....");
	28	fflush(stdout);
	29
	30	flint_randinit(state);
	31	arb_init(m);
	32	arb_init(s);
	33	arb_init(mp);
	34	arb_init(sp);
	35	arb_init(tmp);
	36	rand = _arb_vec_init(N);
	37	prec = 299;
	38
	39	for (iter = 0; iter < N; iter++)
	40	{
	41	arb_urandom(rand + iter, state, prec);
	42	arb_add(m, m, rand + iter, prec);
	43	}
	44
	45	arb_div_si(m, m, N, prec);
	46
	47	for (iter = 0; iter < N; iter++)
	48	{
	49	arb_sub(tmp, rand + iter, m, prec);
	50	arb_sqr(tmp, tmp, prec);
	51	arb_add(s, s, tmp, prec);
	52	}
	53
	54	arb_div_si(s, s, N, prec);
	55
	56	/* one percent deviation */
	57	arb_set_str(mp, "0.5 +/- 0.005", prec);
	58	arb_set_str(sp, "0.083333 +/- 0.00083", prec);
	59
	60	if (!arb_contains(mp, m))
	61	{
	62	flint_printf("FAIL: mean\n\n");
	63	flint_printf("m = "); arb_printd(m, 15); flint_printf("\n\n");
	64	flint_abort();
	65	}
	66
	67	if (!arb_contains(sp, s))
	68	{
	69	flint_printf("FAIL: variance\n\n");
	70	flint_printf("s = "); arb_printd(s, 15); flint_printf("\n\n");
	71	flint_abort();
	72	}
	73
	74	_arb_vec_clear(rand, N);
	75	arb_clear(m);
	76	arb_clear(s);
	77	arb_clear(mp);
	78	arb_clear(sp);
	79	arb_clear(tmp);
	80	flint_randclear(state);
	81	flint_cleanup();
	82	flint_printf("PASS\n");
	83	return EXIT_SUCCESS;
	84	}

+36

-0

arb/urandom.c less more

	0	/*
	1	Copyright (C) 2021 Albin Ahlbäck
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb.h"
	12
	13	void
	14	arb_urandom(arb_t x, flint_rand_t state, slong bits)
	15	{
	16	slong prec = bits;
	17	fmpz_t n;
	18	fmpz_t t;
	19
	20	prec += 128;
	21
	22	fmpz_init(n);
	23	fmpz_one(n);
	24	fmpz_mul_2exp(n, n, (ulong) prec);
	25
	26	fmpz_init(t);
	27	fmpz_randm(t, state, n);
	28
	29	arb_set_round_fmpz(x, t, bits);
	30	arb_mul_2exp_si(x, x, -prec);
	31
	32	fmpz_clear(n);
	33	fmpz_clear(t);
	34	}
	35

+1

-1

arb/version.c less more

10	10
11	11	#include "arb.h"
12	12
13		const char * arb_version = "2.19.0";
	13	const char * arb_version = "2.20.0";

+27

-2

arb.h less more

26	26	#endif
27	27
28	28	#define __ARB_VERSION 2
29		#define __ARB_VERSION_MINOR 19
	29	#define __ARB_VERSION_MINOR 20
30	30	#define __ARB_VERSION_PATCHLEVEL 0
31		#define ARB_VERSION "2.19.0"
	31	#define ARB_VERSION "2.20.0"
32	32	#define __ARB_RELEASE (__ARB_VERSION * 10000 + \
33	33	__ARB_VERSION_MINOR * 100 + \
34	34	__ARB_VERSION_PATCHLEVEL)

332	332
333	333	void arb_randtest_special(arb_t x, flint_rand_t state, slong prec, slong mag_bits);
334	334
	335	void arb_urandom(arb_t x, flint_rand_t state, slong prec);
	336
335	337	void arb_add_error_arf(arb_t x, const arf_t err);
336	338
337	339	void arb_add_error_2exp_si(arb_t x, slong err);

431	433
432	434	void arb_approx_dot(arb_t res, const arb_t initial, int subtract,
433	435	arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec);
	436
	437	void arb_dot_ui(arb_t res, const arb_t initial, int subtract,
	438	arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec);
	439	void arb_dot_si(arb_t res, const arb_t initial, int subtract,
	440	arb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec);
	441	void arb_dot_uiui(arb_t res, const arb_t initial, int subtract,
	442	arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec);
	443	void arb_dot_siui(arb_t res, const arb_t initial, int subtract,
	444	arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec);
	445	void arb_dot_fmpz(arb_t res, const arb_t initial, int subtract,
	446	arb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec);
434	447
435	448	void arb_div(arb_t z, const arb_t x, const arb_t y, slong prec);
436	449	void arb_div_arf(arb_t z, const arb_t x, const arf_t y, slong prec);

639	652	}
640	653
641	654	ARB_INLINE void
	655	_arb_vec_printn(arb_srcptr vec, slong len, slong ndigits, ulong flags)
	656	{
	657	slong i;
	658	for (i = 0; i < len; i++)
	659	{
	660	arb_printn(vec + i, ndigits, flags);
	661	if (i < len - 1)
	662	flint_printf(", ");
	663	}
	664	}
	665
	666	ARB_INLINE void
642	667	_arb_vec_zero(arb_ptr A, slong n)
643	668	{
644	669	slong i;

+32

-1

arb_fmpz_poly/complex_roots.c less more

11	11	#include "arb_fmpz_poly.h"
12	12	#include "flint/profiler.h"
13	13
14		int check_accuracy(acb_ptr vec, slong len, slong prec)
	14	static int check_accuracy(acb_ptr vec, slong len, slong prec)
15	15	{
16	16	slong i;
17	17

19	19	{
20	20	if (acb_rel_accuracy_bits(vec + i) < prec)
21	21	return 0;
	22	}
	23
	24	return 1;
	25	}
	26
	27	static int check_isolation(acb_srcptr roots, slong len)
	28	{
	29	slong i, j;
	30
	31	for (i = 0; i < len; i++)
	32	{
	33	if (arf_sgn(arb_midref(acb_imagref(roots + i))) >= 0)
	34	{
	35	for (j = i + 1; j < len; j++)
	36	{
	37	if (arf_sgn(arb_midref(acb_imagref(roots + j))) >= 0)
	38	{
	39	if (acb_overlaps(roots + i, roots + j))
	40	return 0;
	41	}
	42	}
	43	}
22	44	}
23	45
24	46	return 1;

148	170	arb_zero(acb_imagref(roots + i));
149	171	}
150	172
	173	if (!check_isolation(roots, deg))
	174	{
	175	/* extremely unlikely */
	176	if (flags & ARB_FMPZ_POLY_ROOTS_VERBOSE)
	177	flint_printf("isolation failure!\n");
	178
	179	continue;
	180	}
	181
151	182	if (flags & ARB_FMPZ_POLY_ROOTS_VERBOSE)
152	183	flint_printf("done!\n");
153	184

+5

-6

arb_fmpz_poly/evaluate_acb_rectangular.c less more

14	14	_arb_fmpz_poly_evaluate_acb_rectangular(acb_t y, const fmpz * poly,
15	15	slong len, const acb_t x, slong prec)
16	16	{
17		slong i, j, m, r;
	17	slong i, m, r;
18	18	acb_ptr xs;
19	19	acb_t s, t, c;
20	20

35	35	_acb_vec_set_powers(xs, x, m + 1, prec);
36	36
37	37	acb_set_fmpz(y, poly + (r - 1) * m);
38		for (j = 1; (r - 1) * m + j < len; j++)
39		acb_addmul_fmpz(y, xs + j, poly + (r - 1) * m + j, prec);
	38	acb_dot_fmpz(y, y, 0, xs + 1, 1,
	39	poly + (r - 1) * m + 1, 1, len - (r - 1) * m - 1, prec);
40	40
41	41	for (i = r - 2; i >= 0; i--)
42	42	{
43	43	acb_set_fmpz(s, poly + i * m);
44		for (j = 1; j < m; j++)
45		acb_addmul_fmpz(s, xs + j, poly + i * m + j, prec);
46
	44	acb_dot_fmpz(s, s, 0, xs + 1, 1,
	45	poly + i * m + 1, 1, m - 1, prec);
47	46	acb_mul(y, y, xs + m, prec);
48	47	acb_add(y, y, s, prec);
49	48	}

+15

-8

arb_fmpz_poly/evaluate_arb.c less more

14	14	_arb_fmpz_poly_evaluate_arb(arb_t res, const fmpz * f, slong len,
15	15	const arb_t x, slong prec)
16	16	{
17		if ((prec >= 1024) && (len >= 5 + 20000 / prec))
	17	if (len >= 6 && len >= 5 + 2500 / (FLINT_MAX(prec, 64) + 64))
18	18	{
19		slong fbits;
	19	/* todo: improve this tuning? */
	20	if (prec > 1024)
	21	{
	22	slong fbits;
	23	fbits = _fmpz_vec_max_bits(f, len);
	24	fbits = FLINT_ABS(fbits);
20	25
21		fbits = _fmpz_vec_max_bits(f, len);
	26	if (fbits > prec / 2)
	27	{
	28	_arb_fmpz_poly_evaluate_arb_horner(res, f, len, x, prec);
	29	return;
	30	}
	31	}
22	32
23		if (fbits <= prec / 2)
24		{
25		_arb_fmpz_poly_evaluate_arb_rectangular(res, f, len, x, prec);
26		return;
27		}
	33	_arb_fmpz_poly_evaluate_arb_rectangular(res, f, len, x, prec);
	34	return;
28	35	}
29	36
30	37	_arb_fmpz_poly_evaluate_arb_horner(res, f, len, x, prec);

+5

-6

arb_fmpz_poly/evaluate_arb_rectangular.c less more

14	14	_arb_fmpz_poly_evaluate_arb_rectangular(arb_t y, const fmpz * poly,
15	15	slong len, const arb_t x, slong prec)
16	16	{
17		slong i, j, m, r;
	17	slong i, m, r;
18	18	arb_ptr xs;
19	19	arb_t s, t, c;
20	20

35	35	_arb_vec_set_powers(xs, x, m + 1, prec);
36	36
37	37	arb_set_fmpz(y, poly + (r - 1) * m);
38		for (j = 1; (r - 1) * m + j < len; j++)
39		arb_addmul_fmpz(y, xs + j, poly + (r - 1) * m + j, prec);
	38	arb_dot_fmpz(y, y, 0, xs + 1, 1,
	39	poly + (r - 1) * m + 1, 1, len - (r - 1) * m - 1, prec);
40	40
41	41	for (i = r - 2; i >= 0; i--)
42	42	{
43	43	arb_set_fmpz(s, poly + i * m);
44		for (j = 1; j < m; j++)
45		arb_addmul_fmpz(s, xs + j, poly + i * m + j, prec);
46
	44	arb_dot_fmpz(s, s, 0, xs + 1, 1,
	45	poly + i * m + 1, 1, m - 1, prec);
47	46	arb_mul(y, y, xs + m, prec);
48	47	arb_add(y, y, s, prec);
49	48	}

+18

-1

arb_fmpz_poly/test/t-complex_roots.c less more

19	19	arb_poly_t rpoly;
20	20	arb_t lead;
21	21
22		slong i, num_real, num_upper, deg;
	22	slong i, j, num_real, num_upper, deg;
23	23
24	24	deg = fmpz_poly_degree(poly);
25	25

61	61	flint_printf("rpoly:\n");
62	62	arb_poly_printd(rpoly, 30); flint_printf("\n\n");
63	63	flint_abort();
	64	}
	65
	66	for (i = 0; i < deg; i++)
	67	{
	68	for (j = i + 1; j < deg; j++)
	69	{
	70	if (acb_overlaps(roots + i, roots + j))
	71	{
	72	flint_printf("FAIL! (isolation)\n");
	73	flint_printf("deg = %wd, num_real = %wd, num_upper = %wd\n\n", deg, num_real, num_upper);
	74	for (i = 0; i < deg; i++)
	75	{
	76	acb_printn(roots + i, 30, 0);
	77	flint_printf("\n");
	78	}
	79	}
	80	}
64	81	}
65	82
66	83	_arb_vec_clear(real, num_real);

+53

-0

arb_hypgeom/rising_ui.c less more

	0	/*
	1	Copyright (C) 2014 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	void
	14	arb_hypgeom_rising_ui(arb_t y, const arb_t x, ulong n, slong prec)
	15	{
	16	if (n < FLINT_MAX(prec, 100))
	17	{
	18	arb_hypgeom_rising_ui_rec(y, x, n, prec);
	19	}
	20	else
	21	{
	22	arb_t t;
	23	arb_init(t);
	24	arb_add_ui(t, x, n, prec);
	25	arb_gamma(t, t, prec);
	26	arb_rgamma(y, x, prec);
	27	arb_mul(y, y, t, prec);
	28	arb_clear(t);
	29	}
	30	}
	31
	32	void
	33	arb_hypgeom_rising(arb_t y, const arb_t x, const arb_t n, slong prec)
	34	{
	35	if (arb_is_int(n) && arf_sgn(arb_midref(n)) >= 0 &&
	36	arf_cmpabs_ui(arb_midref(n), FLINT_MAX(prec, 100)) < 0)
	37	{
	38	arb_hypgeom_rising_ui_rec(y, x,
	39	arf_get_si(arb_midref(n), ARF_RND_DOWN), prec);
	40	}
	41	else
	42	{
	43	arb_t t;
	44	arb_init(t);
	45	arb_add(t, x, n, prec);
	46	arb_gamma(t, t, prec);
	47	arb_rgamma(y, x, prec);
	48	arb_mul(y, y, t, prec);
	49	arb_clear(t);
	50	}
	51	}
	52

+69

-0

arb_hypgeom/rising_ui_bs.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	static void
	14	bsplit(arb_t y, const arb_t x, ulong a, ulong b, slong prec)
	15	{
	16	if (b - a <= 16)
	17	{
	18	if (a == 0)
	19	{
	20	arb_hypgeom_rising_ui_forward(y, x, b, prec);
	21	}
	22	else
	23	{
	24	arb_add_ui(y, x, a, prec);
	25	arb_hypgeom_rising_ui_forward(y, y, b - a, prec);
	26	}
	27	}
	28	else
	29	{
	30	arb_t t, u;
	31	ulong m = a + (b - a) / 2;
	32
	33	arb_init(t);
	34	arb_init(u);
	35
	36	bsplit(t, x, a, m, prec);
	37	bsplit(u, x, m, b, prec);
	38
	39	arb_mul(y, t, u, prec);
	40
	41	arb_clear(t);
	42	arb_clear(u);
	43	}
	44	}
	45
	46	void
	47	arb_hypgeom_rising_ui_bs(arb_t res, const arb_t x, ulong n, slong prec)
	48	{
	49	if (n <= 1)
	50	{
	51	if (n == 0)
	52	arb_one(res);
	53	else
	54	arb_set_round(res, x, prec);
	55	return;
	56	}
	57
	58	{
	59	arb_t t;
	60	slong wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
	61
	62	arb_init(t);
	63	bsplit(t, x, 0, n, wp);
	64	arb_set_round(res, t, prec);
	65	arb_clear(t);
	66	}
	67	}
	68

+83

-0

arb_hypgeom/rising_ui_forward.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	int
	14	_arf_increment_fast(arf_t x, slong prec)
	15	{
	16	if (arf_sgn(x) > 0)
	17	{
	18	mp_limb_t hi, v, cy;
	19	mp_ptr xptr;
	20	mp_size_t xn;
	21	slong xexp;
	22
	23	xexp = ARF_EXP(x);
	24
	25	if (xexp >= 1 && xexp <= FLINT_BITS - 1)
	26	{
	27	ARF_GET_MPN_READONLY(xptr, xn, x);
	28
	29	hi = xptr[xn - 1];
	30	v = hi + (UWORD(1) << (FLINT_BITS - xexp));
	31	cy = v < hi;
	32
	33	if (cy == 0)
	34	{
	35	xptr[xn - 1] = v;
	36	return 0;
	37	}
	38	}
	39	}
	40
	41	return arf_add_ui(x, x, 1, prec, ARF_RND_DOWN);
	42	}
	43
	44	void
	45	_arb_increment_fast(arb_t x, slong prec)
	46	{
	47	if (_arf_increment_fast(arb_midref(x), prec))
	48	arf_mag_add_ulp(arb_radref(x), arb_radref(x), arb_midref(x), prec);
	49	}
	50
	51	void
	52	arb_hypgeom_rising_ui_forward(arb_t res, const arb_t x, ulong n, slong prec)
	53	{
	54	arb_t t;
	55	ulong k;
	56	slong wp;
	57
	58	if (n <= 1)
	59	{
	60	if (n == 0)
	61	arb_one(res);
	62	else
	63	arb_set_round(res, x, prec);
	64	return;
	65	}
	66
	67	wp = prec + FLINT_BIT_COUNT(n);
	68
	69	arb_init(t);
	70
	71	arb_add_ui(t, x, 1, wp);
	72	arb_mul(res, x, t, (n == 2) ? prec : wp);
	73
	74	for (k = 2; k < n; k++)
	75	{
	76	_arb_increment_fast(t, wp);
	77	arb_mul(res, res, t, k == (n - 1) ? prec : wp);
	78	}
	79
	80	arb_clear(t);
	81	}
	82

+44

-0

arb_hypgeom/rising_ui_jet.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	void
	14	arb_hypgeom_rising_ui_jet(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
	15	{
	16	if (n <= 7)
	17	{
	18	arb_hypgeom_rising_ui_jet_powsum(res, x, n, len, prec);
	19	}
	20	else if (len == 2)
	21	{
	22	if (n <= 30 \|\| arb_bits(x) >= prec / 128)
	23	arb_hypgeom_rising_ui_jet_rs(res, x, n, 0, len, prec);
	24	else
	25	arb_hypgeom_rising_ui_jet_bs(res, x, n, len, prec);
	26	}
	27	else
	28	{
	29	if (n <= 20 \|\| (n <= 200 && prec > 400 * n && arb_bits(x) >= prec / 4))
	30	{
	31	arb_hypgeom_rising_ui_jet_powsum(res, x, n, len, prec);
	32	}
	33	else if (len >= 64 \|\| (arb_bits(x) + 1 < prec / 1024 && n >= 32))
	34	{
	35	arb_hypgeom_rising_ui_jet_bs(res, x, n, len, prec);
	36	}
	37	else
	38	{
	39	arb_hypgeom_rising_ui_jet_rs(res, x, n, 0, len, prec);
	40	}
	41	}
	42	}
	43

+84

-0

arb_hypgeom/rising_ui_jet_bs.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	static void
	14	bsplit(arb_ptr res, const arb_t x, ulong a, ulong b, slong trunc, slong prec)
	15	{
	16	trunc = FLINT_MIN(trunc, b - a + 1);
	17
	18	if (b - a <= 12)
	19	{
	20	if (a == 0)
	21	{
	22	arb_hypgeom_rising_ui_jet_powsum(res, x, b - a, FLINT_MIN(trunc, b - a + 1), prec);
	23	}
	24	else
	25	{
	26	arb_t t;
	27	arb_init(t);
	28	arb_add_ui(t, x, a, prec);
	29	arb_hypgeom_rising_ui_jet_powsum(res, t, b - a, FLINT_MIN(trunc, b - a + 1), prec);
	30	arb_clear(t);
	31	}
	32	}
	33	else
	34	{
	35	arb_ptr L, R;
	36	slong len1, len2;
	37
	38	slong m = a + (b - a) / 2;
	39
	40	len1 = poly_pow_length(2, m - a, trunc);
	41	len2 = poly_pow_length(2, b - m, trunc);
	42
	43	L = _arb_vec_init(len1 + len2);
	44	R = L + len1;
	45
	46	bsplit(L, x, a, m, trunc, prec);
	47	bsplit(R, x, m, b, trunc, prec);
	48
	49	_arb_poly_mullow(res, L, len1, R, len2,
	50	FLINT_MIN(trunc, len1 + len2 - 1), prec);
	51
	52	_arb_vec_clear(L, len1 + len2);
	53	}
	54	}
	55
	56	void
	57	arb_hypgeom_rising_ui_jet_bs(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
	58	{
	59	if (len == 0)
	60	return;
	61
	62	if (len > n + 1)
	63	{
	64	_arb_vec_zero(res + n + 1, len - n - 1);
	65	len = n + 1;
	66	}
	67
	68	if (len == n + 1)
	69	{
	70	arb_one(res + n);
	71	len = n;
	72	}
	73
	74	if (n <= 1)
	75	{
	76	if (n == 1)
	77	arb_set_round(res, x, prec);
	78	return;
	79	}
	80
	81	bsplit(res, x, 0, n, len, prec);
	82	}
	83

+159

-0

arb_hypgeom/rising_ui_jet_powsum.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	void
	14	arb_hypgeom_rising_ui_jet_powsum(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
	15	{
	16	slong i, j, k, wp;
	17	arb_ptr xpow;
	18	TMP_INIT;
	19
	20	if (len == 0)
	21	return;
	22
	23	if (len > n + 1)
	24	{
	25	_arb_vec_zero(res + n + 1, len - n - 1);
	26	len = n + 1;
	27	}
	28
	29	if (len == n + 1)
	30	{
	31	arb_one(res + n);
	32	len = n;
	33	}
	34
	35	if (n <= 1)
	36	{
	37	if (n == 1)
	38	arb_set_round(res, x, prec);
	39	return;
	40	}
	41
	42	if (len == 1)
	43	{
	44	arb_hypgeom_rising_ui_rs(res, x, n, 0, prec);
	45	return;
	46	}
	47
	48	if (n == 2)
	49	{
	50	arb_mul_2exp_si(res + 1, x, 1);
	51	arb_add_ui(res + 1, res + 1, 1, prec);
	52	arb_mul(res, x, x, prec + 4);
	53	arb_add(res, res, x, prec);
	54	return;
	55	}
	56
	57	if (n <= 12 \|\| (FLINT_BITS == 64 && n <= 20))
	58	{
	59	mp_ptr c;
	60	TMP_START;
	61
	62	wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
	63	c = TMP_ALLOC(sizeof(mp_limb_t) * (n + 1) * len);
	64
	65	_nmod_vec_zero(c, (n + 1) * len);
	66
	67	c[0] = 0;
	68	c[1] = 1;
	69	c[(n + 1) + 0] = 1;
	70
	71	for (i = 2; i <= n; i++)
	72	{
	73	for (j = FLINT_MIN(len - 1, i); j >= 0; j--)
	74	{
	75	slong ln, pos;
	76
	77	ln = i + 1 - j;
	78	pos = (n + 1) * j;
	79	if (i == j)
	80	{
	81	c[pos] = 1;
	82	}
	83	else
	84	{
	85	c[pos + ln - 1] = c[pos + ln - 2];
	86	for (k = ln - 2; k >= 1; k--)
	87	c[pos + k] = c[pos + k] * (i - 1) + c[pos + k - 1];
	88	c[pos + 0] *= (i - 1);
	89	if (j != 0)
	90	for (k = ln - 1; k >= 0; k--)
	91	c[pos + k] += c[pos - (n + 1) + k];
	92	}
	93	}
	94	}
	95
	96	xpow = _arb_vec_init(n + 1);
	97	_arb_vec_set_powers(xpow, x, n + 1, wp);
	98
	99	arb_dot_ui(res, NULL, 0, xpow + 1, 1, c + 1, 1, n, prec);
	100
	101	for (i = 1; i < len; i++)
	102	arb_dot_ui(res + i, NULL, 0, xpow, 1, c + (n + 1) * i, 1, n + 1 - i, prec);
	103
	104	_arb_vec_clear(xpow, n + 1);
	105	TMP_END;
	106	}
	107	else
	108	{
	109	fmpz * c;
	110
	111	wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
	112	c = _fmpz_vec_init((n + 1) * len);
	113
	114	fmpz_one(c + 1);
	115	fmpz_one(c + (n + 1) + 0);
	116
	117	for (i = 2; i <= n; i++)
	118	{
	119	for (j = FLINT_MIN(len - 1, i); j >= 0; j--)
	120	{
	121	slong ln, pos;
	122
	123	ln = i + 1 - j;
	124	pos = (n + 1) * j;
	125	if (i == j)
	126	{
	127	fmpz_one(c + pos);
	128	}
	129	else
	130	{
	131	fmpz_set(c + pos + ln - 1, c + pos + ln - 2);
	132	for (k = ln - 2; k >= 1; k--)
	133	{
	134	fmpz_mul_ui(c + pos + k, c + pos + k, i - 1);
	135	fmpz_add(c + pos + k, c + pos + k, c + pos + k - 1);
	136	}
	137
	138	fmpz_mul_ui(c + pos + 0, c + pos + 0, i - 1);
	139	if (j != 0)
	140	for (k = ln - 1; k >= 0; k--)
	141	fmpz_add(c + pos + k, c + pos + k, c + pos - (n + 1) + k);
	142	}
	143	}
	144	}
	145
	146	xpow = _arb_vec_init(n + 1);
	147	_arb_vec_set_powers(xpow, x, n + 1, wp);
	148
	149	arb_dot_fmpz(res, NULL, 0, xpow + 1, 1, c + 1, 1, n, prec);
	150
	151	for (i = 1; i < len; i++)
	152	arb_dot_fmpz(res + i, NULL, 0, xpow, 1, c + (n + 1) * i, 1, n + 1 - i, prec);
	153
	154	_arb_vec_clear(xpow, n + 1);
	155	_fmpz_vec_clear(c, (n + 1) * len);
	156	}
	157	}
	158

+173

-0

arb_hypgeom/rising_ui_jet_rs.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	void
	14	arb_hypgeom_rising_ui_jet_rs(arb_ptr res, const arb_t x, ulong n, ulong m, slong len, slong prec)
	15	{
	16	slong i, j, k, l, m0, xmlen, tlen, ulen, climbs, climbs_max, wp;
	17	arb_ptr tmp, xpow;
	18	arb_ptr t, u;
	19	mp_ptr c;
	20	TMP_INIT;
	21
	22	if (len == 0)
	23	return;
	24
	25	if (len > n + 1)
	26	{
	27	_arb_vec_zero(res + n + 1, len - n - 1);
	28	len = n + 1;
	29	}
	30
	31	if (len == n + 1)
	32	{
	33	arb_one(res + n);
	34	len = n;
	35	}
	36
	37	if (n <= 1)
	38	{
	39	if (n == 1)
	40	arb_set_round(res, x, prec);
	41	return;
	42	}
	43
	44	if (len == 1)
	45	{
	46	arb_hypgeom_rising_ui_rs(res, x, n, m, prec);
	47	return;
	48	}
	49
	50	TMP_START;
	51
	52	if (m == 0)
	53	{
	54	if (n <= 7)
	55	m = n;
	56	else if (n <= 12)
	57	m = (n + 1) / 2;
	58	else if (n <= 32)
	59	m = (n + 2) / 3;
	60	else
	61	{
	62	m0 = n_sqrt(n);
	63	m = 8 + 0.5 * pow(prec, 0.4);
	64	m = FLINT_MIN(m, m0);
	65	m = FLINT_MIN(m, 64);
	66	}
	67	}
	68
	69	wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
	70
	71	climbs_max = FLINT_BIT_COUNT(n - 1) * m;
	72	c = TMP_ALLOC(sizeof(mp_limb_t) * climbs_max * m);
	73
	74	/* length of (x+t)^m */
	75	xmlen = FLINT_MIN(len, m + 1);
	76
	77	tmp = _arb_vec_init(2 * len + (m + 1) * xmlen);
	78	t = tmp;
	79	u = tmp + len;
	80	xpow = tmp + 2 * len;
	81
	82	_arb_vec_set_powers(xpow, x, m + 1, wp);
	83
	84	tlen = 1;
	85
	86	/* First derivatives */
	87	for (i = 1; i <= m; i++)
	88	arb_mul_ui(xpow + (m + 1) + i, xpow + i - 1, i, wp);
	89
	90	/* Higher derivatives if we need them */
	91	if (len >= 3)
	92	{
	93	fmpz * f = _fmpz_vec_init(len);
	94
	95	fmpz_one(f + 0);
	96	fmpz_one(f + 1);
	97
	98	for (i = 2; i <= m; i++)
	99	{
	100	for (j = FLINT_MIN(xmlen - 1, i + 1); j >= 1; j--)
	101	fmpz_add(f + j, f + j, f + j - 1);
	102
	103	for (j = 2; j < FLINT_MIN(xmlen, i + 1); j++)
	104	arb_mul_fmpz(xpow + (m + 1) * j + i, xpow + i - j, f + j, wp);
	105	}
	106
	107	_fmpz_vec_clear(f, len);
	108	}
	109
	110	for (k = 0; k < n; k += m)
	111	{
	112	l = FLINT_MIN(m, n - k);
	113	climbs = FLINT_BIT_COUNT(k + l - 1) * l;
	114	climbs = (climbs + FLINT_BITS - 1) / FLINT_BITS;
	115
	116	ulen = FLINT_MIN(len, l + 1);
	117
	118	if (l == 1)
	119	{
	120	arb_add_ui(u, x, k, wp);
	121	arb_one(u + 1);
	122	}
	123	else
	124	{
	125	if (climbs == 1)
	126	{
	127	_arb_hypgeom_rising_coeffs_1(c, k, l);
	128	for (j = 0; j < ulen; j++)
	129	arb_dot_ui(u + j, xpow + (m + 1) * j + l, 0, xpow + (m + 1) * j + j, 1, c + j, 1, l - j, wp);
	130	}
	131	else if (climbs == 2)
	132	{
	133	_arb_hypgeom_rising_coeffs_2(c, k, l);
	134	for (j = 0; j < ulen; j++)
	135	arb_dot_uiui(u + j, xpow + (m + 1) * j + l, 0, xpow + (m + 1) * j + j, 1, c + 2 * j, 1, l - j, wp);
	136	}
	137	else
	138	{
	139	fmpz * f = (fmpz *) c;
	140
	141	for (i = 0; i < l; i++)
	142	fmpz_init(f + i);
	143
	144	_arb_hypgeom_rising_coeffs_fmpz(f, k, l);
	145
	146	for (j = 0; j < ulen; j++)
	147	arb_dot_fmpz(u + j, xpow + (m + 1) * j + l, 0, xpow + (m + 1) * j + j, 1, f + j, 1, l - j, wp);
	148
	149	for (i = 0; i < l; i++)
	150	fmpz_clear(f + i);
	151	}
	152	}
	153
	154	if (k == 0)
	155	{
	156	tlen = ulen;
	157	_arb_vec_swap(t, u, ulen);
	158	}
	159	else
	160	{
	161	_arb_poly_mullow(res, t, tlen, u, ulen, FLINT_MIN(len, tlen + ulen - 1), wp);
	162	tlen = FLINT_MIN(len, tlen + ulen - 1);
	163	_arb_vec_swap(t, res, tlen);
	164	}
	165	}
	166
	167	_arb_vec_set_round(res, t, len, prec);
	168
	169	_arb_vec_clear(tmp, 2 * len + (m + 1) * xmlen);
	170	TMP_END;
	171	}
	172

+46

-0

arb_hypgeom/rising_ui_rec.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	void
	14	arb_hypgeom_rising_ui_rec(arb_t res, const arb_t x, ulong n, slong prec)
	15	{
	16	if (n <= 1)
	17	{
	18	if (n == 0)
	19	arb_one(res);
	20	else
	21	arb_set_round(res, x, prec);
	22	return;
	23	}
	24
	25	if (n == 2 && prec <= 1024)
	26	{
	27	if (res != x)
	28	arb_set(res, x);
	29	arb_addmul(res, x, x, prec);
	30	return;
	31	}
	32
	33	if ((prec < 512 && n <= 20) \|\| (n <= FLINT_MIN(80, 6000 / prec)))
	34	{
	35	arb_hypgeom_rising_ui_forward(res, x, n, prec);
	36	}
	37	else
	38	{
	39	if (n >= 20 && arb_bits(x) < prec / 8)
	40	arb_hypgeom_rising_ui_bs(res, x, n, prec);
	41	else
	42	arb_hypgeom_rising_ui_rs(res, x, n, 0, prec);
	43	}
	44	}
	45

+198

-0

arb_hypgeom/rising_ui_rs.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	void
	14	_arb_hypgeom_rising_coeffs_1(ulong * c, ulong k, slong l)
	15	{
	16	slong i, j;
	17	ulong d;
	18
	19	if (l < 2) flint_abort();
	20
	21	c[0] = k * (k + 1);
	22	c[1] = 2 * k + 1;
	23
	24	for (i = 2; i < l; i++)
	25	{
	26	d = k + i;
	27
	28	c[i] = c[i - 1] + d;
	29
	30	for (j = i - 1; j >= 1; j--)
	31	c[j] = c[j] * d + c[j - 1];
	32
	33	c[0] = c[0] * d;
	34	}
	35	}
	36
	37	void
	38	_arb_hypgeom_rising_coeffs_2(ulong * c, ulong k, slong l)
	39	{
	40	slong i, j;
	41	ulong d;
	42	mp_limb_t hi, lo;
	43
	44	if (l < 2) flint_abort();
	45
	46	umul_ppmm(c[1], c[0], k, k + 1);
	47
	48	c[2] = 2 * k + 1;
	49	c[3] = 0;
	50
	51	for (i = 2; i < l; i++)
	52	{
	53	d = k + i;
	54
	55	add_ssaaaa(c[2 * i + 1], c[2 * i], c[2 * i - 1], c[2 * i - 2], 0, d);
	56
	57	for (j = i - 1; j >= 1; j--)
	58	{
	59	umul_ppmm(hi, lo, c[2 * j], d);
	60	hi += c[2 * j + 1] * d;
	61	add_ssaaaa(c[2 * j + 1], c[2 * j], hi, lo, c[2 * j - 1], c[2 * j - 2]);
	62	}
	63
	64	umul_ppmm(hi, lo, c[0], d);
	65	c[0] = lo;
	66	c[1] = c[1] * d + hi;
	67	}
	68	}
	69
	70	void
	71	_arb_hypgeom_rising_coeffs_fmpz(fmpz * c, ulong k, slong l)
	72	{
	73	slong i, j;
	74	ulong d;
	75
	76	if (l < 2) flint_abort();
	77
	78	fmpz_set_ui(c + 0, k);
	79	fmpz_mul_ui(c + 0, c + 0, k + 1);
	80	fmpz_set_ui(c + 1, 2 * k + 1);
	81
	82	for (i = 2; i < l; i++)
	83	{
	84	d = k + i;
	85
	86	fmpz_add_ui(c + i, c + i - 1, d);
	87
	88	for (j = i - 1; j >= 1; j--)
	89	{
	90	fmpz_mul_ui(c + j, c + j, d);
	91	fmpz_add(c + j, c + j, c + j - 1);
	92	}
	93
	94	fmpz_mul_ui(c + 0, c + 0, d);
	95	}
	96	}
	97
	98	void
	99	arb_hypgeom_rising_ui_rs(arb_t res, const arb_t x, ulong n, ulong m, slong prec)
	100	{
	101	slong i, k, l, m0, climbs, climbs_max, wp;
	102	arb_ptr xpow;
	103	arb_t t, u;
	104	mp_ptr c;
	105	TMP_INIT;
	106
	107	if (n <= 1)
	108	{
	109	if (n == 0)
	110	arb_one(res);
	111	else
	112	arb_set_round(res, x, prec);
	113	return;
	114	}
	115
	116	TMP_START;
	117
	118	if (m == 0)
	119	{
	120	if (n <= 6)
	121	m = 1 + (prec >= 1024);
	122	else if (n <= 16)
	123	m = 4;
	124	else if (n <= 50)
	125	m = 6;
	126	else
	127	{
	128	m0 = n_sqrt(n);
	129	m = 8 + 0.2 * pow(FLINT_MAX(0, prec - 4096), 0.4);
	130	m = FLINT_MIN(m, m0);
	131	m = FLINT_MIN(m, 60);
	132	}
	133	}
	134
	135	wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));
	136
	137	climbs_max = FLINT_BIT_COUNT(n - 1) * m;
	138	c = TMP_ALLOC(sizeof(mp_limb_t) * climbs_max * m);
	139
	140	xpow = _arb_vec_init(m + 1);
	141	_arb_vec_set_powers(xpow, x, m + 1, wp);
	142	arb_init(t);
	143	arb_init(u);
	144
	145	for (k = 0; k < n; k += m)
	146	{
	147	l = FLINT_MIN(m, n - k);
	148	climbs = FLINT_BIT_COUNT(k + l - 1) * l;
	149	climbs = (climbs + FLINT_BITS - 1) / FLINT_BITS;
	150
	151	/* assumes l >= 2 */
	152	if (l == 1)
	153	{
	154	arb_add_ui(u, x, k, wp);
	155	}
	156	else
	157	{
	158	if (climbs == 1)
	159	{
	160	_arb_hypgeom_rising_coeffs_1(c, k, l);
	161	arb_dot_ui(u, xpow + l, 0, xpow, 1, c, 1, l, wp);
	162	}
	163	else if (climbs == 2)
	164	{
	165	_arb_hypgeom_rising_coeffs_2(c, k, l);
	166	arb_dot_uiui(u, xpow + l, 0, xpow, 1, c, 1, l, wp);
	167	}
	168	else
	169	{
	170	fmpz * f = (fmpz *) c;
	171
	172	for (i = 0; i < l; i++)
	173	fmpz_init(f + i);
	174
	175	_arb_hypgeom_rising_coeffs_fmpz(f, k, l);
	176
	177	arb_dot_fmpz(u, xpow + l, 0, xpow, 1, f, 1, l, wp);
	178
	179	for (i = 0; i < l; i++)
	180	fmpz_clear(f + i);
	181	}
	182	}
	183
	184	if (k == 0)
	185	arb_swap(t, u);
	186	else
	187	arb_mul(t, t, u, wp);
	188	}
	189
	190	arb_set_round(res, t, prec);
	191
	192	arb_clear(t);
	193	arb_clear(u);
	194	_arb_vec_clear(xpow, m + 1);
	195	TMP_END;
	196	}
	197

+108

-0

arb_hypgeom/test/t-rising_ui.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	void
	14	rising_algorithm(arb_t res, const arb_t x, ulong n, ulong m, slong prec, int alg, int alias)
	15	{
	16	if (alias)
	17	{
	18	arb_set(res, x);
	19	rising_algorithm(res, res, n, m, prec, alg, 0);
	20	return;
	21	}
	22
	23	if (alg == 0)
	24	arb_hypgeom_rising_ui_rs(res, x, n, m, prec);
	25	else if (alg == 1)
	26	arb_hypgeom_rising_ui_forward(res, x, n, prec);
	27	else if (alg == 2)
	28	arb_hypgeom_rising_ui_bs(res, x, n, prec);
	29	else if (alg == 3)
	30	arb_hypgeom_rising_ui_rec(res, x, n, prec);
	31	else
	32	arb_hypgeom_rising_ui(res, x, n, prec);
	33	}
	34
	35	int main()
	36	{
	37	slong iter;
	38	flint_rand_t state;
	39
	40	flint_printf("rising_ui....");
	41	fflush(stdout);
	42
	43	flint_randinit(state);
	44
	45	for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
	46	{
	47	arb_t x, xk, y, ya, yb, yayb;
	48	ulong k, n, m1, m2, m3;
	49	slong prec;
	50	int alg1, alg2, alg3, alias1, alias2, alias3;
	51
	52	prec = 2 + n_randint(state, 200);
	53	k = n_randint(state, 10);
	54	n = n_randint(state, 50);
	55	m1 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n + k, 1));
	56	m2 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(k, 1));
	57	m3 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n, 1));
	58	alg1 = n_randint(state, 5);
	59	alg2 = n_randint(state, 5);
	60	alg3 = n_randint(state, 5);
	61	alias1 = n_randint(state, 2);
	62	alias2 = n_randint(state, 2);
	63	alias3 = n_randint(state, 2);
	64
	65	if (n_randint(state, 100) == 0)
	66	n += 100;
	67
	68	arb_init(x);
	69	arb_init(xk);
	70	arb_init(y);
	71	arb_init(ya);
	72	arb_init(yb);
	73	arb_init(yayb);
	74
	75	arb_randtest(x, state, prec, 10 + n_randint(state, 200));
	76	arb_add_ui(xk, x, k, prec);
	77
	78	rising_algorithm(y, x, n + k, m1, prec, alg1, alias1);
	79	rising_algorithm(ya, x, k, m2, prec, alg2, alias2);
	80	rising_algorithm(yb, xk, n, m3, prec, alg3, alias3);
	81	arb_mul(yayb, ya, yb, prec);
	82
	83	if (!arb_overlaps(y, yayb))
	84	{
	85	flint_printf("FAIL\n\n");
	86	flint_printf("k = %wu, n = %wu, m1 = %wu, m2 = %wu, m3 = %wu\n\n", k, n, m1, m2, m3);
	87	flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
	88	flint_printf("y = "); arb_printn(y, 100, 0); flint_printf("\n\n");
	89	flint_printf("ya = "); arb_printn(ya, 100, 0); flint_printf("\n\n");
	90	flint_printf("yb = "); arb_printn(yb, 100, 0); flint_printf("\n\n");
	91	flint_printf("yayb = "); arb_printn(yayb, 100, 0); flint_printf("\n\n");
	92	flint_abort();
	93	}
	94
	95	arb_clear(x);
	96	arb_clear(xk);
	97	arb_clear(y);
	98	arb_clear(ya);
	99	arb_clear(yb);
	100	arb_clear(yayb);
	101	}
	102
	103	flint_randclear(state);
	104	flint_cleanup();
	105	flint_printf("PASS\n");
	106	return EXIT_SUCCESS;
	107	}

+95

-0

arb_hypgeom/test/t-rising_ui_jet.c less more

	0	/*
	1	Copyright (C) 2021 Fredrik Johansson
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arb_hypgeom.h"
	12
	13	void
	14	rising_algorithm(arb_ptr res, const arb_t x, ulong n, ulong m, slong len, slong prec, int alg)
	15	{
	16	if (alg == 0)
	17	arb_hypgeom_rising_ui_jet_powsum(res, x, n, len, prec);
	18	else if (alg == 1)
	19	arb_hypgeom_rising_ui_jet_rs(res, x, n, m, len, prec);
	20	else if (alg == 2)
	21	arb_hypgeom_rising_ui_jet_bs(res, x, n, len, prec);
	22	else
	23	arb_hypgeom_rising_ui_jet(res, x, n, len, prec);
	24	}
	25
	26	int main()
	27	{
	28	slong iter;
	29	flint_rand_t state;
	30
	31	flint_printf("rising_ui_jet....");
	32	fflush(stdout);
	33
	34	flint_randinit(state);
	35
	36	for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
	37	{
	38	arb_t x, xk;
	39	arb_ptr y, ya, yb, yayb;
	40	ulong k, n, m1, m2, m3, len;
	41	slong prec;
	42	int alg1, alg2, alg3;
	43
	44	prec = 2 + n_randint(state, 200);
	45	len = 1 + n_randint(state, 6);
	46	k = n_randint(state, 10);
	47	n = n_randint(state, 50);
	48	m1 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n + k, 1));
	49	m2 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(k, 1));
	50	m3 = n_randint(state, 2) ? 0 : 1 + n_randint(state, FLINT_MAX(n, 1));
	51	alg1 = n_randint(state, 4);
	52	alg2 = n_randint(state, 4);
	53	alg3 = n_randint(state, 4);
	54
	55	arb_init(x);
	56	arb_init(xk);
	57	y = _arb_vec_init(len);
	58	ya = _arb_vec_init(len);
	59	yb = _arb_vec_init(len);
	60	yayb = _arb_vec_init(len);
	61
	62	arb_randtest(x, state, prec, 10);
	63	arb_add_ui(xk, x, k, prec);
	64
	65	rising_algorithm(y, x, n + k, m1, len, prec, alg1);
	66	rising_algorithm(ya, x, k, m2, len, prec, alg2);
	67	rising_algorithm(yb, xk, n, m3, len, prec, alg3);
	68	_arb_poly_mullow(yayb, ya, len, yb, len, len, prec);
	69
	70	if (!_arb_poly_overlaps(y, len, yayb, len))
	71	{
	72	flint_printf("FAIL\n\n");
	73	flint_printf("len = %wd, k = %wu, n = %wu, m1 = %wu, m2 = %wu, m3 = %wu\n\n", len, k, n, m1, m2, m3);
	74	flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
	75	flint_printf("y = "); _arb_vec_printn(y, len, 100, 0); flint_printf("\n\n");
	76	flint_printf("ya = "); _arb_vec_printn(ya, len, 100, 0); flint_printf("\n\n");
	77	flint_printf("yb = "); _arb_vec_printn(yb, len, 100, 0); flint_printf("\n\n");
	78	flint_printf("yayb = "); _arb_vec_printn(yayb, len, 100, 0); flint_printf("\n\n");
	79	flint_abort();
	80	}
	81
	82	arb_clear(x);
	83	arb_clear(xk);
	84	_arb_vec_clear(y, len);
	85	_arb_vec_clear(ya, len);
	86	_arb_vec_clear(yb, len);
	87	_arb_vec_clear(yayb, len);
	88	}
	89
	90	flint_randclear(state);
	91	flint_cleanup();
	92	flint_printf("PASS\n");
	93	return EXIT_SUCCESS;
	94	}

+16

-0

arb_hypgeom.h less more

17	17	#ifdef __cplusplus
18	18	extern "C" {
19	19	#endif
	20
	21	void _arb_hypgeom_rising_coeffs_1(ulong * c, ulong k, slong l);
	22	void _arb_hypgeom_rising_coeffs_2(ulong * c, ulong k, slong l);
	23	void _arb_hypgeom_rising_coeffs_fmpz(fmpz * c, ulong k, slong l);
	24
	25	void arb_hypgeom_rising_ui_forward(arb_t res, const arb_t x, ulong n, slong prec);
	26	void arb_hypgeom_rising_ui_rs(arb_t res, const arb_t x, ulong n, ulong m, slong prec);
	27	void arb_hypgeom_rising_ui_bs(arb_t res, const arb_t x, ulong n, slong prec);
	28	void arb_hypgeom_rising_ui_rec(arb_t res, const arb_t x, ulong n, slong prec);
	29	void arb_hypgeom_rising_ui(arb_t y, const arb_t x, ulong n, slong prec);
	30	void arb_hypgeom_rising(arb_t y, const arb_t x, const arb_t n, slong prec);
	31
	32	void arb_hypgeom_rising_ui_jet_powsum(arb_ptr res, const arb_t x, ulong n, slong len, slong prec);
	33	void arb_hypgeom_rising_ui_jet_rs(arb_ptr res, const arb_t x, ulong n, ulong m, slong len, slong prec);
	34	void arb_hypgeom_rising_ui_jet_bs(arb_ptr res, const arb_t x, ulong n, slong len, slong prec);
	35	void arb_hypgeom_rising_ui_jet(arb_ptr res, const arb_t x, ulong n, slong len, slong prec);
20	36
21	37	void arb_hypgeom_pfq(arb_t res, arb_srcptr a, slong p, arb_srcptr b, slong q,
22	38	const arb_t z, int regularized, slong prec);

+1

-1

arb_mat/approx_mul.c less more

30	30	arb_mat_t T;
31	31	arb_mat_init(T, ar, bc);
32	32	arb_mat_approx_mul_classical(T, A, B, prec);
33		arb_mat_swap(T, C);
	33	arb_mat_swap_entrywise(T, C);
34	34	arb_mat_clear(T);
35	35	return;
36	36	}

+1

-1

arb_mat/mul_block.c less more

237	237	arb_mat_t T;
238	238	arb_mat_init(T, M, P);
239	239	arb_mat_mul_block(T, A, B, prec);
240		arb_mat_swap(T, C);
	240	arb_mat_swap_entrywise(T, C);
241	241	arb_mat_clear(T);
242	242	return;
243	243	}

+1

-1

arb_mat/mul_classical.c less more

44	44	arb_mat_t T;
45	45	arb_mat_init(T, ar, bc);
46	46	arb_mat_mul_classical(T, A, B, prec);
47		arb_mat_swap(T, C);
	47	arb_mat_swap_entrywise(T, C);
48	48	arb_mat_clear(T);
49	49	return;
50	50	}

+1

-1

arb_mat/mul_threaded.c less more

86	86	arb_mat_t T;
87	87	arb_mat_init(T, ar, bc);
88	88	arb_mat_mul_threaded(T, A, B, prec);
89		arb_mat_swap(T, C);
	89	arb_mat_swap_entrywise(T, C);
90	90	arb_mat_clear(T);
91	91	return;
92	92	}

+10

-0

arb_mat.h less more

61	61	arb_mat_struct t = *mat1;
62	62	mat1 = mat2;
63	63	*mat2 = t;
	64	}
	65
	66	ARB_MAT_INLINE void
	67	arb_mat_swap_entrywise(arb_mat_t mat1, arb_mat_t mat2)
	68	{
	69	slong i, j;
	70
	71	for (i = 0; i < arb_mat_nrows(mat1); i++)
	72	for (j = 0; j < arb_mat_ncols(mat1); j++)
	73	arb_swap(arb_mat_entry(mat2, i, j), arb_mat_entry(mat1, i, j));
64	74	}
65	75
66	76	/* Window matrices */

+36

-0

arf/urandom.c less more

	0	/*
	1	Copyright (C) 2021 Albin Ahlbäck
	2
	3	This file is part of Arb.
	4
	5	Arb is free software: you can redistribute it and/or modify it under
	6	the terms of the GNU Lesser General Public License (LGPL) as published
	7	by the Free Software Foundation; either version 2.1 of the License, or
	8	(at your option) any later version. See <http://www.gnu.org/licenses/>.
	9	*/
	10
	11	#include "arf.h"
	12
	13	void
	14	arf_urandom(arf_t x, flint_rand_t state, slong bits, arf_rnd_t rnd)
	15	{
	16	slong prec = bits;
	17	fmpz_t n;
	18	fmpz_t t;
	19
	20	prec += 128;
	21
	22	fmpz_init(n);
	23	fmpz_one(n);
	24	fmpz_mul_2exp(n, n, (ulong) prec);
	25
	26	fmpz_init(t);
	27	fmpz_randm(t, state, n);
	28
	29	arf_set_round_fmpz(x, t, bits, rnd);
	30	arf_mul_2exp_si(x, x, -prec);
	31
	32	fmpz_clear(n);
	33	fmpz_clear(t);
	34	}
	35

+3

-1

arf.h less more

104	104	#define ARF_XSIZE(x) ((x)->size)
105	105
106	106	/* Construct size field value from size in limbs and sign bit. */
107		#define ARF_MAKE_XSIZE(size, sgnbit) ((((mp_size_t) size) << 1) \| sgnbit)
	107	#define ARF_MAKE_XSIZE(size, sgnbit) ((((mp_size_t) size) << 1) \| (sgnbit))
108	108
109	109	/* The limb size, and the sign bit. */
110	110	#define ARF_SIZE(x) (ARF_XSIZE(x) >> 1)

795	795	void arf_randtest_not_zero(arf_t x, flint_rand_t state, slong bits, slong mag_bits);
796	796
797	797	void arf_randtest_special(arf_t x, flint_rand_t state, slong bits, slong mag_bits);
	798
	799	void arf_urandom(arf_t x, flint_rand_t state, slong bits, arf_rnd_t rnd);
798	800
799	801	#define MUL_MPFR_MIN_LIMBS 25
800	802	#define MUL_MPFR_MAX_LIMBS 10000

+1

-1

configure less more

9	9	# arb => soname
10	10	# 2.7.0 => 0.0.0
11	11	ARB_MAJOR=2
12		ARB_MINOR=10
	12	ARB_MINOR=11
13	13	ARB_PATCH=0
14	14
15	15	PREFIX="/usr/local"

+14

-5

doc/source/acb.rst less more

521	521	-------------------------------------------------------------------------------
522	522
523	523	.. function:: void acb_dot_precise(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)
524
525		.. function:: void acb_dot_simple(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)
526
527		.. function:: void acb_dot(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)
	524	void acb_dot_simple(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)
	525	void acb_dot(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec)
528	526
529	527	Computes the dot product of the vectors x and y, setting
530	528	res to `s + (-1)^{subtract} \sum_{i=0}^{len-1} x_i y_i`.

563	561	The radii of the inputs are ignored (only the midpoints are read)
564	562	and only the midpoint of the output is written.
565	563
	564	.. function:: void acb_dot_ui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	565	void acb_dot_si(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec)
	566	void acb_dot_uiui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	567	void acb_dot_siui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	568	void acb_dot_fmpz(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec)
	569
	570	Equivalent to :func:`acb_dot`, but with integers in the array y.
	571	The uiui and siui versions take an array of double-limb integers
	572	as input; the siui version assumes that these represent signed
	573	integers in two's complement form.
	574
566	575	Mathematical constants
567	576	-------------------------------------------------------------------------------
568	577

773	782
774	783	Computes `\operatorname{sech}(z) = 1 / \cosh(z)`.
775	784
776		.. function:: void acb_csch(acb_t res, const arb_t z, slong prec)
	785	.. function:: void acb_csch(acb_t res, const acb_t z, slong prec)
777	786
778	787	Computes `\operatorname{csch}(z) = 1 / \sinh(z)`.
779	788

+4

-0

doc/source/acb_elliptic.rst less more

298	298	\frac{\theta_4^2(z,\tau)}{\theta_1^2(z,\tau)} -
299	299	\frac{\pi^2}{3} \left[ \theta_2^4(0,\tau) + \theta_3^4(0,\tau)\right].
300	300
	301	.. function:: void acb_elliptic_p_prime(acb_t res, const acb_t z, const acb_t tau, slong prec)
	302
	303	Computes the derivative `\wp'(z, \tau)` of Weierstrass's elliptic function `\wp(z, \tau)`.
	304
301	305	.. function:: void acb_elliptic_p_jet(acb_ptr res, const acb_t z, const acb_t tau, slong len, slong prec)
302	306
303	307	Computes the formal power series `\wp(z + x, \tau) \in \mathbb{C}[[x]]`,

+8

-0

doc/source/acb_hypgeom.rst less more

14	14	In a looser sense, we understand "hypergeometric functions" to be
15	15	linear combinations of generalized hypergeometric functions
16	16	with prefactors that are products of exponentials, powers, and gamma functions.
	17
	18	Rising factorials
	19	-------------------------------------------------------------------------------
	20
	21	.. function:: void acb_hypgeom_rising_ui_rs(acb_t res, const acb_t x, ulong n, ulong m, slong prec)
	22
	23	Computes the rising factorial `(x)_n` using rectangular splitting.
	24	The splitting parameter m can be set to zero to choose automatically.
17	25
18	26	Convergent series
19	27	-------------------------------------------------------------------------------

+32

-9

doc/source/arb.rst less more

155	155
156	156	.. function:: void arb_set_round_fmpz(arb_t y, const fmpz_t x, slong prec)
157	157
158		Sets y to the value of x, rounded to prec bits.
	158	Sets y to the value of x, rounded to prec bits in the direction
	159	towards zero.
159	160
160	161	.. function:: void arb_set_round_fmpz_2exp(arb_t y, const fmpz_t x, const fmpz_t e, slong prec)
161	162
162		Sets y to `x \cdot 2^e`, rounded to prec bits.
	163	Sets y to `x \cdot 2^e`, rounded to prec bits in the direction
	164	towards zero.
163	165
164	166	.. function:: void arb_set_fmpq(arb_t y, const fmpq_t x, slong prec)
165	167
166		Sets y to the rational number x, rounded to prec bits.
	168	Sets y to the rational number x, rounded to prec bits in the direction
	169	towards zero.
167	170
168	171	.. function:: int arb_set_str(arb_t res, const char * inp, slong prec)
169	172

198	201	digits may be printed.
199	202
200	203	If ARB_STR_NO_RADIUS is added to flags, the radius is not
201		included in the output if at least 1 digit of the midpoint
202		can be printed.
	204	included in the output. Unless ARB_STR_MORE is set, the output is
	205	rounded so that the midpoint is correct to 1 ulp. As a special case,
	206	if there are no significant digits after rounding, the result will
	207	be shown as ``0e+n``, meaning that the result is between
	208	``-1e+n`` and ``1e+n`` (following the contract that the output is
	209	correct to within one unit in the only shown digit).
203	210
204	211	By adding a multiple m of ARB_STR_CONDENSE to flags, strings
205	212	of more than three times m consecutive digits are condensed, only

362	369	that representing the endpoints as exact rational numbers would
363	370	cause overflows.
364	371
	372	.. function:: void arb_urandom(arb_t x, flint_rand_t state, slong prec, arf_rnd_t rnd)
	373
	374	Sets x to a uniformly distributed random number in the interval
	375	`[0, 1]`. The method uses rounding from integers to floats, hence the
	376	radius might not be `0`.
	377
365	378	Radius and interval operations
366	379	-------------------------------------------------------------------------------
367	380

852	865	-------------------------------------------------------------------------------
853	866
854	867	.. function:: void arb_dot_precise(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec)
855
856		.. function:: void arb_dot_simple(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec)
857
858		.. function:: void arb_dot(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec)
	868	void arb_dot_simple(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec)
	869	void arb_dot(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec)
859	870
860	871	Computes the dot product of the vectors x and y, setting
861	872	res to `s + (-1)^{subtract} \sum_{i=0}^{len-1} x_i y_i`.

893	904	Computes an approximate dot product without error bounds.
894	905	The radii of the inputs are ignored (only the midpoints are read)
895	906	and only the midpoint of the output is written.
	907
	908	.. function:: void arb_dot_ui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	909	void arb_dot_si(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec)
	910	void arb_dot_uiui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	911	void arb_dot_siui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec)
	912	void arb_dot_fmpz(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec)
	913
	914	Equivalent to :func:`arb_dot`, but with integers in the array y.
	915	The uiui and siui versions take an array of double-limb integers
	916	as input; the siui version assumes that these represent signed
	917	integers in two's complement form.
	918
896	919
897	920	Powers and roots
898	921	-------------------------------------------------------------------------------

+6

-2

doc/source/arb_fmpz_poly.rst less more

65	65	.. function:: void arb_fmpz_poly_complex_roots(acb_ptr roots, const fmpz_poly_t poly, int flags, slong prec)
66	66
67	67	Writes to roots all the real and complex roots of the polynomial poly,
68		computed to prec accurate bits.
	68	computed to at least prec accurate bits.
	69	The root enclosures are guaranteed to be disjoint, so that
	70	all roots are isolated.
	71
69	72	The real roots are written first in ascending order (with
70	73	the imaginary parts set exactly to zero). The following
71	74	nonreal roots are written in arbitrary order, but with conjugate pairs

93	96
94	97	All roots are refined to a relative accuracy of at least prec bits.
95	98	The output values will generally have higher actual precision,
96		depending on the precision used internally by the algorithm.
	99	depending on the precision needed for isolation and the
	100	precision used internally by the algorithm.
97	101
98	102	This implementation should be adequate for general use, but it is not
99	103	currently competitive with state-of-the-art isolation

+59

-7

doc/source/arb_hypgeom.rst less more

14	14
15	15	This module also provides certain functions exclusive to real variables,
16	16	such as functions for computing real roots of common special functions.
	17
	18	Rising factorials
	19	-------------------------------------------------------------------------------
	20
	21	.. function:: void _arb_hypgeom_rising_coeffs_1(ulong * c, ulong k, slong n)
	22	void _arb_hypgeom_rising_coeffs_2(ulong * c, ulong k, slong n)
	23	void _arb_hypgeom_rising_coeffs_fmpz(fmpz * c, ulong k, slong n)
	24
	25	Sets c to the coefficients of the rising factorial polynomial
	26	`(X+k)_n`. The 1 and 2 versions respectively
	27	compute single-word and double-word coefficients, without checking for
	28	overflow, while the fmpz version allows arbitrarily large coefficients.
	29	These functions are mostly intended for internal use; the fmpz version
	30	does not use an asymptotically fast algorithm.
	31	The degree n must be at least 2.
	32
	33	.. function:: void arb_hypgeom_rising_ui_forward(arb_t res, const arb_t x, ulong n, slong prec)
	34	void arb_hypgeom_rising_ui_bs(arb_t res, const arb_t x, ulong n, slong prec)
	35	void arb_hypgeom_rising_ui_rs(arb_t res, const arb_t x, ulong n, ulong m, slong prec)
	36	void arb_hypgeom_rising_ui_rec(arb_t res, const arb_t x, ulong n, slong prec)
	37	void arb_hypgeom_rising_ui(arb_t y, const arb_t x, ulong n, slong prec)
	38	void arb_hypgeom_rising(arb_t y, const arb_t x, const arb_t n, slong prec)
	39
	40	Computes the rising factorial `(x)_n`.
	41
	42	The forward version uses the forward recurrence.
	43	The bs version uses binary splitting.
	44	The rs version uses rectangular splitting. It takes an extra tuning
	45	parameter m which can be set to zero to choose automatically.
	46	The rec version chooses an algorithm automatically, avoiding
	47	use of the gamma function (so that it can be used in the computation
	48	of the gamma function).
	49	The default versions (rising_ui and rising_ui) choose an algorithm
	50	automatically and may additionally fall back on the gamma function.
	51
	52	.. function:: void arb_hypgeom_rising_ui_jet_powsum(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
	53	void arb_hypgeom_rising_ui_jet_bs(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
	54	void arb_hypgeom_rising_ui_jet_rs(arb_ptr res, const arb_t x, ulong n, ulong m, slong len, slong prec)
	55	void arb_hypgeom_rising_ui_jet(arb_ptr res, const arb_t x, ulong n, slong len, slong prec)
	56
	57	Computes the jet of the rising factorial `(x)_n`, truncated to length len.
	58	In other words, constructs the polynomial `(X + x)_n \in \mathbb{R}[X]`,
	59	truncated if `\operatorname{len} < n + 1` (and zero-extended
	60	if `\operatorname{len} > n + 1`).
	61
	62	The powsum version computes the sequence of powers of x and forms integral
	63	linear combinations of these.
	64	The bs version uses binary splitting.
	65	The rs version uses rectangular splitting. It takes an extra tuning
	66	parameter m which can be set to zero to choose automatically.
	67	The default version chooses an algorithm automatically.
	68
	69
	70	Binomial coefficients
	71	-------------------------------------------------------------------------------
	72
	73	.. function:: void arb_hypgeom_central_bin_ui(arb_t res, ulong n, slong prec)
	74
	75	Computes the central binomial coefficient `{2n \choose n}`.
17	76
18	77	Generalized hypergeometric function
19	78	-------------------------------------------------------------------------------

448	507
449	508	Computes the dilogarithm `\operatorname{Li}_2(z)`.
450	509
451		Hypergeometric sequences
452		-------------------------------------------------------------------------------
453
454		.. function:: void arb_hypgeom_central_bin_ui(arb_t res, ulong n, slong prec)
455
456		Computes the central binomial coefficient `{2n \choose n}`.
457

+6

-0

doc/source/arf.rst less more

481	481	Identical to :func:`arf_randtest`, except that the output occasionally
482	482	is set to an infinity or NaN.
483	483
	484	.. function:: void arf_urandom(arf_t res, flint_rand_t state, slong bits, arf_rnd_t rnd)
	485
	486	Sets res to a uniformly distributed random number in the interval
	487	`[0, 1]`. The method uses rounding from integers to floats based on the
	488	rounding mode rnd.
	489
484	490	Input and output
485	491	-------------------------------------------------------------------------------
486	492

+4

-0

doc/source/credits.rst less more

69	69	* Joel Dahne - feedback and improvements for Legendre functions
70	70	* Gianfranco Costamagna - bug reports, Debian testing
71	71	* Julian Rüth - serialization support
	72	* Michael Orlitzky - build system patches
	73	* David Berghaus - aliased window matrix multiplication
	74	* Albin Ahlbäck - uniformly distributed random numbers
	75	* Daniel Schultz - derivative of Weierstrass elliptic function
72	76
73	77	Funding
74	78	-------------------------------------------------------------------------------

+23

-0

doc/source/history.rst less more

11	11	Old versions of the documentation
12	12	-------------------------------------------------------------------------------
13	13
	14	* http://arblib.org/arb-2.20.0.pdf
14	15	* http://arblib.org/arb-2.19.0.pdf
15	16	* http://arblib.org/arb-2.18.1.pdf
16	17	* http://arblib.org/arb-2.18.0.pdf

31	32	* http://arblib.org/arb-2.5.0.pdf
32	33	* http://arblib.org/arb-2.4.0.pdf
33	34	* http://arblib.org/arb-2.3.0.pdf
	35
	36	2021-07-25 -- version 2.20.0
	37	-------------------------------------------------------------------------------
	38
	39	* Flint 2.18 support.
	40	* Change arb_get_str with ARB_STR_NO_RADIUS: [+/- 1.20e-15] now prints as 0e-14.
	41	* Uniformly distributed random number functions arf_urandom, arb_urandom
	42	(contributed by Albin Ahlbäck).
	43	* Use quasilinear algorithm in arb_gamma_fmpq for all small fractions.
	44	* Added derivative of Weierstrass elliptic function (acb_elliptic_p_prime)
	45	(contributed by Daniel Schultz).
	46	* Added dot products with integer coefficients: arb_dot_fmpz, arb_dot_siui,
	47	arb_dot_uiui, arb_dot_si, arb_dot_ui, acb_dot_fmpz, acb_dot_siui,
	48	acb_dot_uiui, acb_dot_si, acb_dot_ui.
	49	* Faster arb_fmpz_poly_evaluate_arb and arb_fmpz_poly_evaluate_acb.
	50	* Explicitly guarantee that roots are isolated in arb_fmpz_poly_complex_roots
	51	(could previously theoretically fail when using the deflation hack).
	52	* Use GNUInstallDirs in CMakeLists.txt to support standard GNU installation
	53	directories (contributed by Michael Orlitzky).
	54	* Fixed bug for aliased multiplication of window matrices (contributed by
	55	David Berghaus).
	56	* Documentation fixes (contributed by Joel Dahne, Hanno Rein).
34	57
35	58	2020-12-06 -- version 2.19.0
36	59	-------------------------------------------------------------------------------