Last active
August 20, 2017 05:48
Star
You must be signed in to star a gist
Julia/exp and log integral functions/Ei function in Julia - Part 3.ipynb
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"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