Let ℍ be the upper half plane, the set of complex real numbers with positive imaginary part. When we measure distances the way we’ve discussed in the last couple posts, the geometry of ℍ is hyperbolic. What is a circle of radius _r_ in ℍ? The same as a circle in any geometry: it’s the set of points a fixed distance _r_ from a center. But when you draw a circle using one metric, it may look very different when viewed from the perspective of another metric. Suppose we put on glasses that gave us a hyperbolic perspective on ℍ, draw a circle of radius _r_ centered at _i_ , then take off the hyperbolic glasses and put on Euclidean glasses. What would our drawing look like? In the previous post we gave several equivalent expressions for the hyperbolic metric. We’ll use the first one here. Here the Fraktur letter ℑ stands for imaginary part. So the set of points in a circle of radius _r_ centered at _i_ is Divide the expression for _d_(_x_ + _iy_ , _i_) by 2, apply sinh, and square. This gives us which is an equation for a Euclidean circle. If we multiply both sides by 4 _y_ and complete the square, we find that the center of the circle is (0, cosh(_r_)) and the radius is sinh(_r_). So to recap, if we put on our hyperbolic glasses and draw a circle, then switch out these glasses for Euclidean glasses, the figure we drew again looks like a circle. To put it another way, a hyperbolic viewer and a Euclidean viewer would agree that a circle has been draw. However, the two viewers would disagree where the center of the circle is located, and they would disagree on the radius. Both could agree that the center is on the imaginary axis, but the hyperbolic viewer would say the imaginary part of the center is 1 and the Euclidean viewer would say it’s cosh(_r_). The hyperbolic observer would say the circle has radius _r_ , but the Euclidean observer would say it has radius sinh(_r_). For small _r_ , they nearly agree. But for large _r_ , the circle appears much larger to the Euclidean observer. Note that if you asked a Euclidean observer to draw a circle of radius 100, centered at (0, 1), he would say that the circle will extend outside of the half plane. A hyperbolic observer would disagree. From his perspective, the real axis is infinitely far away and so he can draw a circle of any radius centered at any point and stay within the half plane.

One common model of the hyperbolic plane is the Poincaré upper half plane ℍ. This is the set of points in the complex plane with positive imaginary part. Straight lines are either vertical, a set of points with constant imaginary part, or arcs of circles centered on the real axis. The real axis is not part of ℍ. From the perspective of hyperbolic geometry these are ideal parts, infinitely far away, and not part of the plane itself. We can define a metric on ℍ as follows. To find the distance between two points _u_ and _v_ , draw a line between the two points, and let _a_ and _b_ be the ideal points at the end of the line. By a line we mean a line as defined in the geometry of ℍ, what we would see from our Euclidean perspective as a half circle or a vertical line. Then the distance between _u_ and _v_ is defined as the absolute value of the log of the cross ratio (_u_ , _v_ ; _a_ , _b_). Cross ratios are unchanged by Möbius transformations, and so Möbius transformations are isometries. Another common model of hyperbolic geometry is the Poincaré disk. We can use the same metric on the Poincaré disk because the Möbius transformation maps the upper half plane to the unit disk. This is very similar to how the Smith chart is created by mapping a grid in the right half plane to the unit disk.

I was reading about Shackleton’s incredible expedition to Antarctica, and the Weddell Sea features prominently. That name sounded familiar, and I was trying to remember where I’d heard of Weddell in math. I figured out that it wasn’t Weddell exactly but Weddle I was thinking of. The Weddell Sea is named after James Weddell (1787–1834). Weddle’s integration rule is named after Thomas Weddle (1817–1853). I wrote about Weddle’s integration rule a couple years ago. Weddle’s rule, also known as Bode’s rule, is as follows. Let’s try this on integrating sin(_x_) from 1 to 2. If we divide the interval [1, 2] into 6 subintervals, _h_ = 1/6. The 8th derivative of sin(_x_) is also sin(_x_), so it is bounded by 1. So we would expect the absolute value of the error to be bounded by 9 / (69 × 1400). Let’s see what happens in practice. import numpy as np x = np.linspace(1, 2, 7) h = (2 - 1)/6 weights = (h/140)*np.array([41, 216, 27, 272, 27, 216, 41]) approx = np.dot(weights, np.sin(x)) exact = np.cos(1) - np.cos(2) print("Error: ", abs(approx - exact) ) print("Expected error: ", 9/(1400*6**9)) Here’s the output: Error: 6.321198009473505e-10 Expected error: 6.379009079626363e-10 ## Related posts * Pushing numerical integration software to its limits * Approximating rapidly diverging integrals * Integration by darts

