xii. Number Theory Miscellany

Smaller number-theoretic tools: intersections of arithmetic progressions, coprime bases, Pythagorean triples, counting square roots of 1 and $k$-th powers, residues of $ai+b$ in a range, Fibonacci numbers, power towers, balanced ternary and NAF, elliptic curves over ${\mathbb{F}}_p$, and Gaussian integers. Dependencies from other notes: extgcd, invMod (GCD), crtMerge, crt, garner, diophantine, countSolutions (Congruences and CRT), floorSum (Floor Sums), power (Modular Arithmetic, used by PowerTower and twoSquares); the simplest fraction in an interval is simplestBetween in Big Integers and Fractions.

Function	Answers	Cost
apIntersection(a1, d1, a2, d2)	$\{a_1+d_{1}k\}\cap \{a_2+d_{2}k\}$ ($k\ge 0$) as first + step * t	$O(\log )$
coprimeBasis(vals)	pairwise coprime numbers that generate every value multiplicatively	$O(k^2\log V)$
forEachPythagorean(N, f)	all primitive $a^2+b^2=c^2$, $c\le N$	$O(N)$
countSquareRootsOfOne(m)	$\#\{x\in [0,m):x^2\equiv 1\}$	$O(\sqrt{m})$
countKthPowers(n, k)	$\#\{a^k\mathrm{mod}n:a\in \mathbb{Z}\}$	$O(\sqrt{n})$
countResiduesInRange(a, b, m, n, lo, hi)	$\#\{i\in [0,n):lo\le (ai+b)\mathrm{mod}m<hi\}$	$O(\log m)$ by floor sums

pair<ll, ll> apIntersection(ll a1, ll d1, ll a2, ll d2) {   // d1, d2 >= 1; returns {first, step}, first = -1 if empty
    auto [x, l] = crtMerge(((a1 % d1) + d1) % d1, d1, ((a2 % d2) + d2) % d2, d2);
    if (x < 0) return {-1, -1};
    ll lo = max(a1, a2);                                   // first common element >= both starting points
    if (x < lo) x += (lo - x + l - 1) / l * l;
    return {x, l};
}
 
vector<ll> coprimeBasis(const vector<ll>& vals) {          // values >= 1
    vector<ll> basis;
    for (ll v : vals) {
        vector<ll> work{v};
        while (!work.empty()) {
            ll x = work.back(); work.pop_back();
            if (x == 1) continue;
            bool split = false;
            for (int i = 0; i < sz(basis); i++) {
                ll g = gcd(basis[i], x);
                if (g == 1) continue;
                ll b = basis[i]; basis.erase(basis.begin() + i);
                work.push_back(g); work.push_back(b / g); work.push_back(x / g);
                split = true; break;
            }
            if (!split) basis.push_back(x);
        }
    }
    sort(all(basis)); basis.erase(unique(all(basis)), basis.end());
    return basis;
}
 
template<class F> void forEachPythagorean(ll N, F f) {   // primitive triples (a, b, c), c <= N
    for (ll m = 2; m * m + 1 <= N; m++)
        for (ll k = 1 + (m % 2); k < m && m * m + k * k <= N; k += 2)
            if (gcd(m, k) == 1) f(m * m - k * k, 2 * m * k, m * m + k * k);
}
 
ll countSquareRootsOfOne(ll m) {
    ll res = 1;
    for (ll p = 2; p * p <= m; p++) {
        if (m % p) continue;
        int e = 0; while (m % p == 0) m /= p, e++;
        res *= p == 2 ? (e == 1 ? 1 : e == 2 ? 2 : 4) : 2;
    }
    if (m > 1) res *= m == 2 ? 1 : 2;
    return res;
}
 
ll countKthPowers(ll n, ll k) {                            // k >= 1; values of a^k mod n, including 0 and non-units
    auto ipow = [](ll b, ll e) { ll r = 1; while (e-- > 0) r *= b; return r; };
    auto f = [&](auto&& self, ll p, ll e) -> ll {          // count modulo p^e
        if (e <= 0) return 1;                              // only 0
        if (p == 2) {
            if (e == 1) return 2;
            ll u = ipow(2, e - 2) / gcd(k, ipow(2, e - 2));
            if (k % 2) u *= 2;
            return u + self(self, 2, e - k);
        }
        ll phi = ipow(p, e) - ipow(p, e - 1);
        return phi / gcd(k, phi) + self(self, p, e - k);
    };
    ll res = 1, x = n;
    for (ll p = 2; p * p <= x; p++) { if (x % p) continue; int e = 0; while (x % p == 0) x /= p, e++; res *= f(f, p, e); }
    if (x > 1) res *= f(f, x, 1);
    return res;
}
 
