code: def mpf_bernoulli(n, prec, rnd=None):
      """Computation of Bernoulli numbers (numerically)"""
      if n < 2:
          if n < 0:
              raise ValueError(
                  "Bernoulli numbers only defined for n >= 0"
              )
          if n == 0:
              return fone
          if n == 1:
              return mpf_neg(fhalf)
      # For odd n > 1, the Bernoulli numbers are zero
      if n & 1:
          return fzero
      # If precision is extremely high, we can save time by computing
      # the Bernoulli number at a lower precision that is sufficient to
      # obtain the exact fraction, round to the exact fraction, and
      # convert the fraction back to an mpf value at the original precision
      if (
          prec > BERNOULLI_PREC_CUTOFF
          and prec > bernoulli_size(n) * 1.1 + 1000
      ):
          p, q = bernfrac(n)
          return from_rational(p, q, prec, rnd or round_floor)
      if n > MAX_BERNOULLI_CACHE:
          return mpf_bernoulli_huge(n, prec, rnd)
      wp = prec + 30
      # Reuse nearby precisions
      wp += 32 - (prec & 31)
      cached = bernoulli_cache.get(wp)
      if cached:
          numbers, state = cached
          if n in numbers:
              if not rnd:
                  return numbers[n]
              return mpf_pos(numbers[n], prec, rnd)
          m, bin, bin1 = state
          if n - m > 10:
              return mpf_bernoulli_huge(n, prec, rnd)
      else:
          if n > 10:
              return mpf_bernoulli_huge(n, prec, rnd)
          numbers = {0: fone}
          m, bin, bin1 = state = [2, MPZ(10), MPZ_ONE]
          bernoulli_cache[wp] = (numbers, state)
      while m <= n:
          # print m
          case = m % 6
          # Accurately estimate size of B_m so we can use
          # fixed point math without using too much precision
          szbm = bernoulli_size(m)
          s = 0
          sexp = max(0, szbm) - wp
          if m < 6:
              a = MPZ_ZERO
          else:
              a = bin1
          for j in xrange(1, m // 6 + 1):
              usign, uman, uexp, ubc = u = numbers[m - 6 * j]
              if usign:
                  uman = -uman
              s += lshift(a * uman, uexp - sexp)
              # Update inner binomial coefficient
              j6 = 6 * j
              a *= (
                  (m - 5 - j6)
                  * (m - 4 - j6)
                  * (m - 3 - j6)
                  * (m - 2 - j6)
                  * (m - 1 - j6)
                  * (m - j6)
              )
              a //= (
                  (4 + j6)
                  * (5 + j6)
                  * (6 + j6)
                  * (7 + j6)
                  * (8 + j6)
                  * (9 + j6)
              )
          if case == 0:
              b = mpf_rdiv_int(m + 3, f3, wp)
          if case == 2:
              b = mpf_rdiv_int(m + 3, f3, wp)
          if case == 4:
              b = mpf_rdiv_int(-m - 3, f6, wp)
          s = from_man_exp(s, sexp, wp)
          b = mpf_div(mpf_sub(b, s, wp), from_int(bin), wp)
          numbers[m] = b
          m += 2
          # Update outer binomial coefficient
          bin = bin * ((m + 2) * (m + 3)) // (m * (m - 1))
          if m > 6:
              bin1 = bin1 * ((2 + m) * (3 + m)) // ((m - 7) * (m - 6))
          state[:] = [m, bin, bin1]
      return numbers[n]
  
imports: import math
  import sys
  from .backend import xrange
  from .backend import MPZ, MPZ_ZERO, MPZ_ONE, MPZ_THREE, gmpy
  from .libintmath import list_primes, ifac, ifac2, moebius
  from .libmpf import round_floor, round_ceiling, round_down, round_up, round_nearest, round_fast, lshift, sqrt_fixed, isqrt_fast, fzero, fone, fnone, fhalf, ftwo, finf, fninf, fnan, from_int, to_int, to_fixed, from_man_exp, from_rational, mpf_pos, mpf_neg, mpf_abs, mpf_add, mpf_sub, mpf_mul, mpf_mul_int, mpf_div, mpf_sqrt, mpf_pow_int, mpf_rdiv_int, mpf_perturb, mpf_le, mpf_lt, mpf_gt, mpf_shift, negative_rnd, reciprocal_rnd, bitcount, to_float, mpf_floor, mpf_sign, ComplexResult
  from .libelefun import constant_memo, def_mpf_constant, mpf_pi, pi_fixed, ln2_fixed, log_int_fixed, mpf_ln2, mpf_exp, mpf_log, mpf_pow, mpf_cosh, mpf_cos_sin, mpf_cosh_sinh, mpf_cos_sin_pi, mpf_cos_pi, mpf_sin_pi, ln_sqrt2pi_fixed, mpf_ln_sqrt2pi, sqrtpi_fixed, mpf_sqrtpi, cos_sin_fixed, exp_fixed
  from .libmpc import mpc_zero, mpc_one, mpc_half, mpc_two, mpc_abs, mpc_shift, mpc_pos, mpc_neg, mpc_add, mpc_sub, mpc_mul, mpc_div, mpc_add_mpf, mpc_mul_mpf, mpc_div_mpf, mpc_mpf_div, mpc_mul_int, mpc_pow_int, mpc_log, mpc_exp, mpc_pow, mpc_cos_pi, mpc_sin_pi, mpc_reciprocal, mpc_square, mpc_sub_mpf
label_desc: Use image node_id 9 for calling a global function with example usage: mpf_bernoulli(n, prec, rnd) and returns: numbers, fzero, from_rational, mpf_bernoulli_huge, fone, mpf_neg, mpf_pos, mpf_bernoulli_huge, mpf_bernoulli_huge, numbers