Problem 226
Blancmange curve vs circle intersection area.
View problem on Project Euler
Performance comparison
| Metric | Our solution | Best known |
|---|
| Time complexity | O(n) | O(n) |
|---|
| Space complexity | O(1) | O(1) |
|---|
| Approach | Flow solution | Numerical iteration |
|---|
| Verdict | Optimal |
|---|
Flow source
# Project Euler 226
# Blancmange curve vs circle intersection area.
extern {
function sqrt(x: f64) -> f64
function floor(x: f64) -> f64
function fabs(x: f64) -> f64
function pow(x: f64, y: f64) -> f64
}
const EPSILON: f64 = 0.00000001
function s_curve(x: f64) -> f64 {
let mut result: f64 = 0.0
let mut n: i32 = 0
while true {
let power: f64 = pow(2.0, n as f64)
let parameter: f64 = power * x
let mut s_val: f64 = parameter - floor(parameter)
if s_val > 0.5 { s_val = 1.0 - s_val }
let add: f64 = s_val / power
result = result + add
if add < EPSILON { return result }
n = n + 1
}
return result
}
function find_intersection(cx: f64, cy: f64, radius: f64, x0: f64, step0: f64) -> f64 {
let mut x: f64 = x0
let mut step: f64 = step0
while true {
let y: f64 = s_curve(x)
let dx: f64 = x - cx
let dy: f64 = y - cy
let distance: f64 = sqrt(dx * dx + dy * dy)
if fabs(distance - radius) < EPSILON { return x }
let mut turn: bool = false
if distance < radius {
if dx > 0.0 && step < 0.0 { turn = true }
if dx < 0.0 && step > 0.0 { turn = true }
} else {
if dx > 0.0 && step > 0.0 { turn = true }
if dx < 0.0 && step < 0.0 { turn = true }
}
if turn { step = 0.0 - step / 2.0 }
x = x + step
}
return x
}
function integrate(cx: f64, cy: f64, radius: f64, from_x: f64, to_x: f64, step: f64) -> f64 {
let mut result: f64 = 0.0
let mut x: f64 = from_x
while x <= to_x {
let upper: f64 = s_curve(x)
let lower: f64 = cy - sqrt(radius * radius - (x - cx) * (x - cx))
result = result + (upper - lower) * step
x = x + step
}
return result
}
function main() -> i32 {
let cx: f64 = 0.25
let cy: f64 = 0.5
let radius: f64 = 0.25
let from_x: f64 = find_intersection(cx, cy, radius, cx, 0.0 - 0.1)
let to_x: f64 = find_intersection(cx, cy, radius, cx, 0.1)
let area: f64 = integrate(cx, cy, radius, from_x, to_x, 0.00001)
printf("%.8f\n", area)
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 s_curve_f64(double x);
double find_intersection_f64_f64_f64_f64_f64(double cx, double cy, double radius, double x0, double step0);
double integrate_f64_f64_f64_f64_f64_f64(double cx, double cy, double radius, double from_x, double to_x, double step);
int32_t main(void);
static const double EPSILON = 0.00000001;
double s_curve_f64(double x) {
double result = 0.0;
int32_t n = 0;
while (1) {
double power = pow(2.0, ((double)(n)));
double parameter = (power * x);
double s_val = (parameter - floor(parameter));
if (s_val > 0.5) {
s_val = (1.0 - s_val);
}
double add = (s_val / power);
result = (result + add);
if (add < EPSILON) {
return result;
}
n = (n + 1);
}
return result;
}
double find_intersection_f64_f64_f64_f64_f64(double cx, double cy, double radius, double x0, double step0) {
double x = x0;
double step = step0;
while (1) {
double y = s_curve_f64(x);
double dx = (x - cx);
double dy = (y - cy);
double distance = sqrt(((dx * dx) + (dy * dy)));
if (fabs((distance - radius)) < EPSILON) {
return x;
}
bool turn = 0;
if (distance < radius) {
if ((dx > 0.0 && step < 0.0)) {
turn = 1;
}
if ((dx < 0.0 && step > 0.0)) {
turn = 1;
}
} else {
if ((dx > 0.0 && step > 0.0)) {
turn = 1;
}
if ((dx < 0.0 && step < 0.0)) {
turn = 1;
}
}
if (turn) {
step = (0.0 - (step / 2.0));
}
x = (x + step);
}
return x;
}
double integrate_f64_f64_f64_f64_f64_f64(double cx, double cy, double radius, double from_x, double to_x, double step) {
double result = 0.0;
double x = from_x;
while (x <= to_x) {
double upper = s_curve_f64(x);
double lower = (cy - sqrt(((radius * radius) - ((x - cx) * (x - cx)))));
result = (result + ((upper - lower) * step));
x = (x + step);
}
return result;
}
int32_t main(void) {
double cx = 0.25;
double cy = 0.5;
double radius = 0.25;
double from_x = find_intersection_f64_f64_f64_f64_f64(cx, cy, radius, cx, (0.0 - 0.1));
double to_x = find_intersection_f64_f64_f64_f64_f64(cx, cy, radius, cx, 0.1);
double area = integrate_f64_f64_f64_f64_f64_f64(cx, cy, radius, from_x, to_x, 0.00001);
printf("%.8f\n", area);
return 0;
}
Generated MLIR
module {
llvm.func @printf(!llvm.ptr, ...) -> i32
llvm.mlir.global internal constant @str_0("%.8f\n\00") {addr_space = 0 : i32} : !llvm.array<6 x i8>
func.func private @sqrt(f64) -> f64
func.func private @floor(f64) -> f64
func.func private @fabs(f64) -> f64
func.func private @pow(f64, f64) -> f64
// Constant: EPSILON
llvm.mlir.global internal constant @EPSILON(0.00000001 : f64) : f64
func.func @s_curve(%arg0: f64) -> f64 {
%0 = arith.constant 0.0 : f32
%1 = arith.extf %0 : f32 to f64
%2 = llvm.mlir.constant(1 : i64) : i64
%3 = llvm.alloca %2 x f64 : (i64) -> !llvm.ptr
llvm.store %1, %3 : f64, !llvm.ptr
%4 = arith.constant 0 : i32
%5 = llvm.mlir.constant(1 : i64) : i64
%6 = llvm.alloca %5 x i32 : (i64) -> !llvm.ptr
llvm.store %4, %6 : i32, !llvm.ptr
cf.br ^bb0
^bb0:
%7 = arith.constant 1 : i1
cf.cond_br %7, ^bb1, ^bb2
^bb1:
%9 = arith.constant 2.0 : f32
%10 = llvm.load %6 : !llvm.ptr -> i32
%11 = arith.sitofp %10 : i32 to f64
%12 = arith.extf %9 : f32 to f64
%8 = func.call @pow(%12, %11) : (f64, f64) -> f64
%13 = arith.mulf %8, %arg0 : f64
%14 = func.call @floor(%13) : (f64) -> f64
%15 = arith.subf %13, %14 : f64
%16 = llvm.mlir.constant(1 : i64) : i64
%17 = llvm.alloca %16 x f64 : (i64) -> !llvm.ptr
llvm.store %15, %17 : f64, !llvm.ptr
%18 = llvm.load %17 : !llvm.ptr -> f64
%19 = arith.constant 0.5 : f32
%21 = arith.extf %19 : f32 to f64
%20 = arith.cmpf ogt, %18, %21 : f64
cf.cond_br %20, ^bb3, ^bb4
^bb3:
%22 = arith.constant 1.0 : f32
%23 = llvm.load %17 : !llvm.ptr -> f64
%25 = arith.extf %22 : f32 to f64
%24 = arith.subf %25, %23 : f64
llvm.store %24, %17 : f64, !llvm.ptr
cf.br ^bb5
^bb4:
cf.br ^bb5
^bb5:
%26 = llvm.load %17 : !llvm.ptr -> f64
%27 = arith.divf %26, %8 : f64
%28 = llvm.load %3 : !llvm.ptr -> f64
%29 = arith.addf %28, %27 : f64
llvm.store %29, %3 : f64, !llvm.ptr
%30 = llvm.mlir.addressof @EPSILON : !llvm.ptr
%31 = llvm.load %30 : !llvm.ptr -> f64
%32 = arith.cmpf olt, %27, %31 : f64
cf.cond_br %32, ^bb6, ^bb7
^bb6:
%33 = llvm.load %3 : !llvm.ptr -> f64
func.return %33 : f64
^bb7:
cf.br ^bb8
^bb8:
%34 = llvm.load %6 : !llvm.ptr -> i32
%35 = arith.constant 1 : i32
%36 = arith.addi %34, %35 : i32
llvm.store %36, %6 : i32, !llvm.ptr
cf.br ^bb0
^bb2:
%37 = llvm.load %3 : !llvm.ptr -> f64
func.return %37 : f64
}
func.func @find_intersection(%arg0: f64, %arg1: f64, %arg2: f64, %arg3: f64, %arg4: f64) -> f64 {
%38 = llvm.mlir.constant(1 : i64) : i64
%39 = llvm.alloca %38 x f64 : (i64) -> !llvm.ptr
llvm.store %arg3, %39 : f64, !llvm.ptr
%40 = llvm.mlir.constant(1 : i64) : i64
%41 = llvm.alloca %40 x f64 : (i64) -> !llvm.ptr
llvm.store %arg4, %41 : f64, !llvm.ptr
cf.br ^bb9
^bb9:
%42 = arith.constant 1 : i1
cf.cond_br %42, ^bb10, ^bb11
^bb10:
%44 = llvm.load %39 : !llvm.ptr -> f64
%43 = func.call @s_curve(%44) : (f64) -> f64
%45 = llvm.load %39 : !llvm.ptr -> f64
%46 = arith.subf %45, %arg0 : f64
%47 = arith.subf %43, %arg1 : f64
%48 = arith.mulf %46, %46 : f64
%49 = arith.mulf %47, %47 : f64
%50 = arith.addf %48, %49 : f64
%51 = math.sqrt %50 : f64
%52 = arith.subf %51, %arg2 : f64
%53 = math.absf %52 : f64
%54 = llvm.mlir.addressof @EPSILON : !llvm.ptr
%55 = llvm.load %54 : !llvm.ptr -> f64
%56 = arith.cmpf olt, %53, %55 : f64
cf.cond_br %56, ^bb12, ^bb13
^bb12:
%57 = llvm.load %39 : !llvm.ptr -> f64
func.return %57 : f64
^bb13:
cf.br ^bb14
^bb14:
%58 = arith.constant 0 : i1
%59 = llvm.mlir.constant(1 : i64) : i64
%60 = llvm.alloca %59 x i1 : (i64) -> !llvm.ptr
llvm.store %58, %60 : i1, !llvm.ptr
%61 = arith.cmpf olt, %51, %arg2 : f64
cf.cond_br %61, ^bb15, ^bb16
^bb15:
%62 = arith.constant 0.0 : f32
%64 = arith.extf %62 : f32 to f64
%63 = arith.cmpf ogt, %46, %64 : f64
%65 = scf.if %63 -> (i1) {
%66 = llvm.load %41 : !llvm.ptr -> f64
%67 = arith.constant 0.0 : f32
%69 = arith.extf %67 : f32 to f64
%68 = arith.cmpf olt, %66, %69 : f64
scf.yield %68 : i1
} else {
%70 = arith.constant false
scf.yield %70 : i1
}
cf.cond_br %65, ^bb18, ^bb19
^bb18:
%71 = arith.constant 1 : i1
llvm.store %71, %60 : i1, !llvm.ptr
cf.br ^bb20
^bb19:
cf.br ^bb20
^bb20:
%72 = arith.constant 0.0 : f32
%74 = arith.extf %72 : f32 to f64
%73 = arith.cmpf olt, %46, %74 : f64
%75 = scf.if %73 -> (i1) {
%76 = llvm.load %41 : !llvm.ptr -> f64
%77 = arith.constant 0.0 : f32
%79 = arith.extf %77 : f32 to f64
%78 = arith.cmpf ogt, %76, %79 : f64
scf.yield %78 : i1
} else {
%80 = arith.constant false
scf.yield %80 : i1
}
cf.cond_br %75, ^bb21, ^bb22
^bb21:
%81 = arith.constant 1 : i1
llvm.store %81, %60 : i1, !llvm.ptr
cf.br ^bb23
^bb22:
cf.br ^bb23
^bb23:
cf.br ^bb17
^bb16:
%82 = arith.constant 0.0 : f32
%84 = arith.extf %82 : f32 to f64
%83 = arith.cmpf ogt, %46, %84 : f64
%85 = scf.if %83 -> (i1) {
%86 = llvm.load %41 : !llvm.ptr -> f64
%87 = arith.constant 0.0 : f32
%89 = arith.extf %87 : f32 to f64
%88 = arith.cmpf ogt, %86, %89 : f64
scf.yield %88 : i1
} else {
%90 = arith.constant false
scf.yield %90 : i1
}
cf.cond_br %85, ^bb24, ^bb25
^bb24:
%91 = arith.constant 1 : i1
llvm.store %91, %60 : i1, !llvm.ptr
cf.br ^bb26
^bb25:
cf.br ^bb26
^bb26:
%92 = arith.constant 0.0 : f32
%94 = arith.extf %92 : f32 to f64
%93 = arith.cmpf olt, %46, %94 : f64
%95 = scf.if %93 -> (i1) {
%96 = llvm.load %41 : !llvm.ptr -> f64
%97 = arith.constant 0.0 : f32
%99 = arith.extf %97 : f32 to f64
%98 = arith.cmpf olt, %96, %99 : f64
scf.yield %98 : i1
} else {
%100 = arith.constant false
scf.yield %100 : i1
}
cf.cond_br %95, ^bb27, ^bb28
^bb27:
%101 = arith.constant 1 : i1
llvm.store %101, %60 : i1, !llvm.ptr
cf.br ^bb29
^bb28:
cf.br ^bb29
^bb29:
cf.br ^bb17
^bb17:
%102 = llvm.load %60 : !llvm.ptr -> i1
cf.cond_br %102, ^bb30, ^bb31
^bb30:
%103 = arith.constant 0.0 : f32
%104 = llvm.load %41 : !llvm.ptr -> f64
%105 = arith.constant 2.0 : f32
%107 = arith.extf %105 : f32 to f64
%106 = arith.divf %104, %107 : f64
%109 = arith.extf %103 : f32 to f64
%108 = arith.subf %109, %106 : f64
llvm.store %108, %41 : f64, !llvm.ptr
cf.br ^bb32
^bb31:
cf.br ^bb32
^bb32:
%110 = llvm.load %39 : !llvm.ptr -> f64
%111 = llvm.load %41 : !llvm.ptr -> f64
%112 = arith.addf %110, %111 : f64
llvm.store %112, %39 : f64, !llvm.ptr
cf.br ^bb9
^bb11:
%113 = llvm.load %39 : !llvm.ptr -> f64
func.return %113 : f64
}
func.func @integrate(%arg0: f64, %arg1: f64, %arg2: f64, %arg3: f64, %arg4: f64, %arg5: f64) -> f64 {
%114 = arith.constant 0.0 : f32
%115 = arith.extf %114 : f32 to f64
%116 = llvm.mlir.constant(1 : i64) : i64
%117 = llvm.alloca %116 x f64 : (i64) -> !llvm.ptr
llvm.store %115, %117 : f64, !llvm.ptr
%118 = llvm.mlir.constant(1 : i64) : i64
%119 = llvm.alloca %118 x f64 : (i64) -> !llvm.ptr
llvm.store %arg3, %119 : f64, !llvm.ptr
cf.br ^bb33
^bb33:
%120 = llvm.load %119 : !llvm.ptr -> f64
%121 = arith.cmpf ole, %120, %arg4 : f64
cf.cond_br %121, ^bb34, ^bb35
^bb34:
%123 = llvm.load %119 : !llvm.ptr -> f64
%122 = func.call @s_curve(%123) : (f64) -> f64
%124 = arith.mulf %arg2, %arg2 : f64
%125 = llvm.load %119 : !llvm.ptr -> f64
%126 = arith.subf %125, %arg0 : f64
%127 = llvm.load %119 : !llvm.ptr -> f64
%128 = arith.subf %127, %arg0 : f64
%129 = arith.mulf %126, %128 : f64
%130 = arith.subf %124, %129 : f64
%131 = math.sqrt %130 : f64
%132 = arith.subf %arg1, %131 : f64
%133 = llvm.load %117 : !llvm.ptr -> f64
%134 = arith.subf %122, %132 : f64
%135 = arith.mulf %134, %arg5 : f64
%136 = arith.addf %133, %135 : f64
llvm.store %136, %117 : f64, !llvm.ptr
%137 = llvm.load %119 : !llvm.ptr -> f64
%138 = arith.addf %137, %arg5 : f64
llvm.store %138, %119 : f64, !llvm.ptr
cf.br ^bb33
^bb35:
%139 = llvm.load %117 : !llvm.ptr -> f64
func.return %139 : f64
}
func.func @main() -> i32 {
%140 = arith.constant 0.25 : f32
%141 = arith.extf %140 : f32 to f64
%142 = arith.constant 0.5 : f32
%143 = arith.extf %142 : f32 to f64
%144 = arith.constant 0.25 : f32
%145 = arith.extf %144 : f32 to f64
%147 = arith.constant 0.0 : f32
%148 = arith.constant 0.1 : f32
%149 = arith.subf %147, %148 : f32
%150 = arith.extf %149 : f32 to f64
%146 = func.call @find_intersection(%141, %143, %145, %141, %150) : (f64, f64, f64, f64, f64) -> f64
%152 = arith.constant 0.1 : f32
%153 = arith.extf %152 : f32 to f64
%151 = func.call @find_intersection(%141, %143, %145, %141, %153) : (f64, f64, f64, f64, f64) -> f64
%155 = arith.constant 0.00001 : f32
%156 = arith.extf %155 : f32 to f64
%154 = func.call @integrate(%141, %143, %145, %146, %151, %156) : (f64, f64, f64, f64, f64, f64) -> f64
%157 = llvm.mlir.addressof @str_0 : !llvm.ptr
%158 = llvm.call @printf(%157, %154) vararg(!llvm.func<i32 (ptr, ...)>) : (!llvm.ptr, f64) -> i32
%159 = arith.constant 0 : i32
func.return %159 : i32
}
}