Browse Source
drop ngspice internal implementation of erfc()
drop ngspice internal implementation of erfc()
which these days is guaranteed to be provided by <math.h>
note,
our own implementation was incorrect anyway.
it evaluated to
erfc_ngspice(x) = erfc(fabs(x))
pre-master-46
rlar
9 years ago
8 changed files with 1 additions and 187 deletions
-
4src/include/ngspice/missing_math.h
-
3src/maths/misc/Makefile.am
-
70src/maths/misc/erfc.c
-
105src/maths/misc/test_erfc.c
-
1visualc/sharedspice.vcxproj
-
3visualc/src/include/ngspice/config.h
-
1visualc/vngspice-fftw.vcxproj
-
1visualc/vngspice.vcxproj
@ -1,70 +0,0 @@ |
|||||
/********** |
|
||||
Copyright 1991 Regents of the University of California. All rights reserved. |
|
||||
Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group |
|
||||
**********/ |
|
||||
|
|
||||
#include "ngspice/ngspice.h" |
|
||||
|
|
||||
#ifndef HAVE_ERFC |
|
||||
|
|
||||
/* erfc computes the erfc(x) the code is from sedan's derfc.f */ |
|
||||
double erfc (double x) |
|
||||
{ |
|
||||
double sqrtPi, n, temp1, xSq, sum1, sum2; |
|
||||
sqrtPi = sqrt( M_PI ); |
|
||||
x = ABS( x ); |
|
||||
n = 1.0; |
|
||||
xSq = 2.0 * x * x; |
|
||||
sum1 = 0.0; |
|
||||
|
|
||||
if ( x > 3.23 ) { |
|
||||
/* asymptotic expansion */ |
|
||||
temp1 = exp( - x * x ) / ( sqrtPi * x ); |
|
||||
sum2 = temp1; |
|
||||
|
|
||||
while ( sum1 != sum2 ) { |
|
||||
sum1 = sum2; |
|
||||
temp1 = -1.0 * ( temp1 / xSq ); |
|
||||
sum2 += temp1; |
|
||||
n += 2.0; |
|
||||
} |
|
||||
return( sum2 ); |
|
||||
} |
|
||||
else { |
|
||||
/* series expansion for small x */ |
|
||||
temp1 = ( 2.0 / sqrtPi ) * exp( - x * x ) * x; |
|
||||
sum2 = temp1; |
|
||||
while ( sum1 != sum2 ) { |
|
||||
n += 2.0; |
|
||||
sum1 = sum2; |
|
||||
temp1 *= xSq / n; |
|
||||
sum2 += temp1; |
|
||||
} |
|
||||
return( 1.0 - sum2 ); |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
/* From C. Hastings, Jr., Approximations for digital computers, |
|
||||
Princeton Univ. Press, 1955. |
|
||||
Approximation accurate to within 1.5E-7 |
|
||||
(making some assumptions about your machine's floating point mechanism) |
|
||||
*/ |
|
||||
double |
|
||||
ierfc(double x) |
|
||||
{ |
|
||||
double t, z; |
|
||||
|
|
||||
t = 1/(1 + 0.3275911*x); |
|
||||
z = 1.061405429; |
|
||||
z = -1.453152027 + t * z; |
|
||||
z = 1.421413741 + t * z; |
|
||||
z = -0.284496736 + t * z; |
|
||||
z = 0.254829592 + t * z; |
|
||||
z = exp(-x*x) * t * z; |
|
||||
|
|
||||
return(z); |
|
||||
} |
|
||||
|
|
||||
#else |
|
||||
int Dummy_Symbol_5; |
|
||||
#endif |
|
||||
@ -1,105 +0,0 @@ |
|||||
/* Paolo Nenzi 2002 - This program tests erfc function |
|
||||
* implementations. |
|
||||
*/ |
|
||||
|
|
||||
|
|
||||
/* |
|
||||
* The situation on erfc functions in spice/cider: |
|
||||
* |
|
||||
* First we have the ierfc in spice, a sort of interpolation, which is |
|
||||
* fast to compute but gives is not so "good" |
|
||||
* Then we have derfc from cider, which is accurate but slow, the code is from sedan's derfc.f . |
|
||||
* Both above are only valid for x > 0.0 |
|
||||
* Then we have glibc/os specific implementation. |
|
||||
* |
|
||||
* Proposal: |
|
||||
* |
|
||||
* Use glibc/os specific implementation as default and then test cider one. |
|
||||
*/ |
|
||||
|
|
||||
|
|
||||
#include <stdio.h> |
|
||||
#include <stdlib.h> |
|
||||
#include <math.h> |
|
||||
#ifdef HAVE_FPU_CTRL |
|
||||
#include <fpu_control.h> |
|
||||
#endif |
|
||||
|
|
||||
|
|
||||
double |
|
||||
derfc(double x) |
|
||||
{ |
|
||||
double sqrtPi, n, temp1, xSq, sum1, sum2; |
|
||||
sqrtPi = sqrt( M_PI ); |
|
||||
x = fabs( x ); /* only x > 0 interested */ |
|
||||
n = 1.0; |
|
||||
xSq = 2.0 * x * x; |
|
||||
sum1 = 0.0; |
|
||||
|
|
||||
if ( x > 3.23 ) { |
|
||||
/* asymptotic expansion */ |
|
||||
temp1 = exp( - x * x ) / ( sqrtPi * x ); |
|
||||
sum2 = temp1; |
|
||||
|
|
||||
while ( sum1 != sum2 ) { |
|
||||
sum1 = sum2; |
|
||||
temp1 = -1.0 * ( temp1 / xSq ); |
|
||||
sum2 += temp1; |
|
||||
n += 2.0; |
|
||||
} |
|
||||
return( sum2 ); |
|
||||
} |
|
||||
else { |
|
||||
/* series expansion for small x */ |
|
||||
temp1 = ( 2.0 / sqrtPi ) * exp( - x * x ) * x; |
|
||||
sum2 = temp1; |
|
||||
while ( sum1 != sum2 ) { |
|
||||
n += 2.0; |
|
||||
sum1 = sum2; |
|
||||
temp1 *= xSq / n; |
|
||||
sum2 += temp1; |
|
||||
} |
|
||||
return( 1.0 - sum2 ); |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
double |
|
||||
ierfc(double x) |
|
||||
{ |
|
||||
double t, z; |
|
||||
|
|
||||
t = 1/(1 + 0.3275911*x); |
|
||||
z = 1.061405429; |
|
||||
z = -1.453152027 + t * z; |
|
||||
z = 1.421413741 + t * z; |
|
||||
z = -0.284496736 + t * z; |
|
||||
z = 0.254829592 + t * z; |
|
||||
z = exp(-x*x) * t * z; |
|
||||
|
|
||||
return(z); |
|
||||
} |
|
||||
|
|
||||
|
|
||||
|
|
||||
int main (void) |
|
||||
{ |
|
||||
double x = -30.0; |
|
||||
double y1= 0.0, y2 = 0.0; |
|
||||
|
|
||||
#ifdef HAVE_FPU_CTRL |
|
||||
fpu_control_t prec; |
|
||||
_FPU_GETCW(prec); |
|
||||
prec &= ~_FPU_EXTENDED; |
|
||||
prec |= _FPU_DOUBLE; |
|
||||
_FPU_SETCW(prec); |
|
||||
#endif |
|
||||
|
|
||||
for (;(x <= 30.0);) |
|
||||
{ |
|
||||
y1 = ierfc(x); |
|
||||
y2 = derfc(x); |
|
||||
printf("x: %f \t ierfc: %e \t derfc: %e \t erfc: %e\n", x, y1, y2, erfc(x) ); |
|
||||
x = x + 1.0; |
|
||||
} |
|
||||
exit(1); |
|
||||
} |
|
||||
Write
Preview
Loading…
Cancel
Save
Reference in new issue