Re: [math-fun] floating point precision and hyperbolic functions
I forgot that math-fun strips out attachments. So below is what I tried to attach in my last email, possibly with some mailer-inserted line breaks. log1p and expm1 are ln1+ and exp-1 below. (In PostScript natural log is ln and exp takes two args--a base and exponent, adding an extra source of error since e can't be expressed exactly; the PostScript trig functions work in degrees.) Notes: log1p(x) = log(1+x) + (x-((1+x)-1))/(1+x) For |x| <=1, we use asinh(x) = (log(1+x/sqrt(x^2+1)) - log(1-x/sqrt(x^2+1)))/2, which is based on atanh(x) = (log(1+x) - log(1-x))/2; while for larger x, we use asinh(x) = sgn(x) log(|x| + sqrt(x^2+1)) PostScript uses single precision floating point; any McClaurin series are truncated accordingly. /sgn {dup 0 ne {dup dup eq {0 gt {1}{-1} ifelse} if} if} bind def /e 2.718281828459045 def /ln2 .6931471805599453 def /log2 {dup ln //ln2 div dup sgn exch abs floor mul % x roundtozero(log2?(x)) exch 2 2 index exp div ln //ln2 div add } bind readonly def /exp-1 {dup abs //ln2 le {dup 9. div 1. add 1 index mul 8. div 1. add 1 index mul 7. div 1. add 1 index mul 6. div 1. add 1 index mul 5. div 1. add 1 index mul 4. div 1. add 1 index mul 3. div 1. add 1 index mul 2. div 1. add mul} {//e exch exp 1 sub} ifelse} bind readonly def /ln1+ {dup 1. add dup ln 3 1 roll dup 3 1 roll 1. sub sub exch div add} bind def /Infinity 2 1000 exp def /-Infinity Infinity neg def /NaN Infinity dup sub def /pi 3.141592653589793 def /asin {dup dup mul 1. exch sub sqrt atan} bind readonly def /acos {dup dup mul 1. exch sub sqrt exch atan} bind readonly def /tan {dup sin exch cos div} bind readonly def /ctn {dup cos exch sin div} bind readonly def /sec {1. exch cos div} bind readonly def /csc {1. exch sin div} bind readonly def /cbrt {dup 0. ne {dup abs .333333333333333333 exp 1 index 0. lt {neg} if dup dup dup add 4 1 roll mul div add 3 div} if} bind readonly def /sinh { dup abs //ln2 le % use power series for better precision { dup dup mul dup 72. div 1. add 1 index mul 42. div 1. add 1 index mul 20. div 1. add mul 6. div 1. add mul } { //e 1 index exp exch //e exch neg exp sub .5 mul } ifelse} bind readonly def /cosh {//e 1 index exp exch e exch neg exp add .5 mul} bind readonly def /cosh-1 {dup abs //ln2 le % use power series for better precision { dup mul dup 56. div 1. add 1 index mul 30. div 1. add 1 index mul 12. div 1. add mul 2. div } {//e 1 index exp exch e exch neg exp add .5 mul 1. sub} ifelse} bind readonly def /tanh {dup sinh exch cosh dup //Infinity lt {div} {pop sgn} ifelse} bind readonly def /asinh { dup dup mul dup 1. le { 1. add sqrt div dup ln1+ exch neg ln1+ sub .5 mul } { dup //Infinity lt { 1. add sqrt 1 index abs add ln } % s^2 finite { pop dup abs ln //ln2 add } ifelse exch 0. lt {neg} if } ifelse } bind readonly def /acosh { dup 1. ge { dup dup mul dup //Infinity lt { 1. sub sqrt add ln} { pop ln //ln2 add } ifelse } {/acosh cvx errordict /rangecheck get exec} ifelse } bind readonly def /acosh1+ { dup 0. ge { dup dup 2. add mul dup //Infinity lt { sqrt add ln1+ } { pop ln1+ //ln2 add } ifelse } {/acosh1+ cvx errordict /rangecheck get exec} ifelse} bind readonly def /atanh { dup abs dup 1. lt { pop dup ln1+ exch neg ln1+ sub .5 mul } { 1. eq { //Infinity mul } % +-infinity { /atanh cvx errordict /rangecheck get exec } % error ifelse } ifelse} bind readonly def
participants (1)
-
Michael Speciner