I mentioned in the previous post that the harmonic numbers _H_ _n_ are never integers for _n_ > 1. In the spirit of that post, I’d like to find the value of _n_ such that _H_ _n_ is closest to a given integer _m_. We have two problems to solve. First, how do we accurately and efficiently compute harmonic numbers? For small _n_ we can directly implement the definition. For large _n_ , the direct approach would be slow and would accumulate floating point error. But in that case we could use the asymptotic approximation from this post. As is often the case, the direct approach gets worse as _n_ increases, but the asymptotic approximation gets better as _n_ increases. Here γ is the Euler-Mascheroni constant. The second problem to solve is how to find the value of _n_ so that _H_ _n_ comes closest to _m_ without trying too many possible values of _n_? We can discard the higher order terms above and see that _n_ is roughly exp(_m_ − γ). Here’s the code. import numpy as np gamma = 0.57721566490153286 def H(n): if n < 1000: return sum([1/k for k in range(1, n+1)]) else: n = float(n) return np.log(n) + gamma + 1/(2*n) - 1/(12*n**3) # return n such that H_n is closest harmonic number to m def nearest_harmonic_number(m): if m == 1: return 1 guess = int(np.exp(m - gamma)) if H(guess) < m: i = guess while H(guess) < m: guess += 1 j = guess else: j = guess while H(guess) > m: guess -= 1 i = guess x = np.array([abs(H(k) - m) for k in range(i, j+1)]) return i + np.argmin(x) We can use this, for example, to find the closest harmonic number to 10. >>> nearest_harmonic_number(10) 12366 >>> H(12366) 9.99996214846655 I wrote the code with integer values of _m_ in mind, but the code works fine with real numbers. For example, we could find the harmonic number closest to √20. >>> nearest_harmonic_number(20**0.5) 49 >>> H(49)**2 20.063280462918804 ## Related posts * Harmonic _e_ * A strange take on the harmonic series * Computing γ

The last couple posts have been about group pairings, specifically Tate pairings as they’re used in cryptography. This post will show that RSA encryption can be seen as a special case of pairing-based cryptography. The idea comes from Ben Lynn’s 2007 dissertation. Lynn is the “L” in BLS signatures—one of the topics in his disserations—and in BLS elliptic curves. A pairing is a bilinear mapping from two groups to a third group _e_ : _G_ 1 × _G_ 2 → _G_ T. Here bilinear means that if _P_ is an element of _G_ 1 and _Q_ is an element of _G_ 2, and _a_ and _b_ are nonnegative integers, then _e_(_aP_ , _bQ_) = _e_(_P_ , _Q_)_ab_. There are more criteria for a pairing to be useful in cryptography, but we won’t need those for this post. Ben Lynn’s disseration mentions that exponentiation is a special case of pairing if you let _G_ 1 and _G_ T be the multiplicative group of the integers mod _r_ and let _G_ 2 be the additive group of integers mod (_r_ − 1). Then you can define a pairing by _e_(_g_ , _a_) = _g_ _a_. Typically you can’t just write down a simple expression for a pairing, but in this case you can. RSA encryption corresponds to _r_ = _pq_ where _p_ and _q_ are large primes. The product _pq_ is made public but the factorization into _p_ and _q_ is held secret. A message [1] is encrypted by exponentiation mod _n_ where the exponent is the public key. In Lynn’s notation, the message is _g_ and the public key is _a_. The security of RSA encryption depends on the fact that you can’t recover _g_ from _g_ _a_ mod _n_ unless you know a trapdoor, the factorization of _n_ [2]. This is true of pairings more generally: it is not practical to recover the inputs to a pairing from the output unless you know a trapdoor. ## Related posts * A quiet change to RSA * My interview with John Tate * More posts on cryptography [1] In practice, RSA isn’t used to encrypt entire messages. Instead, it is used to encrypt a key for a symmetric encryption algorithm such as AES, and that key is used to encrypt the message. This is done for efficiency. [2] Or, more specifically, a private key that can easily be computed if you know the factorization of _n_. It’s conceivable that breaking RSA encryption is easier than factoring, but so far that does not appear to be the case.

