Mathematician: MQL recipes

Recipes

Critical points of a function of two variables
A toy RSA key, end to end
Spectrum and invariants of a symmetric matrix
Cube roots of a complex number
The Basel problem to 30 digits
Gamma and Beta: classic identities checked numerically
Hyperbolic functions and the Gudermannian
Elementary functions of a complex variable
Machine epsilon and catastrophic cancellation
Winning move in the game of Nim

executed on MentDB

Critical points of a function of two variables

Find the critical points of f(x, y) = x^3 - 3x + y^2 - 2y and classify them with the second derivative test (sign of the Hessian determinant D = fxx*fyy - fxy^2).

Pipeline

01 · derivePartial derivatives fx, fymath diff→02 · solveSolve grad f = 0 exactlycas solve→03 · classifyEvaluate D at each pointmath eval

Script

admin

#1. Gradient of f = x^3 - 3x + y^2 - 2y;
-> "[f]" "x^3-3*x+y^2-2*y";
-> "[fx]" (math diff [f] "x");
-> "[fy]" (math diff [f] "y");
#2. Critical points: solve grad f = 0;
-> "[pts]" (cas solve (concat "[\"" [fx] " = 0\", \"" [fy] " = 0\"]") "[\"x\", \"y\"]");
#3. Hessian test D = fxx*fyy - fxy^2;
-> "[D]" (concat "(" (math diff [f] "x" 2) ")*(" (math diff [f] "y" 2) ")-(" (math diff [fx] "y") ")^2");
concat "critical points " [pts] "; D = " [D] "; D(1,1) = " (math eval [D] "{\"x\": 1, \"y\": 1}") " so a local minimum; D(-1,1) = " (math eval [D] "{\"x\": -1, \"y\": 1}") " so a saddle point"

Result

mentdb

critical points [{"x":"-1","y":"1"},{"x":"1","y":"1"}]; D = (6*x)*(2)-(0)^2; D(1,1) = 12 so a local minimum; D(-1,1) = -12 so a saddle point

Go further

Find the minimum numerically with num minimize and a starting point {"x": 0.5, "y": 0}
Expand f around the minimum with cas series to see the quadratic form

executed on MentDB

A toy RSA key, end to end

Build a small RSA key from two primes (public exponent 65537), compute the private exponent with the Carmichael function, then encrypt and decrypt the message 42.

Pipeline

01 · primesPrimes p and qarith next_prime→02 · lambdalambda(n) = lcm(p-1, q-1)arith lcm→03 · keyPrivate key d = e^-1 mod lambdaarith invmod→04 · cipherEncrypt and decryptarith powmod

Script

admin

#1. Two primes and the modulus;
-> "[p]" (arith next_prime 1000);
-> "[q]" (arith next_prime 2000);
-> "[n]" (math calc "[p]*[q]");
#2. Carmichael function and private key;
-> "[lambda]" (arith lcm (math calc "[p]-1") (math calc "[q]-1"));
-> "[d]" (arith invmod 65537 [lambda]);
#3. Encrypt then decrypt the message 42;
-> "[c]" (arith powmod 42 65537 [n]);
concat "n=" [n] " d=" [d] " cipher=" [c] " decrypted=" (arith powmod [c] [d] [n])

Result

mentdb

n=2021027 d=122561 cipher=512387 decrypted=42

Go further

Factor n back with arith factors to see why real keys use 1024-bit primes
Compare lambda(n) with Euler's totient arith phi

executed on MentDB

Spectrum and invariants of a symmetric matrix

Compute the eigenvalues of a symmetric 3x3 matrix, check that their product is the determinant and their sum the trace, and decide if the matrix is positive definite and well conditioned.

Pipeline

01 · eigenEigenvaluesmatrix eigen→02 · detDeterminantmatrix det→03 · traceTrace and summatrix trace→04 · condCondition numbermatrix cond

Script

admin

#1. A symmetric stiffness matrix;
-> "[A]" "[[4, 1, 2], [1, 3, 0], [2, 0, 5]]";
#2. Eigenvalues;
-> "[ev]" (json load select "rx_eig" (matrix eigen [A]) /values);
-> "[l1]" (json select "rx_eig" /values[0]);
-> "[l2]" (json select "rx_eig" /values[1]);
-> "[l3]" (json select "rx_eig" /values[2]);
#3. Check the invariants: product = det, sum = trace;
concat "product of the eigenvalues = " (math round_sig (math calc "[l1]*[l2]*[l3]") 10) ", det = " (matrix det [A]) "; sum = " (stat sum [ev]) ", trace = " (matrix trace [A]) "; smallest eigenvalue " [l3] " > 0: positive definite; cond = " (matrix cond [A])

