Skip to content

Instantly share code, notes, and snippets.

@genkuroki
Last active August 20, 2017 05:48
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save genkuroki/dabc2fb61ceb80e40fbc189676754fa2 to your computer and use it in GitHub Desktop.
Julia/exp and log integral functions/Ei function in Julia - Part 3.ipynb
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"metadata": {},
"cell_type": "markdown",
"source": "# Ei function in Julia - Part 3\n\n**Comparison between Julia versions and gcc verions of the Ei function**\n\nGen KUROKI\n\n2017-08-20\n\nConclusion: Julia versions are on par with gcc versions."
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "versioninfo()",
"execution_count": 1,
"outputs": [
{
"output_type": "stream",
"text": "Julia Version 0.6.0\nCommit 903644385b* (2017-06-19 13:05 UTC)\nPlatform Info:\n OS: Windows (x86_64-w64-mingw32)\n CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz\n WORD_SIZE: 64\n BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)\n LAPACK: libopenblas64_\n LIBM: libopenlibm\n LLVM: libLLVM-3.9.1 (ORCJIT, haswell)\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "run(`gcc --version`)",
"execution_count": 2,
"outputs": [
{
"output_type": "stream",
"text": "gcc (Rev5, Built by MSYS2 project) 5.3.0\r\nCopyright (C) 2015 Free Software Foundation, Inc.\r\nThis is free software; see the source for copying conditions. There is NO\r\nwarranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.\r\n\r\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "## GNU MPFR eint()\n\nsetprecision(256)\nfunction Ei_MPFR(x::BigFloat)\n z = BigFloat()\n ccall((:mpfr_eint, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), \n &z, &x, Base.MPFR.ROUNDING_MODE[])\n return z\nend\nEi_MPFR(big\"0.1\")",
"execution_count": 3,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 3,
"data": {
"text/plain": "-1.622812813969276674965682999227475254223371617672384809534366523239861144920695"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@show Ei_MPFR(big\"-10.0\")\n@show Ei_MPFR(big\"-3.0\")\n@show Ei_MPFR(big\"0.0\")\n@show Ei_MPFR(big\"5.0\")\n@show Ei_MPFR(big\"20.0\")\n@show Ei_MPFR(big\"100.0\");",
"execution_count": 4,
"outputs": [
{
"output_type": "stream",
"text": "Ei_MPFR(@big_str(\"-10.0\")) = BigFloat(NaN, 256)\nEi_MPFR(@big_str(\"-3.0\")) = BigFloat(NaN, 256)\nEi_MPFR(@big_str(\"0.0\")) = BigFloat(-Inf, 256)\nEi_MPFR(@big_str(\"5.0\")) = 4.018527535580317745509142179379586709541908739919593043418287525164264481392828e+01\nEi_MPFR(@big_str(\"20.0\")) = 2.561565266405658882048112080409807182938269835774034260637412340609842955216052e+07\nEi_MPFR(@big_str(\"100.0\")) = 2.715552744853879821914014642310825410295793934191620986249744933654559905245318e+41\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "-log(eps())",
"execution_count": 5,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 5,
"data": {
"text/plain": "36.04365338911715"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# Julia version of\n# Numerical Recipes in C\n# http://www.aip.de/groups/soe/local/numres/bookcpdf/c6-3.pdf\n\nfunction Ei_numres(x::Float64)\n const N = 100\n const g = Float64(eulergamma)\n if x <= 0.0\n error(\"Argument must be positive in Ei\")\n end\n if x < eps()\n return log(x)+g\n end\n if x <= -log(eps())\n s = 0.0\n fact = 1.0\n k = 0.0\n for k in 1:N\n fact *= x/k\n term= fact/k\n s += term\n if term < eps()*s\n break\n end\n end\n if (k >= N)\n error(\"Series failed in Ei\")\n end\n return log(x) + g + s\n else\n s = 0.0\n term = 1.0\n prev = term\n k = 0.0\n for k in 1:N\n term *= k/x\n if term < eps()\n break\n end\n if term < prev\n s += term\n else\n s -= prev\n break\n end\n prev = term\n end\n return exp(x)*(1.0 + s)/x\n end \nend\n\nEi_julia_numres(x::Float64) = Ei_numres(x)\n\nEi_numres(0.1)",
"execution_count": 6,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 6,
"data": {
"text/plain": "-1.6228128139692766"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@show Ei_numres(5.0)\n@show Ei_numres(20.0)\n@show Ei_numres(100.0);",
"execution_count": 7,
"outputs": [
{
"output_type": "stream",
"text": "Ei_numres(5.0) = 40.18527535580317\nEi_numres(20.0) = 2.56156526640566e7\nEi_numres(100.0) = 2.7155527448538802e41\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "BigFloat(γ)",
"execution_count": 8,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 8,
"data": {
"text/plain": "5.772156649015328606065120900824024310421593359399235988057672348848677267776685e-01"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# double (not float) version of\n# Numerical Recipes in C\n# http://www.aip.de/groups/soe/local/numres/bookcpdf/c6-3.pdf\n\nC_code = raw\"\"\"\n#include <stdio.h>\n#include <stddef.h>\n#include <stdlib.h>\n#include <math.h>\n#include <float.h>\n#define EULER 0.577215664901532860606512\n#define MAXIT 100\n#define FPMIN 1.0e-30\n#define EPS (10.0 * DBL_EPSILON)\n\nvoid nrerror(char error_text[])\n{\n fprintf(stderr,\"Numerical Recipes run-time error...\\n\");\n fprintf(stderr,\"%s\\n\",error_text);\n fprintf(stderr,\"...now exiting to system...\\n\");\n exit(1);\n}\n\ndouble ei(double x)\n{\n void nrerror(char error_text[]);\n int k;\n double fact,prev,sum,term;\n if (x <= 0.0) nrerror(\"Bad argument in ei\");\n if (x < FPMIN) return log(x)+EULER;\n if (x <= -log(EPS)) {\n sum=0.0;\n fact=1.0;\n for (k=1;k<=MAXIT;k++) {\n fact *= x/k;\n term=fact/k;\n sum += term;\n if (term < EPS*sum) break;\n }\n if (k > MAXIT) nrerror(\"Series failed in ei\");\n return sum+log(x)+EULER;\n } else {\n sum=0.0;\n term=1.0;\n for (k=1;k<=MAXIT;k++) {\n prev=term;\n term *= k/x;\n if (term < EPS) break;\n if (term < prev) sum += term;\n else {\n sum -= prev;\n break;\n }\n }\n return exp(x)*(1.0+sum)/x;\n }\n}\n\"\"\"\n\nfilename = tempname()\nfilenamedl = filename * \".\" * Libdl.dlext\n\nopen(`gcc -Wall -O3 -march=native -xc -shared -o $filenamedl -`, \"w\") do f\n print(f, C_code)\nend\n\nrun(`ls -l $filenamedl`)",
"execution_count": 9,
"outputs": [
{
"output_type": "stream",
"text": "-rwxr-xr-x 1 genkuroki 197610 123491 Aug 20 14:43 C:\\Users\\GENKUR~1\\AppData\\Local\\Temp\\jl_A218.tmp.dll\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "const Ei_numres_lib = filename\nEi_gcc_numres(x::Float64) = ccall((:ei, Ei_numres_lib), Float64, (Float64,), x)",
"execution_count": 10,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 10,
"data": {
"text/plain": "Ei_gcc_numres (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@show Ei_gcc_numres(5.0)\n@show Ei_gcc_numres(20.0)\n@show Ei_gcc_numres(100.0);",
"execution_count": 11,
"outputs": [
{
"output_type": "stream",
"text": "Ei_gcc_numres(5.0) = 40.18527535580317\nEi_gcc_numres(20.0) = 2.5615652664056588e7\nEi_gcc_numres(100.0) = 2.7155527448538787e41\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# Julia Float64 (not long double) version of\n# http://www.mymathlib.com/c_source/functions/exponential_integrals/exponential_integral_Ei.c\n\nfunction Ei_mml(x :: Float64)\n if x < -5.0\n return Ei_mml_Continued_Fraction(x)\n elseif x == 0.0\n return -typemax(Float64)\n elseif x < 6.8\n return Ei_mml_Power_Series(x)\n elseif x < 50.0\n return Ei_mml_Argument_Addition_Series(x)\n else\n return Ei_mml_Continued_Fraction(x)\n end\nend\n\nEi_julia_mml(x) = Ei_mml(x)\n\nfunction Ei_mml_Continued_Fraction(x :: Float64)\n Am1 = 1.0\n A0 = 0.0\n Bm1 = 0.0\n B0 = 1.0\n a = exp(x)\n b = -x + 1.0\n Ap1 = b * A0 + a * Am1\n Bp1 = b * B0 + a * Bm1\n j = 1;\n\n a = 1.0\n while abs(Ap1 * B0 - A0 * Bp1) > eps() * abs(A0 * Bp1)\n if abs(Bp1) > 1.0\n Am1 = A0 / Bp1\n A0 = Ap1 / Bp1\n Bm1 = B0 / Bp1\n B0 = 1.0\n else\n Am1 = A0\n A0 = Ap1\n Bm1 = B0\n B0 = Bp1\n end\n a = -j * j\n b += 2.0\n Ap1 = b * A0 + a * Am1\n Bp1 = b * B0 + a * Bm1\n j += 1\n end\n return -Ap1 / Bp1\nend\n\nfunction Ei_mml_Power_Series(x :: Float64)\n xn = -x\n Sn = -x\n Sm1 = 0.0\n hsum = 1.0\n g = Float64(eulergamma)\n y = 1.0\n fact = 1.0\n \n if x == 0.0\n return -typemax(Float64)\n end\n \n while abs(Sn - Sm1) > eps() * abs(Sm1)\n Sm1 = Sn\n y += 1.0\n xn *= -x\n fact *= y\n hsum += (1.0 / y)\n Sn += hsum * xn / fact\n end\n return g + log(abs(x)) - exp(x) * Sn\nend\n\nfunction Ei_mml_Argument_Addition_Series(x :: Float64)\n ei = [\n 1.915047433355013959531e2, 4.403798995348382689974e2,\n 1.037878290717089587658e3, 2.492228976241877759138e3,\n 6.071406374098611507965e3, 1.495953266639752885229e4,\n 3.719768849068903560439e4, 9.319251363396537129882e4,\n 2.349558524907683035782e5, 5.955609986708370018502e5,\n 1.516637894042516884433e6, 3.877904330597443502996e6,\n 9.950907251046844760026e6, 2.561565266405658882048e7,\n 6.612718635548492136250e7, 1.711446713003636684975e8,\n 4.439663698302712208698e8, 1.154115391849182948287e9,\n 3.005950906525548689841e9, 7.842940991898186370453e9,\n 2.049649711988081236484e10, 5.364511859231469415605e10,\n 1.405991957584069047340e11, 3.689732094072741970640e11,\n 9.694555759683939661662e11, 2.550043566357786926147e12,\n 6.714640184076497558707e12, 1.769803724411626854310e13,\n 4.669055014466159544500e13, 1.232852079912097685431e14,\n 3.257988998672263996790e14, 8.616388199965786544948e14,\n 2.280446200301902595341e15, 6.039718263611241578359e15,\n 1.600664914324504111070e16, 4.244796092136850759368e16,\n 1.126348290166966760275e17, 2.990444718632336675058e17,\n 7.943916035704453771510e17, 2.111342388647824195000e18,\n 5.614329680810343111535e18, 1.493630213112993142255e19,\n 3.975442747903744836007e19, 1.058563689713169096306e20\n ]\n k = Int(floor(x + 0.5))\n j = 0\n xx = Float64(k)\n dx = x - xx\n xxj = xx\n edx = exp(dx)\n Sm = 1.0\n Sn = (edx - 1.0) / xxj\n term = typemax(Float64)\n fact = 1.0\n dxj = 1.0\n\n while abs(term) > eps() * abs(Sn)\n j += 1\n fact *= Float64(j)\n xxj *= xx\n dxj *= -dx\n Sm += (dxj / fact)\n term = ( fact * (edx * Sm - 1.0) ) / xxj\n Sn += term\n end\n \n #return ei[k-7] + Sn * expl(xx);\n return ei[k-7+1] + Sn * exp(xx) \nend\n\n@show Ei_mml(-10.0)\n@show Ei_mml(-3.0)\n@show Ei_mml(0.0)\n@show Ei_mml(5.0)\n@show Ei_mml(20.0)\n@show Ei_mml(100.0);",
"execution_count": 12,
"outputs": [
{
"output_type": "stream",
"text": "Ei_mml(-10.0) = -4.156968929685325e-6\nEi_mml(-3.0) = -0.013048381094195927\nEi_mml(0.0) = -Inf\nEi_mml(5.0) = 40.1852753558037\nEi_mml(20.0) = 2.5615652664056588e7\nEi_mml(100.0) = 2.7155527448538795e41\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# http://www.mymathlib.com/c_source/functions/exponential_integrals/exponential_integral_Ei.c\n#\n# long double version\n\nC_code = raw\"\"\"\n////////////////////////////////////////////////////////////////////////////////\n// File: exponential_integral_Ei.c //\n// Routine(s): //\n// Exponential_Integral_Ei //\n// xExponential_Integral_Ei //\n////////////////////////////////////////////////////////////////////////////////\n\n#include <math.h> // required for fabsl(), expl() and logl() \n#include <float.h> // required for LDBL_EPSILON, DBL_MAX\n\n// Internally Defined Routines //\ndouble Exponential_Integral_Ei( double x );\nlong double xExponential_Integral_Ei( long double x );\n\nstatic long double Continued_Fraction_Ei( long double x );\nstatic long double Power_Series_Ei( long double x );\nstatic long double Argument_Addition_Series_Ei( long double x);\n\n\n// Internally Defined Constants //\nstatic const long double epsilon = 10.0 * LDBL_EPSILON;\n\n////////////////////////////////////////////////////////////////////////////////\n// double Exponential_Integral_Ei( double x ) //\n// //\n// Description: //\n// The exponential integral Ei(x) is the integral with integrand //\n// exp(t) / t //\n// where the integral extends from -inf to x. //\n// Note that there is a singularity at t = 0. Therefore for x > 0, the //\n// integral is defined to be the Cauchy principal value: //\n// lim { I[-inf, -eta] exp(-t) dt / t + I[eta, x] exp(-t) dt / t } //\n// in which the limit is taken as eta > 0 approaches 0 and I[a,b] //\n// denotes the integral from a to b. //\n// //\n// Arguments: //\n// double x The argument of the exponential integral Ei(). //\n// //\n// Return Value: //\n// The value of the exponential integral Ei evaluated at x. //\n// If x = 0.0, then Ei is -inf and -DBL_MAX is returned. //\n// //\n// Example: //\n// double y, x; //\n// //\n// ( code to initialize x ) //\n// //\n// y = Exponential_Integral_Ei( x ); //\n////////////////////////////////////////////////////////////////////////////////\ndouble Exponential_Integral_Ei( double x )\n{\n return (double) xExponential_Integral_Ei( (long double) x);\n}\n\n\n////////////////////////////////////////////////////////////////////////////////\n// long double xExponential_Integral_Ei( long double x ) //\n// //\n// Description: //\n// The exponential integral Ei(x) is the integral with integrand //\n// exp(t) / t //\n// where the integral extends from -inf to x. //\n// Note that there is a singularity at t = 0. Therefore for x > 0, the //\n// integral is defined to be the Cauchy principal value: //\n// lim { I[-inf, -eta] exp(-t) dt / t + I[eta, x] exp(-t) dt / t } //\n// in which the limit is taken as eta > 0 approaches 0 and I[a,b] //\n// denotes the integral from a to b. //\n// //\n// Arguments: //\n// long double x The argument of the exponential integral Ei(). //\n// //\n// Return Value: //\n// The value of the exponential integral Ei evaluated at x. //\n// If x = 0.0, then Ei is -inf and -DBL_MAX is returned. //\n// //\n// Example: //\n// long double y, x; //\n// //\n// ( code to initialize x ) //\n// //\n// y = xExponential_Integral_Ei( x ); //\n////////////////////////////////////////////////////////////////////////////////\n\nlong double xExponential_Integral_Ei( long double x )\n{\n if ( x < -5.0L ) return Continued_Fraction_Ei(x);\n if ( x == 0.0L ) return -DBL_MAX;\n if ( x < 6.8L ) return Power_Series_Ei(x);\n if ( x < 50.0L ) return Argument_Addition_Series_Ei(x);\n return Continued_Fraction_Ei(x);\n}\n\n////////////////////////////////////////////////////////////////////////////////\n// static long double Continued_Fraction_Ei( long double x ) //\n// //\n// Description: //\n// For x < -5 or x > 50, the continued fraction representation of Ei //\n// converges fairly rapidly. //\n// //\n// The continued fraction expansion of Ei(x) is: //\n// Ei(x) = -exp(x) { 1/(-x+1-) 1/(-x+3-) 4/(-x+5-) 9/(-x+7-) ... }. //\n// //\n// //\n// Arguments: //\n// long double x //\n// The argument of the exponential integral Ei(). //\n// //\n// Return Value: //\n// The value of the exponential integral Ei evaluated at x. //\n////////////////////////////////////////////////////////////////////////////////\n\nstatic long double Continued_Fraction_Ei( long double x )\n{\n long double Am1 = 1.0L;\n long double A0 = 0.0L;\n long double Bm1 = 0.0L;\n long double B0 = 1.0L;\n long double a = expl(x);\n long double b = -x + 1.0L;\n long double Ap1 = b * A0 + a * Am1;\n long double Bp1 = b * B0 + a * Bm1;\n int j = 1;\n\n a = 1.0L;\n while ( fabsl(Ap1 * B0 - A0 * Bp1) > epsilon * fabsl(A0 * Bp1) ) {\n if ( fabsl(Bp1) > 1.0L) {\n Am1 = A0 / Bp1;\n A0 = Ap1 / Bp1;\n Bm1 = B0 / Bp1;\n B0 = 1.0L;\n } else {\n Am1 = A0;\n A0 = Ap1;\n Bm1 = B0;\n B0 = Bp1;\n }\n a = -j * j;\n b += 2.0L;\n Ap1 = b * A0 + a * Am1;\n Bp1 = b * B0 + a * Bm1;\n j += 1;\n }\n return (-Ap1 / Bp1);\n}\n\n\n////////////////////////////////////////////////////////////////////////////////\n// static long double Power_Series_Ei( long double x ) //\n// //\n// Description: //\n// For -5 < x < 6.8, the power series representation for //\n// (Ei(x) - gamma - ln|x|)/exp(x) is used, where gamma is Euler's gamma //\n// constant. //\n// Note that for x = 0.0, Ei is -inf. In which case -DBL_MAX is //\n// returned. //\n// //\n// The power series expansion of (Ei(x) - gamma - ln|x|) / exp(x) is //\n// - Sum(1 + 1/2 + ... + 1/j) (-x)^j / j!, where the Sum extends //\n// from j = 1 to inf. //\n// //\n// Arguments: //\n// long double x //\n// The argument of the exponential integral Ei(). //\n// //\n// Return Value: //\n// The value of the exponential integral Ei evaluated at x. //\n////////////////////////////////////////////////////////////////////////////////\n\nstatic long double Power_Series_Ei( long double x )\n{ \n long double xn = -x;\n long double Sn = -x;\n long double Sm1 = 0.0L;\n long double hsum = 1.0L;\n long double g = 0.5772156649015328606065121L;\n long double y = 1.0L;\n long double factorial = 1.0L;\n \n if ( x == 0.0L ) return (long double) -DBL_MAX;\n \n while ( fabsl(Sn - Sm1) > epsilon * fabsl(Sm1) ) {\n Sm1 = Sn;\n y += 1.0L;\n xn *= (-x);\n factorial *= y;\n hsum += (1.0 / y);\n Sn += hsum * xn / factorial;\n }\n return (g + logl(fabsl(x)) - expl(x) * Sn);\n}\n\n\n////////////////////////////////////////////////////////////////////////////////\n// static long double Argument_Addition_Series_Ei(long double x) //\n// //\n// Description: //\n// For 6.8 < x < 50.0, the argument addition series is used to calculate //\n// Ei. //\n// //\n// The argument addition series for Ei(x) is: //\n// Ei(x+dx) = Ei(x) + exp(x) Sum j! [exp(j) expj(-dx) - 1] / x^(j+1), //\n// where the Sum extends from j = 0 to inf, |x| > |dx| and expj(y) is //\n// the exponential polynomial expj(y) = Sum y^k / k!, the Sum extending //\n// from k = 0 to k = j. //\n// //\n// Arguments: //\n// long double x //\n// The argument of the exponential integral Ei(). //\n// //\n// Return Value: //\n// The value of the exponential integral Ei evaluated at x. //\n////////////////////////////////////////////////////////////////////////////////\nstatic long double Argument_Addition_Series_Ei(long double x)\n{\n static long double ei[] = {\n 1.915047433355013959531e2L, 4.403798995348382689974e2L,\n 1.037878290717089587658e3L, 2.492228976241877759138e3L,\n 6.071406374098611507965e3L, 1.495953266639752885229e4L,\n 3.719768849068903560439e4L, 9.319251363396537129882e4L,\n 2.349558524907683035782e5L, 5.955609986708370018502e5L,\n 1.516637894042516884433e6L, 3.877904330597443502996e6L,\n 9.950907251046844760026e6L, 2.561565266405658882048e7L,\n 6.612718635548492136250e7L, 1.711446713003636684975e8L,\n 4.439663698302712208698e8L, 1.154115391849182948287e9L,\n 3.005950906525548689841e9L, 7.842940991898186370453e9L,\n 2.049649711988081236484e10L, 5.364511859231469415605e10L,\n 1.405991957584069047340e11L, 3.689732094072741970640e11L,\n 9.694555759683939661662e11L, 2.550043566357786926147e12L,\n 6.714640184076497558707e12L, 1.769803724411626854310e13L,\n 4.669055014466159544500e13L, 1.232852079912097685431e14L,\n 3.257988998672263996790e14L, 8.616388199965786544948e14L,\n 2.280446200301902595341e15L, 6.039718263611241578359e15L,\n 1.600664914324504111070e16L, 4.244796092136850759368e16L,\n 1.126348290166966760275e17L, 2.990444718632336675058e17L,\n 7.943916035704453771510e17L, 2.111342388647824195000e18L,\n 5.614329680810343111535e18L, 1.493630213112993142255e19L,\n 3.975442747903744836007e19L, 1.058563689713169096306e20L\n };\n int k = (int) (x + 0.5);\n int j = 0;\n long double xx = (long double) k;\n long double dx = x - xx;\n long double xxj = xx;\n long double edx = expl(dx);\n long double Sm = 1.0L;\n long double Sn = (edx - 1.0L) / xxj;\n long double term = DBL_MAX;\n long double factorial = 1.0L;\n long double dxj = 1.0L;\n\n while (fabsl(term) > epsilon * fabsl(Sn) ) {\n j++;\n factorial *= (long double) j;\n xxj *= xx;\n dxj *= (-dx);\n Sm += (dxj / factorial);\n term = ( factorial * (edx * Sm - 1.0L) ) / xxj;\n Sn += term;\n }\n \n return ei[k-7] + Sn * expl(xx); \n}\n\"\"\"\n\nfilename = tempname()\nfilenamedl = filename * \".\" * Libdl.dlext\n\nopen(`gcc -Wall -O3 -march=native -xc -shared -o $filenamedl -`, \"w\") do f\n print(f, C_code)\nend\n\nrun(`ls -l $filenamedl`)",
"execution_count": 13,
"outputs": [
{
"output_type": "stream",
"text": "-rwxr-xr-x 1 genkuroki 197610 124409 Aug 20 14:43 C:\\Users\\GENKUR~1\\AppData\\Local\\Temp\\jl_A3BE.tmp.dll\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "const Ei_mml_ld_lib = filename\nEi_gcc_mml_ld(x::Float64) = ccall((:Exponential_Integral_Ei, Ei_mml_ld_lib), Float64, (Float64,), x)",
"execution_count": 14,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 14,
"data": {
"text/plain": "Ei_gcc_mml_ld (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@show Ei_gcc_mml_ld(-10.0)\n@show Ei_gcc_mml_ld(-3.0)\n@show Ei_gcc_mml_ld(0.0)\n@show Ei_gcc_mml_ld(5.0)\n@show Ei_gcc_mml_ld(20.0)\n@show Ei_gcc_mml_ld(100.0);",
"execution_count": 15,
"outputs": [
{
"output_type": "stream",
"text": "Ei_gcc_mml_ld(-10.0) = -4.156968929685325e-6\nEi_gcc_mml_ld(-3.0) = -0.013048381094197037\nEi_gcc_mml_ld(0.0) = -1.7976931348623157e308\nEi_gcc_mml_ld(5.0) = 40.18527535580318\nEi_gcc_mml_ld(20.0) = 2.5615652664056588e7\nEi_gcc_mml_ld(100.0) = 2.71555274485388e41\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# double (not long double) version of \n# http://www.mymathlib.com/c_source/functions/exponential_integrals/exponential_integral_Ei.c\n\nC_code = raw\"\"\"\n////////////////////////////////////////////////////////////////////////////////\n// File: exponential_integral_Ei.c //\n// Routine(s): //\n// Exponential_Integral_Ei //\n////////////////////////////////////////////////////////////////////////////////\n\n#include <math.h> // required for fabsl(), expl() and logl() \n#include <float.h> // required for DBL_EPSILON, DBL_MAX\n\n// Internally Defined Routines //\ndouble Exponential_Integral_Ei( double x );\n\nstatic double Continued_Fraction_Ei( double x );\nstatic double Power_Series_Ei( double x );\nstatic double Argument_Addition_Series_Ei( double x);\n\n\n// Internally Defined Constants //\nstatic const double epsilon = 10.0 * DBL_EPSILON;\n\n////////////////////////////////////////////////////////////////////////////////\n// double Exponential_Integral_Ei( double x ) //\n// //\n// Description: //\n// The exponential integral Ei(x) is the integral with integrand //\n// exp(t) / t //\n// where the integral extends from -inf to x. //\n// Note that there is a singularity at t = 0. Therefore for x > 0, the //\n// integral is defined to be the Cauchy principal value: //\n// lim { I[-inf, -eta] exp(-t) dt / t + I[eta, x] exp(-t) dt / t } //\n// in which the limit is taken as eta > 0 approaches 0 and I[a,b] //\n// denotes the integral from a to b. //\n// //\n// Arguments: //\n// double x The argument of the exponential integral Ei(). //\n// //\n// Return Value: //\n// The value of the exponential integral Ei evaluated at x. //\n// If x = 0.0, then Ei is -inf and -DBL_MAX is returned. //\n// //\n// Example: //\n// double y, x; //\n// //\n// ( code to initialize x ) //\n// //\n// y = Exponential_Integral_Ei( x ); //\n////////////////////////////////////////////////////////////////////////////////\n\ndouble Exponential_Integral_Ei( double x )\n{\n if ( x < -5.0L ) return Continued_Fraction_Ei(x);\n if ( x == 0.0L ) return -DBL_MAX;\n if ( x < 6.8L ) return Power_Series_Ei(x);\n if ( x < 50.0L ) return Argument_Addition_Series_Ei(x);\n return Continued_Fraction_Ei(x);\n}\n\n////////////////////////////////////////////////////////////////////////////////\n// static double Continued_Fraction_Ei( double x ) //\n// //\n// Description: //\n// For x < -5 or x > 50, the continued fraction representation of Ei //\n// converges fairly rapidly. //\n// //\n// The continued fraction expansion of Ei(x) is: //\n// Ei(x) = -exp(x) { 1/(-x+1-) 1/(-x+3-) 4/(-x+5-) 9/(-x+7-) ... }. //\n// //\n// //\n// Arguments: //\n// double x //\n// The argument of the exponential integral Ei(). //\n// //\n// Return Value: //\n// The value of the exponential integral Ei evaluated at x. //\n////////////////////////////////////////////////////////////////////////////////\n\nstatic double Continued_Fraction_Ei( double x )\n{\n double Am1 = 1.0;\n double A0 = 0.0;\n double Bm1 = 0.0;\n double B0 = 1.0;\n double a = expl(x);\n double b = -x + 1.0;\n double Ap1 = b * A0 + a * Am1;\n double Bp1 = b * B0 + a * Bm1;\n int j = 1;\n\n a = 1.0;\n while ( fabsl(Ap1 * B0 - A0 * Bp1) > epsilon * fabsl(A0 * Bp1) ) {\n if ( fabsl(Bp1) > 1.0L) {\n Am1 = A0 / Bp1;\n A0 = Ap1 / Bp1;\n Bm1 = B0 / Bp1;\n B0 = 1.0L;\n } else {\n Am1 = A0;\n A0 = Ap1;\n Bm1 = B0;\n B0 = Bp1;\n }\n a = -j * j;\n b += 2.0L;\n Ap1 = b * A0 + a * Am1;\n Bp1 = b * B0 + a * Bm1;\n j += 1;\n }\n return (-Ap1 / Bp1);\n}\n\n\n////////////////////////////////////////////////////////////////////////////////\n// static double Power_Series_Ei( double x ) //\n// //\n// Description: //\n// For -5 < x < 6.8, the power series representation for //\n// (Ei(x) - gamma - ln|x|)/exp(x) is used, where gamma is Euler's gamma //\n// constant. //\n// Note that for x = 0.0, Ei is -inf. In which case -DBL_MAX is //\n// returned. //\n// //\n// The power series expansion of (Ei(x) - gamma - ln|x|) / exp(x) is //\n// - Sum(1 + 1/2 + ... + 1/j) (-x)^j / j!, where the Sum extends //\n// from j = 1 to inf. //\n// //\n// Arguments: //\n// double x //\n// The argument of the exponential integral Ei(). //\n// //\n// Return Value: //\n// The value of the exponential integral Ei evaluated at x. //\n////////////////////////////////////////////////////////////////////////////////\n\nstatic double Power_Series_Ei( double x )\n{ \n double xn = -x;\n double Sn = -x;\n double Sm1 = 0.0;\n double hsum = 1.0;\n double g = 0.5772156649015328606065121;\n double y = 1.0;\n double factorial = 1.0;\n \n if ( x == 0.0 ) return (double) -DBL_MAX;\n \n while ( fabsl(Sn - Sm1) > epsilon * fabsl(Sm1) ) {\n Sm1 = Sn;\n y += 1.0;\n xn *= (-x);\n factorial *= y;\n hsum += (1.0 / y);\n Sn += hsum * xn / factorial;\n }\n return (g + logl(fabsl(x)) - expl(x) * Sn);\n}\n\n\n////////////////////////////////////////////////////////////////////////////////\n// static double Argument_Addition_Series_Ei(double x) //\n// //\n// Description: //\n// For 6.8 < x < 50.0, the argument addition series is used to calculate //\n// Ei. //\n// //\n// The argument addition series for Ei(x) is: //\n// Ei(x+dx) = Ei(x) + exp(x) Sum j! [exp(j) expj(-dx) - 1] / x^(j+1), //\n// where the Sum extends from j = 0 to inf, |x| > |dx| and expj(y) is //\n// the exponential polynomial expj(y) = Sum y^k / k!, the Sum extending //\n// from k = 0 to k = j. //\n// //\n// Arguments: //\n// double x //\n// The argument of the exponential integral Ei(). //\n// //\n// Return Value: //\n// The value of the exponential integral Ei evaluated at x. //\n////////////////////////////////////////////////////////////////////////////////\nstatic double Argument_Addition_Series_Ei(double x)\n{\n static double ei[] = {\n 1.915047433355013959531e2, 4.403798995348382689974e2,\n 1.037878290717089587658e3, 2.492228976241877759138e3,\n 6.071406374098611507965e3, 1.495953266639752885229e4,\n 3.719768849068903560439e4, 9.319251363396537129882e4,\n 2.349558524907683035782e5, 5.955609986708370018502e5,\n 1.516637894042516884433e6, 3.877904330597443502996e6,\n 9.950907251046844760026e6, 2.561565266405658882048e7,\n 6.612718635548492136250e7, 1.711446713003636684975e8,\n 4.439663698302712208698e8, 1.154115391849182948287e9,\n 3.005950906525548689841e9, 7.842940991898186370453e9,\n 2.049649711988081236484e10, 5.364511859231469415605e10,\n 1.405991957584069047340e11, 3.689732094072741970640e11,\n 9.694555759683939661662e11, 2.550043566357786926147e12,\n 6.714640184076497558707e12, 1.769803724411626854310e13,\n 4.669055014466159544500e13, 1.232852079912097685431e14,\n 3.257988998672263996790e14, 8.616388199965786544948e14,\n 2.280446200301902595341e15, 6.039718263611241578359e15,\n 1.600664914324504111070e16, 4.244796092136850759368e16,\n 1.126348290166966760275e17, 2.990444718632336675058e17,\n 7.943916035704453771510e17, 2.111342388647824195000e18,\n 5.614329680810343111535e18, 1.493630213112993142255e19,\n 3.975442747903744836007e19, 1.058563689713169096306e20\n };\n int k = (int) (x + 0.5);\n int j = 0;\n double xx = (double) k;\n double dx = x - xx;\n double xxj = xx;\n double edx = expl(dx);\n double Sm = 1.0;\n double Sn = (edx - 1.0) / xxj;\n double term = DBL_MAX;\n double factorial = 1.0;\n double dxj = 1.0;\n\n while (fabsl(term) > epsilon * fabsl(Sn) ) {\n j++;\n factorial *= (double) j;\n xxj *= xx;\n dxj *= (-dx);\n Sm += (dxj / factorial);\n term = ( factorial * (edx * Sm - 1.0L) ) / xxj;\n Sn += term;\n }\n \n return ei[k-7] + Sn * expl(xx); \n}\n\"\"\"\n\nfilename = tempname()\nfilenamedl = filename * \".\" * Libdl.dlext\n\nopen(`gcc -Wall -O3 -march=native -xc -shared -o $filenamedl -`, \"w\") do f\n print(f, C_code)\nend\n\nrun(`ls -l $filenamedl`)",
"execution_count": 16,
"outputs": [
{
"output_type": "stream",
"text": "-rwxr-xr-x 1 genkuroki 197610 124878 Aug 20 14:43 C:\\Users\\GENKUR~1\\AppData\\Local\\Temp\\jl_A47B.tmp.dll\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "const Ei_mml_lib = filename\nEi_gcc_mml(x::Float64) = ccall((:Exponential_Integral_Ei, Ei_mml_lib), Float64, (Float64,), x)",
"execution_count": 17,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 17,
"data": {
"text/plain": "Ei_gcc_mml (generic function with 1 method)"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@show Ei_gcc_mml(-10.0)\n@show Ei_gcc_mml(-3.0)\n@show Ei_gcc_mml(0.0)\n@show Ei_gcc_mml(5.0)\n@show Ei_gcc_mml(20.0)\n@show Ei_gcc_mml(100.0);",
"execution_count": 18,
"outputs": [
{
"output_type": "stream",
"text": "Ei_gcc_mml(-10.0) = -4.1569689296853255e-6\nEi_gcc_mml(-3.0) = -0.013048381094195725\nEi_gcc_mml(0.0) = -1.7976931348623157e308\nEi_gcc_mml(5.0) = 40.1852753558037\nEi_gcc_mml(20.0) = 2.5615652664056588e7\nEi_gcc_mml(100.0) = 2.7155527448538795e41\n",
"name": "stdout"
}
]
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "using BenchmarkTools",
"execution_count": 19,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "## test data\n\nxx = Array{Float64}(1000)\nxx[1:300] = rand(300)\nxx[301:600] = 10*rand(300)\nxx[601:1000] = 100*rand(400)\nxxx = BigFloat.(xx);",
"execution_count": 20,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "a = Ei_MPFR.(xxx)\n@benchmark Ei_MPFR.(xxx)",
"execution_count": 21,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 21,
"data": {
"text/plain": "BenchmarkTools.Trial: \n memory estimate: 14.07 MiB\n allocs estimate: 239558\n --------------\n minimum time: 111.744 ms (0.00% GC)\n median time: 115.471 ms (0.39% GC)\n mean time: 115.425 ms (0.27% GC)\n maximum time: 119.996 ms (0.47% GC)\n --------------\n samples: 44\n evals/sample: 1"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "b1 = Ei_julia_numres.(xx)\n@benchmark Ei_julia_numres.(xx)",
"execution_count": 22,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 22,
"data": {
"text/plain": "BenchmarkTools.Trial: \n memory estimate: 8.95 KiB\n allocs estimate: 24\n --------------\n minimum time: 243.063 μs (0.00% GC)\n median time: 246.795 μs (0.00% GC)\n mean time: 256.155 μs (0.45% GC)\n maximum time: 3.290 ms (89.70% GC)\n --------------\n samples: 10000\n evals/sample: 1"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "b2 = Ei_gcc_numres.(xx)\n@benchmark Ei_gcc_numres.(xx)",
"execution_count": 23,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 23,
"data": {
"text/plain": "BenchmarkTools.Trial: \n memory estimate: 8.95 KiB\n allocs estimate: 24\n --------------\n minimum time: 224.868 μs (0.00% GC)\n median time: 234.199 μs (0.00% GC)\n mean time: 242.584 μs (0.48% GC)\n maximum time: 3.267 ms (91.82% GC)\n --------------\n samples: 10000\n evals/sample: 1"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "c1 = Ei_julia_mml.(xx)\n@benchmark Ei_julia_mml.(xx)",
"execution_count": 24,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 24,
"data": {
"text/plain": "BenchmarkTools.Trial: \n memory estimate: 126.20 KiB\n allocs estimate: 292\n --------------\n minimum time: 216.470 μs (0.00% GC)\n median time: 224.868 μs (0.00% GC)\n mean time: 237.090 μs (1.83% GC)\n maximum time: 1.300 ms (78.79% GC)\n --------------\n samples: 10000\n evals/sample: 1"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "c2 = Ei_gcc_mml.(xx)\n@benchmark Ei_gcc_mml.(xx)",
"execution_count": 25,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 25,
"data": {
"text/plain": "BenchmarkTools.Trial: \n memory estimate: 8.95 KiB\n allocs estimate: 24\n --------------\n minimum time: 229.533 μs (0.00% GC)\n median time: 241.664 μs (0.00% GC)\n mean time: 253.340 μs (0.48% GC)\n maximum time: 3.361 ms (91.57% GC)\n --------------\n samples: 10000\n evals/sample: 1"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# long double version\nc3 = Ei_gcc_mml_ld.(xx)\n@benchmark Ei_gcc_mml_ld.(xx)",
"execution_count": 26,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 26,
"data": {
"text/plain": "BenchmarkTools.Trial: \n memory estimate: 8.95 KiB\n allocs estimate: 24\n --------------\n minimum time: 297.647 μs (0.00% GC)\n median time: 318.175 μs (0.00% GC)\n mean time: 327.904 μs (0.36% GC)\n maximum time: 3.353 ms (89.36% GC)\n --------------\n samples: 10000\n evals/sample: 1"
},
"metadata": {}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "## relative errors\n\nrelerr(x,y) = maximum(abs.((x-y)./x))\n#relerr(x,y) = mean(abs.((x-y)./x))\n\n@show relerr(a,b1)\n@show relerr(a,b2)\n@show relerr(b1,b2)\n@show relerr(a,c1)\n@show relerr(a,c2)\n@show relerr(a,c3)\n@show relerr(c1,c2)\n@show relerr(c1,c3)\n@show relerr(c2,c3);",
"execution_count": 27,
"outputs": [
{
"output_type": "stream",
"text": "relerr(a, b1) = 1.005585819477779495751362591355570077067842218286070327489762875086146251181161e-14\nrelerr(a, b2) = 9.96278881620033203826728818697194707938857824660676443884115173339962824546806e-15\nrelerr(b1, b2) = 1.0505710210083766e-14\nrelerr(a, c1) = 6.420235995405747815736081076464784038727911952451185886415447765919815420617500e-13\nrelerr(a, c2) = 6.420235995405747815736081076464784038727911952451185886415447765919815420617500e-13\nrelerr(a, c3) = 2.524608923095217767621910416084252002116510232757506109669689148582591931051396e-16\nrelerr(c1, c2) = 7.38682749146511e-15\nrelerr(c1, c3) = 6.422294833842673e-13\nrelerr(c2, c3) = 6.422294833842673e-13\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true,
"collapsed": true
},
"cell_type": "code",
"source": "",
"execution_count": null,
"outputs": []
}
],
"metadata": {
"_draft": {
"nbviewer_url": "https://gist.github.com/dabc2fb61ceb80e40fbc189676754fa2"
},
"gist": {
"id": "dabc2fb61ceb80e40fbc189676754fa2",
"data": {
"description": "Julia/exp and log integral functions/Ei function in Julia - Part 3.ipynb",
"public": true
}
},
"kernelspec": {
"name": "julia-0.6",
"display_name": "Julia 0.6.0",
"language": "julia"
},
"language_info": {
"file_extension": ".jl",
"name": "julia",
"mimetype": "application/julia",
"version": "0.6.0"
},
"toc": {
"threshold": 4,
"number_sections": true,
"toc_cell": false,
"toc_window_display": false,
"toc_section_display": "block",
"sideBar": true,
"navigate_menu": true,
"moveMenuLeft": true,
"widenNotebook": false,
"colors": {
"hover_highlight": "#DAA520",
"selected_highlight": "#FFD700",
"running_highlight": "#FF0000",
"wrapper_background": "#FFFFFF",
"sidebar_border": "#EEEEEE",
"navigate_text": "#333333",
"navigate_num": "#000000"
},
"nav_menu": {
"height": "12px",
"width": "252px"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment