/* Copyright (C) 2016 Fredrik Johansson This file is part of Arb. Arb is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. See . */ #include "acb_dirichlet.h" /* Laurent expansions at s = 1 of first 10 principal L-functions */ /* with mpmath: chis = [[1],[0,1],[0,1,1],[0,1,0,1],[0,1,1,1,1],[0,1,0,0,0,1],[0,1,1,1,1,1,1], [0,1,0,1,0,1,0,1],[0,1,1,0,1,1,0,1,1],[0,1,0,1,0,0,0,1,0,1]] mp.dps = 40 for chi in chis: phi = chi.count(1); q = len(chi) L = lambda s: dirichlet(s, chi) - phi/((s-1)*q) c0 = taylor(L, 1, 0, method="quad") c1 = taylor(L, 1, 5, singular=True)[1:] for c in c0 + c1: print nstr(c, 20) + ",", print */ #define TESTQ 10 #define TESTLEN 6 static const double laurent_data[TESTQ][TESTLEN] = { {0.57721566490153286061, 0.072815845483676724861, -0.0048451815964361592423, -0.00034230573671722431103, 0.000096890419394470835728, -6.6110318108421891813e-6}, {0.63518142273073908501, 0.11634237461305384831, -0.018765738937942729408, 0.00061334298434914532242, 0.00042338142025747308027, -0.00010545096447379519004}, {0.7510145394903918042, 0.058764477744540050414, -0.019011359100973296683, 0.0056382252365739175151, -0.0009550480622176659462, 0.000021808301216554848718}, {0.63518142273073908501, 0.11634237461305384831, -0.018765738937942729408, 0.00061334298434914532242, 0.00042338142025747308027, -0.00010545096447379519004}, {0.78366011440804636341, -0.014977808062405260803, 0.0090104707969118845102, 0.003603799084856807634, -0.0029351216034181476022, 0.00093077685173004747355}, {0.60655632993184433857, 0.2095885418562151802, -0.060844893711330538429, 0.0068080382961291386117, 0.0022236616427578346453, -0.0013581825996235430782}, {0.77274344835207292411, -0.047596894381510269689, 0.035406039531261788462, -0.0054159870134630085898, -0.0019749752308692423114, 0.0014492998471928196325}, {0.63518142273073908501, 0.11634237461305384831, -0.018765738937942729408, 0.00061334298434914532242, 0.00042338142025747308027, -0.00010545096447379519004}, {0.7510145394903918042, 0.058764477744540050414, -0.019011359100973296683, 0.0056382252365739175151, -0.0009550480622176659462, 0.000021808301216554848718}, {0.66908892942800130547, 0.16801639259476784034, -0.072611999814034642781, 0.024624650443138705595, -0.004951850872731033514, -0.00020178815459414925709} }; int main() { slong iter; flint_rand_t state; flint_printf("l_jet...."); fflush(stdout); flint_randinit(state); /* test Laurent series at s = 1 */ { acb_t s, t; dirichlet_group_t G; dirichlet_char_t chi; acb_ptr vec; ulong q; slong i; acb_init(s); acb_init(t); vec = _acb_vec_init(TESTLEN); acb_one(s); for (q = 1; q <= TESTQ; q++) { dirichlet_group_init(G, q); dirichlet_char_init(chi, G); acb_dirichlet_l_jet(vec, s, G, chi, 1, TESTLEN, 100); for (i = 0; i < TESTLEN; i++) { acb_set_d(t, laurent_data[q - 1][i]); mag_set_d(arb_radref(acb_realref(t)), fabs(laurent_data[q - 1][i]) * 1e-14); if (!acb_overlaps(vec + i, t)) { flint_printf("FAIL: Laurent series\n\n"); flint_printf("q = %wu i = %wd\n\n", q, i); flint_printf("r1 = "); acb_printn(vec + i, 50, 0); flint_printf("\n\n"); flint_printf("r2 = "); acb_printn(t, 50, 0); flint_printf("\n\n"); flint_abort(); } } dirichlet_char_clear(chi); dirichlet_group_clear(G); } acb_clear(s); acb_clear(t); _acb_vec_clear(vec, TESTLEN); } /* test self-consistency */ for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) { acb_t s; dirichlet_group_t G; dirichlet_char_t chi; acb_ptr vec1, vec2; slong len1, len2; slong prec1, prec2; int deflate1, deflate2; ulong q, k; slong i; len1 = n_randint(state, 5); len2 = n_randint(state, 5); prec1 = 2 + n_randint(state, 100); prec2 = 2 + n_randint(state, 100); deflate1 = n_randint(state, 2); deflate2 = n_randint(state, 2); q = 1 + n_randint(state, 20); k = n_randint(state, n_euler_phi(q)); dirichlet_group_init(G, q); dirichlet_char_init(chi, G); dirichlet_char_index(chi, G, k); acb_init(s); vec1 = _acb_vec_init(len1); vec2 = _acb_vec_init(len2); if (n_randint(state, 4) == 0) acb_one(s); else acb_randtest(s, state, 2 + n_randint(state, 200), 2); acb_dirichlet_l_jet(vec1, s, G, chi, deflate1, len1, prec1); acb_dirichlet_l_jet(vec2, s, G, chi, deflate2, len2, prec2); if (deflate1 != deflate2 && dirichlet_char_is_principal(G, chi)) { /* add or subtract phi(q)/((s+x-1)q) */ acb_t t, u; acb_init(t); acb_init(u); acb_set_ui(t, n_euler_phi(q)); acb_div_ui(t, t, q, prec1); acb_sub_ui(u, s, 1, prec1); for (i = 0; i < len1; i++) { acb_div(t, t, u, prec1); if (deflate1) acb_add(vec1 + i, vec1 + i, t, prec1); else acb_sub(vec1 + i, vec1 + i, t, prec1); acb_neg(t, t); } acb_clear(t); acb_clear(u); } for (i = 0; i < FLINT_MIN(len1, len2); i++) { if (!acb_overlaps(vec1 + i, vec2 + i)) { flint_printf("FAIL: overlap\n\n"); flint_printf("iter = %wd q = %wu k = %wu i = %wd\n\n", iter, q, k, i); flint_printf("s = "); acb_printn(s, 50, 0); flint_printf("\n\n"); flint_printf("r1 = "); acb_printn(vec1 + i, 50, 0); flint_printf("\n\n"); flint_printf("r2 = "); acb_printn(vec2 + i, 50, 0); flint_printf("\n\n"); flint_abort(); } } dirichlet_char_clear(chi); dirichlet_group_clear(G); acb_clear(s); _acb_vec_clear(vec1, len1); _acb_vec_clear(vec2, len2); } flint_randclear(state); flint_cleanup(); flint_printf("PASS\n"); return EXIT_SUCCESS; }