Result

mentdb

product of the eigenvalues = 43.00000000, det = 43; sum = 12, trace = 12; smallest eigenvalue 1.854897308799577 > 0: positive definite; cond = 3.595389920856739

Go further

Confirm positive definiteness with matrix cholesky
Get the singular values with matrix svd

executed on MentDB

Cube roots of a complex number

Solve z^3 = 8i: get the three roots, write the principal root in polar form, and check z^3 and the fact that the roots sum to zero.

Pipeline

01 · rootsThe n-th rootscomplex roots→02 · polarModulus and argumentcomplex to_polar→03 · checkRaise to the power 3complex pow→04 · sumSum of the rootscomplex add

Script

admin

#1. The three cube roots of 8i;
-> "[roots]" (complex roots "8i" 3);
#2. Principal root and its polar form;
-> "[z]" (json load select "rx_cr" [roots] /[0]);
-> "[r]" (json load select "rx_cp" (complex to_polar [z]) /r);
-> "[deg]" (math rad_to_deg (json select "rx_cp" /theta));
#3. Check z^3 and the sum of the roots;
concat "roots " [roots] "; principal root " [z] " = " [r] " at " (math round [deg] 6) " deg; z^3 = " (complex pow [z] 3) "; sum of the roots = " (complex add (complex add (json select "rx_cr" /[0]) (json select "rx_cr" /[1])) (json select "rx_cr" /[2]))

Result

mentdb

roots ["1.7320508075688774+i","-1.7320508075688774+i","-2i"]; principal root 1.7320508075688774+i = 2 at 30 deg; z^3 = 8.000000000000002i; sum of the roots = 0

Go further

Get the same roots from the polynomial z^3 - 8i with num poly_roots
Use math calc with i directly: (1.7320508075688772+i)^3

executed on MentDB

The Basel problem to 30 digits

Euler's sum of 1/k^2 equals pi^2/6. The script gets the exact value from the CAS, the same constant from the zeta function, pi^2/6 to 30 significant digits with exact decimals, and checks gamma(1/2) = sqrt(pi) and the 1/N tail of the partial sums.

Pipeline

01 · exactSymbolic sumcas eval→02 · zetazeta(2)math eval→03 · decimalspi^2/6 to 30 digitsbigdec pi→04 · checkCompare to 16 digitsbigdec compare→05 · tailPartial sum and tailnum sum

Script

admin

#1. Exact value of the Basel sum with the CAS;
-> "[exact]" (cas eval "Sum(1/k^2, {k, 1, Infinity})");
#2. The same constant is zeta of 2;
-> "[z2]" (math eval "zeta(2)");
#3. Pi squared over 6 with 30 significant digits;
-> "[basel]" (bigdec div (bigdec pow (bigdec pi 40) 2 40) 6 30);
#4. Do both values agree to 16 digits;
-> "[same]" (bigdec compare (bigdec precision [basel] 16) (bigdec precision [z2] 16));
#5. Same pi gives gamma of 1/2 = sqrt of pi;
-> "[g12]" (bigdec sqrt (bigdec pi 40) 30);
#6. Slow convergence: after N terms the tail is close to 1/N;
-> "[s]" (num sum "1/k^2" "k" 1 1000);
-> "[tail]" (math calc "1000*([z2]-[s])");
concat "{\"exact\": \"" [exact] "\", \"zeta2\": " [z2] ", \"pi2_over_6\": \"" [basel] "\", \"compare_16_digits\": " [same] ", \"sqrt_pi\": \"" [g12] "\", \"gamma_half\": " (math eval "gamma(0.5)") ", \"N_times_tail_N1000\": " [tail] "}"

Result

mentdb

{"exact": "Pi^2/6", "zeta2": 1.6449340668482264, "pi2_over_6": "1.64493406684822643647241516665", "compare_16_digits": 0, "sqrt_pi": "1.77245385090551602729816748334", "gamma_half": 1.772453850905516, "N_times_tail_N1000": 0.9995001666667225}

Go further

Try math eval "zeta(4)" against pi^4/90, again with bigdec pow.
bigdec sqrt and bigdec div keep any number of digits; bigdec precision rounds to significant digits.

executed on MentDB

Gamma and Beta: classic identities checked numerically

Check Gamma(1/2)^2 = pi, B(1/2, 1/2) = pi, the reflection formula Gamma(x)Gamma(1-x) = pi/sin(pi x), and how close Stirling's formula gets to ln 100!.

Pipeline