ll countResiduesInRange(ll a, ll b, ll m, ll n, ll lo, ll hi) {   // 0 <= lo <= hi <= m
    a = ((a % m) + m) % m; b = ((b % m) + m) % m;
    auto below = [&](ll t) -> ll {                         // #i with (a i + b) mod m < t
        if (t <= 0) return 0;
        if (t >= m) return n;
        return n - (floorSum(n, m, a, b + m - t) - floorSum(n, m, a, b));
    };
    return below(hi) - below(lo);
}

Proof: the miscellany functions

AP intersection. $x$ lies in both progressions (ignoring $k\ge 0$) iff $x\equiv a_1(\mathrm{mod}{d}_1)$ and $x\equiv a_2(\mathrm{mod}{d}_2)$, which by the generalized CRT is empty or one class modulo $\mathrm{lcm}(d_1,d_2)$; the smallest member $\ge \max (a_1,a_2)$ is the first common element with $k\ge 0$ in both.

Coprime basis. Invariant: the product of all numbers in basis and work (with multiplicity) generates, multiplicatively, every processed value, and basis is pairwise coprime. Splitting $b$ and $x$ with $g=\gcd (b,x)>1$ into $g,b/g,x/g$ keeps generation ($b=g\cdot \frac{b}{g}$, $x=g\cdot \frac{x}{g}$) and replaces two numbers of product $bx$ by three of product $\frac{bx}{g}<bx$ (ones are dropped), so the product of everything in basis and work strictly decreases and the process terminates; an element is appended only when it is coprime to the whole basis.

Pythagorean triples (Euclid). For a primitive triple, exactly one of $a,b$ is even (squares are $0,1\mathrm{mod}4$, so two odd legs would make $c^2\equiv 2$); with $b$ even, $\frac{b^2}{4}=\frac{c-a}{2}\cdot \frac{c+a}{2}$ with coprime factors, so both are squares $k^2,m^2$: $a=m^2-k^2$, $b=2mk$, $c=m^2+k^2$ with $m>k$, $\gcd (m,k)=1$, $m-k$ odd; conversely these are primitive. Each triple arises once.

Square roots of 1. By the CRT, count per prime power. For odd $p^e$, $(x-1)(x+1)\equiv 0$ and $p$ cannot divide both factors (their difference is 2), so $x\equiv \pm 1$: 2 roots. Modulo $2^e$: 1 root for $e=1$, 2 for $e=2$ ($\pm 1$), and 4 for $e\ge 3$ ($\pm 1$, $2^{e-1}\pm 1$), since both factors are even and one must be divisible by $2^{e-1}$.

$k$-th powers modulo $p^e$. Write $a=p^{j}u$ with $u$ a unit. Units: for odd $p$ the unit group is cyclic of order $\phi =p^{e-1}(p-1)$, and the image of $x\mapsto x^k$ has $\frac{\phi }{\gcd (k,\phi )}$ elements. For $p=2$, $e\ge 2$, the units are $\{\pm 1\}\times ⟨5⟩\cong {\mathbb{Z}}_2\times {\mathbb{Z}}_{2^{e-2}}$: the image has $\frac{2^{e-2}}{\gcd (k,2^{e-2})}$ elements times 2 if $k$ is odd (the $\pm$ part survives). Multiples of $p$: $a=pb$ gives $a^k=p^k{b}^k$, and multiplication by $p^k$ is injective from ${\mathbb{Z}}_{p^{e-k}}$ into ${\mathbb{Z}}_{p^e}$, so these values correspond one-to-one to the $k$-th powers modulo $p^{e-k}$ (all residues there, including 0): the recursion self(p, e - k), with $p^{e-k}\le 1$ leaving only 0. By the CRT the counts multiply over prime powers. Checked by brute force for all $n\le 300$ and $k\le 8$.

