Problem 620

Planetary Gears: G(500). Radii scale with circumference; with c = s+p+q the p-planet triangle (D, s+q, s+p) and the q-planet triangle are the same triangle with the two sides swapped, so the two meshing phases satisfy A_p + A_q = s + (p+q)/2 identically and the whole assembly reduces to one condition: A(D) = b*u + a*(pi - w) must be a multiple of 1/2 (a = (s+p)/2pi, b = (s+q)/2pi, u and w the triangle angles at the two gear centres). A increases monotonically from 0 as the centre distance D grows from q-p side-by-side degeneracy (excluded: the two equal planets would coincide) up to the 1cm-gap bound D = (p+q)/2pi - 1, so g = floor(2 * A(Dmax)). Verified: g(16,5,5,6) = 9, G(16) = 9, G(20) = 205.

View problem on Project Euler

Performance comparison

Flow source

# Project Euler 620
# Planetary Gears: G(500).
#
# Radii scale with circumference; with c = s+p+q the p-planet triangle
# (D, s+q, s+p) and the q-planet triangle are the same triangle with the
# two sides swapped, so the two meshing phases satisfy
# A_p + A_q = s + (p+q)/2 identically and the whole assembly reduces to
# one condition: A(D) = b*u + a*(pi - w) must be a multiple of 1/2
# (a = (s+p)/2pi, b = (s+q)/2pi, u and w the triangle angles at the two
# gear centres).  A increases monotonically from 0 as the centre distance
# D grows from q-p side-by-side degeneracy (excluded: the two equal
# planets would coincide) up to the 1cm-gap bound D = (p+q)/2pi - 1, so
#   g = floor(2 * A(Dmax)).
# Verified: g(16,5,5,6) = 9, G(16) = 9, G(20) = 205.

extern {
    function calloc(n: i64, size: i64) -> ptr<void>
    function free(p: ptr<void>) -> void
    function acos(x: f64) -> f64
}

const N: i64 = 500
const PI: f64 = 3.14159265358979323846

function g_val(s: i64, p: i64, q: i64) -> i64 {
    let a: f64 = ((s + p) as f64) / (2.0 * PI)
    let b: f64 = ((s + q) as f64) / (2.0 * PI)
    let Dlo: f64 = b - a
    let D: f64 = ((p + q) as f64) / (2.0 * PI) - 1.0
    if D <= Dlo {
        return 0
    }
    let mut cu: f64 = (D * D + b * b - a * a) / (2.0 * D * b)
    let mut cw: f64 = (D * D + a * a - b * b) / (2.0 * D * a)
    if cu > 1.0 {
        cu = 1.0
    }
    if cu < 0.0 - 1.0 {
        cu = 0.0 - 1.0
    }
    if cw > 1.0 {
        cw = 1.0
    }
    if cw < 0.0 - 1.0 {
        cw = 0.0 - 1.0
    }
    let u: f64 = acos(cu)
    let w: f64 = acos(cw)
    let A: f64 = b * u + a * (PI - w)
    return (2.0 * A + 0.000000000001) as i64
}

function main() -> i32 {
    let mut tot: i64 = 0
    let mut s: i64 = 5
    while s + 11 <= N {
        let mut p: i64 = 5
        while s + 2 * p + 1 <= N {
            let mut q: i64 = p + 1
            while s + p + q <= N {
                tot = tot + g_val(s, p, q)
                q = q + 1
            }
            p = p + 1
        }
        s = s + 1
    }
    printf("%lld\n", tot)
    return 0
}

Generated C

#include <stdint.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/* Flow runtime helpers */
typedef struct flow_temp_node { struct flow_temp_node* next; } flow_temp_node;
static flow_temp_node* flow_temp_head = NULL;
static int flow_temp_atexit_set = 0;
__attribute__((unused)) static void flow_temp_free_all(void) {
    while (flow_temp_head) {
        flow_temp_node* n = flow_temp_head;
        flow_temp_head = n->next;
        free(n);
    }
}
__attribute__((unused)) static void* flow_temp_alloc(size_t nbytes) {
    flow_temp_node* node = (flow_temp_node*)malloc(sizeof(flow_temp_node) + nbytes);
    if (!node) return NULL;
    node->next = flow_temp_head;
    flow_temp_head = node;
    if (!flow_temp_atexit_set) {
        flow_temp_atexit_set = 1;
        atexit(flow_temp_free_all);
    }
    return (void*)(node + 1);
}
#ifndef FLOW_DIAG
#define FLOW_DIAG(msg) fprintf(stderr, "%s", (msg))
#endif
#ifndef FLOW_LOG
#define FLOW_LOG(fmt, ...) printf(fmt, __VA_ARGS__)
#endif
#ifndef FLOW_LOG_EMPTY
#define FLOW_LOG_EMPTY(fmt) printf(fmt)
#endif
static char* flow_strcat(const char* a, const char* b) {
    size_t la = strlen(a ? a : ""), lb = strlen(b ? b : "");
    char* r = (char*)flow_temp_alloc(la + lb + 1);
    if (!r) return NULL;
    if (la) memcpy(r, a, la);
    if (lb) memcpy(r + la, b, lb);
    r[la + lb] = '\0';
    return r;
}

