def ref_metric__hatted_quantities(SymPySimplifyExpressions=True):
    enable_rfm_precompute = False
    if (par.parval_from_str((thismodule + '::enable_rfm_precompute')) == 'True'):
        enable_rfm_precompute = True
    DIM = par.parval_from_str('grid::DIM')
    global ReU, ReDD, ghatDD
    ReU = ixp.zerorank1()
    ReDD = ixp.zerorank2()
    ghatDD = ixp.zerorank2()
    if (not enable_rfm_precompute):
        for i in range(DIM):
            scalefactor_orthog[i] = sp.sympify(scalefactor_orthog[i])
            ghatDD[i][i] = (scalefactor_orthog[i] ** 2)
            ReU[i] = (1 / scalefactor_orthog[i])
            for j in range(DIM):
                ReDD[i][j] = (scalefactor_orthog[i] * scalefactor_orthog[j])
    else:
        for i in range(DIM):
            scalefactor_orthog_funcform[i] = sp.sympify(scalefactor_orthog_funcform[i])
            ghatDD[i][i] = (scalefactor_orthog_funcform[i] ** 2)
            ReU[i] = (1 / scalefactor_orthog_funcform[i])
            for j in range(DIM):
                ReDD[i][j] = (scalefactor_orthog_funcform[i] * scalefactor_orthog_funcform[j])
    global ghatUU
    global detgammahat
    ghatUU = ixp.zerorank2()
    (ghatUU, detgammahat) = ixp.symm_matrix_inverter3x3(ghatDD)
    for i in range(DIM):
        for j in range(DIM):
            if (ReDD[i][j] != ReDD[j][i]):
                print(((((((((((('Error: ReDD[' + str(i)) + '][') + str(j)) + '] != ReDD[') + str(j)) + '][') + str(i)) + ': ') + str(ReDD[i][j])) + '!=') + str(ReDD[j][i])))
                sys.exit(1)
            if (ghatDD[i][j] != ghatDD[j][i]):
                print(((((((((((('Error: ghatDD[' + str(i)) + '][') + str(j)) + '] != ghatDD[') + str(j)) + '][') + str(i)) + ': ') + str(ghatDD[i][j])) + '!=') + str(ghatDD[j][i])))
                sys.exit(1)
            if (ghatUU[i][j] != ghatUU[j][i]):
                print(((((((((((('Error: ghatUU[' + str(i)) + '][') + str(j)) + '] != ghatUU[') + str(j)) + '][') + str(i)) + ': ') + str(ghatUU[i][j])) + '!=') + str(ghatUU[j][i])))
                sys.exit(1)
    global detgammahatdD, detgammahatdDD
    detgammahatdD = ixp.zerorank1(DIM)
    detgammahatdDD = ixp.zerorank2(DIM)
    for i in range(DIM):
        detgammahatdD[i] = sp.diff(detgammahat, xx[i])
        for j in range(DIM):
            detgammahatdDD[i][j] = sp.diff(detgammahatdD[i], xx[j])
    global ReUdD, ReUdDD
    ReUdD = ixp.zerorank2(DIM)
    ReUdDD = ixp.zerorank3(DIM)
    for i in range(DIM):
        for j in range(DIM):
            ReUdD[i][j] = sp.diff(ReU[i], xx[j])
            for k in range(DIM):
                ReUdDD[i][j][k] = sp.diff(ReUdD[i][j], xx[k])
    global ReDDdD, ReDDdDD
    ReDDdD = ixp.zerorank3(DIM)
    ReDDdDD = ixp.zerorank4(DIM)
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                ReDDdD[i][j][k] = sp.diff(ReDD[i][j], xx[k])
                for l in range(DIM):
                    ReDDdDD[i][j][k][l] = sp.diff(ReDDdD[i][j][k], xx[l])
    global ghatDDdD, ghatDDdDD
    ghatDDdD = ixp.zerorank3(DIM)
    ghatDDdDD = ixp.zerorank4(DIM)
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                if (SymPySimplifyExpressions == True):
                    ghatDDdD[i][j][k] = sp.simplify(sp.diff(ghatDD[i][j], xx[k]))
                else:
                    ghatDDdD[i][j][k] = sp.diff(ghatDD[i][j], xx[k])
                for l in range(DIM):
                    ghatDDdDD[i][j][k][l] = sp.diff(ghatDDdD[i][j][k], xx[l])
    global GammahatUDD
    GammahatUDD = ixp.zerorank3(DIM)
    for i in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                for m in range(DIM):
                    GammahatUDD[i][k][l] += ((sp.Rational(1, 2) * ghatUU[i][m]) * ((ghatDDdD[m][k][l] + ghatDDdD[m][l][k]) - ghatDDdD[k][l][m]))
    global GammahatUDDdD
    GammahatUDDdD = ixp.zerorank4(DIM)
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    GammahatUDDdD[i][j][k][l] = sp.diff(GammahatUDD[i][j][k], xx[l])
    if (not enable_rfm_precompute):
        return

    def make_replacements(expr):
        sympy_version = sp.__version__.replace('rc', '...').replace('b', '...')
        sympy_version_decimal = float((int(sympy_version.split('.')[0]) + (int(sympy_version.split('.')[1]) / 10.0)))
        is_old_sympy_version = (sympy_version_decimal < 1.2)
        rule = {}
        for item in sp.preorder_traversal(expr):
            if (item.func == sp.Derivative):
                strfunc = str(item.args[0]).split('_funcform(', 1)[0]
                if is_old_sympy_version:
                    (var, order) = (str(item.args[1])[2:], (len(item.args) - 1))
                else:
                    (var, order) = (str(item.args[1][0])[2:], item.args[1][1])
                oper = (('__D' + ('D' * (order - 1))) + (var * order))
                rule[item] = sp.sympify((strfunc + oper))
        expr = expr.xreplace(rule)
        rule = {}
        for item in sp.preorder_traversal(expr):
            if ('_funcform' in str(item.func)):
                strfunc = str(item.func).split('_funcform', 1)[0]
                rule[item] = sp.sympify(strfunc)
        return expr.xreplace(rule)
    detgammahat = make_replacements(detgammahat)
    for i in range(DIM):
        ReU[i] = make_replacements(ReU[i])
        detgammahatdD[i] = make_replacements(detgammahatdD[i])
        for j in range(DIM):
            ReDD[i][j] = make_replacements(ReDD[i][j])
            ReUdD[i][j] = make_replacements(ReUdD[i][j])
            ghatDD[i][j] = make_replacements(ghatDD[i][j])
            ghatUU[i][j] = make_replacements(ghatUU[i][j])
            detgammahatdDD[i][j] = make_replacements(detgammahatdDD[i][j])
            for k in range(DIM):
                ReDDdD[i][j][k] = make_replacements(ReDDdD[i][j][k])
                ReUdDD[i][j][k] = make_replacements(ReUdDD[i][j][k])
                ghatDDdD[i][j][k] = make_replacements(ghatDDdD[i][j][k])
                GammahatUDD[i][j][k] = make_replacements(GammahatUDD[i][j][k])
                for l in range(DIM):
                    ReDDdDD[i][j][k][l] = make_replacements(ReDDdDD[i][j][k][l])
                    ghatDDdDD[i][j][k][l] = make_replacements(ghatDDdDD[i][j][k][l])
                    GammahatUDDdD[i][j][k][l] = make_replacements(GammahatUDDdD[i][j][k][l])
    freevars = []
    freevars.extend(detgammahat.free_symbols)
    for i in range(DIM):
        freevars.extend(ReU[i].free_symbols)
        freevars.extend(detgammahatdD[i].free_symbols)
        for j in range(DIM):
            freevars.extend(ReDD[i][j].free_symbols)
            freevars.extend(ReUdD[i][j].free_symbols)
            freevars.extend(ghatDD[i][j].free_symbols)
            freevars.extend(ghatUU[i][j].free_symbols)
            freevars.extend(detgammahatdDD[i][j].free_symbols)
            for k in range(DIM):
                freevars.extend(ReDDdD[i][j][k].free_symbols)
                freevars.extend(ReUdDD[i][j][k].free_symbols)
                freevars.extend(ghatDDdD[i][j][k].free_symbols)
                freevars.extend(GammahatUDD[i][j][k].free_symbols)
                for l in range(DIM):
                    freevars.extend(ReDDdDD[i][j][k][l].free_symbols)
                    freevars.extend(ghatDDdDD[i][j][k][l].free_symbols)
                    freevars.extend(GammahatUDDdD[i][j][k][l].free_symbols)
    freevars_uniq = superfast_uniq(freevars)
    freevars_uniq_xx_indep = []
    for i in range(len(freevars_uniq)):
        freevars_uniq_xx_indep.append(freevars_uniq[i])
    freevars_uniq_vals = []
    for (i, var) in enumerate(freevars_uniq):
        basename = str(var).split('__')[0].replace('_funcform', '')
        derivatv = ''
        if ('__' in str(var)):
            derivatv = str(var).split('__')[1].replace('_funcform', '')
        if (basename == 'f0_of_xx0'):
            basefunc = f0_of_xx0
        elif (basename == 'f1_of_xx1'):
            basefunc = f1_of_xx1
        elif (basename == 'f2_of_xx0_xx1'):
            basefunc = f2_of_xx0_xx1
        elif (basename == 'f3_of_xx0'):
            basefunc = f3_of_xx0
        elif (basename == 'f4_of_xx2'):
            basefunc = f4_of_xx2
        else:
            print((('Error: function inside ' + str(var)) + ' undefined.'))
            sys.exit(1)
        diff_result = basefunc
        if (derivatv == ''):
            pass
        else:
            derivorder = derivatv.replace('d', '').replace('D', '').replace('0', '0 ').replace('1', '1 ').replace('2', '2 ').split(' ')
            for derivdirn in derivorder:
                if (derivdirn != ''):
                    derivwrt = xx[int(derivdirn)]
                    diff_result = sp.diff(diff_result, derivwrt)
        freevars_uniq_vals.append(diff_result)
        frees_uniq = superfast_uniq(diff_result.free_symbols)
        has_xx_dependence = False
        for dirn in range(3):
            if (gri.xx[dirn] in frees_uniq):
                has_xx_dependence = True
        if (not has_xx_dependence):
            freevars_uniq_xx_indep[i] = diff_result
    for (varidx, freevar) in enumerate(freevars_uniq):
        detgammahat = detgammahat.subs(freevar, freevars_uniq_xx_indep[varidx])
        for i in range(DIM):
            ReU[i] = ReU[i].subs(freevar, freevars_uniq_xx_indep[varidx])
            detgammahatdD[i] = detgammahatdD[i].subs(freevar, freevars_uniq_xx_indep[varidx])
            for j in range(DIM):
                ReDD[i][j] = ReDD[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                ReUdD[i][j] = ReUdD[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                ghatDD[i][j] = ghatDD[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                ghatUU[i][j] = ghatUU[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                detgammahatdDD[i][j] = detgammahatdDD[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                for k in range(DIM):
                    ReDDdD[i][j][k] = ReDDdD[i][j][k].subs(freevar, freevars_uniq_xx_indep[varidx])
                    ReUdDD[i][j][k] = ReUdDD[i][j][k].subs(freevar, freevars_uniq_xx_indep[varidx])
                    ghatDDdD[i][j][k] = ghatDDdD[i][j][k].subs(freevar, freevars_uniq_xx_indep[varidx])
                    GammahatUDD[i][j][k] = GammahatUDD[i][j][k].subs(freevar, freevars_uniq_xx_indep[varidx])
                    for l in range(DIM):
                        ReDDdDD[i][j][k][l] = ReDDdDD[i][j][k][l].subs(freevar, freevars_uniq_xx_indep[varidx])
                        ghatDDdDD[i][j][k][l] = ghatDDdDD[i][j][k][l].subs(freevar, freevars_uniq_xx_indep[varidx])
                        GammahatUDDdD[i][j][k][l] = GammahatUDDdD[i][j][k][l].subs(freevar, freevars_uniq_xx_indep[varidx])
    struct_str = 'typedef struct __rfmstruct__ {\n'
    define_str = ''
    malloc_str = ''
    if (par.parval_from_str('reference_metric::rfm_precompute_to_Cfunctions_and_NRPy_basic_defines') == 'False'):
        malloc_str = 'rfm_struct rfmstruct;\n'
    freemm_str = ''
    readvr_str = ['', '', '']
    readvr_SIMD_outer_str = ['', '', '']
    readvr_SIMD_inner_str = ['', '', '']
    if (len(freevars_uniq_xx_indep) > 0):
        (freevars_uniq_xx_indep, freevars_uniq_vals) = (list(x) for x in zip(*sorted(zip(freevars_uniq_xx_indep, freevars_uniq_vals), key=str)))
    which_freevar = 0
    for expr in freevars_uniq_vals:
        if ('_of_xx' in str(freevars_uniq_xx_indep[which_freevar])):
            frees = expr.free_symbols
            frees_uniq = superfast_uniq(frees)
            xx_list = []
            malloc_size = 1
            for i in range(3):
                if (gri.xx[i] in frees_uniq):
                    xx_list.append(gri.xx[i])
                    malloc_size *= gri.Nxx_plus_2NGHOSTS[i]
            struct_str += (('\tREAL *restrict ' + str(freevars_uniq_xx_indep[which_freevar])) + ';\n')
            malloc_str += (((('rfmstruct.' + str(freevars_uniq_xx_indep[which_freevar])) + ' = (REAL *)malloc(sizeof(REAL)*') + str(malloc_size)) + ');\n')
            freemm_str += (('free(rfmstruct.' + str(freevars_uniq_xx_indep[which_freevar])) + ');\n')
            output_define_and_readvr = False
            for dirn in range(3):
                if ((gri.xx[dirn] in frees_uniq) and (not (gri.xx[((dirn + 1) % 3)] in frees_uniq)) and (not (gri.xx[((dirn + 2) % 3)] in frees_uniq))):
                    define_str += (((((((('for(int i' + str(dirn)) + '=0;i') + str(dirn)) + '<Nxx_plus_2NGHOSTS') + str(dirn)) + ';i') + str(dirn)) + '++) {\n')
                    define_str += (((((('  const REAL xx' + str(dirn)) + ' = xx[') + str(dirn)) + '][i') + str(dirn)) + '];\n')
                    define_str += (((((('  rfmstruct.' + str(freevars_uniq_xx_indep[which_freevar])) + '[i') + str(dirn)) + '] = ') + str(sp.ccode(freevars_uniq_vals[which_freevar]))) + ';\n')
                    define_str += '}\n\n'
                    readvr_str[dirn] += (((((('const REAL ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i') + str(dirn)) + '];\n')
                    readvr_SIMD_outer_str[dirn] += (((((('const double NOSIMD' + str(freevars_uniq_xx_indep[which_freevar])) + ' = rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i') + str(dirn)) + ']; ')
                    readvr_SIMD_outer_str[dirn] += (((('const REAL_SIMD_ARRAY ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = ConstSIMD(NOSIMD') + str(freevars_uniq_xx_indep[which_freevar])) + ');\n')
                    readvr_SIMD_inner_str[dirn] += (((((('const REAL_SIMD_ARRAY ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = ReadSIMD(&rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i') + str(dirn)) + ']);\n')
                    output_define_and_readvr = True
            if ((not output_define_and_readvr) and (gri.xx[0] in frees_uniq) and (gri.xx[1] in frees_uniq)):
                define_str += (((('\nfor(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++) for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) {\n  const REAL xx0 = xx[0][i0];\n  const REAL xx1 = xx[1][i1];\n  rfmstruct.' + str(freevars_uniq_xx_indep[which_freevar])) + '[i0 + Nxx_plus_2NGHOSTS0*i1] = ') + str(sp.ccode(freevars_uniq_vals[which_freevar]))) + ';\n}\n\n')
                readvr_str[0] += (((('const REAL ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i0 + Nxx_plus_2NGHOSTS0*i1];\n')
                readvr_SIMD_outer_str[0] += (((('const double NOSIMD' + str(freevars_uniq_xx_indep[which_freevar])) + ' = rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i0 + Nxx_plus_2NGHOSTS0*i1]; ')
                readvr_SIMD_outer_str[0] += (((('const REAL_SIMD_ARRAY ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = ConstSIMD(NOSIMD') + str(freevars_uniq_xx_indep[which_freevar])) + ');\n')
                readvr_SIMD_inner_str[0] += (((('const REAL_SIMD_ARRAY ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = ReadSIMD(&rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i0 + Nxx_plus_2NGHOSTS0*i1]);\n')
                output_define_and_readvr = True
            if (not output_define_and_readvr):
                print((('ERROR: Could not figure out the (xx0,xx1,xx2) dependency within the expression for ' + str(freevars_uniq_xx_indep[which_freevar])) + ':'))
                print(str(freevars_uniq_vals[which_freevar]))
                sys.exit(1)
        which_freevar += 1
    struct_str += '} rfm_struct;\n'
    outdir = par.parval_from_str((thismodule + '::rfm_precompute_Ccode_outdir'))
    if (par.parval_from_str((thismodule + '::rfm_precompute_to_Cfunctions_and_NRPy_basic_defines')) == 'False'):
        with open(os.path.join(outdir, 'rfm_struct__declare.h'), 'w') as file:
            file.write(struct_str)
        with open(os.path.join(outdir, 'rfm_struct__malloc.h'), 'w') as file:
            file.write(malloc_str)
        with open(os.path.join(outdir, 'rfm_struct__define.h'), 'w') as file:
            file.write(define_str)
        with open(os.path.join(outdir, 'rfm_struct__freemem.h'), 'w') as file:
            file.write(freemm_str)
    else:
        global NRPy_basic_defines_str
        NRPy_basic_defines_str = struct_str
        global rfm_struct__malloc, rfm_struct__define, rfm_struct__freemem
        rfm_struct__malloc = malloc_str
        rfm_struct__define = define_str
        rfm_struct__freemem = freemm_str
    for i in range(3):
        with open(os.path.join(outdir, (('rfm_struct__read' + str(i)) + '.h')), 'w') as file:
            file.write(readvr_str[i])
        with open(os.path.join(outdir, (('rfm_struct__SIMD_outer_read' + str(i)) + '.h')), 'w') as file:
            file.write(readvr_SIMD_outer_str[i])
        with open(os.path.join(outdir, (('rfm_struct__SIMD_inner_read' + str(i)) + '.h')), 'w') as file:
            file.write(readvr_SIMD_inner_str[i])
def ref_metric__hatted_quantities(SymPySimplifyExpressions=True):
    enable_rfm_precompute = False
    if (par.parval_from_str((thismodule + '::enable_rfm_precompute')) == 'True'):
        enable_rfm_precompute = True
    DIM = par.parval_from_str('grid::DIM')
    global ReU, ReDD, ghatDD
    ReU = ixp.zerorank1()
    ReDD = ixp.zerorank2()
    ghatDD = ixp.zerorank2()
    if (not enable_rfm_precompute):
        for i in range(DIM):
            scalefactor_orthog[i] = sp.sympify(scalefactor_orthog[i])
            ghatDD[i][i] = (scalefactor_orthog[i] ** 2)
            ReU[i] = (1 / scalefactor_orthog[i])
            for j in range(DIM):
                ReDD[i][j] = (scalefactor_orthog[i] * scalefactor_orthog[j])
    else:
        for i in range(DIM):
            scalefactor_orthog_funcform[i] = sp.sympify(scalefactor_orthog_funcform[i])
            ghatDD[i][i] = (scalefactor_orthog_funcform[i] ** 2)
            ReU[i] = (1 / scalefactor_orthog_funcform[i])
            for j in range(DIM):
                ReDD[i][j] = (scalefactor_orthog_funcform[i] * scalefactor_orthog_funcform[j])
    global ghatUU
    global detgammahat
    ghatUU = ixp.zerorank2()
    (ghatUU, detgammahat) = ixp.symm_matrix_inverter3x3(ghatDD)
    for i in range(DIM):
        for j in range(DIM):
            if (ReDD[i][j] != ReDD[j][i]):
                print(((((((((((('Error: ReDD[' + str(i)) + '][') + str(j)) + '] != ReDD[') + str(j)) + '][') + str(i)) + ': ') + str(ReDD[i][j])) + '!=') + str(ReDD[j][i])))
                sys.exit(1)
            if (ghatDD[i][j] != ghatDD[j][i]):
                print(((((((((((('Error: ghatDD[' + str(i)) + '][') + str(j)) + '] != ghatDD[') + str(j)) + '][') + str(i)) + ': ') + str(ghatDD[i][j])) + '!=') + str(ghatDD[j][i])))
                sys.exit(1)
            if (ghatUU[i][j] != ghatUU[j][i]):
                print(((((((((((('Error: ghatUU[' + str(i)) + '][') + str(j)) + '] != ghatUU[') + str(j)) + '][') + str(i)) + ': ') + str(ghatUU[i][j])) + '!=') + str(ghatUU[j][i])))
                sys.exit(1)
    global detgammahatdD, detgammahatdDD
    detgammahatdD = ixp.zerorank1(DIM)
    detgammahatdDD = ixp.zerorank2(DIM)
    for i in range(DIM):
        detgammahatdD[i] = sp.diff(detgammahat, xx[i])
        for j in range(DIM):
            detgammahatdDD[i][j] = sp.diff(detgammahatdD[i], xx[j])
    global ReUdD, ReUdDD
    ReUdD = ixp.zerorank2(DIM)
    ReUdDD = ixp.zerorank3(DIM)
    for i in range(DIM):
        for j in range(DIM):
            ReUdD[i][j] = sp.diff(ReU[i], xx[j])
            for k in range(DIM):
                ReUdDD[i][j][k] = sp.diff(ReUdD[i][j], xx[k])
    global ReDDdD, ReDDdDD
    ReDDdD = ixp.zerorank3(DIM)
    ReDDdDD = ixp.zerorank4(DIM)
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                ReDDdD[i][j][k] = sp.diff(ReDD[i][j], xx[k])
                for l in range(DIM):
                    ReDDdDD[i][j][k][l] = sp.diff(ReDDdD[i][j][k], xx[l])
    global ghatDDdD, ghatDDdDD
    ghatDDdD = ixp.zerorank3(DIM)
    ghatDDdDD = ixp.zerorank4(DIM)
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                if (SymPySimplifyExpressions == True):
                    ghatDDdD[i][j][k] = sp.simplify(sp.diff(ghatDD[i][j], xx[k]))
                else:
                    ghatDDdD[i][j][k] = sp.diff(ghatDD[i][j], xx[k])
                for l in range(DIM):
                    ghatDDdDD[i][j][k][l] = sp.diff(ghatDDdD[i][j][k], xx[l])
    global GammahatUDD
    GammahatUDD = ixp.zerorank3(DIM)
    for i in range(DIM):
        for k in range(DIM):
            for l in range(DIM):
                for m in range(DIM):
                    GammahatUDD[i][k][l] += ((sp.Rational(1, 2) * ghatUU[i][m]) * ((ghatDDdD[m][k][l] + ghatDDdD[m][l][k]) - ghatDDdD[k][l][m]))
    global GammahatUDDdD
    GammahatUDDdD = ixp.zerorank4(DIM)
    for i in range(DIM):
        for j in range(DIM):
            for k in range(DIM):
                for l in range(DIM):
                    GammahatUDDdD[i][j][k][l] = sp.diff(GammahatUDD[i][j][k], xx[l])
    if (not enable_rfm_precompute):
        return

    def make_replacements(expr):
        sympy_version = sp.__version__.replace('rc', '...').replace('b', '...')
        sympy_version_decimal = float((int(sympy_version.split('.')[0]) + (int(sympy_version.split('.')[1]) / 10.0)))
        is_old_sympy_version = (sympy_version_decimal < 1.2)
        rule = {}
        for item in sp.preorder_traversal(expr):
            if (item.func == sp.Derivative):
                strfunc = str(item.args[0]).split('_funcform(', 1)[0]
                if is_old_sympy_version:
                    (var, order) = (str(item.args[1])[2:], (len(item.args) - 1))
                else:
                    (var, order) = (str(item.args[1][0])[2:], item.args[1][1])
                oper = (('__D' + ('D' * (order - 1))) + (var * order))
                rule[item] = sp.sympify((strfunc + oper))
        expr = expr.xreplace(rule)
        rule = {}
        for item in sp.preorder_traversal(expr):
            if ('_funcform' in str(item.func)):
                strfunc = str(item.func).split('_funcform', 1)[0]
                rule[item] = sp.sympify(strfunc)
        return expr.xreplace(rule)
    detgammahat = make_replacements(detgammahat)
    for i in range(DIM):
        ReU[i] = make_replacements(ReU[i])
        detgammahatdD[i] = make_replacements(detgammahatdD[i])
        for j in range(DIM):
            ReDD[i][j] = make_replacements(ReDD[i][j])
            ReUdD[i][j] = make_replacements(ReUdD[i][j])
            ghatDD[i][j] = make_replacements(ghatDD[i][j])
            ghatUU[i][j] = make_replacements(ghatUU[i][j])
            detgammahatdDD[i][j] = make_replacements(detgammahatdDD[i][j])
            for k in range(DIM):
                ReDDdD[i][j][k] = make_replacements(ReDDdD[i][j][k])
                ReUdDD[i][j][k] = make_replacements(ReUdDD[i][j][k])
                ghatDDdD[i][j][k] = make_replacements(ghatDDdD[i][j][k])
                GammahatUDD[i][j][k] = make_replacements(GammahatUDD[i][j][k])
                for l in range(DIM):
                    ReDDdDD[i][j][k][l] = make_replacements(ReDDdDD[i][j][k][l])
                    ghatDDdDD[i][j][k][l] = make_replacements(ghatDDdDD[i][j][k][l])
                    GammahatUDDdD[i][j][k][l] = make_replacements(GammahatUDDdD[i][j][k][l])
    freevars = []
    freevars.extend(detgammahat.free_symbols)
    for i in range(DIM):
        freevars.extend(ReU[i].free_symbols)
        freevars.extend(detgammahatdD[i].free_symbols)
        for j in range(DIM):
            freevars.extend(ReDD[i][j].free_symbols)
            freevars.extend(ReUdD[i][j].free_symbols)
            freevars.extend(ghatDD[i][j].free_symbols)
            freevars.extend(ghatUU[i][j].free_symbols)
            freevars.extend(detgammahatdDD[i][j].free_symbols)
            for k in range(DIM):
                freevars.extend(ReDDdD[i][j][k].free_symbols)
                freevars.extend(ReUdDD[i][j][k].free_symbols)
                freevars.extend(ghatDDdD[i][j][k].free_symbols)
                freevars.extend(GammahatUDD[i][j][k].free_symbols)
                for l in range(DIM):
                    freevars.extend(ReDDdDD[i][j][k][l].free_symbols)
                    freevars.extend(ghatDDdDD[i][j][k][l].free_symbols)
                    freevars.extend(GammahatUDDdD[i][j][k][l].free_symbols)
    freevars_uniq = superfast_uniq(freevars)
    freevars_uniq_xx_indep = []
    for i in range(len(freevars_uniq)):
        freevars_uniq_xx_indep.append(freevars_uniq[i])
    freevars_uniq_vals = []
    for (i, var) in enumerate(freevars_uniq):
        basename = str(var).split('__')[0].replace('_funcform', )
        derivatv = 
        if ('__' in str(var)):
            derivatv = str(var).split('__')[1].replace('_funcform', )
        if (basename == 'f0_of_xx0'):
            basefunc = f0_of_xx0
        elif (basename == 'f1_of_xx1'):
            basefunc = f1_of_xx1
        elif (basename == 'f2_of_xx0_xx1'):
            basefunc = f2_of_xx0_xx1
        elif (basename == 'f3_of_xx0'):
            basefunc = f3_of_xx0
        elif (basename == 'f4_of_xx2'):
            basefunc = f4_of_xx2
        else:
            print((('Error: function inside ' + str(var)) + ' undefined.'))
            sys.exit(1)
        diff_result = basefunc
        if (derivatv == ):
            pass
        else:
            derivorder = derivatv.replace('d', ).replace('D', ).replace('0', '0 ').replace('1', '1 ').replace('2', '2 ').split(' ')
            for derivdirn in derivorder:
                if (derivdirn != ):
                    derivwrt = xx[int(derivdirn)]
                    diff_result = sp.diff(diff_result, derivwrt)
        freevars_uniq_vals.append(diff_result)
        frees_uniq = superfast_uniq(diff_result.free_symbols)
        has_xx_dependence = False
        for dirn in range(3):
            if (gri.xx[dirn] in frees_uniq):
                has_xx_dependence = True
        if (not has_xx_dependence):
            freevars_uniq_xx_indep[i] = diff_result
    for (varidx, freevar) in enumerate(freevars_uniq):
        detgammahat = detgammahat.subs(freevar, freevars_uniq_xx_indep[varidx])
        for i in range(DIM):
            ReU[i] = ReU[i].subs(freevar, freevars_uniq_xx_indep[varidx])
            detgammahatdD[i] = detgammahatdD[i].subs(freevar, freevars_uniq_xx_indep[varidx])
            for j in range(DIM):
                ReDD[i][j] = ReDD[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                ReUdD[i][j] = ReUdD[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                ghatDD[i][j] = ghatDD[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                ghatUU[i][j] = ghatUU[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                detgammahatdDD[i][j] = detgammahatdDD[i][j].subs(freevar, freevars_uniq_xx_indep[varidx])
                for k in range(DIM):
                    ReDDdD[i][j][k] = ReDDdD[i][j][k].subs(freevar, freevars_uniq_xx_indep[varidx])
                    ReUdDD[i][j][k] = ReUdDD[i][j][k].subs(freevar, freevars_uniq_xx_indep[varidx])
                    ghatDDdD[i][j][k] = ghatDDdD[i][j][k].subs(freevar, freevars_uniq_xx_indep[varidx])
                    GammahatUDD[i][j][k] = GammahatUDD[i][j][k].subs(freevar, freevars_uniq_xx_indep[varidx])
                    for l in range(DIM):
                        ReDDdDD[i][j][k][l] = ReDDdDD[i][j][k][l].subs(freevar, freevars_uniq_xx_indep[varidx])
                        ghatDDdDD[i][j][k][l] = ghatDDdDD[i][j][k][l].subs(freevar, freevars_uniq_xx_indep[varidx])
                        GammahatUDDdD[i][j][k][l] = GammahatUDDdD[i][j][k][l].subs(freevar, freevars_uniq_xx_indep[varidx])
    struct_str = 'typedef struct __rfmstruct__ {\n'
    define_str = 
    malloc_str = 
    if (par.parval_from_str('reference_metric::rfm_precompute_to_Cfunctions_and_NRPy_basic_defines') == 'False'):
        malloc_str = 'rfm_struct rfmstruct;\n'
    freemm_str = 
    readvr_str = [, , ]
    readvr_SIMD_outer_str = [, , ]
    readvr_SIMD_inner_str = [, , ]
    if (len(freevars_uniq_xx_indep) > 0):
        (freevars_uniq_xx_indep, freevars_uniq_vals) = (list(x) for x in zip(*sorted(zip(freevars_uniq_xx_indep, freevars_uniq_vals), key=str)))
    which_freevar = 0
    for expr in freevars_uniq_vals:
        if ('_of_xx' in str(freevars_uniq_xx_indep[which_freevar])):
            frees = expr.free_symbols
            frees_uniq = superfast_uniq(frees)
            xx_list = []
            malloc_size = 1
            for i in range(3):
                if (gri.xx[i] in frees_uniq):
                    xx_list.append(gri.xx[i])
                    malloc_size *= gri.Nxx_plus_2NGHOSTS[i]
            struct_str += (('\tREAL *restrict ' + str(freevars_uniq_xx_indep[which_freevar])) + ';\n')
            malloc_str += (((('rfmstruct.' + str(freevars_uniq_xx_indep[which_freevar])) + ' = (REAL *)malloc(sizeof(REAL)*') + str(malloc_size)) + ');\n')
            freemm_str += (('free(rfmstruct.' + str(freevars_uniq_xx_indep[which_freevar])) + ');\n')
            output_define_and_readvr = False
            for dirn in range(3):
                if ((gri.xx[dirn] in frees_uniq) and (not (gri.xx[((dirn + 1) % 3)] in frees_uniq)) and (not (gri.xx[((dirn + 2) % 3)] in frees_uniq))):
                    define_str += (((((((('for(int i' + str(dirn)) + '=0;i') + str(dirn)) + '<Nxx_plus_2NGHOSTS') + str(dirn)) + ';i') + str(dirn)) + '++) {\n')
                    define_str += (((((('  const REAL xx' + str(dirn)) + ' = xx[') + str(dirn)) + '][i') + str(dirn)) + '];\n')
                    define_str += (((((('  rfmstruct.' + str(freevars_uniq_xx_indep[which_freevar])) + '[i') + str(dirn)) + '] = ') + str(sp.ccode(freevars_uniq_vals[which_freevar]))) + ';\n')
                    define_str += '}\n\n'
                    readvr_str[dirn] += (((((('const REAL ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i') + str(dirn)) + '];\n')
                    readvr_SIMD_outer_str[dirn] += (((((('const double NOSIMD' + str(freevars_uniq_xx_indep[which_freevar])) + ' = rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i') + str(dirn)) + ']; ')
                    readvr_SIMD_outer_str[dirn] += (((('const REAL_SIMD_ARRAY ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = ConstSIMD(NOSIMD') + str(freevars_uniq_xx_indep[which_freevar])) + ');\n')
                    readvr_SIMD_inner_str[dirn] += (((((('const REAL_SIMD_ARRAY ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = ReadSIMD(&rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i') + str(dirn)) + ']);\n')
                    output_define_and_readvr = True
            if ((not output_define_and_readvr) and (gri.xx[0] in frees_uniq) and (gri.xx[1] in frees_uniq)):
                define_str += (((('\nfor(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++) for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) {\n  const REAL xx0 = xx[0][i0];\n  const REAL xx1 = xx[1][i1];\n  rfmstruct.' + str(freevars_uniq_xx_indep[which_freevar])) + '[i0 + Nxx_plus_2NGHOSTS0*i1] = ') + str(sp.ccode(freevars_uniq_vals[which_freevar]))) + ';\n}\n\n')
                readvr_str[0] += (((('const REAL ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i0 + Nxx_plus_2NGHOSTS0*i1];\n')
                readvr_SIMD_outer_str[0] += (((('const double NOSIMD' + str(freevars_uniq_xx_indep[which_freevar])) + ' = rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i0 + Nxx_plus_2NGHOSTS0*i1]; ')
                readvr_SIMD_outer_str[0] += (((('const REAL_SIMD_ARRAY ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = ConstSIMD(NOSIMD') + str(freevars_uniq_xx_indep[which_freevar])) + ');\n')
                readvr_SIMD_inner_str[0] += (((('const REAL_SIMD_ARRAY ' + str(freevars_uniq_xx_indep[which_freevar])) + ' = ReadSIMD(&rfmstruct->') + str(freevars_uniq_xx_indep[which_freevar])) + '[i0 + Nxx_plus_2NGHOSTS0*i1]);\n')
                output_define_and_readvr = True
            if (not output_define_and_readvr):
                print((('ERROR: Could not figure out the (xx0,xx1,xx2) dependency within the expression for ' + str(freevars_uniq_xx_indep[which_freevar])) + ':'))
                print(str(freevars_uniq_vals[which_freevar]))
                sys.exit(1)
        which_freevar += 1
    struct_str += '} rfm_struct;\n'
    outdir = par.parval_from_str((thismodule + '::rfm_precompute_Ccode_outdir'))
    if (par.parval_from_str((thismodule + '::rfm_precompute_to_Cfunctions_and_NRPy_basic_defines')) == 'False'):
        with open(os.path.join(outdir, 'rfm_struct__declare.h'), 'w') as file:
            file.write(struct_str)
        with open(os.path.join(outdir, 'rfm_struct__malloc.h'), 'w') as file:
            file.write(malloc_str)
        with open(os.path.join(outdir, 'rfm_struct__define.h'), 'w') as file:
            file.write(define_str)
        with open(os.path.join(outdir, 'rfm_struct__freemem.h'), 'w') as file:
            file.write(freemm_str)
    else:
        global NRPy_basic_defines_str
        NRPy_basic_defines_str = struct_str
        global rfm_struct__malloc, rfm_struct__define, rfm_struct__freemem
        rfm_struct__malloc = malloc_str
        rfm_struct__define = define_str
        rfm_struct__freemem = freemm_str
    for i in range(3):
        with open(os.path.join(outdir, (('rfm_struct__read' + str(i)) + '.h')), 'w') as file:
            file.write(readvr_str[i])
        with open(os.path.join(outdir, (('rfm_struct__SIMD_outer_read' + str(i)) + '.h')), 'w') as file:
            file.write(readvr_SIMD_outer_str[i])
        with open(os.path.join(outdir, (('rfm_struct__SIMD_inner_read' + str(i)) + '.h')), 'w') as file:
            file.write(readvr_SIMD_inner_str[i])