approximate imaginary-quadratic class numbers using analytic class number formula

src/sage/rings/number_field/order.py

@@ -120,6 +120,82 @@ def quadratic_order_class_number(disc):
     return ZZ(h)
 
 
+def quadratic_order_approximate_class_number(disc, *, bound=10**4):
+    r"""
+    Return *an approximation of* the class number of
+    the quadratic order of given discriminant.
+
+    Currently only implemented for maximal orders
+    in imaginary-quadratic fields.
+
+    EXAMPLES::
+
+        sage: from sage.rings.number_field.order import quadratic_order_approximate_class_number
+        sage: QuadraticField(-419).class_number()
+        9
+        sage: quadratic_order_approximate_class_number(-419)  # rel tol .01
+        9.01653836091712
+
+    ::
+
+        sage: from sage.rings.number_field.order import quadratic_order_approximate_class_number
+        sage: d = 100000000000031
+        sage: QuadraticField(-d).class_number(proof=False)
+        14414435
+        sage: round(quadratic_order_approximate_class_number(-d))  # rel tol .01
+        14407657
+        sage: round(quadratic_order_approximate_class_number(-d, bound=10**6))  # rel tol .01
+        14413626
+
+    Test it against the exact class number computed for the CSIDH-512 prime (source: https://eprint.iacr.org/2019/498.pdf)::
+
+        sage: from sage.rings.number_field.order import quadratic_order_approximate_class_number
+        sage: p = 4 * prod(primes(3,374)) * 587 - 1
+        sage: hreal = 84884147409828091725676728670213067387206838101828807864190286991865870575397
+        sage: assert not hreal * BQFClassGroup(-p).random_element()
+        sage: h = round(quadratic_order_approximate_class_number(-p, bound=10**3)); h  # rel tol .01
+        85020334529027955331134025285584937708277525762952749243936367281020276898911
+        sage: RR(h / hreal)  # abs tol .01
+        1.00160438813790
+        sage: h = round(quadratic_order_approximate_class_number(-p, bound=10**4)); h  # rel tol .01
+        84396416322932187013685015858232028818143153418602750494380687212063263982976
+        sage: RR(h / hreal)  # abs tol .01
+        0.994254155790231
+        sage: h = round(quadratic_order_approximate_class_number(-p, bound=10**5)); h  # rel tol .01
+        84823787383264935642168065590216697209124465727576867303448690101681967969348
+        sage: RR(h / hreal)  # abs tol .01
+        0.999288912848806
+        sage: h = round(quadratic_order_approximate_class_number(-p, bound=10**6)); h  # rel tol .01, long time (2s)
+        84884627342209070883738394179700676127184729325091471467872632417652836154863
+        sage: RR(h / hreal)  # abs tol .01, long time (2s)
+        1.00000565396951
+
+    ALGORITHM: Finite approximation of the infinite product given by
+    the analytic class number formula, using primes up to ``bound``.
+    """
+    disc = ZZ(disc)
+    if disc >= 0:
+        raise NotImplementedError('only imaginary-quadratic fields supported')
+    if not disc.is_fundamental_discriminant():
+        raise NotImplementedError('only fundamental discriminants supported')
+
+    from sage.rings.real_mpfr import RealField
+    from sage.arith.misc import primes, kronecker_symbol
+    from sage.symbolic.constants import pi
+
+    w = 6 if disc == -3 else 4 if disc == -4 else 2
+    RR = RealField(max(53, disc.bit_length()))  # wild guess!
+
+    # compute numerator and denominator separately for speed
+    L1 = L2 = RR(1)
+    for ell in primes(bound):
+        L1 *= ell
+        L2 *= ell - kronecker_symbol(disc, ell)
+    L = L1 / L2
+
+    return RR(w * abs(disc).sqrt() * L / (2 * pi))
+
+
 class OrderFactory(UniqueFactory):
     r"""
     Abstract base class for factories creating orders, such as
@@ -1091,6 +1167,15 @@ def class_number(self, proof=None):
         r"""
         Return the class number of this order.
 
+        .. NOTE::
+
+            For some applications (e.g., in algorithms for computing
+            class groups) it is required to merely *approximate* the
+            class number. The function
+            :func:`quadratic_order_approximate_class_number`
+            can be used to compute such an approximation (currently
+            restricted to maximal imaginary-quadratic orders).
+
         EXAMPLES::
 
             sage: ZZ[2^(1/3)].class_number()                                            # needs sage.symbolic