1 /* Approximation of Dickman's rho function.
2
3 Copyright 2010, 2012 Paul Zimmermann
4
5 This file is part of CADO-NFS.
6
7 CADO-NFS is free software; you can redistribute it and/or modify it under the
8 terms of the GNU Lesser General Public License as published by the Free
9 Software Foundation; either version 2.1 of the License, or (at your option)
10 any later version.
11
12 CADO-NFS is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
14 A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
15 details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with CADO-NFS; see the file COPYING. If not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 */
21
22 #include "cado.h" // IWYU pragma: keep
23 #include <stdio.h>
24 #include <math.h>
25 #include "rho.h"
26
27 /* return an approximation of rho(x) with relative error <= 2.1e-6 (about)
28 for x <= 15, and an absolute error <= 3.5e-21 for x > 15. */
29 double
dickman_rho(double x)30 dickman_rho (double x)
31 {
32 if (x < 0)
33 {
34 return 0.0;
35 }
36
37 if (x <= 1.0)
38 return 1.0; /* no error */
39
40 if (x <= 2.0) /* 1 < x <= 2 */
41 {
42 x = 2 * x - 3.0; /* -1 <= x <= 1 */
43 /* the following approximation was computed with Sage:
44 f1 = dickman_rho.power_series(1, 64)
45 then using Maple:
46 p:=numapprox[minimax](f1,x=-1..1,10,1/f1).
47 It gives a relative error less than 1.2e-9
48 (numapprox[infnorm](f1/p-1,x=-1..1))
49 In this interval, rho(x) = 1 - log(x), but the math library log()
50 function is slower than the (relatively short) power series here. */
51 return .59453489206592253066+(-.33333334045356381396+(.55555549124525324224e-1+(-.12345536626584334599e-1+(.30864445307049269724e-2+(-.82383544119364408582e-3+(.22867526556051420719e-3+(-.63554707267886054080e-4+(.18727717457043009514e-4+(-.73223168705152723341e-5+.21206262864513086463e-5*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
52 }
53
54 if (x <= 3.0) /* 2 < x <= 3 */
55 {
56 x = 2 * x - 5.0; /* -1 <= x <= 1 */
57 /* the following approximation was computed with Sage:
58 f2 = dickman_rho.power_series(2, 64)
59 then using Maple:
60 p:=numapprox[minimax](f2,x=-1..1,10,1/f2).
61 It gives a relative error less than 8.9e-10
62 (numapprox[infnorm](f2/p-1,x=-1..1)) */
63 return .13031956190159661490+(-.11890697945147816983+(.45224027972316855319e-1+(-.97335515428611864000e-2+(.20773419306178919493e-2+(-.45596184871832967815e-3+(.10336375145598596368e-3+(-.23948090479838959902e-4+(.58842414590359757174e-5+(-.17744830412711835918e-5+.42395344408760226490e-6*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
64 }
65
66 if (x <= 4.0) /* 3 < x <= 4 */
67 {
68 x = 2 * x - 7.0; /* -1 <= x <= 1 */
69 /* the following approximation was computed with Sage:
70 f3 = dickman_rho.power_series(3, 64)
71 then using Maple:
72 p:=numapprox[minimax](f3,x=-1..1,10,1/f3).
73 It gives a relative error less than 7.7e-10
74 (numapprox[infnorm](f3/p-1,x=-1..1)) */
75 return .16229593252987041396e-1+(-.18617080355889960242e-1+(.98231465619823138710e-2+(-.30890609038683816355e-2+(.67860233545741575724e-3+(-.13691972815251836202e-3+(.27142792736320271161e-4+(-.54020882619058856352e-5+(.11184852221470961276e-5+(-.26822513121459348050e-6+.53524419961799178404e-7*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
76 }
77
78 if (x <= 5.0) /* 4 < x <= 5 */
79 {
80 x = 2 * x - 9.0; /* -1 <= x <= 1 */
81 /* the following approximation was computed with Sage:
82 f4 = dickman_rho.power_series(4, 64)
83 then using Maple:
84 p:=numapprox[minimax](f4,x=-1..1,10,1/f4).
85 It gives a relative error less than 8.0e-10
86 (numapprox[infnorm](f4/p-1,x=-1..1)) */
87 return .13701177421087314589e-2+(-.18032881447340571269e-2+(.11344648599460666353e-2+(-.44785452702335769551e-3+(.12312893763747451118e-3+(-.26025855424244979420e-4+(.49439561514862370128e-5+(-.89908455922138763242e-6+(.16408096470458054605e-6+(-.32859809253342863007e-7+.55954809519625267389e-8*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
88 }
89
90 if (x <= 6.0) /* 5 < x <= 6 */
91 {
92 x = 2 * x - 11.0; /* -1 <= x <= 1 */
93 /* the following approximation was computed with Sage:
94 f5 = dickman_rho.power_series(5, 64)
95 then using Maple:
96 p:=numapprox[minimax](f5,x=-1..1,10,1/f5).
97 It gives a relative error less than 1.1e-9
98 (numapprox[infnorm](f5/p-1,x=-1..1)) */
99 return .86018611204769414547e-4+(-.12455615869107014099e-3+(.87629281627689306328e-4+(-.39688578340310188782e-4+(.12884593599298603240e-4+(-.31758435723373651519e-5+(.63480088103905159169e-6+(-.11347189099214240765e-6+(.19390866486609035589e-7+(-.34493644515794744033e-8+.52005397653635760845e-9*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
100 }
101
102 if (x <= 7.0) /* 6 < x <= 7 */
103 {
104 x = 2 * x - 13.0; /* -1 <= x <= 1 */
105 /* the following approximation was computed with Sage:
106 f6 = dickman_rho.power_series(6, 64)
107 then using Maple:
108 p:=numapprox[minimax](f6,x=-1..1,10,1/f6).
109 It gives a relative error less than 1.6e-9
110 (numapprox[infnorm](f6/p-1,x=-1..1)) */
111 return .42503555236283394611e-5+(-.66168162622648216578e-5+(.50451140426428793943e-5+(-.25056278110193975867e-5+(.90780066639285274491e-6+(-.25409423435688445924e-6+(.56994106166489666850e-7+(-.10719440462879689904e-7+(.18244728209885020524e-8+(-.30691539036795813350e-9+.42848520313019466802e-10*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
112 }
113
114 if (x <= 8.0) /* 7 < x <= 8 */
115 {
116 x = 2 * x - 15.0; /* -1 <= x <= 1 */
117 /* the following approximation was computed with Sage:
118 f7 = dickman_rho.power_series(7, 64)
119 then using Maple:
120 p:=numapprox[minimax](f7,x=-1..1,10,1/f7).
121 It gives a relative error less than 2.6e-9
122 (numapprox[infnorm](f7/p-1,x=-1..1)) */
123 return .17178674964103956208e-6+(-.28335703551339176870e-6+(.23000575057461623263e-6+(-.12233608702949464058e-6+(.47877500071532198696e-7+(-.14657787264925438894e-7+(.36368896326425931590e-8+(-.74970784803203153691e-9+(.13390834297490308260e-9+(-.22532330111335516039e-10+.30448478436493554124e-11*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
124 }
125
126 if (x <= 9.0) /* 8 < x <= 9 */
127 {
128 x = 2 * x - 17.0; /* -1 <= x <= 1 */
129 /* the following approximation was computed with Sage:
130 f8 = dickman_rho.power_series(8, 64)
131 then using Maple:
132 p:=numapprox[minimax](f8,x=-1..1,10,1/f8).
133 It gives a relative error less than 4.4e-9
134 (numapprox[infnorm](f8/p-1,x=-1..1)) */
135 return .58405695885954012372e-8+(-.10105102932563601939e-7+(.86312378192120507636e-8+(-.48483948355429166070e-8+(.20129738917787451004e-8+(-.65800969013673567377e-9+(.17591641970556292848e-9+(-.39380340623747158286e-10+(.75914903051312238224e-11+(-.13345077002021028171e-11+.18138420114766605833e-12*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
136 }
137
138 if (x <= 10.0) /* 9 < x <= 10 */
139 {
140 x = 2 * x - 19.0; /* -1 <= x <= 1 */
141 /* the following approximation was computed with Sage:
142 f9 = dickman_rho.power_series(9, 64)
143 then using Maple:
144 p:=numapprox[minimax](f9,x=-1..1,10,1/f9).
145 It gives a relative error less than 7.5e-9
146 (numapprox[infnorm](f9/p-1,x=-1..1)) */
147 return .17063527516986966493e-9+(-.30739839907061929947e-9+(.27401311519313991860e-9+(-.16103964872970941164e-9+(.70152211692303181302e-10+(-.24143747383189013980e-10+(.68287507297199575436e-11+(-.16282751120203734207e-11+(.33676191519721199042e-12+(-.63207992345353337190e-13+.88821884863349801119e-14*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
148 }
149
150 if (x <= 11.0) /* 10 < x <= 11 */
151 {
152 x = 2 * x - 21.0; /* -1 <= x <= 1 */
153 /* the following approximation was computed with Sage:
154 f10 = dickman_rho.power_series(10, 64)
155 then using Maple:
156 Digits:=20:
157 p:=numapprox[minimax](f10,x=-1..1,10,1/f10);
158 It gives a relative error less than 1.2e-8
159 (numapprox[infnorm](f10/p-1,x=-1..1)) */
160 return .43559526824795217852e-11+(-.81254892374415436871e-11+(.75124702775657186758e-11+(-.45879032842537328635e-11+(.20810226978994617073e-11+(-.74742298793568571039e-12+(.22118242218483402875e-12+(-.55382747996458896000e-13+(.12115887170289387926e-13+(-.24203055184086403014e-14+.35552802442946562758e-15*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
161 }
162
163 if (x <= 12.0) /* 11 < x <= 12 */
164 {
165 x = 2 * x - 23.0; /* -1 <= x <= 1 */
166 /* the following approximation was computed with Sage:
167 f11 = dickman_rho.power_series(11, 128)
168 then using Maple:
169 Digits:=20:
170 p:=numapprox[minimax](f11,x=-1..1,10,1/f11);
171 It gives a relative error less than 1.9e-8
172 (numapprox[infnorm](f11/p-1,x=-1..1)) */
173 return .98476422850126247236e-13+(-.18938924316092442150e-12+(.18075811456452207443e-12+(-.11411568337724330786e-12+(.53590747173713771259e-13+(-.19960825800354613327e-13+(.61359753086709747219e-14+(-.15991438833241030242e-14+(.36598455644175594635e-15+(-.76985111329665577357e-16+.11768795877242596285e-16*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
174 }
175
176 if (x <= 13.0) /* 12 < x <= 13 */
177 {
178 x = 2 * x - 25.0; /* -1 <= x <= 1 */
179 /* the following approximation was computed with Sage:
180 f12 = dickman_rho.power_series(12, 80)
181 then using Maple:
182 Digits:=20:
183 p:=numapprox[minimax](f12,x=-1..1,10,1/f12);
184 It gives a relative error less than 2.7e-8
185 (numapprox[infnorm](f12/p-1,x=-1..1)) */
186 return .19934633853036819003e-14+(-.39390567960226078379e-14+(.38665626810909948382e-14+(-.25132152019284365490e-14+(.12165850034414917031e-14+(-.46768237317409899068e-15+(.14856141016839901209e-15+(-.40061895880500065473e-16+(.95227478488277343487e-17+(-.20899104913735568600e-17+.32985921210508602052e-18*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
187 }
188
189 if (x <= 14.0) /* 13 < x <= 14 */
190 {
191 x = 2 * x - 27.0; /* -1 <= x <= 1 */
192 /* the following approximation was computed with Sage:
193 f13 = dickman_rho.power_series(13, 80)
194 then using Maple:
195 Digits:=20:
196 p:=numapprox[minimax](f13,x=-1..1,10,1/f13);
197 It gives a relative error less than 3.7e-8
198 (numapprox[infnorm](f13/p-1,x=-1..1)) */
199 return .36468388101796492224e-16+(-.73831973563371575758e-16+(.74312671625965740134e-16+(-.49570201041373844974e-16+(.24648267921024808625e-16+(-.97426560588907822003e-17+(.31850597161723004412e-17+(-.88485037859050841047e-18+(.21736499712587315716e-18+(-.49450337815428832609e-19+.80095079043879918725e-20*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
200 }
201
202 if (x <= 15.0) /* 14 < x <= 15 */
203 {
204 x = 2 * x - 29.0; /* -1 <= x <= 1 */
205 /* the following approximation was computed with Sage:
206 f14 = dickman_rho.power_series(14, 128)
207 then using Maple:
208 Digits:=20:
209 p:=numapprox[minimax](f14,x=-1..1,10,1/f14);
210 It gives a relative error less than 4.9e-8
211 (numapprox[infnorm](f14/p-1,x=-1..1)) */
212 return .60765098932992206115e-18+(-.12575305192897762077e-17+(.12946447870858888856e-17+(-.88393035793877424780e-18+(.45020754487871940709e-18+(-.18241969586978266280e-18+(.61176643040672211241e-19+(-.17448851858831783538e-19+(.44129951768989290506e-20+(-.10359223715080577635e-20+.17147161493537144558e-21*x)*x)*x)*x)*x)*x)*x)*x)*x)*x;
213 }
214
215 /* for x > 15, Dickman's function differs from the approximation
216 below by at most 3.5e-21 */
217 {
218 static int count = 0;
219 if (count ++ == 0)
220 fprintf (stderr, "# Warning: Dickman rho is imprecise for x > 15\n");
221 }
222 return 0.0332357434363490 * pow (x, -x);
223 }
224
225 #ifndef M_EULER
226 #define M_EULER 0.57721566490153286060651209008240243104
227 #endif
228
229 /* Return the probability that a random integer close to N is N^(1/x)-smooth */
230 double
dickman_rho_local(const double x,const double N)231 dickman_rho_local (const double x, const double N)
232 {
233 return dickman_rho(x) - M_EULER * dickman_rho(x - 1.) / log(N);
234 }
235