rho.c

/* Dickman's rho function (to compute probability of success of ecm).

Copyright 2004, 2005, 2006, 2008, 2009, 2010, 2011, 2012, 2013
Alexander Kruppa, Paul Zimmermann.

This file is part of the ECM Library.

The ECM Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The ECM Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the ECM Library; see the file COPYING.LIB.  If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */

/* define TESTDRIVE to compile rho as a stand-alone program, in which case
   you need to have libgsl installed */

#include "config.h"
#if defined(TESTDRIVE)
#define _ISOC99_SOURCE 1
#endif
#if defined(DEBUG_NUMINTEGRATE) || defined(TESTDRIVE)
# include <stdio.h>
#endif
#include <stdlib.h>
#include <math.h>
#if defined(TESTDRIVE)
#include <string.h>
#include <primesieve.h>
#endif
#if defined(TESTDRIVE)
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_integration.h>
#endif
#include "ecm-impl.h"

/* For Suyama's curves, we have a known torsion factor of 12 = 2^2*3^1, and
   an average extra exponent of 1/2 for 2, and 1/3 for 3 due to the probability
   that the group order divided by 12 is divisible by 2 or 3, thus on average
   we should have 2^2.5*3^1.333 ~ 24.5, however experimentally we have
   2^3.323*3^1.687 ~ 63.9 (see Alexander Kruppa's thesis, Table 5.1 page 96,
   row sigma=2, http://tel.archives-ouvertes.fr/tel-00477005/en/).
   The exp(ECM_EXTRA_SMOOTHNESS) value takes into account the extra
   smoothness with respect to a random number. */
#ifndef ECM_EXTRA_SMOOTHNESS
#define ECM_EXTRA_SMOOTHNESS 3.134
#endif

#define M_PI_SQR   9.869604401089358619 /* Pi^2 */
#define M_PI_SQR_6 1.644934066848226436 /* Pi^2/6 */
/* gsl_math.h defines M_EULER */
#ifndef M_EULER
#define M_EULER    0.577215664901532861
#endif
#define M_EULER_1   0.422784335098467139 /* 1 - Euler */

#ifndef MAX
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#endif
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))