<gh_stars>10-100
#!/usr/bin/python

import numpy as np
import scipy.special as spc
import matplotlib.pyplot as plt

'''

Taylor erf(x)/x
C = 2/sqrt(pi)
C +
- x^2/3  
+ x^4/5  
- x^6/21  
+ x^8/108  
- x^10/660  
+ x^12/4680  
- x^14/37800  
+ x^16/342720  
- x^18/3447360  
+ x^20/38102400  
- x^22/459043200  
+ x^24/5987520000  
(Taylor series)


Taylor x/erf(x)
C = sqrt(pi)/2 
C + 
+ x^2  1/6 
+ x^4  1/180 
- x^6  11/3780 
+ x^8  29/226800 
+ x^10 547/7484400 
- x^12 216341/20432412000 
- x^14 79979/61297236000 
+ x^16 56151197/125046361440000 
+ x^18 3251221/9978699642912000 
- x^20 232514321849/16464854410804800000 
+ x^22 110360356871/87390381103502400000 
+ x^24 185940275380001/563906295520600032000000 
+ O(x^25)

Interesting Functions

exp(-x^2)
(exp(exp(-x^2))-1)/(e-1)
(exp(-exp(x^2)))*e
exp(-1/(1-x^2))*e   # Bump Function https://en.wikipedia.org/wiki/Bump_function

exp(-abs(x))
exp(-exp(abs(x)))*e

Approximate Gaussian like this
1-x^2(1.0-0.43*x^2+0.072*x^4)

'''

def erfx_10( x ):
    xx = x*x
    #even =   0.9850156203557686  +xx*(-0.02758911857985105  +xx*(-0.00125751573583003  +xx*(0.0005522596605522894  +xx*(-4.933984495902465e-05 +xx*(2.886081672787144e-06 ) ) ) ) ) 
    #odd  =  -0.13893174458211105 +xx*(-0.007637751932154861 +xx*(0.0006635698078029163 +xx*(-3.030234586967172e-05 +xx*(-5.95301392208869e-08  +xx*(-3.048663627380765e-07 ) ) ) ) )
    even =   0.9850156202961753  +xx*(-0.02756061032579559  +xx*(-0.00188409579491924  +xx*(-0.003098629936170076 +xx*(-0.001348858853909826  +xx*(-3.98946569988845e-05 ) ) ) ) ) 
    odd  =  -0.13893350387140332 +xx*(-0.007664292475021448 +xx*( 0.003046826535877866 +xx*( 0.002879338499080343 +xx*( 0.0003260490382458129 +xx*( 1.97093650414204e-06 ) ) ) ) )
    t = even + x*odd
    t*=t # ^2
    t*=t # ^4
    t*=t # ^8
    y = 1./( t + x )
    # === x>5.5
    mask = x>4.5
    y[mask]=1./x[mask]
    return y

def erfx_8( x ):
    xx = x*x
    # === x<1
    #x4*xx*xx
    #even =  1.1283791670955126   +xx*(  0.11283791670955126  +xx*(  0.005223977625442188  +xx*(  0.00012055332981789664 +xx*(  1.6462114365889246e-06 +xx*(  1.4807192815879218e-08 );
    #odd  = -0.37612638903183754  +xx*( -0.026866170645131252 +xx*( -0.0008548327023450852 +xx*( -1.492565035840625e-05  +xx*( -1.6365844691234924e-07 +xx*( -1.2290555301717926e-09 );
    #y = even + xx*odd
    #y = 1.1283791388189353 +xx*(-0.3761241556350163 +xx*(0.1128095820183334 +xx*(-0.02673482602072087 +xx*(0.004946345338720848 +xx*(-0.0005754051770340968 ) ) ) ) ) 
    y = 1.1283791662308296 +xx*(-0.3761262972953429 +xx*(0.1128363404233098 +xx*(-0.02685603827999912 +xx*(0.005192885862299865 +xx*(-0.0008053004722300972 +xx*(8.004020068129447e-05 ) ) ) ) ) )
    # === x>1
    #even =   0.98067267667069    +xx*(-0.08183565358823473 +xx*(-0.05303116762372068 +xx*(-0.004743027700272434  ) ) ) 
    #odd  =  -0.11542445820973404 +xx*( 0.06129092918652679 +xx*( 0.02344302411223819 +xx*( 0.0003677467102954218 ) ) ) 
    #even =  0.9878239907607781  +xx*(0.03719497985018601 +xx*(0.1415245895895332  +xx*(0.06052603691072229  +xx*(0.004745181744943059   +xx*(5.244735766994635e-05 ) ) ) ) ) 
    #odd  = -0.15922365218901394 +xx*(-0.1277842231558644 +xx*(-0.1122972554364347 +xx*(-0.02115819731384116 +xx*(-0.0006620172597403533 +xx*(-1.80492148568745e-06 ) ) ) ) )
    even =   0.9903386741213333  +xx*( 0.08180278811069948 +xx*( 0.219787883285348  +xx*( 0.0893543139653664  +xx*( 0.0071698531450102   +xx*( 8.644883946761633e-05 ) ) ) ) ) 
    odd  =  -0.17511814497584813 +xx*(-0.2010794452848663  +xx*(-0.1692686167813105 +xx*(-0.03129254573733003 +xx*(-0.001037968593234627 +xx*(-3.164137211658646e-06 ) ) ) ) ) 
    t = even + x*odd
    t*=t # ^2
    t*=t # ^4
    t*=t # ^8
    mask = x>1.
    y[mask] = 1./( t + x )[mask]
    # === x>5.5
    mask = x>4.5
    y[mask]=1./x[mask]
    return y

