Back
Download
See GEN_MU for a Sieve based method of generating values of the Mobius
mu-function, where mu(n) = 0 if n has a repeated factor and otherwise
mu(n)= (-1)^k where n is the product of k distinct primes.
GEN_MU 30
0 1 _1 _1 0 _1 1 _1 0 0 1 _1 0 _1 1 1 0 _1 0 1 0 1 1 1 0 0 1 0 0 1
Note: The function SIEVE is found in math.d4f.
Z<-ZETA.D4F
|class description
SQRT2<-#sqrt 2
LOG2 <- #ln 2
$CS_LUT<-" ", #av 176+i3 | Chaos Cascade. Lookup table
$CS_CNT<-1500
$PM_GAIN<-60.0
Z<-ADD_ALIAS ZETA.D4S
"ZETA.D4S"
edzeta|WW.ED"zeta.d4f ",ARGS
ldzeta|#copy"zeta.d4f" & ZETA.D4F
plotzeta|H<-PLOTZETA ARGS
zeta_info|3 #sed"X:d;$zeta.inf" & Z<-WW.FSCR 3
zbars|DS_PLOT 0, (?145), 0.5 7
mzbars|(GEN_MU 600) DS_SCREEN 0, (0.1*?2500), 0.5 5+?3
paz|Z<-P_AZETA ARGS
"zeta.inf">>EOF
These functions are based on the Zeta function of Number Theory.
DS_PLOT Animation of adding terms in a divergent series.
SZETA Attempt to calculate zeta on the critical line.
Very inaccurate.
Use alternating series a-zeta for calculations.
plot |a-zeta(0.5+it)| for various intervals as 400x500 xbm image.
use Newton-Raphson to seek zeros.
track phase changes. clock animation.
EOF
Z<-DS_PLOT BOX;LN;N;S;T;X;Y;Z0;WW;SI;TXT;EV.ALX
|Animate convergence of a dirichlet series in a box
|Plot contours of z+=u(x,y,t)
EV.ALX<-"10 #wput 0 0,#screen & GO"
$GAIN<-1.0 #ifndef "$GAIN"
CNT<-25000 #ifndef "$CS_CNT"
BOX<-0 0 1 1 #ifndef "BOX"
WW<-0 0,SZ<-#screen
N<-r$CS_LUT |Lookup table
|Map screen to unit square
Z<-yiSZ
X<-BOX[1]+BOX[3]* Z[1]%1.0*SZ[1]
Y<-BOX[0]+BOX[2]* m Z[0]%1.0*SZ[0]
T<-1 & Z0<-SZr0.00
AA: ->(T=CNT)/XX
LN<-#ln T
S<-(2 o X*LN) * #exp -0.5+Y*LN
Z<-Z0<-Z0+S & ->(T>CNT)/XX
Z<-"N" #av ($GAIN*Z) #mod 48.0
TXT<-SZ r $CS_LUT[Z #mod N]
TXT #wput WW
SI<-,Z #mod 16
CC: (SZrSI) #wput WW
WAIT 1
T++ & ->AA
XX: Z<-""
Z<-VEC DS_SCREEN BOX;LN;N;S;T;X;Y;Z0;WW;SI;TXT;VT;EV.ALX
|Animate convergence of a dirichlet series in a box
|Plot contours of z+=u(x,y,t)
EV.ALX<-"10 #wput 0 0,#screen & GO"
$GAIN<-1.0 #ifndef "$GAIN"
VEC<-(25000 r 1) #ifndef "VEC"
BOX<-0 0 1 1 #ifndef "BOX"
$CS_LUT<-" ", #av 176+i3 | Chaos Cascade. Lookup table
CNT<-rVEC
WW<-0 0,SZ<-#screen
N<-r$CS_LUT |Lookup table
|Map screen to unit square
Z<-yiSZ
X<-BOX[1]+BOX[3]* Z[1]%1.0*SZ[1]
Y<-BOX[0]+BOX[2]* m Z[0]%1.0*SZ[0]
T<-1 & Z0<-SZr0.00
AA: ->(T=CNT)/XX
->(0=VT<-VEC[T])/DD
LN<-#ln T
S<-(2 o X*LN) * #exp -0.5+Y*LN
Z<-Z0<-Z0+VT*S & ->(T>CNT)/XX
Z<-"N" #av ($GAIN*Z) #mod 48.0
TXT<-SZ r $CS_LUT[Z #mod N]
TXT #wput WW
SI<-,Z #mod 16
CC: (SZrSI) #wput WW
WAIT 1
DD: T++ & ->AA
XX: Z<-""
Z<-GEN_MU NMAX;J;L;P;P2;PL
|Generate table of Mobius function
Z<-NMAXr1
PL<-2,SIEVE (NMAX+1)%2
Z[1]<-1 & J<-0
AA: ->(NMAX < 2*P<-PL[J++])/BB
L<-(NMAX+P2)%P2<-P*P & Z[P2*iL]<-0
L<-P*i(NMAX+P)%P & Z[L]<--Z[L]
->AA
BB:Z[P]<-_1
Z<-ZETA S
|Zeta(s=x+iy)
|Z(s)= 1/1-2**(1-s) * sum (-1)**(n-1)/n**s
Z<-SZETA T;D;T2;X;Y
|Zeta (0.5+iT)
D<-3-2*2 o T2<-T * LOG2
X<-1-SQRT2* 2 o T2
Y<-SQRT2 * 1 o T2
Z<-25000 ALT_CSUM T
Z<-Z C_MULT (X,Y)%D
Z<-A C_MULT B;X;Y
|Multiply two vectors of complex numbers
X<-A*B & Y<-A*mB
Z<-(X[:0]-X[:1]),Y[:0]+Y[:1]
Z<-CNT ALT_CSUM T;N;XD;XL
|Sum (-1)**n-1 * 1/n**0.5+it
CNT<-10000 #ifndef "CNT"
XL<-#ln N<-1+iCNT
XD<-(CNT r 1 _1) * #sqrt N
Z<-((CNT,2) r 2 1) o T*2/yXL
Z<-+/Z%2/yXD
Z<-CNT IZETA T;N;XD;XL;XN
| 1%ZETA(0.5+iT)
CNT<-10000 #ifndef"CNT"
->(0<#nc"MU_SEQ")/A0
MU_SEQ<-1 d GEN_MU 1+CNT
A0: XL<- #ln N<-1+iCNT
XD<-#sqrt N
XN<-((CNT,2)r 2 1) o T*2/yXL
Z<-CNT PLOTZETA STR;A;B;C21;DT;N;T;LN;XD
|for T=a step dt to b: print alt_zeta(0.5+iT) : next T
N<-rZ<-#fi STR<-"" #ifndef "STR"
A<-Z[0] & DT<-(0.01,Z[1])[N>1] & B<-((A+5)*N<3)+Z[2]
N<-1+iCNT<-10000 #ifndef "CNT"
T<-A & Z<-3 #sed":d"
C21<-CNT/:2 1 & LN<-2/y #ln N & XD<-2/y(CNT r 1 _1)* #sqrt N
AA: ->(T>B)/XX
Z<-+/(C21 o T*LN)%XD
STR<-8 4 z T, Z, #sqrt +/Z*Z
STR
3 #print STR
T<-T+DT & ->AA
XX: Z<-3 #sed ":cm"
Z<-M PM.PLOT YX;I;J;N
|Make ascii plot of points in size M
|Join YX[0 1], YX[1 2] ... etc.
Z<-(*/M) r " "
Z[MbYX]<-"O"
N<-1tr:YX & I<-0
AA: ->(I=N-1)/XX
T<-#kjvec -/YX[I+1 0]
J<-(+\T)+(rT)rYX[I]
Z[MbJ]<-"O"
I++ & ->AA
XX: Z<-MrZ
Z<-P_AZETA STR;N;IY;YX
|Make .xpm plot from values
Z<-PLOTZETA STR
IY<-"N" #av $PM_GAIN * #fi ," ", 0 24 d Z
YX<-(yIY),yiN<-rIY
Z<-#theta (400,N) PM.PLOT YX
Z<-"azeta" PM2XBM Z
Z<-FN PM2XBM X;NB;I;M;T
|Make .xbm file of an byte matrix
FN<-"bitmap" #ifndef "FN"
M<-r:X & NB<-(7+M[1])%8 & I<-0
7 #sed":d"
#print "#define ",FN, "_width ",zM[1]
#print "#define ",FN, "_height ",zM[0]
#print "static char ",FN,"_bits[]={"
AA: ->(I=1tM)/XX
T<-,m(NB,8)r(8*NB)t " " #ne X[I:]
T<-1m,((NB,3)r",0x"),(NB,2)r"0123456789abcdef"[(4r2)b(2*NB,2)rT]
#print T
I++ & ->AA
XX: #print "};"
Z<-WW.FSCR 7
Back to the Top