01 · constantpi as a stored constantconstant math_pi→02 · gammaGamma, Beta and log-Gammamath gamma→03 · betaEuler Beta integralmath beta→04 · stirlingln 100! vs Stirlingmath lngamma

Script

admin

#1. pi from the constant, its value only;
-> "[pi]" (string regexp_extract (constant math_pi) "[-0-9.]+E?[-0-9]*" 0);
-> "[g]" (math gamma 0.5);
-> "[b]" (math beta 0.5 0.5);
#2. Reflection formula at x = 0.3;
-> "[gg]" (math calc (concat (math gamma 0.3) "*" (math gamma 0.7)));
-> "[refl]" (math calc "[pi]/sin([pi]*0.3)");
#3. Stirling for 100! with lngamma;
-> "[ln]" (math lngamma 101);
-> "[st]" (math calc "0.5*ln(2*[pi]*100)+100*ln(100)-100");
concat "Gamma(1/2)^2 = " (math calc "[g]^2") " vs pi = " [pi] "; B(1/2,1/2) = " [b] "; Gamma(0.3)Gamma(0.7) = " [gg] " vs pi/sin(0.3pi) = " [refl] "; ln 100! = " [ln] ", Stirling = " [st] ", relative error " (math round_sig (math calc "([ln]-[st])/[ln]") 3) "; e = " (math e)

Result

mentdb

Gamma(1/2)^2 = 3.1415926535897936 vs pi = 3.1415926535897931; B(1/2,1/2) = 3.1415926535897936; Gamma(0.3)Gamma(0.7) = 3.8832220774509336 vs pi/sin(0.3pi) = 3.8832220774509327; ln 100! = 363.73937555556347, Stirling = 363.7385422250079, relative error 0.00000229; e = 2.718281828459045

Go further

Add the 1/(12n) correction term to Stirling and watch the error drop by three orders of magnitude.
Check the duplication formula Gamma(x)Gamma(x+1/2) = 2^(1-2x) sqrt(pi) Gamma(2x) with math gamma.

executed on MentDB

Hyperbolic functions and the Gudermannian

At x = 0.8, verify cosh^2 - sinh^2 = 1, the inverse functions asinh, acosh, atanh, the closed form atanh(1/2) = ln(3)/2, and compute the Gudermannian gd(x) = atan(sinh x) = 2 atan(tanh(x/2)).

Pipeline

01 · directsinh, cosh, tanhmath sinh→02 · inverseRound trip through the inversesmath asinh→03 · closedatanh as a logarithmmath atanh→04 · gdGudermannian two waysmath atan

Script

admin

-> "[x]" 0.8;
-> "[s]" (math sinh [x]);
-> "[c]" (math cosh [x]);
-> "[t]" (math tanh [x]);
#Gudermannian, two equivalent formulas;
-> "[gd1]" (math atan [s]);
-> "[gd2]" (math calc (concat "2*atan(" (math tanh (math calc "[x]/2")) ")"));
concat "cosh^2 - sinh^2 = " (math calc "[c]^2-[s]^2") "; asinh sinh x = " (math asinh [s]) ", acosh cosh x = " (math acosh [c]) ", atanh tanh x = " (math atanh [t]) "; atanh 0.5 = " (math atanh 0.5) " = ln 3 / 2 = " (math calc "ln(3)/2") "; gd(0.8) = " [gd1] " = " [gd2]

Result

mentdb

cosh^2 - sinh^2 = 1.0000000000000002; asinh sinh x = 0.8000000000000002, acosh cosh x = 0.8000000000000002, atanh tanh x = 0.8000000000000003; atanh 0.5 = 0.5493061443340549 = ln 3 / 2 = 0.5493061443340548; gd(0.8) = 0.7262048227415289 = 0.7262048227415289

Go further

The Gudermannian links the Mercator projection to latitude: y = gd^-1(phi) = asinh(tan phi).
Tabulate the identities on a grid with num table to see the rounding error grow with x.

executed on MentDB

Elementary functions of a complex variable

For z = 1 + 2i, check sin^2 z + cos^2 z = 1, tan(atan z) = z, sinh(iz) = i sin z, compute the principal log of i, the real number i^i = exp(-pi/2), cosh(i pi) = -1 and acos 2, which is purely imaginary.

Pipeline

01 · trigsin and cos of zcomplex sin→02 · logPrincipal logarithm and i^icomplex log→03 · inverseatan then tancomplex atan→04 · hyperbolicsinh, cosh, acoscomplex sinh

Script

admin