def erfx(x,n=12):
    C = 2./np.sqrt(np.pi)
    coefs = [
        +C,
        -C/3,  
        +C/(5*2),  
        -C/(21*2),  
        +C/(108*2),  
        -C/(660*2),  
        +C/(4680*2),  
        -C/(37800*2),  
        +C/(342720*2),  
        -C/(3447360*2),  
        +C/(38102400*2),  
        -C/(459043200*2),  
        +C/(5987520000*2)
    ]
    print coefs
    xx  = x*x
    xxn = xx.copy()
    y = coefs[0] + xx*0
    for c in coefs[1:n]:
        y   += c * xxn
        xxn *= xx
    return y  
    #return np.polyval(coefs[0:n:-1],x**2)


# --- erf(x)/x using only even power (r^2 = norm2() ... now need for sqrt() )
# x from [0.0 ... 2.0]
# 1.1157459916800467 +xx*(-0.300015116925183 +xx*(0.03810680424304402 ) ) 
# 1.1262494925748456 +xx*(-0.3551787973499269 +xx*(0.07944958267393389 +xx*(-0.007572833000122792 ) ) ) 
# 1.1280654295948946 +xx*(-0.3715352604737077 +xx*(0.1019253904806975 +xx*(-0.0173041498495804 +xx*(0.001302114871465865 ) ) ) )
# 1.1283383276699155 +xx*(-0.3752922093037583 +xx*(0.1100609758703164 +xx*(-0.02340094259354313 +xx*(0.003151285196038148 +xx*(-0.0001950095184811707 ) ) ) ) )
# 1.1283786692048163 +xx*(-0.3761089561355552 +xx*(0.1127378973259718 +xx*(-0.02664608686425424 +xx*(0.004982434860751369 +xx*(-0.0007062918753425185 +xx*(6.654103432737757e-05 +xx*(-3.024474296897672e-06 ) ) ) ) ) ) ) 
# 1.1283791196079755 +xx*(-0.3761243083255276 +xx*(0.1128229500233847 +xx*(-0.02682457598467265 +xx*(0.005165555881937873 +xx*(-0.0008079368750211671 +xx*(9.769929458012389e-05 +xx*(-7.984652771242384e-06 +xx*(3.200439628869844e-07 ) ) ) ) ) ) ) )
# ---- 
# x from [2.0 ... 4.0]
# 0.6585801071176838 +xx*(-0.05043913967216926 +xx*(0.001598154897464572 ) ) 
# 0.7650646717484043 +xx*(-0.08915317524088477 +xx*(0.005858085149332657 +xx*(-0.0001444584146456651 ) ) )
# 0.8535159884584169 +xx*(-0.1324706729650264 +xx*(0.01326519653646158 +xx*(-0.0006726625793671091 +xx*(1.337289131772557e-05 ) ) ) ) 
# 0.9266007111474646 +xx*(-0.1774442808198486 +xx*(0.02373217835463652 +xx*(-0.001828598530278551 +xx*(7.425326687938126e-05 +xx*(-1.229824158621717e-06 ) ) ) ) ) 
# 1.0312476394960062 +xx*(-0.260692477348949 +xx*(0.05073748368315148 +xx*(-0.006465093541651888 +xx*(0.0005301118907990944 +xx*(-2.694648346015247e-05 +xx*(7.719652935084254e-07 +xx*(-9.51949258309881e-09 ) ) ) ) ) ) )
# 1.0651711067725838 +xx*(-0.29431909718531 +xx*(0.0648373478287063 +xx*(-0.009733992287520466 +xx*(0.000988982142813379 +xx*(-6.694610752440582e-05 +xx*(2.890141347286707e-06 +xx*(-7.193459185608764e-08 +xx*(7.849719763086682e-10 ) ) ) ) ) ) ) )


def erf_4(x):
    p = 1. + x*( 0.278393 + x*( 0.230389 + x*(0.000972 + x*0.078108 )))
    p=p*p; p=p*p;
    return 1. - 1./p;

def erf_6(x):
    p = 1. + x*( 0.0705230784 + x*( 0.0422820123 + x*( 0.0092705272 + x*( 0.0001520143 + x*( 0.0002765672 + x*0.0000430638 )))))
    p=p*p; p=p*p; p=p*p; p=p*p;
    return 1. - 1./p;

if __name__ == "__main__":

    #xs = np.arange( -1.0, 6.0, 0.05 )
    xs = np.arange( -1.0, 8.0, 0.05 )
    x = xs.copy()
    #n = 10
    #colors = plt.cm.jet(np.linspace(0.,1.,n+1))

    y_ref = spc.erf( x )/x
    #y4    = erf_4  ( x )/x
    #y6    = erf_6  ( x )/x
    y0    = 1/x
    yn  = erfx(x)
    y8  = erfx_8 (x)
    y10 = erfx_10(x)

    plt.plot( x, y_ref,'k',  lw=3, label="ref")
    #plt.plot( x, y4   ,':', label="y4")
    #plt.plot( x, y6   ,':', label="y6")
    plt.plot( x, y0   ,':', label="y0")
    plt.plot( x, yn   ,':', label="yn")
    plt.plot( x, y8   ,':', label="y8")
    plt.plot( x, y10   ,':', label="y10")
    
    #plt.plot( x, abs(y4-y_ref), '--',label="err4" )
    #plt.plot( x, abs(y6-y_ref), '--',label="err6" )
    plt.plot( x, abs(y0-y_ref), '--',label="err0" )
    plt.plot( x, abs(yn-y_ref), '--',label="errN" )
    plt.plot( x, abs(y8 -y_ref), '--',label="err8" )
    plt.plot( x, abs(y10-y_ref), '--',label="err10" )

    plt.yscale('log')
    plt.ylim(1e-16,1.5)
    plt.legend()
    plt.grid()
    plt.show()