## Elliptic curve Diffie-Hellman Given a point _P_ on an elliptic curve _E_ , and a random number _a_ , _aP_ means to add _P_ to itself _a_ times, using the addition on _E_. The point _aP_ can be computed efficiently, even if _a_ is a very large number [1]. However, if _E_ has a large number of points, and if _a_ is chosen at random from a large range, then it is not practical to compute _a_ given _P_ and _aP_. This is the elliptic curve version of the discrete logarithm problem, and its presumed difficulty is the basis of the security of Diffie-Hellman key exchange. ## Two-party Diffie-Hellman With two-party Diffie-Hellman key exchange, two parties, Alice and Bob, generate random private keys _a_ and _b_ respectively. They agree on a point _P_ on an elliptic curve _E_. Alice computes _aP_ and sends it to Bob. Simultaneously Bob computes _bP_ and sends it to Alice. Then Alice can compute _a_(_bP_) = (_ab_)_P_ and Bob can compute _b_(_aP_) = (_ba_)_P_ = (_ab_)_P_. Then both Alice and Bob know a shared secret, the point (_ab_)_P_ on _E_ , but neither party has revealed a private key. ## Three-party Diffie-Hellman You could extend the approach above to three parties, say adding Carol, but this would require extra communication: Alice could send (_ab_)_P_ to Carol, which she could multiply by her private key _c_ to obtain _abcP_. Similarly, everyone else could arrive at _abcP_. Each person has to do a computation, send and receive a message, do another computation, and send an receive another message. Joux [2] came up with a way to do Diffie-Hellman key exchange with three people and only one round of sending and receiving messages. The set up uses a pairing _e_( , ) of two elliptic curve subgroups, _G_ 1 and _G_ 2, as in the previous post. Fix generators _P_ ∈ _G_ 1 and _Q_ ∈ _G_ 2. Each party multiplies _P_ and _Q_ by their private key and sends the results to the other two parties. Alice receives _bP_ from Bob and _cQ_ from Carol. This is enough for her to compute _e_(_bP_ , _cQ_)_a_ = _e_(_P_ , _Q_)_abc_. Similarly, Bob receives _aP_ from Alice and _cQ_ from Carol, enabling him to compute _e_(_aP_ , _cQ_)_b_ = _e_(_P_ , _Q_)_abc_. And finally, Carol receives _aP_ from Alice and _bQ_ from Bob, enabling her to compute _e_(_aP_ , _bQ_)_c_ = _e_(_P_ , _Q_)_abc_. So all three parties can compute the shared secret _e_(_P_ , _Q_)_abc_. but no party knows the other parties’ private keys. ## Footnotes [1] If you want to multiply a point by 2100, for example, you don’t carry out 2100 additions; you carry out 100 doublings. Of course not every positive integer is a power of 2, but every positive integer is the _sum_ of powers of 2, i.e. it can be written in binary. So as you’re doing your doublings, sum the terms that correspond to 1s in the binary representation of the number you’re multiplying by. [2] Antoine Joux. A One Round Protocol for Tripartite Diffie–Hellman. Journal of Cryptology (2004) 17: 263–276.