-> "[z]" "1+2i";
-> "[s]" (complex sin [z]);
-> "[c]" (complex cos [z]);
-> "[one]" (complex add (complex mul [s] [s]) (complex mul [c] [c]));
#i^i = exp of i log i;
-> "[ii]" (complex exp (complex mul "i" (complex log "i")));
-> "[tt]" (complex tan (complex atan [z]));
-> "[sh]" (complex sinh (complex mul "i" [z]));
concat "sin^2+cos^2 = " [one] "; log i = " (complex log "i") ", i^i = " [ii] "; tan atan z = " [tt] "; sinh iz = " [sh] " = i sin z = " (complex mul "i" [s]) "; cosh i pi = " (complex cosh "3.141592653589793i") "; acos 2 = " (complex acos 2)

Result

mentdb

sin^2+cos^2 = 1.0000000000000027; log i = 1.5707963267948966i, i^i = 0.20787957635076193; tan atan z = 1.0000000000000002+2i; sinh iz = -1.959601041421606+3.165778513216168i = i sin z = -1.959601041421606+3.165778513216168i; cosh i pi = -1; acos 2 = 1.3169578969248164i

Go further

Map a circle through complex log and watch the branch cut on the negative real axis.
Compare complex tanh with -i tan(iz).

executed on MentDB

Machine epsilon and catastrophic cancellation

Measure the spacing of doubles (ulp) at 1, 10^6 and 10^16, then show why exp(x) - 1 and ln(1 + x) must be computed with expm1 and log1p when x = 10^-12.

Pipeline

01 · ulpSpacing of the doublesmath ulp→02 · naiveexp then subtractmath log→03 · stableStable expm1math expm1→04 · stableStable log1pmath log1p

Script

admin

-> "[x]" 1E-12;
#1. Naive formulas lose 4 digits;
-> "[naive]" (math calc "exp([x])-1");
-> "[lnaive]" (math log (math calc "1+[x]"));
#2. Stable library functions;
-> "[good]" (math expm1 [x]);
-> "[lgood]" (math log1p [x]);
concat "ulp 1 = " (math ulp 1) ", ulp 1E6 = " (math ulp 1000000) ", ulp 1E16 = " (math ulp 1E16) "; exp x - 1: naive " [naive] " vs expm1 " [good] "; ln 1+x: naive " [lnaive] " vs log1p " [lgood] "; relative error naive " (math round_sig (math calc "abs([naive]-[good])/[good]") 2)

Result

mentdb

ulp 1 = 2.220446049250313E-16, ulp 1E6 = 1.1641532182693481E-10, ulp 1E16 = 2; exp x - 1: naive 1.000088900582341E-12 vs expm1 1.0000000000005E-12; ln 1+x: naive 1.000088900581841E-12 vs log1p 9.999999999995E-13; relative error naive 0.000089

Go further

The relative error of the naive formula is about ulp(1)/x: try x = 1E-8 and 1E-15.
Use bigdec functions to get a reference value with 40 digits.

executed on MentDB

Winning move in the game of Nim

Bouton's theorem: a Nim position is lost for the player to move when the XOR of the heap sizes is 0. With heaps 5, 3, 4, compute the nim-sum, the target size of each heap, and a few bitwise facts (Gray code of 13).

Pipeline

01 · sortHeaps in ordermath order_int→02 · nimsumNim-sum = XOR of the heapsmath bit_xor→03 · bitsAND / OR of heap sizesmath bit_and→04 · grayGray code n XOR n/2math trunc

Script

admin

#Nim with heaps 5 3 4;
-> "[h]" (math order_int "[5, 3, 4]");
-> "[s]" (math bit_xor (math bit_xor 3 4) 5);
#A heap h can move to h XOR s when it is smaller;
-> "[m3]" (math bit_xor 3 [s]);
-> "[m4]" (math bit_xor 4 [s]);
-> "[m5]" (math bit_xor 5 [s]);
concat "heaps " [h] ", nim-sum = " [s] "; targets 3->" [m3] " 4->" [m4] " 5->" [m5] " so take 2 from the heap of 3; 3 AND 5 = " (math bit_and 3 5) ", 3 OR 4 = " (math bit_or 3 4) ", Gray code of 13 = " (math bit_xor 13 (math trunc (math calc "13/2")))

Result

mentdb

heaps [3,4,5], nim-sum = 2; targets 3->1 4->6 5->7 so take 2 from the heap of 3; 3 AND 5 = 1, 3 OR 4 = 7, Gray code of 13 = 11

Go further

Only the heaps whose size has the leading bit of the nim-sum set give a legal winning move.
Sprague-Grundy: the same XOR rule combines any impartial games.