blob: d68b2e17d6ae9e32d7945f94458633ec2ccb92c3 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
|
% -*- Mode: Postscript -*-
% pst-math.pro --- PostScript header file pst-math.pro
%
% Author : Christophe JORSSEN
% Author : Herbert Voß <hvoss@tug.org>
% Last Mod : $Date: 2018/12/16 $
% Version : 0.65 $
%
/PI 3.14159265359 def
/ENeperian 2.71828182846 def
%
/DegToRad {PI mul 180 div} bind def
/RadToDeg {180 mul PI div} bind def
%
/COS {RadToDeg cos} bind def
/SIN {RadToDeg sin} bind def
/TAN {dup SIN exch COS Div} bind def
/tan {dup sin exch cos Div} bind def
/ATAN {neg -1 atan 180 sub DegToRad} bind def
/ACOS {dup dup mul neg 1 add sqrt exch atan DegToRad} bind def
/acos {dup dup mul neg 1 add sqrt exch atan} bind def
/ASEC {1 exch Div ACOS} bind def
/ASIN {neg dup dup mul neg 1 add sqrt neg atan 180 sub DegToRad} bind def
/asin {neg dup dup mul neg 1 add sqrt neg atan 180 sub} bind def
/ACSC {1 exch Div ASIN} bind def
%
/EXP {ENeperian exch exp} bind def
%
/COSH {dup EXP exch neg EXP add 2 div} bind def
/SINH {dup EXP exch neg EXP sub 2 div} bind def
/TANH {dup SINH exch COSH div} bind def
/ACOSH {dup dup mul 1 sub sqrt add ln} bind def
/ASINH {dup dup mul 1 add sqrt add ln} bind def
/ATANH {dup 1 add exch neg 1 add Div ln 2 div} bind def
%
%/SINC {dup SIN exch Div} bind def
/SINC { dup 0 eq { pop 1 } { dup SIN exch div } ifelse } bind def
/GAUSS {dup mul 2 mul dup 4 -2 roll sub dup mul exch div neg EXP exch PI mul sqrt div} bind def
%
/GAMMA { 2 dict begin % hv 2007-08-30
/z exch def
1.000000000190015 % p(0)
0 1 5 { % on stack is 0 1 2 3 4 5
dup % n-1 n-1
[ 76.18009172947146
-86.50532032941677
24.0140982483091
-1.231739572450155
0.1208650973866179E-2
-0.5395239384953E-5 ] exch get exch % p(n) n-1
1 add z add div % p(n)/(z+n)
add % build the sum
} for
Pi 2 mul sqrt z div mul
z 5.5 add z 0.5 add Power mul ENeperian z 5.5 add neg Power mul % Power is controlled exp
end } bind def
%
/GAMMALN {dup dup dup 5.5 add dup ln 3 -1 roll .5 add mul sub neg 1.000000000190015
0 1 5 {
[76.18009172947146 -86.50532032941677 24.0140982483091 -1.231739572450155
.1208650973866179E-2 -.5395239384953E-5 2.5066282746310005] exch get
4 -1 roll 1 add dup 5 1 roll div add} for
4 -1 roll div 2.5066282746310005 mul ln add exch pop} bind def
/BETA {2 copy add GAMMALN neg exch GAMMALN 3 -1 roll GAMMALN EXP} bind def
%
/HORNER {aload length
dup 2 add -1 roll
exch 1 sub {
dup 4 1 roll
mul add exch
} repeat
pop
} bind def
%
/BESSEL_J0 {dup abs 8 lt {
dup mul dup [57568490574 -13362590354 651619640.7 -11214424.18 77392.33017 -184.9052456] HORNER
exch [57568490411 1029532985 9494680.718 59272.64853 267.8532712 1] HORNER
Div}
{abs dup .636619772 exch div sqrt exch dup .785398164 sub exch 8 exch div dup dup mul dup
[1 -1.098628627E-2 .2734510407E-4 -.2073370639E-5 .2093887211E-6] HORNER
3 index COS mul
exch [-.1562499995E-1 .1430488765E-3 -.6911147651E-5 .7621095161E-6 -.934945152E-7] HORNER
4 -1 roll SIN mul 3 -1 roll mul neg add mul}
ifelse} bind def
%
/BESSEL_Y0 {dup 8 lt {
dup dup mul dup [-2957821389 7062834065 -512359803.6 10879881.29 -86327.92757 228.4622733] HORNER
exch [40076544269 745249964.8 7189466.438 47447.26470 226.1030244 1] HORNER
Div exch dup ln exch BESSEL_J0 .636619772 mul mul add}
{dup .636619772 exch div sqrt exch dup .785398164 sub exch 8 exch div dup dup mul dup
[1 -.1098628627E-2 .2734510407E-4 -.2073370639E-5 .2093887211E-6] HORNER
3 index SIN mul
exch [-.1562499995E-1 .1430488765E-3 -.6911147651E-5 .7621095161E-6 -.934945152E-7] HORNER
4 -1 roll COS mul 3 -1 roll mul add mul}
ifelse} bind def
%
/BESSEL_J1 {dup abs 8 lt {
dup dup mul dup 3 -2 roll [72362614232 -7895059235 242396853.1 -2972611.439 15704.48260 -30.16036606] HORNER mul
exch [144725228442 2300535178 18583304.74 99447.43394 376.9991397 1] HORNER
Div}
{dup abs dup .636619772 exch div sqrt exch dup 2.356194491 sub exch 8 exch div dup dup mul dup
[1 .183105E-2 -.3516396496E-4 .2457520174E-5 -.240337019E-6] HORNER
3 index COS mul
exch [.04687499995 6.2002690873E-3 .8449199096E-5 -.88228987E-6 .105787412E-6] HORNER
4 -1 roll SIN mul 3 -1 roll mul neg add mul exch dup abs Div mul}
ifelse} bind def
%
/BESSEL_Y1 {dup 8 lt {
dup dup dup mul dup [-.4900604943E13 .1275274390E13 -.5153428139E11 .7349264551E9 -.4237922726E7 .8511937935E4] HORNER
exch [.2499580570E14 .4244419664E12 .3733650367E10 .2245904002E8 .1020426050E6 .3549632885E3 1] HORNER
Div mul exch dup dup ln exch BESSEL_J1 mul exch 1 exch div sub .636619772 mul add}
{dup .636619772 exch div sqrt exch dup 2.356194491 sub exch 8 exch div dup dup mul dup
[1 .183105E-2 -.3516396496E-4 .2457520174E-5 -.240337019E-6] HORNER
3 index SIN mul
exch [.04687499995 -.2002690873E-3 .8449199096E-5 6.88228987E-6 .105787412E-6] HORNER
4 -1 roll COS mul 3 -1 roll mul add mul}
ifelse} bind def
%
% En cours...
/BESSEL_Yn {dup 0 eq {pop BESSEL_Y0}{dup 1 eq {pop BESSEL_Y1}{
exch dup BESSEL_Y0 exch dup BESSEL_Y1 exch 2 exch Div {
mul 3 -1 roll mul 2 index sub pstack} for
} ifelse } ifelse } bind def
%
/SIMPSON { 1 dict begin %% on stack a b var f ierr Dominik Rodriguez
3 index 5 index sub % compute h
1 % a b var f ierr h n
4 index 7 index def 3 index exec % a b var f ierr h n f(a)
5 index 7 index def 4 index exec add % a b var f ierr h n f(a)+f(b)
5 index 8 index 4 index 2 div add def 4 index exec % a b var f ierr h n f(a)+f(b) f(a+h/2)
exch 1 index 4 mul add 0 % a b var f ierr h n old=f(a+h/2) Estim=f(a)+f(b)+4f(a+h/2) NbLoop
{ % a b var f ierr h n old Estim NbLoop
5 -1 roll 2 div dup 6 1 roll % h<-h/2
5 -1 roll 2 mul 5 1 roll % n<-2n
% a b var f ierr h n old Estim NbLoop h
2 div 10 index add 0 % a b var f ierr h n old Estim NbLoop a+h/2 Cumul
5 index {
1 index 10 index exch def 8 index exec add exch 6 index add exch
} repeat % a b var f ierr h n old Estim NbLoop a+nh/2 Cumul
exch pop % a b var f ierr h n old Estim NbLoop New
2 index 1 index 4 mul 6 -1 roll 2 mul sub sub % a b var f ierr h n Estim NbLoop New Diff
4 -1 roll 2 mul 1 index sub 4 1 roll % a b var f ierr h n Estim NbLoop New Diff
exch 4 1 roll % a b var f ierr h n old Estim NbLoop Diff
5 index 6 div mul abs 6 index lt { exit } if
1 add dup 9 eq { exit } if
} loop % a b var f ierr h n old Estim NbLoop
exch 5 -1 roll 6 div mul mark 10 2 roll cleartomark
end
} def
% ------------------------------------ math stuff ----------------------------------
%
% Matrix A in arrays of rows A[[row1][row2]...]
% with [row1]=[a11 a12 ... b1]
% returns on stack solution vector X=[x1 x2 ... xn]
/SolveLinEqSystem { % on stack matrix M=[A,b] (A*x=b)
10 dict begin % hold all ocal
/A exch def
/Rows A length def % Rows = number of rows
/Cols A 0 get length def % Cols = number of columns
/Index [ 0 1 Rows 1 sub { } for ] def % Index = [0 1 2 ... Rows-1]
/col 0 def
/row 0 def
/PR Rows array def % PR[c] = pivot row for row row
{ % starts the loop, find pivot entry in row r
col Cols ge row Rows ge or { exit } if % col < Cols and row < Rows else exit
/pRow row def % pRow = pivot row
/max A row get col get abs def % get A[row[col]], first A[0,0]
row 1 add 1 Rows 1 sub { % starts for loop 1 1 Rows-1
/j exch def % index counter
/x A j get col get abs def % get A[j[r]]
x max gt { % x>max, then save position
/pRow j def
/max x def
} if
} for % now we have the row with biggest A[0,1]
% with pRow = the pivot row
max 0 gt { % swap entries pRow and row in i
/tmp Index row get def
Index row Index pRow get put
Index pRow tmp put % and columns pRow and row in A
/tmp A row get def
A row A pRow get put
A pRow tmp put % pivot
/row0 A row get def % the pivoting row
/p0 row0 col get def % the pivot value
row 1 add 1 Rows 1 sub { % start for loop
/j exch def
/c1 A j get def
/p c1 col get p0 div def
c1 col p put % subtract (p1/p0)*row[i] from row[j]
col 1 add 1 Cols 1 sub { % start for loop
/i exch def
c1 dup i exch % c1 i c1
i get row0 i get p mul sub put
} for
} for
PR row col put
/col col 1 add def
/row row 1 add def
}{ % all zero entries
/row row 1 add def % continue loop with same row
} ifelse
} loop
/X A def % solution vector
A Rows 1 sub get dup
Cols 1 sub get exch
Cols 2 sub get div
X Rows 1 sub 3 -1 roll put % X[n]
Rows 2 sub -1 0 { % for loop to calculate X[i]
/xi exch def % current index
A xi get % i-th row
/Axi exch def
/sum 0 def
Cols 2 sub -1 xi 1 add {
/n exch def
/sum sum Axi n get X n get mul add def
} for
Axi Cols 1 sub get % b=Axi[Cols-1]
sum sub % b-sum
Axi xi get div % b-sum / Axi[xi]
X xi 3 -1 roll put % X[xi]
} for
X
end
} def
%
/c@_0 2.515517 def
/c@_1 0.802853 def
/c@_2 0.010328 def
/d@_1 1.432788 def
/d@_2 0.189269 def
/d@_3 0.001308 def
/norminv {
5 dict begin
neg 1 add 1 exch div ln 2 mul sqrt
/t exch def
/t2 t dup mul def
/t3 t2 t mul def
c@_0 c@_1 t mul add c@_2 t2 mul add 1 d@_1 t mul add
d@_2 t2 mul add d@_3 t3 mul add div neg t add
end
} def
%end{norminv Michael Sharpe}
%
%
% END pst-math.pro
|