Let _p_ be a prime number. Then the integers mod _p_ form a finite field. The number of elements in a finite field must be a power of a prime, i.e. the order _q_ = _p_ _n_ for some _n_. When _n_ > 1, we can take the elements of our field to be polynomials of degree _n_ − 1 with coefficients in the integers mod _p_. Addition works just as you’d expect addition to work, adding coefficients mod _p_ , but multiplication is a little more complicated. You multiply field elements by multiplying their polynomial representatives, but then you divide by an irreducible polynomial and take the remainder. When _n_ = 2, for some _p_ you can define the field by adding an imaginary unit. ## When you can and cannot adjoin an _i_ For some finite fields of order _p_ , you can construct a field of order _p_ ² by joining an element _i_ to the field, very much the way you form the complex numbers from the real numbers. For example, you can create a field with 49 elements by taking pairs of (_a_ , _b_) of integers mod 7 and multiplying them as if they were _a_ + _bi_. So (_a_ , _b_) * (_c_ , _d_) = (_ac_ − _bd_ , _ad_ + _bc_). This is equivalent to choosing the polynomial _x_ ² + 1 as your irreducible polynomial and following every polynomial multiplication by taking the remainder modulo _x_ ² + 1. This works for a field with 49 elements, but not for a field of 25 elements. That’s because over the integers mod 5 the polynomial _x_ ² + 1 already has a root. Two of them in fact: _x_ = 2 or _x_ = 3. So you could say that mod 5, _i_ = 2. Or _i_ = 3 if you prefer. You can still form a field of 25 elements by taking pairs of elements from a field of 5 elements, but you have to choose a different polynomial as your irreducible polynomial because _x_ ² + 1 is **not** irreducible because _x_ ² + 1 = (_x_ − 2)(_x_ + 2) when working over the integers mod 5. You could use _x_ ² + _x_ + 1 as your irreducible polynomial. To prove that this polynomial is irreducible mod 5, plug in the numbers 0, 1, 2, 3, and 4 and confirm that none of them make the polynomial equal 0. In general, you can create a field of order _p_ ² by adjoining an element _i_ if and only if _p_ = 3 mod 4. Next we’ll look at an example of making a very large finite field even larger by adding an imaginary element. ## Example from Ethereum The Ethereum virtual machine has support for a pairing—more on that in a future post—of two elliptic curves, **BN254** and **alt_bn128**. The BN254 curve is defined by _y_ ² = _x_ ³ + 3 over the field _F p_, the integers mod _p_ , where _p_ = 21888242871839275222246405745257275088696311157297823662689037894645226208583. The curve alt_bn128 is defined by _y_ ² = _x_ ³ + 3/(9 + _i_) over the field _F p_[_i_], i.e. the field _F p_, with an element _i_ adjoined. Note the that last two digits of _p_ are 83, and so _p_ is congruent to 3 mod 4. ## Special point on curve The Ethereum documentation (EIP-197) singles out a particular point (_x_ , _y_) on alt_bn128: _x_ = _a_ + _bi_ _y _=_c_ + _di___ where _a_ = 10857046999023057135944570762232829481370756359578518086990519993285655852781 _b_ = 11559732032986387107991004021392285783925812861821192530917403151452391805634 _c_ = 8495653923123431417604973247489272438418190587263600148770280649306958101930 _d_ = 4082367875863433681332203403145435568316851327593401208105741076214120093531. We will show that this point is on the curve as an exercise in working in the field _F p_[_i_]. We’ll write Python code from scratch, not using any libraries, so all the details will be explicit. def add(pair0, pair1, p): a, b = pair0 c, d = pair1 return ((a + c) % p, (b + d) % p) def mult(pair0, pair1, p): a, b = pair0 c, d = pair1 return ((a*c - b*d) % p, (b*c + a*d) % p) p = 21888242871839275222246405745257275088696311157297823662689037894645226208583 a = 10857046999023057135944570762232829481370756359578518086990519993285655852781 b = 11559732032986387107991004021392285783925812861821192530917403151452391805634 c = 8495653923123431417604973247489272438418190587263600148770280649306958101930 d = 4082367875863433681332203403145435568316851327593401208105741076214120093531 # Find (e, f) such that (e, f)*(9, 1) = (1, 0). # 9e - f = 1 # e + 9f = 0 # Multiply first equation by 9 and add. e = (9*pow(82, -1, p)) % p f = (-e*pow(9, -1, p)) % p prod = mult((e, f), (9, 1), p) assert(prod[0] == 1 and prod[1] == 0) y2 = mult((c, d), (c, d), p) x3 = mult((a, b), mult((a, b), (a, b), p), p) rhs = add(x3, mult((3, 0), (e, f), p), p) assert(y2[0] == rhs[0]) assert(y2[1] == rhs[1]) ## Related posts * Finite fields * Elliptic curve in Ethereum’s consensus layer * Misleading plots of elliptic curves

