%% $Id: pst-func.pro 689 2023-04-03 06:46:27Z herbert $ %% %% This is file `pst-func.pro', %% %% IMPORTANT NOTICE: %% %% Package `pst-func' %% %% Herbert Voss %% %% This program can be redistributed and/or modified under the terms %% of the LaTeX Project Public License Distributed from CTAN archives %% in directory macros/latex/base/lppl.txt. %% %% DESCRIPTION: %% `pst-func' is a PSTricks package to plot special math functions %% %% %% version 0.19 / 2023-04-02 Herbert Voss % /tx@FuncDict 100 dict def tx@FuncDict begin % /eps1 1.0e-05 def /eps2 1.0e-04 def /eps5 1.0e-05 def /eps8 1.0e-08 def % /PiHalf 1.57079632679489661925640 def /CEuler 0.5772156649 def % Euler-Mascheroni constant % /factorial { % n on stack, returns n! dup 32 gt { pop 1e32 } { dup 0 eq { 1 }{ dup 1 gt { dup 1 sub factorial mul } if } ifelse } ifelse } def % /MoverN { % m n on stack, returns the binomial coefficient m over n 2 dict begin /n exch def /m exch def n 0 eq { 1 }{ m n eq { 1 }{ m factorial n factorial m n sub factorial mul div } ifelse } ifelse end } def % /Pascal [ [ 1 ] % 0 [ 1 1 ] % 1 [ 1 2 1 ] % 2 [ 1 3 3 1 ] % 3 [ 1 4 6 4 1 ] % 4 [ 1 5 10 10 5 1 ] % 5 [ 1 6 15 20 15 6 1 ] % 6 [ 1 7 21 35 35 21 7 1 ] % 7 [ 1 8 28 56 70 56 28 8 1 ] % 8 [ 1 9 36 84 126 126 84 36 9 1 ] % 9 ] def % /GetBezierCoor { % t on stack 10 dict begin % hold all local /t ED /t1 1 t sub def % t1=1-t /Coeff Pascal BezierType get def % get the coefficients 0 0 % initial values for x y BezierType -1 0 { % BezierType,...,2,1,0 /I ED % I=BezierType,...,2,1,0 /J BezierType I sub def % J=0,1,2,...,BezierType /T t I exp Coeff J get mul def % coeff(J)*t^I /T1 t1 J exp def % t1^J Points I dup add 1 add get % y(2*I+1) T mul T1 mul add % the y coordinate exch % y x Points I dup add get % x(2*I) T mul T1 mul add % the x coordinate exch % x y } for % x y on stack end } def /BezierCurve { % on stack [ coors psk@plotpoints BezierType % 10 dict begin /BezierType ED % 2,3,4,5,6,... 1 exch div /epsilon ED % step for Bezier =1/plotpoints ] % [ yi xi ... y3 x3 y2 x2 y1 x1 y0 x0] ps@ReverseOrderOfPoints % [y0 x0 y1 x1 ... yi xi] /Points ED % save Points array epsilon GetBezierCoor % next Bezier point Points 0 get Points 1 get % starting point ArrowA lineto epsilon epsilon 1 epsilon sub { % on stack is the loop variable GetBezierCoor lineto } for 1 epsilon sub GetBezierCoor 1 GetBezierCoor ArrowB lineto moveto % end } def /Bernstein { % on stack tStart tEnd plotpoints i n 12 dict begin % hold all local /envelope ED % plot envelope? /n ED /i ED /ni n i sub def /epsilon ED % step=1/plotpoints /tEnd ED /tStart ED % % B_{i,n}(t)=\binom{n}{i}t^i(1-t)^{n-i} (Bernstein) % f_n(x)=\frac{1}{\sqrt{\pi n\cdot x(1-x)}} (envelope) % n i MoverN /noveri ED % \binom{n}{i} [ % for the array of points tStart epsilon tEnd { dup dup /t ED % leave one on stack neg 1 add /t1 ED % t1=1-t envelope { t t1 mul 4 mul PiHalf mul n mul sqrt 1 exch Div } % envelope { noveri t i exp mul t1 ni exp mul } ifelse % t f(t) ScreenCoor % convert to screen coor } for end false /Lineto /lineto load def Line } def %% /Si { % integral sin from 0 to x (arg on stack) 10 dict begin % hold all local dup 0 eq { pop 0 } { /arg exch def % x /arg2 arg dup mul def /Sum arg def % /sign -1 def /I 3 def /Frac arg2 arg mul 6 div def { % a sequence of x - x^3/(3*3!) + x^5/(5*5!) -...+... Frac I div sign mul dup abs eps5 lt { pop exit } if Sum add /Sum exch def /sign sign neg def /I I 2 add def Frac arg2 mul I 1 sub I mul div /Frac ED % arg I Power dup abs 1e30 gt { pop exit } if % I factorial div I div sign mul % dup abs eps8 lt { pop exit } if % Sum add /Sum exch def % /sign sign neg def % /I I 2 add def } loop Sum } ifelse end } def % /si { % integral sin from x to infty -> si(x)=Si(x)-pi/2 Si PiHalf sub } def % /Ci { % integral cosin from x to infty (arg on stack) 10 dict begin % hold all local abs /arg exch def arg 0 eq { 0 } { /arg2 arg dup mul def /Sum CEuler arg ln add def /sign -1 def /I 2 def /Frac arg2 2 div def % first fraction { Frac I div sign mul dup abs eps5 lt { pop exit } if Sum add /Sum exch def /sign sign neg def /I I 2 add def Frac arg2 mul I 1 sub I mul div /Frac ED } loop Sum } ifelse end } def % /ci { % integral cosin from x to infty -> ci(x)=-Ci(x)+ln(x)+CEuler dup Ci neg exch abs ln add CEuler add } def % /MaxIter 255 def /func { coeff Derivation FuncValue } def /func' { coeff Derivation 1 add FuncValue } def /func'' { coeff Derivation 2 add FuncValue } def % /NewtonMehrfach {% the start value must be on top of the stack /Nx exch def /Iter 0 def { /Iter Iter 1 add def Nx func /F exch def % f(Nx) F abs eps2 lt { exit } if Nx func' /FS exch def % f'(Nx) FS 0 eq { /FS 1.0e-06 def } if Nx func'' /F2S exch def % f''(Nx) 1.0 1.0 F F2S mul FS dup mul div sub div /J exch def J F mul FS div /Diff exch def /Nx Nx Diff sub def Diff abs eps1 lt Iter MaxIter gt or { exit } if } loop Nx % the returned value ist the zero point } def /Steffensen {% the start value must be on top of the stack /y0 exch def % the start value /Iter 0 def /MaxIter 200 def { pstack y0 func /F exch def F abs eps2 lt { exit } if y0 F sub /Phi exch def Phi func /F2 exch def F2 abs eps2 le { exit }{ Phi y0 sub dup mul Phi F2 sub 2 Phi mul sub y0 add Div /Diff exch def y0 Diff sub /y0 exch def Diff abs eps1 le { exit } if } ifelse /Iter Iter 1 add def Iter MaxIter gt { exit } if } loop y0 28 mul % the returned value ist the zero point 0 3 0 360 arc gsave 0 0 1 setrgbcolor fill grestore 1 setlinewidth stroke } def % /Horner {% x [coeff] must be on top of the stack aload length dup 2 add -1 roll exch 1 sub { dup 4 1 roll mul add exch } repeat pop % the y value is on top of the stack } def % /FuncValue {% x [coeff] Derivation must be on top of the stack { aload % a0 a1 a2 ... a(n-1) [array] length % a0 a1 a2 ... a(n-1) n 1 sub /grad exch def % a0 a1 a2 ... a(n-1) grad -1 1 { % for n=grad step -1 until 1 /n exch def % Laufvariable speichern n % a0 a1 a2 ... a(n-1) n mul % a0 a1 a2 ... a(n-1)*n grad 1 add % a0 a1 a2 ... a(n-1)*n grad+1 1 roll % an*na0 a1 a2 ... a(n-2) } for pop % loesche a0 grad array astore % [ a1 a2 ... a(n-2)] } repeat Horner } def % /FindZeros { % dxN dxZ must be on top of the stack (x0..x1 the intervall) => [] 12 dict begin /dxZ exch def /dxN exch def /pstZeros [] def x0 dxZ x1 { % suche Nullstellen /xWert exch def xWert NewtonMehrfach %xWert Steffensen /xNull exch def pstZeros aload length /Laenge exch def % now test if value is a new one Laenge 0 eq { xNull 1 } { /newZero true def Laenge { xNull sub abs dxN lt { /newZero false def } if } repeat pstZeros aload pop newZero { xNull Laenge 1 add } { Laenge } ifelse } ifelse array astore /pstZeros exch def } for pstZeros % the end array is now on the stack end } def % /Simpson { % on stack must be a b M useXVal --- simple version --- % /SFunc must be defined /useX ED % for algebraic functions which uses f(x) /M ED /b ED /a ED /h b a sub M 2 mul div def /s1 0 def /s2 0 def 1 1 M { /simp_k exch def /xVal simp_k 2 mul 1 sub h mul a add def /s1 s1 xVal useX { /x exch def } if SFunc add def } for 1 1 M 1 sub { /simp_k exch def /xVal simp_k 2 mul h mul a add def /s2 s2 xVal useX { /x exch def } if SFunc add def } for /I a useX { /x exch def } if SFunc b useX { /x exch def } if SFunc add s1 4 mul add s2 2 mul add 3 div h mul def } def % % /LogGamma { 5 dict begin % z on stack /z ED /sum 0 def /logGamma_k 1 def { z logGamma_k div dup 1 add ln sub dup abs eps8 lt { pop exit } if sum add /sum exch def /logGamma_k logGamma_k 1 add def } loop sum z ln sub CEuler z mul sub end } def % /ChebyshevT { 5 dict begin % z on stack /xtmp exch def /n exch def 0 0 1 n .5 mul floor { /cheb_k exch def xtmp xtmp mul 1 sub cheb_k exp xtmp n 2 cheb_k mul sub exp mul n 2 cheb_k mul MoverN mul add } for end } def % /ChebyshevU { 5 dict begin % z on stack /xtmp exch def /n exch def 0 0 1 n .5 mul floor { /cheb_k exch def xtmp xtmp mul 1 sub cheb_k exp xtmp n 2 cheb_k mul sub exp mul n 1 add 2 cheb_k mul 1 add MoverN mul add } for end } def % /vasicek{ %density=sqrt((1-R2)/R2)*exp(1/2*(norminv(x)2 - (1/sqrt(R2)*((sqrt(1-R2)*norminv(x)-norminv(pd)))2)) 2 dict begin /pd where { pop }{ /pd 0.22 def } ifelse % element of (0,1) probability of default of portfolio /R2 where { pop }{ /R2 0.11 def } ifelse % element of (0,1) R_Squared of portfolio dup % x x norminv % x norminv(x) dup mul % x norminv(x)^2 exch % norminv(x)2 x norminv % norminv(x)2 norminv(x) 1 R2 sub sqrt mul % norminv(x)2 sqrt(1-R2)*norminv(x) pd norminv sub % norminv(x)2 sqrt(1-R2)*norminv(x)-norminv(pd) R2 sqrt div % norminv(x)2 1/sqrt(R2)*(sqrt(1-R2)*norminv(x)-norminv(pd)) dup mul % norminv(x)2 (1/sqrt(R2)*(sqrt(1-R2)*norminv(x)-norminv(pd)))2 sub % norminv(x)2 -(1/sqrt(R2)*(sqrt(1-R2)*norminv(x)-norminv(pd)))2 2 div % 1/2*(norminv(x)2 -(1/sqrt(R2)*(sqrt(1-R2)*norminv(x)-norminv(pd)))2) ENeperian exch exp % exp(1/2*(norminv(x)2 -(1/sqrt(R2)*(sqrt(1-R2)*norminv(x)-norminv(pd)))2) 1 R2 sub % exp(1/2*(norminv(x)2 -(1/sqrt(R2)*(sqrt(1-R2)*norminv(x)-norminv(pd)))2) 1-R2 R2 div % exp(1/2*(norminv(x)2 -(1/sqrt(R2)*(sqrt(1-R2)*norminv(x)-norminv(pd)))2) (1-R2)/R2 sqrt % exp(1/2*(norminv(x)2 -(1/sqrt(R2)*(sqrt(1-R2)*norminv(x)-norminv(pd)))2) sqrt((1-R2)/R2) mul % sqrt((1-R2)/R2)*exp(1/2*(norminv(x)2 - (1/sqrt(R2)*((sqrt(1-R2)*norminv(x)-norminv(pd)))2)) end } def %end{vasicek density} % % https://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html /ConfHyperFunc { % Confluent Hypergeometric Funtion of the First Kind % on stack must be a b z 10 dict begin /z ED /b ED /a ED /sum 1 def /k 0 def /nom 1 def /denom 1 def { nom a k add mul z mul /nom ED denom b k add mul k 1 add mul /denom ED nom denom div /diff ED sum diff add /sum ED diff 1e-5 lt { exit } { k 1 add /k ED } ifelse k 1000 ge { (Error ConfHyperFunc: k =1000) == exit } if } loop sum end } def % end % /arraySum { % on stack the array 5 dict begin dup length 0 eq {} { /xArray exch def /xSum 0 def /i 0 def xArray length { /xSum xSum xArray i get add def /i i 1 add def } repeat } ifelse xSum end } def % /arrayProd { % on stack the array 5 dict begin dup length 0 eq {} { /xArray exch def /xProd 1 def /i 0 def xArray length { /xProd xProd xArray i get mul def /i i 1 add def } repeat } ifelse xProd end } def % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % subroutines for complex numbers, given as an array [a b] % which is a+bi = Real+i Imag %% Global defined % /cxadd { % [a1 b1] [a2 b2] = [a1+a2 b1+b2] dup 0 get % [a1 b1] [a2 b2] a2 3 -1 roll % [a2 b2] a2 [a1 b1] dup 0 get % [a2 b2] a2 [a1 b1] a1 3 -1 roll % [a2 b2] [a1 b1] a1 a2 add % [a2 b2] [a1 b1] a1+a2 3 1 roll % a1+a2 [a2 b2] [a1 b1] 1 get % a1+a2 [a2 b2] b1 exch 1 get % a1+a2 b1 b2 add 2 array astore } def % /cxneg { % [a b] dup 1 get % [a b] b exch 0 get % b a neg exch neg % -a -b 2 array astore } def % /cxsub { cxneg cxadd } def % same as negative addition % % [a1 b1][a2 b2] = [a1a2-b1b2 a1b2+b1a2] = [a3 b3] /cxmul { % [a1 b1] [a2 b2] dup 0 get % [a1 b1] [a2 b2] a2 exch 1 get % [a1 b1] a2 b2 3 -1 roll % a2 b2 [a1 b1] dup 0 get % a2 b2 [a1 b1] a1 exch 1 get % a2 b2 a1 b1 dup % a2 b2 a1 b1 b1 5 -1 roll dup % b2 a1 b1 b1 a2 a2 3 1 roll mul % b2 a1 b1 a2 b1a2 5 -2 roll dup % b1 a2 b1a2 b2 a1 a1 3 -1 roll dup % b1 a2 b1a2 a1 a1 b2 b2 3 1 roll mul % b1 a2 b1a2 a1 b2 a1b2 4 -1 roll add % b1 a2 a1 b2 b3 4 2 roll mul % b1 b2 b3 a1a2 4 2 roll mul sub % b3 a3 exch 2 array astore } def % % [a b]^2 = [a^2-b^2 2ab] = [a2 b2] /cxsqr { % [a b] square root dup 0 get exch 1 get % a b dup dup mul % a b b^2 3 -1 roll % b b^2 a dup dup mul % b b^2 a a^2 3 -1 roll sub % b a a2 3 1 roll mul 2 mul % a2 b2 2 array astore } def % /cxsqrt { % [a b] % dup cxnorm sqrt /r exch def % cxarg 2 div RadtoDeg dup cos r mul exch sin r mul cxmake2 cxlog % log[a b] 2 cxrdiv % log[a b]/2 aload pop exch % b a 2.781 exch exp % b exp(a) exch cxconv exch % [Re +iIm] exp(a) cxrmul % } def % /cxarg { % [a b] aload pop % a b exch atan % arctan b/a DegtoRad % arg(z)=atan(b/a) } def % % log[a b] = [a^2-b^2 2ab] = [a2 b2] /cxlog { % [a b] dup % [a b][a b] cxnorm % [a b] |z| log % [a b] log|z| exch % log|z|[a b] cxarg % log|z| Theta cxmake2 % [log|z| Theta] } def % % square of magnitude of complex number /cxnorm2 { % [a b] dup 0 get exch 1 get % a b dup mul % a b^2 exch dup mul add % a^2+b^2 } def % /cxnorm { % [a b] cxnorm2 sqrt } def % /cxconj { % conjugent complex dup 0 get exch 1 get % a b neg 2 array astore % [a -b] } def % /cxre { 0 get } def % real value /cxim { 1 get } def % imag value % % 1/[a b] = ([a -b]/(a^2+b^2) /cxrecip { % [a b] dup cxnorm2 exch % n2 [a b] dup 0 get exch 1 get % n2 a b 3 -1 roll % a b n2 dup % a b n2 n2 4 -1 roll exch div % b n2 a/n2 3 1 roll div % a/n2 b/n2 neg 2 array astore } def % /cxmake1 { 0 2 array astore } def % make a complex number, real given /cxmake2 { 2 array astore } def % dito, both given % /cxdiv { cxrecip cxmul } def % % multiplikation by a real number /cxrmul { % [a b] r exch aload pop % r a b 3 -1 roll dup % a b r r 3 1 roll mul % a r b*r 3 1 roll mul % b*r a*r exch 2 array astore % [a*r b*r] } def % % division by a real number /cxrdiv { % [a b] r 1 exch div % [a b] 1/r cxrmul } def % % exp(i theta) = cos(theta)+i sin(theta) polar<->cartesian /cxconv { % theta RadtoDeg dup sin exch cos cxmake2 } def % % cxexp z^k with k as a natural number /cxexp { % z k 3 dict begin dup 0 eq { pop pop [1 0] }{ /k ED /z ED /sol [1 0] def k { sol z cxmul /sol ED } repeat sol } ifelse end } def %