Fixed-point iteration brute-force finder.

zinder.zig

//! A brute-force algorithm for finding a fixed-point iteration formula for an arbitrary constant. It can be used to
//! find the Babylonian method for finding the square root of a number, or to find a way to divide using addition,
//! subtraction and multiplication. More complex constants like pi are already too much for this algorithm, so overall
//! it's pretty useless.
//!
//! Tested with Zig 0.11.0-dev.4191+1bf16b172.

const std = @import("std");

fn neg(a: f64) f64 {
    return -a;
}

fn rec(a: f64) f64 {
    return 1 / a;
}

fn inc(a: f64) f64 {
    return a + 1;
}

fn dec(a: f64) f64 {
    return a - 1;
}

fn double(a: f64) f64 {
    return a * 2;
}

fn halve(a: f64) f64 {
    return a / 2;
}

fn square(a: f64) f64 {
    return a * a;
}

fn alt(a: f64) f64 {
    return if (@mod(a, 2) == 0) 1 else -1;
}

fn sqrt(a: f64) f64 {
    return std.math.sqrt(a);
}

fn add(a: f64, b: f64) f64 {
    return a + b;
}

fn sub(a: f64, b: f64) f64 {
    return a - b;
}

fn mul(a: f64, b: f64) f64 {
    return a * b;
}

fn div(a: f64, b: f64) f64 {
    return a / b;
}

fn pow(a: f64, b: f64) f64 {
    return std.math.pow(f64, a, b);
}

// Modify the following 3 lines to determine what is and isn't allowed in the fixed-point iteration formulas. The more
// stuff there is, the slower the search will end up being.
const consts = [_]f64{ 1, 2 };
const unary_funcs = [_]*const fn (f64) f64{ neg, rec, double, halve, inc, dec };
const binary_funcs = [_]*const fn (f64, f64) f64{ add, sub, mul, div };

const Op = union(enum) {
    push_prev,
    push_iter,
    push_const: u8,
    unary: u8,
    binary: u8,
};

const stack_size = std.math.maxInt(u8);
fn run(ops: []const Op, init: f64, goal: f64, epsilon: f64, max_iter: usize) bool {
    var stack: [stack_size]f64 = undefined;
    var stack_idx: u8 = stack_size;
    var prev = init;
    for (0..max_iter) |iter| {
        var i: u8 = 0;
        while (i < ops.len) : (i += 1) {
            switch (ops[i]) {
                .push_prev => {
                    stack_idx -= 1;
                    stack[stack_idx] = prev;
                },
                .push_iter => {
                    stack_idx -= 1;
                    stack[stack_idx] = @floatFromInt(iter);
                },
                .push_const => |j| if (consts.len != 0) {
                    stack_idx -= 1;
                    stack[stack_idx] = consts[j];
                },
                .unary => |j| if (unary_funcs.len != 0) {
                    stack[stack_idx] = unary_funcs[j](stack[stack_idx]);
                },
                .binary => |j| if (binary_funcs.len != 0) {
                    stack_idx += 1;
                    stack[stack_idx] = binary_funcs[j](stack[stack_idx], stack[stack_idx - 1]);
                },
            }
        }
        const curr = stack[stack_idx];
        if (@fabs(curr - goal) > @fabs(prev - goal)) {
            return false;
        } else if (@fabs(curr - goal) < epsilon) {
            return true;
        }
        prev = curr;
        stack_idx += 1;
    }
    return false;
}

fn findImpl(ops: []Op, init: f64, goal: f64, items: u8, curr_depth: u8, max_depth: u8, epsilon: f64, max_iter: usize) bool {
    if (curr_depth == max_depth) {
        return run(ops, init, goal, epsilon, max_iter);
    }

    for (std.enums.values(@typeInfo(Op).Union.tag_type.?)) |op_type| {
        const depth_left = max_depth - curr_depth;
        switch (op_type) {
            .push_prev, .push_iter => if (depth_left > items) {
                ops[curr_depth] = if (op_type == .push_prev) .push_prev else .push_iter;
                const found = findImpl(ops, init, goal, items + 1, curr_depth + 1, max_depth, epsilon, max_iter);
                if (found) return true;
            },
            .push_const => if (depth_left >= items + 1) {
                var i: u8 = 0;
                while (i < consts.len) : (i += 1) {
                    ops[curr_depth] = .{ .push_const = i };
                    const found = findImpl(ops, init, goal, items + 1, curr_depth + 1, max_depth, epsilon, max_iter);
                    if (found) return true;
                }
            },
            .unary => if (items >= 1 and depth_left >= items) {
                var i: u8 = 0;
                while (i < unary_funcs.len) : (i += 1) {
                    ops[curr_depth] = .{ .unary = i };
                    const found = findImpl(ops, init, goal, items, curr_depth + 1, max_depth, epsilon, max_iter);
                    if (found) return true;
                }
            },
            .binary => if (items >= 2 and depth_left >= items - 1) {
                var i: u8 = 0;
                while (i < binary_funcs.len) : (i += 1) {
                    ops[curr_depth] = .{ .binary = i };
                    const found = findImpl(ops, init, goal, items - 1, curr_depth + 1, max_depth, epsilon, max_iter);
                    if (found) return true;
                }
            },
        }
    }
    return false;
}

fn find(ops: []Op, init: f64, goal: f64, max_depth: u8, epsilon: f64, max_iter: usize) ?u8 {
    var depth: u8 = 1;
    while (depth <= max_depth) : (depth += 1) {
        const found = findImpl(ops[0..depth], init, goal, 0, 0, depth, epsilon, max_iter);
        if (found) return depth;
    }
    return null;
}

pub fn main() !void {
    var ops: [stack_size]Op = undefined;
    var timer = try std.time.Timer.start();

    // Modify the line below to change the number that the fixed-point iteration should converge to, as well as how fast
    // it should converge to it and what number it starts with.
    const found = find(&ops, 2, std.math.sqrt2, 10, 1e-14, 10);

    const time = timer.read();
    if (found) |depth| {
        for (ops[0..depth]) |op| {
            std.debug.print("{any}\n", .{op});
        }
    } else {
        std.debug.print("Not found.\n", .{});
    }
    std.debug.print("Took {d} seconds\n", .{@as(f64, @floatFromInt(time)) / std.time.ns_per_s});
}