Suppose you have $100 to invest in two independent assets, _A_ and _B_ , and you want to minimize volatility. Suppose _A_ is more volatile than _B_. Then putting all your money on _A_ would be the worst thing to do, but putting all your money on _B_ would not be the best thing to do. The optimal allocation would be some mix of _A_ and _B_ , with more (but not all) going to _B_. We will formalize this problem and determine the optimal allocation, then generalize the problem to more assets. ## Two variables Let _X_ and _Y_ be two independent random variables with finite variance and assume at least one of _X_ and _Y_ is not constant. We want to find _t_ that minimizes subject to the constraint 0 ≤ _t_ ≤ 1. Because _X_ and _Y_ are independent, Taking the derivative with respect to _t_ and setting it to zero shows that So the smaller the variance on _Y_ , the less we allocate to _X_. If _Y_ is constant, we allocate nothing to _X_ and go all in on _Y_. If _X_ and _Y_ have equal variance, we allocate an equal amount to each. If _X_ has twice the variance of _Y_ , we allocate 1/3 to _X_ and 2/3 to _Y_. ## Muliple variables Now suppose we have _n_ independent random variables _X_ _i_ for _i_ running from 1 to _n_ , and at least one of the variables is not constant. Then we want to minimize subject to the constraint and all _t_ _i_ non-negative. We can solve this optimization problem with Lagrange multipliers and find that for all 1 ≤ _i_ , _j_ ≤ _n_. These (_n_ − 1) equations along with the constraint that all the _t_ _i_ sum to 1 give us a system of equations whose solution is Incidentally, the denominator has a name: the (_n_ − 1)st elementary symmetric polynomial in _n_ variables. Here the variables are the variances. ## Related posts * Symmetric funcions and U-statistics * Regular solids and Monte Carlo integration

The area of a triangle can be computed directly from the lengths of its sides via Heron’s formula. Here _s_ is the semiperimeter, _s_ = (_a_ + _b_ + _c_)/2. Is there an analogous formula for spherical triangles? It’s not obvious there should be, but there is a formula by Simon Antoine Jean L’Huilier (1750–1840). Here we denote area by _S_ for surface area, rather than _A_ because in the context of spherical trigonometry _A_ usually denotes the angle opposite side _a_. Now tan θ ≈ θ for small θ, and so L’Huilier’s formula reduces to Heron’s formula for small triangles. Imagine the Earth as a sphere of radius 1 and take a spherical triangle with one vertex at the north pole and two vertices on the equator 90° longitude apart. Such a triangle takes of 1/8 of the Earth’s surface area of 4π, so the area should be π/2. You can verify that L’Huilier’s formula gives the correct area. It’s not a proof, but it’s a good sanity check that L’Huilier’s formula is correct for small triangles and for at least one big triangle.

The cross ratio of four points _A_ , _B_ , _C_ , _D_ is defined by where _XY_ denotes the length of the line segment from _X_ to _Y_. The idea of a cross ratio goes back at least as far as Pappus of Alexandria (c. 290 – c. 350 AD). Numerous theorems from geometry are stated in terms of the cross ratio. For example, the cross ratio of four points is unchanged under a projective transformation. ## Complex numbers The cross ratio of four (extended [1]) complex numbers is defined by The absolute value of the complex cross ratio is the cross ratio of the four numbers as points in a plane. The cross ratio is invariant under Möbius transformations, i.e. if _T_ is any Möbius transformation, then This is connected to the invariance of the cross ratio in geometry: Möbius transformations are projective transformations on a complex projective line. (More on that here.) If we fix the first three arguments but leave the last argument variable, then is the unique Möbius transformation mapping _z_ 1, _z_ 2, and _z_ 3 to ∞, 0, and 1 respectively. ## The anharmonic group Suppose (_a_ , _b_ ; _c_ , _d_) = λ ≠ 1. Then there are 4! = 24 permutations of the arguments and 6 corresponding cross ratios: Viewed as functions of λ, these six functions form a group, generated by This group is called the anharmonic group. Four numbers are said to be in harmonic relation if their cross ratio is 1, so the requirement that λ ≠ 1 says that the four numbers are anharmonic. The six elements of the group can be written as ## Hypergeometric transformations When I was looking at the six possible cross ratios for permutations of the arguments, I thought about where I’d seen them before: the linear transformation formulas for hypergeometric functions. These are, for example, equations 15.3.3 through 15.3.9 in A&S. They relate the hypergeometric function _F_(_a_ , _b_ ; _c_ ; _z_) to similar functions where the argument _z_ is replaced with one of the elements of the anharmonic group. I’ve written about these transformations before here. For example, There are deep relationships between hypergeometric functions and projective geometry, so I assume there’s an elegant explanation for the similarity between the transformation formulas and the anharmonic group, though I can’t say right now what it is. ## Related posts * Projective duality * Finite projective planes * Area of the unit disk after a Möbius transformation [1] For completeness we need to include a point at infinity. If one of the _z_ equals ∞ then the terms involving ∞ are dropped from the definition of the cross ratio.