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:
|
function PolEvalX(x: extended; const a: array of extended; n: integer): extended;
var
i: integer;
s: extended;
begin
if n<=0 then begin
PolEvalX := 0.0;
exit;
end;
s := a[n-1];
for i:=n-2 downto 0 do s := s*x + a[i];
PolEvalX := s;
end;
function cterm(n: integer; z: extended): extended;
const
n0=22;
c0: array[0..n0-1] of extended = (
+0.38268343236508977173,
+0.43724046807752044936,
+0.13237657548034352332,
-0.01360502604767418865,
-0.01356762197010358089,
-0.00162372532314446528,
+0.00029705353733379691,
+0.00007943300879521470,
+0.00000046556124614505,
-0.00000143272516309551,
-0.00000010354847112313,
+0.00000001235792708386,
+0.00000000178810838580,
-0.00000000003391414390,
-0.00000000001632663390,
-0.00000000000037851093,
+0.00000000000009327423,
+0.00000000000000522184,
-0.00000000000000033507,
-0.00000000000000003412,
+0.00000000000000000058,
+0.00000000000000000015);
n1=23;
c1: array[0..n1-1] of extended = (
-0.02682510262837534703,
+0.01378477342635185305,
+0.03849125048223508223,
+0.00987106629906207647,
-0.00331075976085840433,
-0.00146478085779541508,
-0.00001320794062487696,
+0.00005922748701847141,
+0.00000598024258537345,
-0.00000096413224561698,
-0.00000018334733722714,
+0.00000000446708756272,
+0.00000000270963508218,
+0.00000000007785288654,
-0.00000000002343762601,
-0.00000000000158301728,
+0.00000000000012119942,
+0.00000000000001458378,
-0.00000000000000028786,
-0.00000000000000008663,
-0.00000000000000000084,
+0.00000000000000000036,
+0.00000000000000000001);
n2=24;
c2: array[0..n2-1] of extended = (
+0.00518854283029316849,
+0.00030946583880634746,
-0.01133594107822937338,
+0.00223304574195814477,
+0.00519663740886233021,
+0.00034399144076208337,
-0.00059106484274705828,
-0.00010229972547935857,
+0.00002088839221699276,
+0.00000592766549309654,
-0.00000016423838362436,
-0.00000015161199700941,
-0.00000000590780369821,
+0.00000000209115148595,
+0.00000000017815649583,
-0.00000000001616407246,
-0.00000000000238069625,
+0.00000000000005398265,
+0.00000000000001975014,
+0.00000000000000023333,
-0.00000000000000011188,
-0.00000000000000000416,
+0.00000000000000000044,
+0.00000000000000000003);
n3=24;
c3: array[0..n3-1] of extended = (
-0.00133971609071945690,
+0.00374421513637939370,
-0.00133031789193214681,
-0.00226546607654717871,
+0.00095484999985067304,
+0.00060100384589636039,
-0.00010128858286776622,
-0.00006865733449299826,
+0.00000059853667915386,
+0.00000333165985123995,
+0.00000021919289102435,
-0.00000007890884245681,
-0.00000000941468508130,
+0.00000000095701162109,
+0.00000000018763137453,
-0.00000000000443783768,
-0.00000000000224267385,
-0.00000000000003627687,
+0.00000000000001763981,
+0.00000000000000079608,
-0.00000000000000009420,
-0.00000000000000000713,
+0.00000000000000000033,
+0.00000000000000000004);
n4=25;
c4: array[0..n4-1] of extended = (
+0.00046483389361763382,
-0.00100566073653404708,
+0.00024044856573725793,
+0.00102830861497023219,
-0.00076578610717556442,
-0.00020365286803084818,
+0.00023212290491068728,
+0.00003260214424386520,
-0.00002557906251794953,
-0.00000410746443891574,
+0.00000117811136403713,
+0.00000024456561422485,
-0.00000002391582476734,
-0.00000000750521420704,
+0.00000000013312279416,
+0.00000000013440626754,
+0.00000000000351377004,
-0.00000000000151915445,
-0.00000000000008915418,
+0.00000000000001119589,
+0.00000000000000105160,
-0.00000000000000005179,
-0.00000000000000000807,
+0.00000000000000000011,
+0.00000000000000000004);
var
q: extended;
begin
q := sqr(z);
case n of
0: cterm := PolEvalX(q,c0,n0);
1: cterm := z*PolEvalX(q,c1,n1);
2: cterm := PolEvalX(q,c2,n2);
3: cterm := z*PolEvalX(q,c3,n3);
else cterm := PolEvalX(q,c4,n4);
end;
end;
function theta(t:extended):extended;
var
x,s: extended;
const
c0 = 1.0/48.0;
c1 = 7.0/5760.0;
c2 = 31.0/80640.0;
c3 = 127.0/430080.0;
begin
x := 1.0/sqr(t);
s := (((c3*x + c2)*x + c1)*x + c0)/t;
theta := 0.5*t*(ln(t/(2*Pi)) - 1.0) - Pi/8 + s;
end;
procedure zetacrit(t: extended; var re,im: extended);
const
m=2;
var
i,n: longint;
axx,axy: extended;
q,y,p,s,r,tt:extended;
begin
q := sqrt(t/(2*Pi));
n := trunc(q);
s := 0.0;
tt := theta(t);
for i:=1 to n do begin
s := s + 1/sqrt(i)*cos(tt - t*ln(i));
end;
s :=2*s;
y := 1.0;
r := 0.0;
p := 2.0*(q-n) - 1.0;
for i:=0 to m do begin
r := r + cterm(i,p)*y;
y := y/q;
end;
r := r/sqrt(q);
if odd(n+1) then r := -r;
sincos(-tt, axy, axx);
re := (s + r)*axx;
im := (s + r)*axy;
end;
|
Gammatesters Hinweisen habe ich die Berechnung hinbekommen und damit das erste kleinere Programm, dass eine ähnliche Grafik erzeugt.