#define __flow_in_arr(arr, val) __extension__ ({ \
    int _found = 0; \
    size_t _n = sizeof(arr)/sizeof((arr)[0]); \
    for (size_t _i = 0; _i < _n; _i++) { \
        if ((arr)[_i] == (val)) { _found = 1; break; } \
    } _found; })

/* Unified fault handler (MISRA #279) — override with -DFLOW_FAULT_HANDLER=fn */
#ifndef FLOW_FAULT_HANDLER
__attribute__((unused)) static inline void flow_fault_handler(const char* msg) {
    fprintf(stderr, "flow: %s\n", msg ? msg : "fault");
    abort();
#if defined(__GNUC__) || defined(__clang__)
    __builtin_unreachable();
#endif
}
#else
#define flow_fault_handler FLOW_FAULT_HANDLER
#endif
#define flow_div_by_zero_handler() flow_fault_handler("division by zero")
#define flow_shift_ub_handler() flow_fault_handler("invalid shift (amount out of range or left-shift of negative)")

#ifndef FLOW_CHECKED_DIV
#define FLOW_CHECKED_DIV(L, R) (((R) != 0) ? ((L) / (R)) : (flow_div_by_zero_handler(), (L) * 0))
#endif
#ifndef FLOW_CHECKED_MOD
#define FLOW_CHECKED_MOD(L, R) (((R) != 0) ? ((L) % (R)) : (flow_div_by_zero_handler(), (L) * 0))
#endif
#ifndef FLOW_CHECKED_SHL
#define FLOW_CHECKED_SHL(L, R) ((((R) >= 0) && ((unsigned long long)(R) < (sizeof(L) * 8ull)) && ((L) >= 0)) ? ((L) << (R)) : (flow_shift_ub_handler(), (L) * 0))
#endif
#ifndef FLOW_CHECKED_SHR
#define FLOW_CHECKED_SHR(L, R) ((((R) >= 0) && ((unsigned long long)(R) < (sizeof(L) * 8ull))) ? ((L) >> (R)) : (flow_shift_ub_handler(), (L) * 0))
#endif

#include <math.h>

void* _ui_state = NULL;

static inline float i32_to_f32(int32_t v) { return (float)v; }

/* Host stub for @gpu kernels (device codegen replaces this). */
static inline int32_t gpu_thread_id(void) { return 0; }

double acos(double x);
int64_t g_val_i64_i64_i64(int64_t s, int64_t p, int64_t q);
int32_t main(void);

static const int64_t N = 500;
static const double PI = 3.14159265358979323846;




int64_t g_val_i64_i64_i64(int64_t s, int64_t p, int64_t q) {
    double a = (((double)((s + p))) / (2.0 * PI));
    double b = (((double)((s + q))) / (2.0 * PI));
    double Dlo = (b - a);
    double D = ((((double)((p + q))) / (2.0 * PI)) - 1.0);
    if (D <= Dlo) {
        return 0;
    }
    double cu = ((((D * D) + (b * b)) - (a * a)) / ((2.0 * D) * b));
    double cw = ((((D * D) + (a * a)) - (b * b)) / ((2.0 * D) * a));
    if (cu > 1.0) {
        cu = 1.0;
    }
    if (cu < (0.0 - 1.0)) {
        cu = (0.0 - 1.0);
    }
    if (cw > 1.0) {
        cw = 1.0;
    }
    if (cw < (0.0 - 1.0)) {
        cw = (0.0 - 1.0);
    }
    double u = acos(cu);
    double w = acos(cw);
    double A = ((b * u) + (a * (PI - w)));
    return ((int64_t)(((2.0 * A) + 0.000000000001)));
}

