Z<-NM.D4F |New math. 2003. _PI <-4*_3o1 _SQRT2 <-#sqrt 2 _LOG2<- #ln 2 $PRIMES<-2,SIEVE 18000 Z<-ADD_ALIAS NM.D4S "NM.D4S" ednm|WW.ED"nm.d4f ",ARGS ldnm|#copy"nm.d4f" & NM.D4F nm_info|3 #sed"X:d;$ nm.inf" & Z<-WW.FSCR 3 bernoulli|Z<-BERNOULLI ARGS mlink|Z<-MLINK ARGS ulam|Z<-ULAM_SPIRAL ARGS nrazeta|Z<-NR_AZETA ARGS nrz|Z<-NR_ZEROS ARGS shuffle|SHUFFLE_CARDS ARGS factors|Z<-FACTORIZE ARGS & #deb 16 0 z Z "nm.inf">>EOF This file contains functions to generate the ULAM spiral. It also contains links to other maths scripts. * ada.d4f Bernoulli numbers. * boxes.d4f Boxes game. * hail.d4f the 'n = (n&1) ? 3*n+1 : n/2;' problem. * kara.d4f see Sloane Integer Sequence Database. * imath.d4f Tutorial on RSA Encryption. * mandel.d4f Mandelbrot set. DOS/WINDOWS * nds.d4f Connected with Dirichlet series. * polya.d4f Tutorial on Polya theory of counting. * power.d4f Voting systems. Banzhaf power index. * qcf.d4f Quadratic continued fractions. Pell's Equation. * rm.d4f Symbolic mathematics. * rtj.d4f Rush Hour Traffic jam solution. * sieve.d4f Old sieve functions from STSCs APL * xpm.d4f Mandelbrot set. Linux. EOF Z<-MLINK STR;EXPR;F;K;M |Select .d4f files EXPR<-"vi" STR<-"-x EXPR<-SHIFT" GETOPTS STR K<-0 & M<-0 #sed"X:d;$nm.inf;g/\*.*\.d4f/;:cm" WHILE _1 #ne K<-(4 8 14 64,K) PMENU M STR<-#deb M[K] & Z<-SHIFT & F<-SHIFT Z<-DO EXPR," ",F WEND Z<-FIB N |Fibonacci sequence, N terms Z<-1 1 WHILE N>rZ;Z<-Z,+/_2tZ;WEND Z<-N P_INV A |Inverse of series a[0]+a[1]x+a[2]x*x.... |assume a[0]=1 Z<-1 & N<-20 #ifndef"N" WHILE N>rZ;Z<-Z,-+/(mZ)*(rZ)t1dA;WEND Z<-X CV Y;I;J;N;T |Convolution or polynomial multiplication Z<-iI<-0 & N<-1+c/J<-,(irX) #outer+ irY & T<-,X #outer* Y WHILE I(D<0)/0 Z<-K<-(_1tX) % (_1tY) ->(D=0)/0 R<-CTRIM X-K*(-RX)tY Z<-CTRIM ((D/0),K)+(1+D)t CTRIM R CDIV Y Z<-CTRIM X |Get rid of excess trailing zeros IF 0=rZ<-(~mf\mX = 0)/X; Z<-0 Z<-V P2STR C;J;N;P |Format a polynomial V<-1 t "x" #ifndef "V" Z<-"" & BREAKIF 0=N<-rC P<-#sstomat "//x",,"/","x","^",zy2+iN-2 J<-(C #ne 0)/iN Z<-(Z #ne " ")/Z<-1d ,"+",(zyC[J]),P[J] Z<-"/+1x/+x/*" #rxsubs Z Z<-"/+-/-/*" #rxsubs Z Z<-"/-1x/-x/*" #rxsubs Z IF V #ne "x"; Z[(Z="x")/irZ]<-V Z<-ULAM N;DZ;T |Ulam spiral to size N |See http://mathworld.wolfram.com/PrimeSpiral.html K<-Z<-0 WHILE N>T<-c/,Z; DZ<-T+1+ic/r:Z ->(A1,A2,A3,A4)[K #mod 4] A1:Z<-Z,ymDZ & ->CC A2:Z<-(mDZ),:Z & ->CC A3:Z<-(yDZ),Z & ->CC A4:Z<-Z,:DZ & ->CC CC:K++ & WEND R<-SIEVE N;J;K;Z |SIEVE OF ERASTOTHENES FOR ODD PRIMES |Get odd primes up to 2*N R<-(1,N-1)/"10" & Z<-iJ<-0 WHILE (J*J)<2*N;Z<-Z,K<-Ri"0";R[K+J*i(N+J-K+1)%J<-1+2*K]<-"1";WEND R[Z]<-"0" & R<-1+2*(R="0")/iN Z<-GEN_MU N;J;L;P;P2;PL |Generate table of Mobius mu function Z<-N r 1 & J<-0 & PL<-2,SIEVE (N+1)%2 WHILE N>=2*P<-PL[J++] Z[P2*i(N+P2)%P2<-P*P]<-0 Z[L]<--Z[L<-P*i(N+P)%P] WEND Z[PL]<-_1 Z<-PRODUCT X;N |Multiple precision product IF 0=N<-rX; X<-"1" IF N<=1; Z<-z X & ->0 N<-N%2 & Z<-(PRODUCT NtX) * PRODUCT NdX Z<-A EA B;C;D;R |Euclidean algorithm with A>B R<-A #mod B & D<-B WHILE ~R #equiv "0"; C<-D & D<-R & R<-C #mod D & WEND Z<-D Z<-XQ_INI |Clear to Accumulate fractions |Global XQ0, XQ1 XQ0<-"0" & XQ1<-"1" & Z<-"0/1" Z<-XQ_ADD X;$0;$1;DEN;NUM;S |Accumulate fractions |Global XQ0, XQ1 S<-"" & Z<-"/" #split X<-#deb X NUM<-(XQ0*$1)+XQ1*$0 & DEN<-XQ1*$1 IF "-"=1tNUM; S<-"-" & NUM <-1dNUM IF ~"1" #equiv Z<-NUM EA DEN; NUM<-NUM%Z & DEN<-DEN%Z XQ0<-S,NUM & XQ1<-DEN & Z<-XQ0,"/",XQ1 Z<-XQ_HS N;I |Harmonic series |same as R_SUM z,1,y1+iN Z<-XQ_INI & I<-0 WHILE I(0 = (rX) f rY)/0 Z<-"/" #split X & NUM<-$0 & DEN<-$1 Z<-"/" #split Y & NUM<-NUM*$0 & DEN<-DEN*$1 IF "-"=1tNUM; S<-"-" & NUM<- 1dNUM IF ~"1" #equiv Z<-NUM EA DEN; NUM<-NUM%Z & DEN<-DEN%Z Z<-S,NUM,"/",DEN Z<-W XQ_SUM X;J;N;T |Sum matrix of signed fractions J<-0 & N<-1tr:X & Z<-XQ_INI -> (0< #nc "W")/BB WHILE J0 BB: N<-N f 1tr:W WHILE J>EOF Bernoulli numbers just come from from long division by power series. They can be used to estimate sums by integrals. The most famous formula in which they arise is Stirling's Approximation for n factorial (n! = 1 x 2 x 3 .... x n-1 x n). They come as coefficients in the expansion:- oo n x --- x ------ = > B[n]--- x --- n! e -1 n=1 The ratio 2n*(2n-1)*B[2n]/B[2n-2] tends to (6*Pi)^2. The magnitude of the numbers B[2n] are given by an approximation to the exact formula :- B[2n] = Zeta(2n)*2*(2n!)/(2*Pi)^2n EOF Z<-BERNOULLI STR;J;ID;N;T;M |Compute Bernoulli numbers as a series of rationals. QR<-0 STR<-"-real QR<-1"GETOPTS "" #ifndef "STR" N<-10 & IF 00 ID<-"" & Z<-"1" & J<-0 WHILE J=J++; S<-(~J #and 1)/"-"; P<-P,",",S,"1/",K; K<-K*zX*1+X<-2*J;WEND P<-#sstomat P CSEC<- #sstomat N RP_INV P I<-1 & Z<-",1/1" & P<-"1" & M2<-"2" |Now get the B[2n] numbers WHILE I<=N P<-P* z V*_1+V<-2*I Q<-"2"* M2-"1" & M2<-M2*"4" V<-(P,"/",Q) Q_MUL CSEC[I] Z<-Z,",",((~I #and 1)/"-"),V I++ & WEND Z<-#sstomat Z "euler_numbers.inf">>EOF Euler numbers occur in the power series for 1/cos(x) = sum E[n]*(x^n)/2n! The radius of convergence of the series pi/2. 2 * #sqrt 2n*(2n-1)*E[n-1]/E[n] tends to pi/2 as n->oo The Euler numbers are all integers. It is easy to prove this. Differentiate the function f(x)=1/cos x n times. The derivative may be expressed as a function (P[n](sin x))/cos(x) and P[n](x) is a polynomial with integer coefficients. The Euler numbers are the constant terms of these polynomials. The polynomials P[n] satisfy the recurrence relationship:- P[n+1](x) = (1-x^2)DP[n](x)+ (n+1)xP[n](x) Here DP[n](x) is the derivative with respect to x. This means the Euler numbers may be calculated by integer addition and multiplication. The Bernoulli numbers which come from the expression x/sin(x) are much harder to calculate. EOF Z<-E2_NUMS N;CSEC;I;K;P;S;V;X |Euler numbers. Expansion for 1/cos x N<-20 #ifndef "N" Z<-"" & I<-0 & K<-"1" WHILE N>=I++; S<-(~I #and 1)/"-";Z<-Z,",",S,"1/",K;K<-K*z(X-1)*X<-2*I;WEND CSEC<-#sstomat N RP_INV #sstomat Z P<-zI<-1 & Z<-",1" WHILE I>EOF "pi = 4 * arctan(1)" 16 10 z 4*_3o1 "zeta(2)=pi*pi/6 si pi=sqrt 6*zeta(2)" 16 10 z #sqrt 6*100000 RZETA 2 "zeta(4)= pi^4/90 so pi = sqrt sqrt 90*zeta(4)" 16 10 z #sqrt #sqrt 90*100000 RZETA 4 EOF Z<-A DS_SUM YX;N;LN;T |Complex Z= sum :a[n]*n^-(x+iy), real coefficients A LN<-#ln 1+iN<-rA Z<-A * #exp -LN*YX[:1] T<-(N/:1 2) o 2/y YX[:0]*LN Z<-_1 1 * +/(2/yZ)*T Z<-CNT CZETA YX;MF |Complex zeta Zeta[s] = A(s)/(1-2^1-s) CNT<-25000 #ifndef "CNT" Z<-(CNT r 1 _1) DS_SUM YX MF<-0 1-(_1 1 * 1 2 o YX[0])*#exp _LOG2*1-YX[1] MF<-(_1 1*MF)%+/MF*MF Z<-Z CMULT MF Z<-A CMULT B |(i a[0]+a[1]) * (i b[0] + b[1]) = i(a0b1+a1b0)+(b1a1-a0b0) Z<-(+/A*mB),--/B*A "NRZ.OPTS" -cnt CNT<-#fi SHIFT -eps EPS<-#fi SHIFT -iter M<-#fi SHIFT -seed SEED<-#fi SHIFT -v QV<-1 Z<-NR_AZETA STR;A;DA;CNT;EPS;LN;N;I;M;CS;SN;SEED;T;DT;TL;QV;VEC |Zeros of real or imaginary parts of sum V[n]*exp(-i t log n) CNT<-10000 & EPS<-0.005 & SEED<-1 & M<-20 & QV<-0 STR<-NRZ.OPTS GETOPTS STR LN<-#ln N<-1+iCNT & VEC<-(CNTr1 _1) * 1%#sqrt N T<-SEED & DT<-0 & I<-0 |Real part WHILE M>I++ CS<-2oTL<-T*LN & SN<-1oTL & A<-+/VEC*CS & DA<--+/VEC*SN*LN IF QV; 10 6 z T, DT, A, DA BREAKIF EPS > #abs DA & BREAKIF EPS > #abs DT<-A%DA T<-T-DT & WEND Z<-T & I<-0 & DT<-0 |Use this root to seek a root for imaginary part WHILE M>I++ CS<-2oTL<-T*LN & SN<-1oTL & A<-+/VEC*SN & DA<-+/VEC*CS*LN IF QV; 10 6 z T, DT, A, DA BREAKIF EPS > #abs DA & BREAKIF EPS > #abs DT<-A%DA T<-T-DT & WEND Z<-Z,T Z<-NR_ZEROS STR;RANGE;STEP;T;V RANGE<-10 50 & STEP<-1 7 #print STR STR<-"-range RANGE<-NARGS 2|-step STEP<-#fi SHIFT" GETOPTS STR T<-RANGE[0] WHILE T #abs -/Z; "seed:",(12 6zT)," zero=",V<-16 8 z Z & #print V T<-T+STEP & WEND Z<-WW.FSCR 7 "ULAM.SC" ulam -size 16384 ulam -size 16384 VAL<-1+GEN_MU 17000 ulam -size 16384 VAL<-1=(i 17000) #mod 7 "ULAM.OPTS" -size N<-#fi SHIFT -o OFILE<-#fi SHIFT -lut LUT<-SHIFT -seq Z<-ULAM_SPIRAL STR;EXP;LUT;N;NMAX;VAL |Render ulam spiral frames N<-30 & OFILE<-"" & LUT<-".123456789ABCDEF" EXP<-"VAL<-N r 0 & P<-2, SIEVE (1+N)%2 & VAL[P]<-1" STR<-ULAM.OPTS GETOPTS STR N<-r,Z<-ULAM N IF STR; EXP<-STR 0 r #do EXP Z<-VAL[Z] IF LUT; Z<-LUT[Z] "ds.inf">>EOF Hardy & Wright. Theorems 287-291 D(s) = sum C[n] n^-s Z(s) C=1 zeta function 1/Z(s) C= MU[n] Mobius function Z(s)*Z(s) C= d[n] Divisor function Z(s-1)/Z(s) C= phi[n] Totient function Z(s)*Z[s-1] C= sigma[n] sum of divisors EOF Z<-DIVISORS N |Run out list of divisors Z<-(0=N #mod Z)/Z<-1+iN Z<-FACTORIZE STR;LST;N;J;P |Factorise a number LST<-$PRIMES N<-#fi STR & Z<-iJ<-0 WHILE N>=P*P<-LST[J++] IF 0=N #mod P; Z<-Z,P & N<-N%P & J-- WEND Z<-Z,(N>1)/N Z<-A D_CV B;D;I;N |Convolution of dirichlet series |z[n] = sum (a[d]*b[n%d]) Z<-0 & I<-1 & N<-f/(rA),rB |maybe N<-(rA)*rB WHILE I<=N Z<-Z,+/A[D]*mB[D<-DIVISORS I] I++ WEND Z<-GEN_PHI N |Series of Euler totient function (phi) up to N Z<-(iN) D_CV GEN_MU N Z<-GEN_SIGMA N;I |Tabulate sum of divisors of N Z<-Nr0 & I<-1 WHILE IP*P<-LST[J++] BREAKIF Z<-0 = N #mod P WEND Z<-1-Z Z<-SCAN_AMI N;PL;X |Scan for amicable numbers in family |Test primes in sequence 3.2^n-1 Z<-i0 & M<-1 & PL<-2,SIEVE 35000 WHILE 00 A<-iN & CNT<-1 WHILE QV;BREAKIF (iN) #equiv A<-A[Z];CNT++;WEND Z<-CNT Z<-S1_NUMS MAX |Stirling nimbers of the first kind |x(x-1)..(x-n+1)=sum{s(n,k)x^k:k=1 to n} Z<-0 1 & J<-0 | S(1,1) WHILE MAX>rZ;1dZ & #print 12 z 1dZ & J-- & Z<-Z CV J,1 & WEND Z<-1dZ Z<-R_SUM STR;A;B;C;D;DEN;K;N;NUM;S;U |Sum rationals stored as 2n pairs p,q for p/q |denominator must be positive | a/b + c/d = (ad+bc)/bd |can be recursive -- divide and conquer N<-1+c/U<-+\" "=STR<-#deb STR IF 1 = N #mod 2; Z<-"0 1" ; ->0 IF N=2; Z<-STR; ->0 ->(N>4)/ADD |Ordinary case A<-SHIFT;B<-SHIFT;C<-SHIFT;D<-SHIFT NUM<-(A*D)+B*C; DEN<-B*D IF DEN #equiv "0"; Z<-1 0; ->0 IF S<-"-"=1tNUM; NUM<-1dNUM IF ~"1" #equiv Z<-NUM EA DEN; NUM<-NUM%Z; DEN<-DEN%Z IF S; NUM<-"-", NUM Z<-NUM," ",DEN ->0 ADD: | Recursion, so split STR K<-(K #mod 2)+K<-N%2 Z<-R_SUM (R_SUM (U=K)/STR Z<-R_FLOOR STR;A;B;S |Floor function of fraction. work for negative values. A<-SHIFT;B<-SHIFT IF S<-"-"=1tA; A<-1dA Z<-A%B IF S; Z<-"-", Z+(~"0" #equiv $REMAINDER)/"1" Z<-R_FRAC STR;S;T | Fractional part of p/q i.e p/q-[p/q] S<-"-"; IF "-"=1 t T<-R_FLOOR STR; S<-"";T<-1dT Z<-R_SUM STR," ",S,T," 1" # John Alen Paulos. Guardian article. # Billboards asked for the first ten digit prime in the # decimal expansion in e. # These functions do this for e and ln(2) # Normally all ten digit primes shouls occur somewhere. Z<-DEC_E N;T;K |Decimal expansion of e to N places. Z<-"0";T<-"1",N/"0"; K<-0 WHILE K0 Z<-3 #sed":cm" 3 #sed":d" IF QY; Z<-yZ #print"sieve -size ", z r Z #print Z Z<-3 XSM2XPM"-o ", OFILE Z<-FB XI_FACTORS N;I;T |Factorize a number by trial division |FB is int vector and N is long |return a vector of length rFB I<-0; Z<-(1+rFB)r0 WHILE I1; Z[rFB]++ Z<-MC_A N;R;T |Make a sequence of values _1 0 1. A[0]=1 Z<-1 R<-3 4 r 1 1 _1 0 1 _1 0 0 1 _1 _1 0 WHILE N >= rZ Z<-,(R[1+*+/Z])[?4] WEND Z<-N DEC_ATAN X;K;S;T |Calculate N digits of arctan 1%X |Compute 4 * atan(1%2)+atan(1%3) N<-50 #ifndef "N" Z<-"0"; T<-("1",N/"0") % zX; K<-0; S<-1 WHILE K<2*N Z<-Z+T % z S*1+2*K S<- -S K++ BREAKIF "0" #equiv T<-T% z X*X WEND