mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
add incgam algorithm for L functions, disabled
This commit is contained in:
Pascal
parent
317707255e
commit
ceae160de3
3 changed files with 184 additions and 41 deletions
|
|
@ -271,6 +271,7 @@ void acb_dirichlet_jacobi_sum(acb_t res, const acb_dirichlet_group_t G, const ac
|
|||
|
||||
void acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
|
||||
void acb_dirichlet_l_vec_hurwitz(acb_ptr res, const acb_t s, const acb_dirichlet_group_t G, slong prec);
|
||||
void acb_dirichlet_l_incgam(acb_t res, const acb_t s, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
|
||||
|
||||
/* Discrete Fourier Transform */
|
||||
|
||||
|
|
|
|||
111
acb_dirichlet/l_incgam.c
Normal file
111
acb_dirichlet/l_incgam.c
Normal file
|
|
@ -0,0 +1,111 @@
|
|||
/*
|
||||
Copyright (C) 2016 Pascal Molin
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_dirichlet.h"
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
ulong
|
||||
_acb_dirichlet_l_incgam_length(const acb_t s, ulong q, slong prec)
|
||||
{
|
||||
#if 0
|
||||
arf_t a;
|
||||
arf_init(a);
|
||||
arb_get_lbound_arf(a, s, 53);
|
||||
a = arf_get_d(at, ARF_RND_DOWN);
|
||||
#else
|
||||
/* as a first approximation for Re(s) <= 1 */
|
||||
return acb_dirichlet_theta_length_d(q, 1., prec);
|
||||
#endif
|
||||
}
|
||||
|
||||
void
|
||||
acb_dirichlet_l_incgam(acb_t res, const acb_t s,
|
||||
const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
|
||||
{
|
||||
slong k, len;
|
||||
acb_ptr a;
|
||||
acb_t s1, s2, args, arg1s, eps;
|
||||
acb_t g1, g2, pq, pqk2;
|
||||
|
||||
len = _acb_dirichlet_l_incgam_length(s, chi->q, prec) + 2;
|
||||
flint_printf("incgam length = %ld\n",len);
|
||||
|
||||
prec += n_clog(len, 2);
|
||||
|
||||
a = _acb_vec_init(len);
|
||||
acb_dirichlet_chi_vec(a, G, chi, len, prec);
|
||||
|
||||
acb_init(args);
|
||||
acb_init(arg1s);
|
||||
acb_init(s1);
|
||||
acb_init(s2);
|
||||
acb_init(eps);
|
||||
acb_init(g1);
|
||||
acb_init(g2);
|
||||
acb_init(pq);
|
||||
acb_init(pqk2);
|
||||
|
||||
acb_set(args, s);
|
||||
acb_one(arg1s);
|
||||
acb_sub(arg1s, arg1s, s, prec);
|
||||
|
||||
if (chi->parity)
|
||||
{
|
||||
acb_add_ui(args, args, 1, prec);
|
||||
acb_add_ui(arg1s, arg1s, 1, prec);
|
||||
|
||||
for (k = 2; k < len; k++)
|
||||
acb_mul_si(a + k, a + k, k, prec);
|
||||
|
||||
}
|
||||
|
||||
acb_mul_2exp_si(args, args, -1);
|
||||
acb_mul_2exp_si(arg1s, arg1s, -1);
|
||||
acb_conj(arg1s, arg1s);
|
||||
|
||||
arb_const_pi(acb_realref(pq), prec);
|
||||
acb_div_ui(pq, pq, G->q, prec);
|
||||
|
||||
acb_dirichlet_root_number(eps, G, chi, prec);
|
||||
|
||||
for (k = 1; k < len; k++)
|
||||
{
|
||||
if (acb_is_zero(a + k))
|
||||
continue;
|
||||
acb_mul_ui(pqk2, pq, k*k, prec);
|
||||
/* FIXME: accuracy can be very bad */
|
||||
acb_hypgeom_gamma_upper(g1, args, pqk2, 2, prec);
|
||||
acb_hypgeom_gamma_upper(g2, arg1s, pqk2, 2, prec);
|
||||
acb_addmul(s1, a + k, g1, prec);
|
||||
acb_addmul(s2, a + k, g2, prec);
|
||||
}
|
||||
|
||||
acb_conj(s2, s2);
|
||||
acb_mul(s2, eps, s2, prec);
|
||||
acb_add(res, s1, s2, prec);
|
||||
|
||||
acb_pow(g1, pq, args, prec);
|
||||
acb_mul(res, res, g1, prec);
|
||||
|
||||
acb_gamma(g1, args, prec);
|
||||
acb_div(res, res, g1, prec);
|
||||
|
||||
acb_clear(args);
|
||||
acb_clear(arg1s);
|
||||
acb_clear(s1);
|
||||
acb_clear(s2);
|
||||
acb_clear(eps);
|
||||
acb_clear(g1);
|
||||
acb_clear(g2);
|
||||
acb_clear(pq);
|
||||
acb_clear(pqk2);
|
||||
_acb_vec_clear(a, len);
|
||||
}
|
||||
|
|
@ -11,9 +11,6 @@
|
|||
|
||||
#include "acb_dirichlet.h"
|
||||
|
||||
#define nq 5
|
||||
#define nx 3
|
||||
|
||||
void
|
||||
test_dft(ulong q)
|
||||
{
|
||||
|
|
@ -85,6 +82,9 @@ test_dft(ulong q)
|
|||
acb_dirichlet_group_clear(G);
|
||||
}
|
||||
|
||||
#define nq 5
|
||||
#define nx 4
|
||||
|
||||
int main()
|
||||
{
|
||||
|
||||
|
|
@ -95,51 +95,61 @@ int main()
|
|||
|
||||
acb_ptr x;
|
||||
/* cannot test at s = 1 with hurwitz */
|
||||
const char * x_r[nx] = { "1", "0.5", "0.5" };
|
||||
const char * x_i[nx] = { "1", "0", "6" };
|
||||
const char * x_r[nx] = { "0.5", "1", "1", "0.5" };
|
||||
const char * x_i[nx] = { "0", "0", "1", "6" };
|
||||
|
||||
acb_t ref, res;
|
||||
/*
|
||||
default(realprecision, 54)
|
||||
X = [ 1 + I, 1/2, 1/2 + 6 * I ]
|
||||
default(realprecision, 100)
|
||||
X = [ 1/2, 1, 1 + I, 1/2 + 6 * I ]
|
||||
C = [Mod(2,3),Mod(4,5),Mod(11,61),Mod(2,91),Mod(3,800)]
|
||||
v = concat([ [lfun(c,x) | x<-X] | c<-C])
|
||||
apply(z->printf("\"%s\",\n",real(z)),v)
|
||||
apply(z->printf("\"%s\",\n",imag(z)),v)
|
||||
apply(z->printf("\"%s\",\n",precision(real(z),70)),v)
|
||||
apply(z->printf("\"%s\",\n",precision(imag(z),70)),v)
|
||||
*/
|
||||
const char * ref_r[nq * nx] = {
|
||||
"0.655527984002548033786648216345221087939439503905627469",
|
||||
"0.480867557696828626181220063235589877776829730832371863",
|
||||
"1.56831301727577320813799211138797101541772722814204757",
|
||||
"0.521271244517346991221550773660594765406476858135844321",
|
||||
"0.231750947504015755883383661760877226427886696409005898",
|
||||
"0.275543455389521803395512886745330595086898302178508437",
|
||||
"0.598221809458540554839300433683735304093606595684903281",
|
||||
"0.489264190003695740292779374874163221990017067040417393",
|
||||
"0.573331076412428980263984182365365715292560207445592018",
|
||||
"0.510279695870740409778738767334484809708615155539404548",
|
||||
"0.635626509594367380604827545000418331455019188562281349",
|
||||
"0.129304857274642475564179442785425797926079767522671163",
|
||||
"1.18088858810025653590356481638012816019876881487868657",
|
||||
"2.17175778983760437737667738885959799183430688287297767",
|
||||
"3.41568550810774629867945639900431994221065497147578087"
|
||||
"0.48086755769682862618122006323558987777682973083237186318155991033493348067023",
|
||||
"0.60459978807807261686469275254738524409468874936424685852329497846270772704212",
|
||||
"0.65552798400254803378664821634522108793943950390562746892997938695272366439543",
|
||||
"1.56831301727577320813799211138797101541772722814204756600556183229594192303169",
|
||||
"0.23175094750401575588338366176087722642788669640900589796611788733984438064542",
|
||||
"0.43040894096400403888943323295060542542457068254028965475700610399256121546113",
|
||||
"0.52127124451734699122155077366059476540647685813584432110209552194219857263213",
|
||||
"0.27554345538952180339551288674533059508689830217850843726588160302316425288119",
|
||||
"0.48926419000369574029277937487416322199001706704041739336662332096055431153634",
|
||||
"0.53487618489463947737652797478014028523937753967835874218217650866316877024685",
|
||||
"0.59822180945854055483930043368373530409360659568490328055967882287521072175030",
|
||||
"0.57333107641242898026398418236536571529256020744559201810788762987710987052691",
|
||||
"0.63562650959436738060482754500041833145501918856228134932017130492895859585857",
|
||||
"0.58238559844425623388521061809515173901502388709897277008305502466245886012927",
|
||||
"0.51027969587074040977873876733448480970861515553940454822311416220939931612454",
|
||||
"0.12930485727464247556417944278542579792607976752267116291530276949569247783180",
|
||||
"2.1717577898376043773766773888595979918343068828729776727671442492863726379365",
|
||||
"1.59480084775553338997043324643146074042666522080647790312953333780585384429160",
|
||||
"1.18088858810025653590356481638012816019876881487868657080290826197336624703523",
|
||||
"3.4156855081077462986794563990043199422106549714757808711688925562811420743252"
|
||||
};
|
||||
const char * ref_i[nq * nx] = {
|
||||
"0.220206044893215842652155131408935133596486560067476474",
|
||||
"0",
|
||||
"-0.969458654385732175077973304161399773236957587792986099",
|
||||
"0.354614573267731219242838516784605303288232150760467423",
|
||||
"0",
|
||||
"-0.995392028773643947872231871832838543767598301887892796",
|
||||
"1.04370497561090171487193145841005574472705644411957863",
|
||||
"-0.108902811943905225853677097712717212629591264759957602",
|
||||
"-0.232114369998608907337769019848201890558327186146689311",
|
||||
"-0.133300066189980774635445078240315148544665020358019145",
|
||||
"0.0119464572932630291870372694406253796888930803905106876",
|
||||
"-0.567660589679294457801153713636532209809112025502518666",
|
||||
"-0.654079942571300523223917775358845549990877148918886474",
|
||||
"0.970337207245832214408380510793679653538607483205616894",
|
||||
"-1.43652482351673593824956935036654893593947145947637807"
|
||||
"0.22020604489321584265215513140893513359648656006747647412735100189927138978498",
|
||||
"-0.96945865438573217507797330416139977323695758779298609887525494651999839519954",
|
||||
"0",
|
||||
"0",
|
||||
"0.35461457326773121924283851678460530328823215076046742298881317219868652838939",
|
||||
"-0.99539202877364394787223187183283854376759830188789279576725217941822771397165",
|
||||
"-0.108902811943905225853677097712717212629591264759957601734509998247516905304233",
|
||||
"0.19365188807114177330261930603607468329782453247898970868809396136721159487113",
|
||||
"1.04370497561090171487193145841005574472705644411957863465786632716182836513300",
|
||||
"-0.23211436999860890733776901984820189055832718614668931112661490548978197438207",
|
||||
"0.0119464572932630291870372694406253796888930803905106876170081736901046822264969",
|
||||
"-0.30766179419579263291559128367018027923130717977616497873373231273162909091033",
|
||||
"-0.13330006618998077463544507824031514854466502035801914491875951675974942110581",
|
||||
"-0.56766058967929445780115371363653220980911202550251866578697696371454503377792",
|
||||
"0.97033720724583221440838051079367965353860748320561689369676918117330599741939",
|
||||
"0.38019821209845078839690749709574556766339551329040418677194351169649156337996",
|
||||
"-0.65407994257130052322391777535884554999087714891888647434188545263517159257326",
|
||||
"-1.43652482351673593824956935036654893593947145947637807266976772120632968400451"
|
||||
};
|
||||
|
||||
flint_printf("l....");
|
||||
|
|
@ -175,13 +185,15 @@ int main()
|
|||
for (j = 0; j < nx; j++)
|
||||
{
|
||||
|
||||
if (arb_set_str(acb_realref(ref), ref_r[i * nx + j], prec - 10) ||
|
||||
arb_set_str(acb_imagref(ref), ref_i[i * nx + j], prec - 10) )
|
||||
if (arb_set_str(acb_realref(ref), ref_r[i * nx + j], prec) ||
|
||||
arb_set_str(acb_imagref(ref), ref_i[i * nx + j], prec) )
|
||||
{
|
||||
flint_printf("error while setting ref <- %s+I*%s\n",
|
||||
ref_r[i * nx + j], ref_i[i * nx + j]);
|
||||
abort();
|
||||
}
|
||||
arb_add_error_2exp_si(acb_realref(ref), -prec+10);
|
||||
arb_add_error_2exp_si(acb_imagref(ref), -prec+10);
|
||||
|
||||
acb_dirichlet_l_hurwitz(res, x + j, G, chi, prec + 10);
|
||||
|
||||
|
|
@ -191,15 +203,34 @@ int main()
|
|||
flint_printf("q = %wu\n", q[i]);
|
||||
flint_printf("m = %wu\n", m[i]);
|
||||
flint_printf("x = ");
|
||||
acb_printd(x, 54);
|
||||
acb_printd(x + j, 54);
|
||||
flint_printf("\nref = ");
|
||||
acb_printd(ref, 54);
|
||||
flint_printf("\nl(chi,x) = ");
|
||||
flint_printf("\nl_hurwitz(chi,x) = ");
|
||||
acb_printd(res, 54);
|
||||
flint_printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
#if INCGAM_FIXED
|
||||
acb_dirichlet_l_incgam(res, x + j, G, chi, prec + 10);
|
||||
|
||||
if (!acb_contains(ref, res))
|
||||
{
|
||||
flint_printf("FAIL:\n\n");
|
||||
flint_printf("q = %wu\n", q[i]);
|
||||
flint_printf("m = %wu\n", m[i]);
|
||||
flint_printf("x = ");
|
||||
acb_printd(x + j, 54);
|
||||
flint_printf("\nref = ");
|
||||
acb_printd(ref, 54);
|
||||
flint_printf("\nl_incgam(chi,x) = ");
|
||||
acb_printd(res, 54);
|
||||
flint_printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
#endif
|
||||
|
||||
}
|
||||
acb_dirichlet_char_clear(chi);
|
||||
acb_dirichlet_group_clear(G);
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue