%% $Id: pst-func.pro 694 2023-04-05 19:17:32Z herbert $
%%
%% This is file `pst-func.pro',
%%
%% IMPORTANT NOTICE:
%%
%% Package `pst-func'
%%
%% Herbert Voss <[email protected]>
%%
%% 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.20 / 2023-04-04  Herbert Voss
%
/tx@FuncDict 100 dict def
tx@FuncDict begin
%
/eps4 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
%
/factorialR {
 dup 1e32 gt { pop 1e32 } {
   dup 0 eq     % check for the argument being 0
     { pop 1 }      % if so, the result is 1
     { dup 1 sub factorial       % call recursively with n - 1
       mul        % multiply the result with n
     } ifelse
   } ifelse
} def
%
%Iterative
/factorial {
 1            % initial value for the product
 1 1          % for's start value and increment
 4 -1 roll    % bring the argument to the top as for's end value
 { mul } for
} 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 eps4 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 eps5 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 eps4 lt { exit } if
   y0 F sub /Phi exch def
   Phi func /F2 exch def
   F2 abs eps4 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 eps5 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
   15 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
%
% https://mathworld.wolfram.com/BetaFunction.html
/BETA {% BETA function  ((p-1)!(q-1)!)/(p+q-1)!
      % on stack must be p q
      2 dict begin
      /q ED
      /p ED
      p 1 sub factorial q 1 sub factorial
      mul
      p q add 1 sub factorial
      div
      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
%

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.