Residues in a range. For an integer $x\ge 0$ and $0<t\le m$: $\lfloor \frac{x+m-t}{m}\rfloor -\lfloor \frac{x}{m}\rfloor$ is 1 if $x\mathrm{mod}m\ge t$ and 0 otherwise. Summing over $x=ai+b$, $i<n$, gives $\#\{(ai+b)\mathrm{mod}m\ge t\}$ as a difference of two floor sums; subtract from $n$. All functions were checked against brute force on random inputs.

Fibonacci numbers

$F_0=0$, $F_1=1$, $F_n=F_{n-1}+F_{n-2}$; $F_n\approx {\phi }^n/\sqrt{5}$ ($\phi =\frac{1+\sqrt{5}}{2}\approx 1.618$); $F_{92}$ is the last that fits in long long, $F_{93}$ the last in unsigned long long.

Method	Time	Use
iteration	$O(n)$	$n\le 10^7$
matrix power ${\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^n=\left[\begin{array}{cc}{F}_{n+1}& F_n\\ F_n& F_{n-1}\end{array}\right]$	$O(\log n)$	$n\le 10^{18}$
fast doubling	$O(\log n)$	same, about 2 times faster, no matrix struct
Binet $F_n=\frac{{\phi }^n-{\psi }^n}{\sqrt{5}}$	$O(1)$	inexact beyond $n\approx 70$

Identity
${\sum }_{i=0}^n{F}_i=F_{n+2}-1$
${\sum }_{i=0}^n{F}_i^2=F_n{F}_{n+1}$
$F_{m+n}=F_m{F}_{n+1}+F_{m-1}{F}_n$	addition formula
$F_{2k}=F_k(2F_{k+1}-F_k)$, $F_{2k+1}=F_k^2+F_{k+1}^2$	doubling
$\gcd (F_m,F_n)=F_{\gcd (m,n)}$, hence $F_n\mid F_{kn}$	strong divisibility: common factors of Fibonacci numbers become gcds of indices
$F_{n-1}{F}_{n+1}-F_n^2=(-1{)}^n$	Cassini
Lucas numbers $L_0=2$, $L_1=1$: $L_n=F_{n-1}+F_{n+1}$, $F_{2n}=F_n{L}_n$

Pisano period $\pi (m)$ ($F_n\mathrm{mod}m$ is periodic)
$\pi (\mathrm{lcm}(m_1,m_2))=\mathrm{lcm}(\pi (m_1),\pi (m_2))$	multiplicative structure
$\pi (m)\le 6m$, with equality iff $m=2\cdot 5^k$	(known theorem)
$\pi (p^k)=p^{k-1}\pi (p)$	conjectured always (Wall), verified far
$p\equiv \pm 1(\mathrm{mod}5)$: $\pi (p)\mid p-1$; $p\equiv \pm 2(\mathrm{mod}5)$: $\pi (p)\mid 2(p+1)$	from $\phi \in {\mathbb{F}}_p$ or ${\mathbb{F}}_{p^2}$
$\pi (10)=60$, $\pi (100)=300$, $\pi (1000)=1500$	so $F_n\mathrm{mod}m=F_{n\mathrm{mod}\pi (m)}\mathrm{mod}m$

Zeckendorf: every positive integer is uniquely a sum of non-consecutive Fibonacci numbers (greedy: take the largest that fits); it underlies Wythoff’s game, Fibonacci coding (self-delimiting codes) and “no two adjacent terms” representations.

Where Fibonacci appears	Count
binary strings of length $n$ with no two adjacent 1s; subsets of $\{\mathrm{1..}n\}$ with no two consecutive	$F_{n+2}$
domino tilings of $2\times n$; stairs climbed 1 or 2 at a time	$F_{n+1}$
worst case of the Euclidean algorithm	consecutive Fibonacci inputs
golden-section search, minimum AVL tree of height $h$	ratio $\phi$, $F_{h+3}-1$ nodes
tribonacci / $k$-bonacci	$O(k^3\log n)$ matrix power, or $O(k\log k\log n)$ Bostan–Mori: Linear Recurrences

pair<ll, ll> fibPair(ll n, ll mod) {                       // {F(n) mod m, F(n+1) mod m}, fast doubling
    if (n == 0) return {0, 1 % mod};
    auto [a, b] = fibPair(n >> 1, mod);
    ll c = (ll)((i128)a * ((2 * b - a + mod) % mod) % mod);  // F(2k)
    ll d = (ll)(((i128)a * a + (i128)b * b) % mod);           // F(2k+1)
    return n & 1 ? pair{d, (c + d) % mod} : pair{c, d};
}
 
vector<ll> zeckendorf(ll n) {                              // n >= 1: non-consecutive Fibonacci terms, decreasing
    vector<ll> f{1, 2}; while (f.back() <= n - f[sz(f) - 2]) f.push_back(f.back() + f[sz(f) - 2]);
    vector<ll> res;
    for (int i = sz(f) - 1; i >= 0 && n > 0; i--) if (f[i] <= n) res.push_back(f[i]), n -= f[i];
    return res;
}

Proof: the doubling formulas, strong divisibility, Cassini, sums and Zeckendorf

Addition and doubling. $M^n=\left[\begin{array}{cc}{F}_{n+1}& F_n\\ F_n& F_{n-1}\end{array}\right]$ by induction; comparing entries of $M^{m+n}=M^m{M}^n$ gives $F_{m+n}=F_m{F}_{n+1}+F_{m-1}{F}_n$. With $m=n=k$: $F_{2k}=F_k{F}_{k+1}+F_{k-1}{F}_k=F_k(2F_{k+1}-F_k)$; with $m=k+1$, $n=k$: $F_{2k+1}=F_{k+1}^2+F_k^2$. Cassini: $\det M^n=(\det M{)}^n=(-1{)}^n$. Sums: ${\sum }_{i\le n}{F}_i=F_{n+2}-1$ by induction ($F_{n+2}-1+F_{n+1}=F_{n+3}-1$); $\sum F_i^2=F_n{F}_{n+1}$ since $F_n{F}_{n+1}-F_{n-1}{F}_n=F_n^2$.

Strong divisibility. Consecutive Fibonacci numbers are coprime (a common divisor of $F_n,F_{n+1}$ divides all earlier terms, down to $F_1=1$). From the addition formula, $F_{m+n}\equiv F_{m-1}{F}_n(\mathrm{mod}{F}_m)$, so $\gcd (F_{m+n},F_m)=\gcd (F_{m-1}{F}_n,F_m)=\gcd (F_n,F_m)$ (as $\gcd (F_{m-1},F_m)=1$). Thus $\gcd (F_a,F_b)$ follows the Euclidean algorithm on the indices and ends at $F_{\gcd (a,b)}$.

Zeckendorf. Existence: greedy takes the largest $F\le n$; the remainder is $<F_{\text{next}}-F=F_{\text{prev}}$, so the next term chosen is not adjacent. Uniqueness: a sum of non-consecutive Fibonacci numbers with largest term $F_j$ is less than $F_{j+1}$ (induction: $F_j+(\text{a sum below }{F}_{j-1})<F_j+F_{j-1}$); so the largest term of any representation of $n$ must be the largest Fibonacci number $\le n$, and induct. Pisano: $F_n\mathrm{mod}m$ is determined by the pair $(F_n,F_{n+1})\mathrm{mod}m$, which is invertible as a recurrence, so the sequence is purely periodic, and the period of the pair modulo $\mathrm{lcm}$ is the lcm of the periods (CRT). The $\pi (m)\le 6m$ bound and the divisibility facts are quoted; checked numerically for $m\le 2000$, together with the fast doubling and Zeckendorf code.

Power towers

$a_l^{a_{l+1}^{{\cdot }^{{\cdot }^{a_r}}}}\mathrm{mod}m$ for any modulus (no coprimality): use $a^x\equiv a^{(x\mathrm{mod}\phi (m))+\phi (m)}(\mathrm{mod}m)$ whenever $x\ge \phi (m)$, and decide “is the exponent at least $\phi (m)$?” from the exact value capped at $2^{62}$, which a tower of more than five numbers $\ge 2$ always exceeds. The modulus chain $m,\phi (m),\phi (\phi (m)),\dots$ reaches 1 within $O(\log m)$ levels. Classic: Codeforces 906D Power Tower. Cost: $O(\log m)$ levels, each $O(\sqrt{m})$ for $\phi$ (memoized) plus $O(\log )$ for the power. (The folklore version that stores ”$x$ if $x<m$, else $x\mathrm{mod}m+m$” can lose the "$\ge m$" information after an exponent has been reduced; carrying the capped exact value avoids that question.)

struct PowerTower {                                        // uses power(a, e, m) from Modular Arithmetic
    static constexpr ll CAP = 1LL << 62;
    vector<ll> a; vector<int> nextOne; map<ll, ll> phiMemo;
    PowerTower(const vector<ll>& a) : a(a), nextOne(sz(a) + 1, sz(a)) {   // a_i >= 1
        for (int i = sz(a) - 1; i >= 0; i--) nextOne[i] = a[i] == 1 ? i : nextOne[i + 1];
    }
    ll phi(ll n) {
        auto it = phiMemo.find(n); if (it != phiMemo.end()) return it->second;
        ll res = n, m = n;
        for (ll p = 2; p * p <= m; p++) if (m % p == 0) { while (m % p == 0) m /= p; res = res / p * (p - 1); }
        if (m > 1) res = res / m * (m - 1);
        return phiMemo[n] = res;
    }
    static ll satPow(ll b, ll e) { ll r = 1; while (e-- > 0 && b > 1) { if (r > CAP / b) return CAP; r *= b; } return r; }
    ll capped(int l, int r) {                              // min(tower a[l..r], CAP)
        int k = min(nextOne[l], r);                        // a 1 at k makes everything above it irrelevant
        if (k - l >= 6) return CAP;                        // at least 6 numbers >= 2: beyond 2^62
        ll v = min(a[k], CAP);
        for (int i = k - 1; i >= l; i--) v = satPow(a[i], min(v, 64LL));
        return v;
    }
    ll tower(int l, int r, ll m) {                         // a[l..r] as a tower, mod m
        if (m == 1) return 0;
        if (l == r) return a[l] % m;
        ll ph = phi(m), ec = capped(l + 1, r);
        ll e = ec < ph ? ec : tower(l + 1, r, ph) + ph;    // exact exponent, or an equivalent one >= phi(m)
        return power(a[l], e, m);
    }
};

Proof: the generalized Euler theorem and the capped tower

Claim: if $x,y\ge \phi (m)$ and $x\equiv y(\mathrm{mod}\phi (m))$, then $a^x\equiv a^y(\mathrm{mod}m)$ for every $a$. By the CRT it suffices modulo each prime power $p^e\parallel m$, using $\phi (p^e)\mid \phi (m)$. If $p\nmid a$, Euler gives $a^{\phi (p^e)}\equiv 1$, and $x-y$ is a multiple of $\phi (p^e)$. If $p\mid a$, then $a^z\equiv 0(\mathrm{mod}{p}^e)$ for every $z\ge e$, and both exponents are $\ge \phi (m)\ge \phi (p^e)=p^{e-1}(p-1)\ge e$.

The code. capped returns the exact tower value when it is below $2^{62}$, and $2^{62}$ otherwise: above a 1 nothing matters ($1^z=1$); six or more numbers $\ge 2$ give at least $2^{2^{2^{2^{2^2}}}}>2^{62}$; otherwise the at most six levels are evaluated from the top with saturation (an exponent $\ge 64$ on a base $\ge 2$ saturates). In tower, if the capped exponent is below $\phi (m)$ it is the exact exponent and the power is exact modulo $m$. Otherwise the true exponent $x$ is $\ge \phi (m)$, the recursive call gives $x\mathrm{mod}\phi (m)$, and $e=(x\mathrm{mod}\phi (m))+\phi (m)$ satisfies the claim’s hypotheses with $y=e$. Since $\phi (m)$ is even for $m>2$ and $\phi$ of an even number is at most half of it, the modulus at least halves every two levels: $O(\log m)$ levels. Checked against exact evaluation on random towers whose true exponent fits in 64 bits, including bases 1.

Balanced ternary and non-adjacent form

Base 3 with digits $\{-1,0,1\}$: every integer, positive or negative, has a unique representation without a sign ($5=9-3-1$, $-5=-9+3+1$, $7=9-3+1$). Conversion: take $n\mathrm{mod}3$; a 2 becomes $-1$ with a carry. Bachet’s balance: weights $1,3,\dots ,3^k$, each on the empty pan ($+1$), with the object ($-1$) or unused ($0$), weigh every integer mass up to $\frac{3^{k+1}-1}{2}$ (Bachet’s Weights). The Soviet Setun computer (1958) used balanced ternary.

Property
representation	unique for every integer
negation	flip every digit; no two’s-complement asymmetry
sign	carried by the leading digit
rounding	truncation gives the nearest value
range of $k$ digits	$[-\frac{3^k-1}{2},\frac{3^k-1}{2}]$, symmetric
digits needed	about $\log_3(2

Problem	Why balanced ternary
weigh any mass with weights $3^k$	the defining example
$n$ as $\pm 3^a\pm 3^b\pm \dots$ with distinct exponents	exactly the representation (always possible, uniquely)
games where each move is $\pm 3^k$	state = digits
minimum number of signed powers of 3 summing to $n$	the number of non-zero digits
ternary search trees, three-way partition	the three-outcome analogy

Non-adjacent form (NAF): base 2 with digits $\{-1,0,1\}$ and no two adjacent non-zero digits; unique, of minimal weight, with on average one third of the digits non-zero (versus one half in binary): exponentiation and elliptic-curve scalar multiplication get about 11% faster (fewer multiplications; subtraction of a point is as cheap as addition). Related systems: binary, Zeckendorf (Fibonacci base, no two adjacent 1s), factorial base (digit $i$ in $[0,i]$; permutation ranking via the Lehmer code), Gray code (consecutive values differ in one bit), negabinary (base $-2$, negatives without a sign): Modular Arithmetic has negative bases and fromBalancedTernary.

vector<int> toBalancedTernary(ll n) {                      // least significant digit first, digits in {-1, 0, 1}
    vector<int> d;
    while (n) {
        int r = int(((n % 3) + 3) % 3);
        if (r == 2) { d.push_back(-1); n = (n + 1) / 3; }   // exact division
        else { d.push_back(r); n = (n - r) / 3; }
    }
    return d;
}
 
vector<int> naf(ll n) {                                    // n >= 0, least significant first, no two adjacent non-zeros
    vector<int> d;
    while (n) {
        if (n & 1) { int z = 2 - int(n % 4); d.push_back(z); n -= z; }   // z = 1 if n = 1 mod 4, -1 if n = 3 mod 4
        else d.push_back(0);
        n >>= 1;
    }
    return d;
}

Proof: balanced ternary and NAF are unique, and NAF density is 1/3

Balanced ternary. The last digit must be $\equiv n(\mathrm{mod}3)$ and lie in $\{-1,0,1\}$, so it is forced ($2\mapsto -1$); then $\frac{n-d_0}{3}$ is an integer of smaller absolute value (for $n\ne 0$) and we recurse: existence and uniqueness by induction on $|n|$. With $k$ digits the extremes are $\pm {\sum }_{i<k}{3}^i=\pm \frac{3^k-1}{2}$, and there are exactly $3^k$ digit strings for the $3^k$ integers in that range. NAF. If $n$ is even the digit is 0. If $n$ is odd the digit $z=\pm 1$ is chosen with $n-z\equiv 0(\mathrm{mod}4)$, so the next digit is 0: no two adjacent non-zeros. Uniqueness: an odd $n$ forces $z\equiv n(\mathrm{mod}2)$ and the non-adjacency forces $\frac{n-z}{2}$ even, which determines $z\mathrm{mod}4$. Density: after a non-zero digit the next is 0, and then each bit is odd with probability $\frac{1}{2}$ independently (for random $n$), so the gap to the next non-zero digit is $1+\mathrm{Geom}(\frac{1}{2})$, mean 3: one non-zero digit per 3 positions. Binary exponentiation then costs about $n$ squarings and $\frac{n}{3}$ multiplications instead of $\frac{n}{2}$: $\frac{1.5-1.33}{1.5}\approx 11\%$ less work. Checked: reconstruction and the digit constraints for all $|n|\le 10^5$.

Elliptic curves over ${\mathbb{F}}_p$

Points of $y^2=x^3+ax+b(\mathrm{mod}p)$ ($p>3$ prime, $4a^3+27b^2\not\equiv 0$) plus the point at infinity $\mathcal{O}$ form an abelian group under the chord-and-tangent law: slope $\frac{y_2-y_1}{x_2-x_1}$, or $\frac{3x^2+a}{2y}$ when doubling; $P+(-P)=\mathcal{O}$. Hasse: $|\#E-(p+1)|\le 2\sqrt{p}$. Scalar multiplication by double-and-add is the one-way function of elliptic-curve cryptography (the discrete logarithm is hard), and Lenstra’s factorization runs the same arithmetic modulo a composite $n$, succeeding when an inversion fails (Theoretical Algorithms). Cost: $O(\log p)$ per addition (one inverse), $O(\log k\log p)$ per multiplication, point counting $O(p\log p)$ by Euler’s criterion.

struct ECPoint {
    ll x = 0, y = 0; bool inf = false;
    static ECPoint infinity() { return {0, 0, true}; }
    bool operator==(const ECPoint& o) const { return inf == o.inf && (inf || (x == o.x && y == o.y)); }
};
 
struct EllipticCurve {
    ll a, b, p;
    EllipticCurve(ll a_, ll b_, ll p_) : a((a_ % p_ + p_) % p_), b((b_ % p_ + p_) % p_), p(p_) {}
    ll mul(ll x, ll y) const { return (ll)((i128)x * y % p); }
    ll pw(ll x, ll e) const { ll r = 1; for (x %= p; e; e >>= 1, x = mul(x, x)) if (e & 1) r = mul(r, x); return r; }
    ll inv(ll x) const { return pw(x, p - 2); }
    bool contains(const ECPoint& P) const { return P.inf || mul(P.y, P.y) == (mul(mul(P.x, P.x), P.x) + mul(a, P.x) + b) % p; }
    ECPoint add(const ECPoint& P, const ECPoint& Q) const {
        if (P.inf) return Q;
        if (Q.inf) return P;
        ll s;
        if (P.x == Q.x) {
            if ((P.y + Q.y) % p == 0) return ECPoint::infinity();
            s = mul((3 * mul(P.x, P.x) + a) % p, inv(2 * P.y % p));
        } else s = mul(((Q.y - P.y) % p + p) % p, inv(((Q.x - P.x) % p + p) % p));
        ll x = ((mul(s, s) - P.x - Q.x) % p + 2 * p) % p;
        ll y = ((mul(s, (P.x - x + p) % p) - P.y) % p + p) % p;
        return {x, y, false};
    }
    ECPoint negate(const ECPoint& P) const { return P.inf ? P : ECPoint{P.x, (p - P.y) % p, false}; }
    ECPoint multiply(ECPoint P, ll k) const { ECPoint R = ECPoint::infinity(); for (; k > 0; k >>= 1, P = add(P, P)) if (k & 1) R = add(R, P); return R; }
    ll countPoints() const {                               // including infinity
        ll cnt = 1;
        for (ll x = 0; x < p; x++) { ll rhs = (mul(mul(x, x), x) + mul(a, x) + b) % p; cnt += rhs == 0 ? 1 : pw(rhs, (p - 1) / 2) == 1 ? 2 : 0; }
        return cnt;
    }
};

Proof: the formulas and the point count

The line through $P,Q$ (or the tangent at $P$) meets the cubic in a third point $R$; substituting $y=sx+c$ into $y^2=x^3+ax+b$ gives a monic cubic in $x$ whose roots sum to $s^2$, so $x_R=s^2-x_P-x_Q$, and $P+Q$ is the reflection of $R$: $y=s(x_P-x_R)-y_P$. The tangent slope comes from implicit differentiation, $2y{y}^{\prime }=3x^2+a$. A vertical line ($x_P=x_Q$, $y_P=-y_Q$) meets the curve at infinity: $P+Q=\mathcal{O}$. Associativity and Hasse’s bound are classical theorems (quoted). Counting: for each $x$, $y^2=r$ has 1 solution if $r=0$, 2 if $r$ is a non-zero square (Euler: $r^{(p-1)/2}=1$) and 0 otherwise. Checked: associativity, commutativity and the identity on random points of random curves, multiply against repeated addition, countPoints against brute force over all $(x,y)$, and Hasse’s bound.

Gaussian integers

$\mathbb{Z}[i]=\{a+bi\}$ with norm $N(a+bi)=a^2+b^2$ is a Euclidean domain (rounding each coordinate of $x/y$ leaves a remainder of norm $\le N(y)/2$), so Euclid’s algorithm computes gcds. Units: $\pm 1,\pm i$. Primes: $1+i$; rational primes $p\equiv 3(\mathrm{mod}4)$; and the two factors $a\pm bi$ of each $p\equiv 1(\mathrm{mod}4)=a^2+b^2$. This gives an $O(\log p)$ two-square decomposition once $t^2\equiv -1(\mathrm{mod}p)$ is known, and the representation count $r_2(n)=4{\sum }_{d\mid n}{\chi }_4(d)$ (with ${\chi }_4(d)=0,1,0,-1$ for $d\equiv 0,1,2,3(\mathrm{mod}4)$).

struct GI { ll a, b; };                                     // a + b i (keep norms below about 4e18)
GI operator-(GI x, GI y) { return {x.a - y.a, x.b - y.b}; }
GI operator*(GI x, GI y) { return {x.a * y.a - x.b * y.b, x.a * y.b + x.b * y.a}; }
ll roundDiv(i128 n, i128 d) { return (ll)(n >= 0 ? (2 * n + d) / (2 * d) : -((-2 * n + d) / (2 * d))); }   // d > 0
GI gquot(GI x, GI y) {                                      // nearest Gaussian integer to x / y
    i128 nn = (i128)y.a * y.a + (i128)y.b * y.b;
    i128 re = (i128)x.a * y.a + (i128)x.b * y.b, im = (i128)x.b * y.a - (i128)x.a * y.b;
    return {roundDiv(re, nn), roundDiv(im, nn)};
}
GI ggcd(GI x, GI y) { while (y.a || y.b) { GI r = x - gquot(x, y) * y; x = y; y = r; } return x; }
pair<ll, ll> twoSquares(ll p) {                             // prime p = 1 (mod 4) -> (a, b) with a^2 + b^2 = p
    ll t = 0;
    for (ll c = 2;; c++) if (power(c, (p - 1) / 2, p) == p - 1) { t = power(c, (p - 1) / 4, p); break; }   // t^2 = -1
    GI g = ggcd({p, 0}, {t, 1});
    return {llabs(g.a), llabs(g.b)};
}

Proof: Euclidean division, two squares, and r_2

Division. Write $x/y=u+vi$ and round to the nearest $q$; then $|x/y-q{|}^2\le (\frac{1}{2}{)}^2+(\frac{1}{2}{)}^2=\frac{1}{2}$, so $N(x-qy)\le \frac{N(y)}{2}$: norms halve each step, $O(\log N)$ iterations, and the last non-zero remainder is a gcd (same argument as for integers). Two squares. With $t^2\equiv -1$, $p\mid t^2+1=(t+i)(t-i)$ but $p$ divides neither factor, so $p$ is not prime in $\mathbb{Z}[i]$: $p=\pi \bar{\pi }$ with $N(\pi )=p$. Some Gaussian prime factor $\pi$ of $p$ divides $t+i$, so $g=\gcd (p,t+i)$ is not a unit, and it is not an associate of $p$ (since $p\nmid t+i$); hence $N(g)$ is a proper divisor $>1$ of $N(p)=p^2$, so $N(g)=p$: $g=a+bi$ with $a^2+b^2=p$. A quadratic non-residue $c$ has $c^{(p-1)/2}\equiv -1$, so $t=c^{(p-1)/4}$ squares to $-1$. $r_2$: representations of $n$ as $a^2+b^2$ correspond to Gaussian integers of norm $n$; counting them through the factorization of $n$ into Gaussian primes gives $4{\prod }_{p\equiv 1(4)}(e_p+1)$ (zero if some $p\equiv 3(\mathrm{mod}4)$ has odd exponent), which equals $4{\sum }_{d\mid n}{\chi }_4(d)$ because ${\chi }_4$ is multiplicative. Checked: twoSquares for every prime $p\equiv 1(\mathrm{mod}4)$ below $10^6$, and $r_2(n)$ against brute force for $n\le 5000$.

Related notes

• Congruences and CRT
• Floor Sums
• Fibonacci Rabbits