int32_t main(void) {
    int64_t tot = 0;
    int64_t s = 5;
    while ((s + 11) <= N) {
        int64_t p = 5;
        while (((s + (2 * p)) + 1) <= N) {
            int64_t q = (p + 1);
            while (((s + p) + q) <= N) {
                tot = (tot + g_val_i64_i64_i64(s, p, q));
                q = (q + 1);
            }
            p = (p + 1);
        }
        s = (s + 1);
    }
    printf("%lld\n", tot);
    return 0;
}

Generated MLIR

module {
  llvm.func @printf(!llvm.ptr, ...) -> i32
  llvm.mlir.global internal constant @str_0("%lld\n\00") {addr_space = 0 : i32} : !llvm.array<6 x i8>
  func.func private @calloc(i64, i64) -> !llvm.ptr
  func.func private @free(!llvm.ptr) -> ()
  func.func private @acos(f64) -> f64
  // Constant: N
  llvm.mlir.global internal constant @N(500 : i64) : i64
  // Constant: PI
  llvm.mlir.global internal constant @PI(3.14159265358979323846 : f64) : f64
  func.func @g_val(%arg0: i64, %arg1: i64, %arg2: i64) -> i64 {
    %0 = arith.addi %arg0, %arg1 : i64
    %1 = arith.sitofp %0 : i64 to f64
    %2 = arith.constant 2.0 : f32
    %3 = llvm.mlir.addressof @PI : !llvm.ptr
    %4 = llvm.load %3 : !llvm.ptr -> f64
    %6 = arith.extf %2 : f32 to f64
    %5 = arith.mulf %6, %4 : f64
    %7 = arith.divf %1, %5 : f64
    %8 = arith.addi %arg0, %arg2 : i64
    %9 = arith.sitofp %8 : i64 to f64
    %10 = arith.constant 2.0 : f32
    %11 = llvm.mlir.addressof @PI : !llvm.ptr
    %12 = llvm.load %11 : !llvm.ptr -> f64
    %14 = arith.extf %10 : f32 to f64
    %13 = arith.mulf %14, %12 : f64
    %15 = arith.divf %9, %13 : f64
    %16 = arith.subf %15, %7 : f64
    %17 = arith.addi %arg1, %arg2 : i64
    %18 = arith.sitofp %17 : i64 to f64
    %19 = arith.constant 2.0 : f32
    %20 = llvm.mlir.addressof @PI : !llvm.ptr
    %21 = llvm.load %20 : !llvm.ptr -> f64
    %23 = arith.extf %19 : f32 to f64
    %22 = arith.mulf %23, %21 : f64
    %24 = arith.divf %18, %22 : f64
    %25 = arith.constant 1.0 : f32
    %27 = arith.extf %25 : f32 to f64
    %26 = arith.subf %24, %27 : f64
    %28 = arith.cmpf ole, %26, %16 : f64
    cf.cond_br %28, ^bb0, ^bb1
    ^bb0:
      %29 = arith.constant 0 : i32
      %30 = arith.extsi %29 : i32 to i64
      func.return %30 : i64
    ^bb1:
      cf.br ^bb2
    ^bb2:
    %31 = arith.mulf %26, %26 : f64
    %32 = arith.mulf %15, %15 : f64
    %33 = arith.addf %31, %32 : f64
    %34 = arith.mulf %7, %7 : f64
    %35 = arith.subf %33, %34 : f64
    %36 = arith.constant 2.0 : f32
    %38 = arith.extf %36 : f32 to f64
    %37 = arith.mulf %38, %26 : f64
    %39 = arith.mulf %37, %15 : f64
    %40 = arith.divf %35, %39 : f64
    %41 = llvm.mlir.constant(1 : i64) : i64
    %42 = llvm.alloca %41 x f64 : (i64) -> !llvm.ptr
    llvm.store %40, %42 : f64, !llvm.ptr
    %43 = arith.mulf %26, %26 : f64
    %44 = arith.mulf %7, %7 : f64
    %45 = arith.addf %43, %44 : f64
    %46 = arith.mulf %15, %15 : f64
    %47 = arith.subf %45, %46 : f64
    %48 = arith.constant 2.0 : f32
    %50 = arith.extf %48 : f32 to f64
    %49 = arith.mulf %50, %26 : f64
    %51 = arith.mulf %49, %7 : f64
    %52 = arith.divf %47, %51 : f64
    %53 = llvm.mlir.constant(1 : i64) : i64
    %54 = llvm.alloca %53 x f64 : (i64) -> !llvm.ptr
    llvm.store %52, %54 : f64, !llvm.ptr
    %55 = llvm.load %42 : !llvm.ptr -> f64
    %56 = arith.constant 1.0 : f32
    %58 = arith.extf %56 : f32 to f64
    %57 = arith.cmpf ogt, %55, %58 : f64
    cf.cond_br %57, ^bb3, ^bb4
    ^bb3:
      %59 = arith.constant 1.0 : f32
      %60 = arith.extf %59 : f32 to f64
      llvm.store %60, %42 : f64, !llvm.ptr
      cf.br ^bb5
    ^bb4:
      cf.br ^bb5
    ^bb5:
    %61 = llvm.load %42 : !llvm.ptr -> f64
    %62 = arith.constant 0.0 : f32
    %63 = arith.constant 1.0 : f32
    %64 = arith.subf %62, %63 : f32
    %66 = arith.extf %64 : f32 to f64
    %65 = arith.cmpf olt, %61, %66 : f64
    cf.cond_br %65, ^bb6, ^bb7
    ^bb6:
      %67 = arith.constant 0.0 : f32
      %68 = arith.constant 1.0 : f32
      %69 = arith.subf %67, %68 : f32
      %70 = arith.extf %69 : f32 to f64
      llvm.store %70, %42 : f64, !llvm.ptr
      cf.br ^bb8
    ^bb7:
      cf.br ^bb8
    ^bb8:
    %71 = llvm.load %54 : !llvm.ptr -> f64
    %72 = arith.constant 1.0 : f32
    %74 = arith.extf %72 : f32 to f64
    %73 = arith.cmpf ogt, %71, %74 : f64
    cf.cond_br %73, ^bb9, ^bb10
    ^bb9:
      %75 = arith.constant 1.0 : f32
      %76 = arith.extf %75 : f32 to f64
      llvm.store %76, %54 : f64, !llvm.ptr
      cf.br ^bb11
    ^bb10:
      cf.br ^bb11
    ^bb11:
    %77 = llvm.load %54 : !llvm.ptr -> f64
    %78 = arith.constant 0.0 : f32
    %79 = arith.constant 1.0 : f32
    %80 = arith.subf %78, %79 : f32
    %82 = arith.extf %80 : f32 to f64
    %81 = arith.cmpf olt, %77, %82 : f64
    cf.cond_br %81, ^bb12, ^bb13
    ^bb12:
      %83 = arith.constant 0.0 : f32
      %84 = arith.constant 1.0 : f32
      %85 = arith.subf %83, %84 : f32
      %86 = arith.extf %85 : f32 to f64
      llvm.store %86, %54 : f64, !llvm.ptr
      cf.br ^bb14
    ^bb13:
      cf.br ^bb14
    ^bb14:
    %88 = llvm.load %42 : !llvm.ptr -> f64
    %87 = func.call @acos(%88) : (f64) -> f64
    %90 = llvm.load %54 : !llvm.ptr -> f64
    %89 = func.call @acos(%90) : (f64) -> f64
    %91 = arith.mulf %15, %87 : f64
    %92 = llvm.mlir.addressof @PI : !llvm.ptr
    %93 = llvm.load %92 : !llvm.ptr -> f64
    %94 = arith.subf %93, %89 : f64
    %95 = arith.mulf %7, %94 : f64
    %96 = arith.addf %91, %95 : f64
    %97 = arith.constant 2.0 : f32
    %99 = arith.extf %97 : f32 to f64
    %98 = arith.mulf %99, %96 : f64
    %100 = arith.constant 0.000000000001 : f32
    %102 = arith.extf %100 : f32 to f64
    %101 = arith.addf %98, %102 : f64
    %103 = arith.fptosi %101 : f64 to i64
    func.return %103 : i64
  }
  func.func @main() -> i32 {
    %104 = arith.constant 0 : i32
    %105 = arith.extsi %104 : i32 to i64
    %106 = llvm.mlir.constant(1 : i64) : i64
    %107 = llvm.alloca %106 x i64 : (i64) -> !llvm.ptr
    llvm.store %105, %107 : i64, !llvm.ptr
    %108 = arith.constant 5 : i32
    %109 = arith.extsi %108 : i32 to i64
    %110 = llvm.mlir.constant(1 : i64) : i64
    %111 = llvm.alloca %110 x i64 : (i64) -> !llvm.ptr
    llvm.store %109, %111 : i64, !llvm.ptr
    cf.br ^bb15
    ^bb15:
    %112 = llvm.load %111 : !llvm.ptr -> i64
    %113 = arith.constant 11 : i32
    %115 = arith.extsi %113 : i32 to i64
    %114 = arith.addi %112, %115 : i64
    %116 = llvm.mlir.addressof @N : !llvm.ptr
    %117 = llvm.load %116 : !llvm.ptr -> i64
    %118 = arith.cmpi sle, %114, %117 : i64
    cf.cond_br %118, ^bb16, ^bb17
    ^bb16:
      %119 = arith.constant 5 : i32
      %120 = arith.extsi %119 : i32 to i64
      %121 = llvm.mlir.constant(1 : i64) : i64
      %122 = llvm.alloca %121 x i64 : (i64) -> !llvm.ptr
      llvm.store %120, %122 : i64, !llvm.ptr
      cf.br ^bb18
      ^bb18:
      %123 = llvm.load %111 : !llvm.ptr -> i64
      %124 = arith.constant 2 : i32
      %125 = llvm.load %122 : !llvm.ptr -> i64
      %127 = arith.extsi %124 : i32 to i64
      %126 = arith.muli %127, %125 : i64
      %128 = arith.addi %123, %126 : i64
      %129 = arith.constant 1 : i32
      %131 = arith.extsi %129 : i32 to i64
      %130 = arith.addi %128, %131 : i64
      %132 = llvm.mlir.addressof @N : !llvm.ptr
      %133 = llvm.load %132 : !llvm.ptr -> i64
      %134 = arith.cmpi sle, %130, %133 : i64
      cf.cond_br %134, ^bb19, ^bb20
      ^bb19:
        %135 = llvm.load %122 : !llvm.ptr -> i64
        %136 = arith.constant 1 : i32
        %138 = arith.extsi %136 : i32 to i64
        %137 = arith.addi %135, %138 : i64
        %139 = llvm.mlir.constant(1 : i64) : i64
        %140 = llvm.alloca %139 x i64 : (i64) -> !llvm.ptr
        llvm.store %137, %140 : i64, !llvm.ptr
        cf.br ^bb21
        ^bb21:
        %141 = llvm.load %111 : !llvm.ptr -> i64
        %142 = llvm.load %122 : !llvm.ptr -> i64
        %143 = arith.addi %141, %142 : i64
        %144 = llvm.load %140 : !llvm.ptr -> i64
        %145 = arith.addi %143, %144 : i64
        %146 = llvm.mlir.addressof @N : !llvm.ptr
        %147 = llvm.load %146 : !llvm.ptr -> i64
        %148 = arith.cmpi sle, %145, %147 : i64
        cf.cond_br %148, ^bb22, ^bb23
        ^bb22:
          %149 = llvm.load %107 : !llvm.ptr -> i64
          %151 = llvm.load %111 : !llvm.ptr -> i64
          %152 = llvm.load %122 : !llvm.ptr -> i64
          %153 = llvm.load %140 : !llvm.ptr -> i64
          %150 = func.call @g_val(%151, %152, %153) : (i64, i64, i64) -> i64
          %154 = arith.addi %149, %150 : i64
          llvm.store %154, %107 : i64, !llvm.ptr
          %155 = llvm.load %140 : !llvm.ptr -> i64
          %156 = arith.constant 1 : i32
          %158 = arith.extsi %156 : i32 to i64
          %157 = arith.addi %155, %158 : i64
          llvm.store %157, %140 : i64, !llvm.ptr
          cf.br ^bb21
        ^bb23:
        %159 = llvm.load %122 : !llvm.ptr -> i64
        %160 = arith.constant 1 : i32
        %162 = arith.extsi %160 : i32 to i64
        %161 = arith.addi %159, %162 : i64
        llvm.store %161, %122 : i64, !llvm.ptr
        cf.br ^bb18
      ^bb20:
      %163 = llvm.load %111 : !llvm.ptr -> i64
      %164 = arith.constant 1 : i32
      %166 = arith.extsi %164 : i32 to i64
      %165 = arith.addi %163, %166 : i64
      llvm.store %165, %111 : i64, !llvm.ptr
      cf.br ^bb15
    ^bb17:
    %167 = llvm.mlir.addressof @str_0 : !llvm.ptr
    %168 = llvm.load %107 : !llvm.ptr -> i64
    %169 = llvm.call @printf(%167, %168) vararg(!llvm.func<i32 (ptr, ...)>) : (!llvm.ptr, i64) -> i32
    %170 = arith.constant 0 : i32
    func.return %170 : i32
  }
}