summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/poly.ml186
-rw-r--r--src/poly.mli82
-rw-r--r--src/polytest.ml51
-rw-r--r--src/sequence.ml35
-rw-r--r--src/sequence.mli13
-rw-r--r--src/util.ml57
-rw-r--r--src/util.mli27
-rw-r--r--src/visualwall.ml203
-rw-r--r--src/wall.ml337
-rw-r--r--src/wall.mli52
-rw-r--r--src/wallgen.ml78
-rw-r--r--src/wallsolve.ml162
12 files changed, 1283 insertions, 0 deletions
diff --git a/src/poly.ml b/src/poly.ml
new file mode 100644
index 0000000..ed32263
--- /dev/null
+++ b/src/poly.ml
@@ -0,0 +1,186 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+module PolyMap = Map.Make(Z)
+
+type t = Z.t PolyMap.t
+
+
+
+let zero =
+ PolyMap.empty
+
+let one =
+ PolyMap.singleton Z.zero Z.one
+
+let singleton ~coeff ~expnt =
+ if coeff = Z.zero
+ then zero
+ else PolyMap.singleton expnt coeff
+
+let ( ~$ ) (a, b) =
+ singleton ~coeff:a ~expnt:b
+
+
+
+let of_int a =
+ PolyMap.singleton Z.zero (Z.of_int a)
+
+let of_zint a =
+ PolyMap.singleton Z.zero a
+
+
+
+let of_arraylist_helper p (coeff,expnt) =
+ if coeff = Z.zero then p
+ else if PolyMap.mem expnt p
+ then PolyMap.add expnt Z.((PolyMap.find expnt p) + coeff) p
+ else PolyMap.add expnt coeff p
+
+let to_list a =
+ List.rev_map (fun (expnt,coeff) -> (coeff,expnt)) (PolyMap.bindings a)
+
+let of_list a =
+ List.fold_left of_arraylist_helper zero a
+
+let to_array a =
+ Array.of_list (to_list a)
+
+let of_array a =
+ Array.fold_left of_arraylist_helper zero a
+
+
+
+let degree p =
+ match (PolyMap.max_binding_opt p) with
+ None -> Z.zero
+ | Some (k,v) -> k
+
+let get_coeff ~expnt p =
+ match (PolyMap.find_opt expnt p) with
+ None -> Z.zero
+ | Some v -> v
+
+
+
+let stringifier (coeff, expnt) =
+ let sign_str = if coeff < Z.zero then "-" else "+" in
+ let coeff_str =
+ if Z.abs coeff = Z.one && expnt <> Z.zero
+ then ""
+ else Z.to_string (Z.abs coeff) in
+ let expnt_str =
+ if expnt = Z.zero
+ then ""
+ else if expnt = Z.one
+ then "x"
+ else "x^" ^ (Z.to_string expnt) in
+ (sign_str, coeff_str ^ expnt_str)
+
+let to_string a =
+ let pieces = List.map stringifier (to_list a)
+ in if List.length pieces = 0
+ then "0"
+ else let first_term_sign = if fst (List.hd pieces) = "-" then "-" else "" in
+ let first_term = first_term_sign ^ (snd (List.hd pieces)) in
+ let term_concat s (x,y) = s ^ " " ^ x ^ " " ^ y in
+ let remaining_terms = List.fold_left term_concat "" (List.tl pieces) in
+ first_term ^ remaining_terms
+
+let print a =
+ print_string (to_string a)
+
+let output channel a =
+ output_string channel (to_string a)
+
+let pp_print format a =
+ let pieces = List.map stringifier (to_list a) in
+ let do_first item =
+ if fst item = "-" then Format.pp_print_string format "-";
+ Format.pp_print_string format (snd item) in
+ let do_rest item =
+ Format.pp_print_string format (" " ^ (fst item));
+ Format.pp_print_space format ();
+ Format.pp_print_string format (snd item) in
+ if List.length pieces = 0
+ then Format.pp_print_string format "0"
+ else (do_first (List.hd pieces); List.iter do_rest (List.tl pieces))
+
+
+
+let neg =
+ PolyMap.map Z.neg
+
+let add a b =
+ let f key x y =
+ if Z.(x + y) = Z.zero
+ then None
+ else Some Z.(x + y)
+ in PolyMap.union f a b
+
+let sub a b =
+ let f key x y =
+ match (x, y) with
+ (None, None) -> None
+ | (Some v, None) -> Some v
+ | (None, Some v) -> Some Z.(~- v)
+ | (Some v1, Some v2) ->
+ if Z.(v1 - v2) = Z.zero
+ then None
+ else Some Z.(v1 - v2)
+ in PolyMap.merge f a b
+
+let single_mul expnt coeff p =
+ let helper expnt2 coeff2 r =
+ PolyMap.add Z.(expnt + expnt2) Z.(coeff * coeff2) r
+ in PolyMap.fold helper p zero
+
+let mul a b =
+ PolyMap.fold (fun k v result -> add result (single_mul k v b)) a zero
+
+let div_rem a b =
+ let rec helper dividend quotient remainder =
+ if dividend = zero
+ then (quotient, remainder)
+ else if degree b > degree dividend
+ then (quotient, add dividend remainder)
+ else let (dsor_expnt, dsor_coeff) = PolyMap.max_binding b in
+ let (dend_expnt, dend_coeff) = PolyMap.max_binding dividend in
+ let (coeff_div, coeff_rem) = Z.div_rem dend_coeff dsor_coeff in
+ let expnt_diff = Z.(dend_expnt - dsor_expnt) in
+ let divisor_adjusted = single_mul expnt_diff coeff_div b in
+ let dividend_adjusted = sub dividend divisor_adjusted in
+ let new_dividend = PolyMap.remove dend_expnt dividend_adjusted in
+ let new_quotient =
+ if coeff_div <> Z.zero
+ then PolyMap.add expnt_diff coeff_div quotient
+ else quotient in
+ let new_remainder =
+ if coeff_rem <> Z.zero
+ then PolyMap.add dend_expnt coeff_rem remainder
+ else remainder in
+ helper new_dividend new_quotient new_remainder
+ in if b = zero
+ then raise Division_by_zero
+ else helper a zero zero
+
+let div a b =
+ fst (div_rem a b)
+
+let rem a b =
+ snd (div_rem a b)
+
+
+
+let ( + ) = add
+
+let ( - ) = sub
+
+let ( * ) = mul
+
+let ( / ) = div
+
+
diff --git a/src/poly.mli b/src/poly.mli
new file mode 100644
index 0000000..4837142
--- /dev/null
+++ b/src/poly.mli
@@ -0,0 +1,82 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+type t
+
+
+
+val zero : t
+
+val one : t
+
+val singleton : coeff:Z.t -> expnt:Z.t -> t
+
+val ( ~$ ) : Z.t * Z.t -> t
+
+
+
+(* These two functions produce polynomials of degree zero. *)
+
+val of_int : int -> t
+
+val of_zint : Z.t -> t
+
+
+
+(* The first values are the coefficients, the second values are the exponents. *)
+(* Values given by to_* functions are in descending exponent order. *)
+
+val to_list : t -> (Z.t * Z.t) list
+
+val of_list : (Z.t * Z.t) list -> t
+
+val to_array : t -> (Z.t * Z.t) array
+
+val of_array : (Z.t * Z.t) array -> t
+
+
+
+val degree : t -> Z.t
+
+val get_coeff : expnt:Z.t -> t -> Z.t
+
+
+
+val to_string : t -> string
+
+val print : t -> unit
+
+val output : Stdlib.out_channel -> t -> unit
+
+val pp_print : Stdlib.Format.formatter -> t -> unit
+
+
+
+val neg : t -> t
+
+val add : t -> t -> t
+
+val sub : t -> t -> t
+
+val mul : t -> t -> t
+
+val div_rem : t -> t -> t * t
+
+val div : t -> t -> t
+
+val rem : t -> t -> t
+
+
+
+val ( + ) : t -> t -> t
+
+val ( - ) : t -> t -> t
+
+val ( * ) : t -> t -> t
+
+val ( / ) : t -> t -> t
+
+
diff --git a/src/polytest.ml b/src/polytest.ml
new file mode 100644
index 0000000..c378279
--- /dev/null
+++ b/src/polytest.ml
@@ -0,0 +1,51 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+let _ =
+ let p1 = Poly.of_list Z.( [(~$1,~$2);(~$7,~$1);(~$3,~$0)] ) in
+ print_string "dividend: "; Poly.print p1; print_newline ();
+
+ let p2 = Poly.of_list Z.( [(~$1,~$1);(~$(-1),~$0)] ) in
+ print_string "divisor: "; Poly.print p2; print_newline ();
+
+ let (p3,p4) = Poly.div_rem p1 p2 in
+ print_string "quotient: "; Poly.print p3; print_newline ();
+ print_string "remainder: "; Poly.print p4; print_newline ();
+
+ let p5 = Poly.add (Poly.mul p3 p2) p4 in
+ print_string "reconstituted: "; Poly.print p5; print_newline ();
+
+ print_newline ();
+
+ let p6 = Poly.of_list Z.( [(~$(-3),~$2);(~$7,~$0)] ) in
+ print_string "dividend: "; Poly.print p6; print_newline ();
+
+ let p7 = Poly.of_list Z.( [(~$2,~$1);(~$2,~$0)] ) in
+ print_string "divisor: "; Poly.print p7; print_newline ();
+
+ let (p8,p9) = Poly.div_rem p6 p7 in
+ print_string "quotient: "; Poly.print p8; print_newline ();
+ print_string "remainder: "; Poly.print p9; print_newline ();
+
+ let p10 = Poly.add (Poly.mul p8 p7) p9 in
+ print_string "reconstituted: "; Poly.print p10; print_newline ();
+
+ print_newline ();
+
+ let p11 = Poly.of_list Z.( [(~$1,~$8);(~$(-1),~$0)] ) in
+ print_string "dividend: "; Poly.print p11; print_newline ();
+
+ let p12 = Poly.of_list Z.( [(~$1,~$1);(~$(-1),~$0)] ) in
+ print_string "divisor: "; Poly.print p12; print_newline ();
+
+ let (p13,p14) = Poly.div_rem p11 p12 in
+ print_string "quotient: "; Poly.print p13; print_newline ();
+ print_string "remainder: "; Poly.print p14; print_newline ();
+
+ let p15 = Poly.add (Poly.mul p13 p12) p14 in
+ print_string "reconstituted: "; Poly.print p15; print_newline ();
+
+
diff --git a/src/sequence.ml b/src/sequence.ml
new file mode 100644
index 0000000..fba7f96
--- /dev/null
+++ b/src/sequence.ml
@@ -0,0 +1,35 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+let fibonacci ~len =
+ let rec fib a b r =
+ if List.length (r) > len - 1 then r
+ else fib b Z.(a + b) (a :: r)
+ in Array.of_list (List.rev (fib Z.zero Z.one []))
+
+
+
+let pagoda ~len =
+ let open Z in
+ let rec b n =
+ if n = zero then zero
+ else (uneven n) / ~$2 mod ~$2
+ and a n =
+ if ~$n = zero then minus_one
+ else b (~$n + one) - b (~$n - one)
+ and uneven n =
+ if is_even n then uneven (n / ~$2) else n
+ in Array.init len a
+
+
+
+let random ~len =
+ Random.self_init ();
+ let f n =
+ Z.of_int ((Random.int 1000000000) - 500000000)
+ in Array.init len f
+
+
diff --git a/src/sequence.mli b/src/sequence.mli
new file mode 100644
index 0000000..a2fefbc
--- /dev/null
+++ b/src/sequence.mli
@@ -0,0 +1,13 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+val fibonacci : len:int -> Z.t array
+
+val pagoda : len:int -> Z.t array
+
+val random : len:int -> Z.t array
+
+
diff --git a/src/util.ml b/src/util.ml
new file mode 100644
index 0000000..93f2c93
--- /dev/null
+++ b/src/util.ml
@@ -0,0 +1,57 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+exception Not_an_integer of string
+
+
+
+let ( @@ ) f g x = f (g x)
+
+let ( |> ) f x = f x
+
+
+
+let call_with_open_input_file ~filename ~func =
+ let ic = open_in filename in
+ let res =
+ try func ic
+ with exn -> close_in ic; raise exn in
+ close_in ic;
+ res
+
+
+
+let call_with_open_output_file ~filename ~func =
+ let oc = open_out filename in
+ let res =
+ try func oc
+ with exn -> close_out oc; raise exn in
+ close_out oc;
+ res
+
+
+
+let read_integer_list args =
+ let convert x =
+ try Z.of_string x
+ with Invalid_argument _ ->
+ raise (Not_an_integer x) in
+ let blank_is_zero x =
+ if x = "" then "0" else x in
+ (Array.map convert) @@
+ Array.of_list @@
+ (List.map (blank_is_zero @@ String.trim)) @@
+ (String.split_on_char ',') |> args
+
+
+
+let read_sequence_file filename =
+ let reader handle =
+ let line = input_line handle in
+ read_integer_list line
+ in call_with_open_input_file ~filename:filename ~func:reader
+
+
diff --git a/src/util.mli b/src/util.mli
new file mode 100644
index 0000000..b1821a4
--- /dev/null
+++ b/src/util.mli
@@ -0,0 +1,27 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+exception Not_an_integer of string
+
+
+
+val ( @@ ) : ('a -> 'b) -> ('c -> 'a) -> 'c -> 'b
+
+val ( |> ) : ('a -> 'b) -> 'a -> 'b
+
+
+
+val call_with_open_input_file : filename:string -> func:(in_channel -> 'a) -> 'a
+
+val call_with_open_output_file : filename:string -> func:(out_channel -> 'a) -> 'a
+
+
+
+val read_integer_list : string -> Z.t array
+
+val read_sequence_file : string -> Z.t array
+
+
diff --git a/src/visualwall.ml b/src/visualwall.ml
new file mode 100644
index 0000000..3c07363
--- /dev/null
+++ b/src/visualwall.ml
@@ -0,0 +1,203 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+let main () =
+ ignore (GtkMain.Main.init ());
+
+ let window = GWindow.window
+ ~title:"Number Wall Viewer"
+ ~width:250 () in
+ ignore (window#connect#destroy ~callback:GMain.quit);
+
+ let main_vbox = GPack.vbox
+ ~packing:window#add () in
+
+ let menubar = GMenu.menu_bar
+ ~packing:(main_vbox#pack ~expand:false) () in
+ let wall_view = GBin.scrolled_window
+ ~hpolicy:`AUTOMATIC
+ ~vpolicy:`AUTOMATIC
+ ~packing:(main_vbox#pack ~expand:true) () in
+ let statusbar = GMisc.statusbar
+ ~packing:(main_vbox#pack ~expand:false) () in
+
+
+
+ let context = statusbar#new_context ~name:"Program status" in
+
+ let set_status ~msg =
+ context#pop ();
+ ignore (context#push msg) in
+
+
+
+ let choose_open_file ~title ?(folder = "") () =
+ let chooser = GWindow.file_chooser_dialog
+ ~action:`OPEN
+ ~title:title () in
+ chooser#add_button_stock `CANCEL `dialog;
+ chooser#add_select_button_stock `OPEN `filechooser;
+ ignore (chooser#set_current_folder folder);
+ let code = chooser#run () in
+ let maybe_result = chooser#filename in
+ chooser#destroy ();
+ match maybe_result with
+ | None -> ""
+ | Some result -> if code = `filechooser then result else "" in
+
+ let choose_save_file ~title ?(folder = "") ?(name = "") () =
+ let chooser = GWindow.file_chooser_dialog
+ ~action:`SAVE
+ ~title:title () in
+ chooser#add_button_stock `CANCEL `dialog;
+ chooser#add_select_button_stock `SAVE `filechooser;
+ chooser#set_current_name name;
+ ignore (chooser#set_current_folder folder);
+ let code = chooser#run () in
+ let maybe_result = chooser#filename in
+ chooser#destroy ();
+ match maybe_result with
+ | None -> ""
+ | Some result -> if code = `filechooser then result else "" in
+
+
+
+ let generate_unused ~basename ~extension =
+ let rec checker num =
+ let testname = basename ^ "-" ^ (string_of_int num) ^ "." ^ extension in
+ if not (Sys.file_exists testname) then testname else checker (num + 1) in
+ checker 1 in
+
+
+
+ let overwrite_check ~filename =
+ if not (Sys.file_exists filename) then true else
+ (let asker = GWindow.dialog
+ ~title:"Confirm Overwrite"
+ ~resizable:false
+ ~width:380 ~height:130 () in
+ let box = GPack.hbox
+ ~packing:asker#vbox#add () in
+ let _ = GMisc.image
+ ~stock:`DIALOG_ERROR
+ ~packing:box#add () in
+ let message = GMisc.label
+ ~text:("Error: " ^ filename ^ " already exists.")
+ ~packing:box#add () in
+ message#set_line_wrap true;
+ asker#add_button "Cancel" `dialog;
+ asker#add_button "Overwrite" `filechooser;
+ let code = asker#run () in
+ asker#destroy ();
+ code = `filechooser) in
+
+
+
+ let new_cb () =
+ print_endline "New..." in
+
+
+
+ let open_cb () =
+ let filename = choose_open_file
+ ~title:"Open Sequence - Number Wall Viewer"
+ ~folder:(Sys.getcwd ()) () in
+ if filename != "" then
+ (set_status ~msg:("Opened " ^ filename);
+ print_endline filename)
+ else
+ (set_status ~msg:"Canceled open") in
+
+ let save_as_cb () =
+ let filename = choose_save_file
+ ~title:"Save Sequence - Number Wall Viewer"
+ ~folder:(Sys.getcwd ())
+ ~name:(generate_unused ~basename:"sequence" ~extension:"txt") () in
+ if filename != "" && overwrite_check ~filename then
+ (set_status ~msg:("Saved sequence to " ^ filename);
+ print_endline filename)
+ else
+ (set_status ~msg:"Canceled save") in
+
+ let export_cb () =
+ let filename = choose_save_file
+ ~title:"Export Image - Number Wall Viewer"
+ ~folder:(Sys.getcwd ())
+ ~name:(generate_unused ~basename:"image" ~extension:"png") () in
+ if filename != "" && overwrite_check ~filename then
+ (set_status ~msg:("Exported image to " ^ filename);
+ print_endline filename)
+ else
+ (set_status ~msg:"Canceled export") in
+
+
+
+ let discard_cb () =
+ print_endline "Discard" in
+
+
+
+ let color_cb () =
+ print_endline "Color..." in
+
+ let modulus_cb () =
+ print_endline "Modulus..." in
+
+
+
+ let about_cb () =
+ let about_window = GWindow.about_dialog
+ ~authors:["Jedidiah Barber"]
+ ~comments:"A simple program to generate and visualise number walls."
+ ~license:"Sunset License v1.0"
+ ~title:"About - Number Wall Viewer"
+ ~modal:true () in
+ let destroy_me button = about_window#destroy () in
+ ignore (about_window#connect#response ~callback:destroy_me);
+ about_window#show () in
+
+
+
+ let file_menu_entries = [
+ `I ("New...", new_cb);
+ `S;
+ `I ("Open...", open_cb);
+ `I ("Save As...", save_as_cb);
+ `I ("Export...", export_cb);
+ `S;
+ `I ("Discard", discard_cb);
+ `I ("Quit", GMain.quit)
+ ] in
+
+ let option_menu_entries = [
+ `I ("Color...", color_cb);
+ `I ("Modulus...", modulus_cb)
+ ] in
+
+ let help_menu_entries = [
+ `I ("About", about_cb)
+ ] in
+
+
+
+ let create_menu ~label ~bar ~entries =
+ let item = GMenu.menu_item ~label ~packing:bar#append () in
+ let menu = GMenu.menu ~packing:item#set_submenu () in
+ GToolbox.build_menu menu ~entries in
+
+ ignore (create_menu ~label:"File" ~bar:menubar ~entries:file_menu_entries);
+ ignore (create_menu ~label:"Options" ~bar:menubar ~entries:option_menu_entries);
+ ignore (create_menu ~label:"Help" ~bar:menubar ~entries:help_menu_entries);
+
+ set_status ~msg:"No wall loaded";
+ window#show ();
+ GMain.Main.main ()
+
+
+
+let _ = main ()
+
+
diff --git a/src/wall.ml b/src/wall.ml
new file mode 100644
index 0000000..94043dc
--- /dev/null
+++ b/src/wall.ml
@@ -0,0 +1,337 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+module type Algebra =
+ sig
+ type t
+ val zero : t
+ val one : t
+ val add : t -> t -> t
+ val sub : t -> t -> t
+ val mul : t -> t -> t
+ val div : t -> t -> t
+ val to_string : t -> string
+ end
+
+
+
+module type S =
+ sig
+ type element
+ type wall
+ val create : dimx:int -> dimy:int -> init:element array -> wall
+ val x_length : wall -> int
+ val y_length : wall -> int
+ val has_value : wall -> x:int -> y:int -> bool
+ val get : wall -> x:int -> y:int -> element
+ val get_opt : wall -> x:int -> y:int -> element option
+ val last_nonzero_row : wall -> int
+ val to_string : wall -> string
+ val print : wall -> unit
+ val output : Stdlib.out_channel -> wall -> unit
+ val pp_print : Stdlib.Format.formatter -> wall -> unit
+ end
+
+
+
+let (let*?) = Option.bind
+
+
+
+module Make (A : Algebra) =
+ struct
+
+ type element = A.t
+ type wall = element option array array
+
+ let x_length wall =
+ Array.length wall
+
+ let y_length wall =
+ Array.length wall.(0)
+
+ let has_value wall ~x ~y =
+ if x >= 0 && x < x_length wall &&
+ y >= 0 && y < y_length wall
+ then match wall.(x).(y) with
+ None -> false
+ | Some _ -> true
+ else false
+
+ let get_opt wall ~x ~y =
+ if x >= 0 && x < x_length wall &&
+ y >= 0 && y < y_length wall
+ then wall.(x).(y)
+ else None
+
+
+
+ let next_nonzero_up wall ~x ~y =
+ let rec walker y1 =
+ let*? v = get_opt wall ~x ~y:y1 in
+ if v <> A.zero then Some y1 else walker (y1-1)
+ in walker y
+
+ let next_nonzero_left wall ~x ~y =
+ let rec walker x1 =
+ let*? v = get_opt wall ~x:x1 ~y in
+ if v <> A.zero then Some x1 else walker (x1-1)
+ in walker x
+
+ let next_nonzero_right wall ~x ~y =
+ let rec walker x1 =
+ let*? v = get_opt wall ~x:x1 ~y in
+ if v <> A.zero then Some x1 else walker (x1+1)
+ in walker x
+
+
+
+ let next_valid_left wall ~x ~y =
+ let rec walker x1 =
+ if x1 < 0
+ then None
+ else if Option.is_some (get_opt wall ~x:x1 ~y)
+ then Some x1
+ else walker (x1-1) in
+ walker x
+
+
+
+ let safe_div a b =
+ if b = A.zero then None else Some (A.div a b)
+
+ let multi_mul items =
+ let f a b = A.mul a b in
+ List.fold_left f A.one items
+
+ let square a =
+ A.mul a a
+
+ let rec power ~base ~num ~denom ~times =
+ if times <= 0
+ then Some base
+ else let*? new_base = safe_div (A.mul base num) denom in
+ power ~base:new_base ~num ~denom ~times:(times-1)
+
+
+
+ (* red x red + green x green = blue squared *)
+ (* target to calculate is red on the bottom of the cross *)
+ (* current cell x two up + one diag up left x one diag up right = one up squared *)
+ let cross_rule wall ~x ~y =
+ let*? green_left = get_opt wall ~x:(x-1) ~y:(y-1) in
+ let*? green_right = get_opt wall ~x:(x+1) ~y:(y-1) in
+ let*? blue = get_opt wall ~x ~y:(y-1) in
+ let*? red_top = get_opt wall ~x ~y:(y-2) in
+ safe_div (A.sub (square blue) (A.mul green_left green_right)) red_top
+
+
+
+ (* outer yellow x inner yellow squared + outer blue x inner blue squared *)
+ (* equals outer grey x inner grey squared + outer red x inner red squared *)
+ (* inner yellow is up, outer yellow is down *)
+ (* inner blue is down, outer blue is up *)
+ (* inner grey is left, outer grey is right *)
+ (* inner red is right, outer red is left *)
+ (* target to calculate is outer yellow *)
+ let long_cross_rule wall ~x ~y =
+ let*? outer_red = get_opt wall ~x:(x-2) ~y:(y-2) in
+ let*? inner_red = get_opt wall ~x:(x+1) ~y:(y-2) in
+ let*? red_product = Some (A.mul outer_red (square inner_red)) in
+ let*? outer_grey = get_opt wall ~x:(x+2) ~y:(y-2) in
+ let*? inner_grey = get_opt wall ~x:(x-1) ~y:(y-2) in
+ let*? grey_product = Some (A.mul outer_grey (square inner_grey)) in
+ let*? grey_red_sum = Some (A.add grey_product red_product) in
+ let*? outer_blue = if y - 4 = -1 then Some A.zero else get_opt wall ~x ~y:(y-4) in
+ let*? inner_blue = get_opt wall ~x ~y:(y-1) in
+ let*? blue_product = Some (A.mul outer_blue (square inner_blue)) in
+ let*? greyred_blue_diff = Some (A.sub grey_red_sum blue_product) in
+ let*? inner_yellow = get_opt wall ~x ~y:(y-3) in
+ safe_div greyred_blue_diff (square inner_yellow)
+
+
+
+ (* All blocks of zeros in the wall come in squares. *)
+ let zero_window wall ~x ~y =
+ let*? top = next_nonzero_up wall ~x ~y:(y-1) in
+ let*? left = next_nonzero_left wall ~x ~y:(top+1) in
+ let*? right = next_nonzero_right wall ~x ~y:(top+1) in
+ if y - 1 > top && right - left > 0 && right - left + 1 > y - 1 - top
+ then Some A.zero
+ else None
+
+
+
+ (* denote ratio of top sequence from left to right as a/b *)
+ (* denote ratio of bottom sequence from right to left as c/d *)
+ (* denote ratio of left sequence from top to bottom as e/f *)
+ (* denote ratio of right sequence from bottom to top as g/h *)
+ (* denote rhs as s *)
+ (* s is 1 when even sized window and -1 when odd sized window *)
+ (* goal is to calculate c/d *)
+ (* then *)
+ (* ((a/b)(c/d))/((e/f)(g/h)) = s *)
+ (* (a/b)(c/d) = s(e/f)(g/h) *)
+ (* (a/b)(c/d) = (seg)/(fh) *)
+ (* c/d = ((seg)/(fh))(b/a) *)
+ (* c/d = (begs)/(afh) *)
+ (* and since that ratio is going from right to left, invert to go left to right *)
+ (* thus to transform bottom left corner number into next number... *)
+ (* multiply by afh then divide by begs *)
+ let horseshoe_rule wall ~x ~y =
+ let*? top = next_nonzero_up wall ~x ~y:(y-1) in
+ let*? left = next_nonzero_left wall ~x ~y:(top+1) in
+ let*? right = next_nonzero_right wall ~x ~y:(top+1) in
+ let*? _ = if y - top <> right - left || right - left <= 2 then None else Some 1 in
+ let*? a = get_opt wall ~x:(left+1) ~y:top in
+ let*? b = get_opt wall ~x:left ~y:top in
+ let*? e = get_opt wall ~x:left ~y:(top+1) in
+ let*? f = get_opt wall ~x:left ~y:top in
+ let*? g = get_opt wall ~x:right ~y:top in
+ let*? h = get_opt wall ~x:right ~y:(top+1) in
+ let*? s = if (right - left - 1) mod 2 = 0 then Some A.one else Some (A.sub A.zero A.one) in
+ let*? prev = next_valid_left wall ~x ~y in
+ power
+ ~base:(Option.get wall.(prev).(y))
+ ~num:(A.mul (A.mul a f) h)
+ ~denom:(A.mul (A.mul (A.mul b e) g) s)
+ ~times:(x - prev)
+
+
+
+ (* denote ratio of top sequence from left to right as a/b *)
+ (* denote ratio of bottom sequence from right to left as c/d *)
+ (* denote ratio of left sequence from top to bottom as e/f *)
+ (* denote ratio of right sequence from bottom to top as g/h *)
+ (* denote sign of grey/red terms as s *)
+ (* goal is to calculate outer yellow *)
+ (* and also to put off division until the very end to ensure integer intermediate results *)
+ (* then *)
+ (* (e/f)(o_bl/i_bl) + s(a/b)(o_gr/i_gr) = (g/h)(o_ye/i_ye) + s(c/d)(o_re/i_re) *)
+ (* (e/f)(o_bl/i_bl) + s(a/b)(o_gr/i_gr) - s(c/d)(o_re/i_re) = (g/h)(o_ye/i_ye) *)
+ (* i_ye((e/f)(o_bl/i_bl) + s(a/b)(o_gr/i_gr) - s(c/d)(o_re/i_re)) = (g/h)o_ye *)
+ (* o_ye = i_ye(h/g)((e/f)(o_bl/i_bl) + s(a/b)(o_gr/i_gr) - s(c/d)(o_re/i_re)) *)
+ (* o_ye = i_ye(h/g)((e*o_bl)/(f*i_bl) + s(a*o_gr)/(b*i_gr) - s(c*o_re)/(d*i_re)) *)
+ (* o_ye = h*i_ye*(e*o_bl*b*i_gr*d*i_re + s*a*o_gr*f*i_bl*d*i_re - s*c*o_re*f*i_bl*b*i_gr) / *)
+ (* (g*f*i_bl*b*i_gr*d*i_re) *)
+ let broken_cross_rule wall ~x ~y =
+ let*? top = next_nonzero_up wall ~x ~y:(y-2) in
+ let*? left = next_nonzero_left wall ~x ~y:(top+1) in
+ let*? right = next_nonzero_right wall ~x ~y:(top+1) in
+ let*? bottom = Some (y - 1) in
+ let*? _ = if bottom - top <> right - left || right - left <= 2 then None else Some 1 in
+ let*? a = get_opt wall ~x:(left+1) ~y:top in
+ let*? b = get_opt wall ~x:left ~y:top in
+ let*? c = get_opt wall ~x:left ~y:bottom in
+ let*? d = get_opt wall ~x:(left+1) ~y:bottom in
+ let*? e = get_opt wall ~x:left ~y:(top+1) in
+ let*? f = get_opt wall ~x:left ~y:top in
+ let*? g = get_opt wall ~x:right ~y:top in
+ let*? h = get_opt wall ~x:right ~y:(top+1) in
+ let*? s = if (right - x) mod 2 = 0 then Some A.one else Some (A.sub A.zero A.one) in
+ let*? inner_yellow = get_opt wall ~x ~y:bottom in
+ let*? inner_blue = get_opt wall ~x:(right-x+left) ~y:top in
+ let*? inner_red = get_opt wall ~x:right ~y:(bottom-right-x) in
+ let*? inner_grey = get_opt wall ~x:left ~y:(right-x+top) in
+ let*? outer_blue = get_opt wall ~x:(right-x+left) ~y:(top-1) in
+ let*? outer_red = get_opt wall ~x:(right+1) ~y:(bottom-right-x) in
+ let*? outer_grey = get_opt wall ~x:(left-1) ~y:(right-x+top) in
+ let*? term_1 = Some (multi_mul [ e; outer_blue; b; inner_grey; d; inner_red]) in
+ let*? term_2 = Some (multi_mul [s; a; outer_grey; f; inner_blue; d; inner_red]) in
+ let*? term_3 = Some (multi_mul [s; c; outer_red; f; inner_blue; b; inner_grey]) in
+ let*? term_sum = Some (A.sub (A.add term_1 term_2) term_3) in
+ let*? numerator = Some (multi_mul [h; inner_yellow; term_sum]) in
+ let*? denominator = Some (multi_mul [g; f; inner_blue; b; inner_grey; d; inner_red]) in
+ safe_div numerator denominator
+
+
+
+ let rec attempt_order wall ~x ~y funcs =
+ if funcs = []
+ then None
+ else let r = (List.hd funcs) wall ~x ~y in
+ if Option.is_some r
+ then r
+ else attempt_order wall ~x ~y (List.tl funcs)
+
+ let create ~dimx ~dimy ~init =
+ let result = Array.make_matrix dimx dimy None in
+ for i = 0 to dimx - 1 do
+ result.(i).(0) <- Some A.one
+ done;
+ for i = 0 to min (dimx - 1) ((Array.length init) - 1) do
+ result.(i).(1) <- Some init.(i)
+ done;
+ try
+ for j = 2 to dimy - 1 do
+ let do_exit = ref true in
+ for i = 0 to dimx - 1 do
+ let v = attempt_order result ~x:i ~y:j
+ [cross_rule; long_cross_rule; zero_window; horseshoe_rule; broken_cross_rule] in
+ result.(i).(j) <- v;
+ if Option.is_some v && Option.get v <> A.zero
+ then do_exit := false
+ done;
+ if !do_exit then raise Exit
+ done;
+ result
+ with Exit ->
+ result
+
+
+
+ let get wall ~x ~y =
+ if x >= 0 && x < x_length wall &&
+ y >= 0 && y < y_length wall
+ then match wall.(x).(y) with
+ None -> raise Not_found
+ | Some v -> v
+ else raise Not_found
+
+ let last_nonzero_row wall =
+ let rec checker x y =
+ if wall.(x).(y) <> Some A.zero && wall.(x).(y) <> None
+ then (if y + 1 = y_length wall then y else checker 0 (y + 1))
+ else (if x + 1 = x_length wall then y - 1 else checker (x + 1) y)
+ in checker 0 0
+
+
+
+ let row_string wall ~y =
+ let result = ref "" in
+ for i = 0 to (x_length wall) - 2 do
+ let v = wall.(i).(y) in
+ if Option.is_some v
+ then result := !result ^ (A.to_string (Option.get v)) ^ ","
+ else result := !result ^ ","
+ done;
+ let v = wall.(x_length wall - 1).(y) in
+ if Option.is_some v
+ then result := !result ^ (A.to_string (Option.get v));
+ !result
+
+ let to_string wall =
+ let result = ref "" in
+ for j = 0 to (y_length wall) - 1 do
+ result := !result ^ (row_string wall ~y:j) ^ (String.make 1 (char_of_int 10))
+ done;
+ !result
+
+ let print wall =
+ print_string (to_string wall)
+
+ let output channel wall =
+ output_string channel (to_string wall)
+
+ let pp_print format wall =
+ for j = 0 to (y_length wall) - 1 do
+ Format.pp_print_string format (row_string wall ~y:j);
+ Format.pp_print_break format 8 0
+ done
+
+ end
+
+
diff --git a/src/wall.mli b/src/wall.mli
new file mode 100644
index 0000000..59d8edb
--- /dev/null
+++ b/src/wall.mli
@@ -0,0 +1,52 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+module type Algebra =
+ sig
+ type t
+
+ val zero : t
+ val one : t
+
+ val add : t -> t -> t
+ val sub : t -> t -> t
+ val mul : t -> t -> t
+ val div : t -> t -> t
+
+ val to_string : t -> string
+ end
+
+
+
+module type S =
+ sig
+ type element
+ type wall
+
+ val create : dimx:int -> dimy:int -> init:element array -> wall
+
+ val x_length : wall -> int
+ val y_length : wall -> int
+
+ val has_value : wall -> x:int -> y:int -> bool
+
+ val get : wall -> x:int -> y:int -> element
+ val get_opt : wall -> x:int -> y:int -> element option
+
+ val last_nonzero_row : wall -> int
+
+ val to_string : wall -> string
+ val print : wall -> unit
+ val output : Stdlib.out_channel -> wall -> unit
+ val pp_print : Stdlib.Format.formatter -> wall -> unit
+ end
+
+
+
+module Make:
+ functor (A : Algebra) -> S with type element = A.t
+
+
diff --git a/src/wallgen.ml b/src/wallgen.ml
new file mode 100644
index 0000000..c27f0eb
--- /dev/null
+++ b/src/wallgen.ml
@@ -0,0 +1,78 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+let usage_msg = "Usage: wallgen [options]"
+let input_seq = ref (Array.make 0 (Z.of_int 0))
+let depth_limit = ref 20
+let use_random = ref false
+
+
+
+let anon_fun item =
+ prerr_endline "Error: Don't know what to do with anonymous argument";
+ exit 2
+
+
+
+let input_immediate items =
+ if Array.length !input_seq > 0
+ then (prerr_endline "Error: Multiple input sequences specified"; exit 2)
+ else try input_seq := Util.read_integer_list items
+ with Util.Not_an_integer s ->
+ prerr_endline ("Error: The argument '" ^ s ^ "' is not an integer");
+ exit 2
+
+let input_from_file filename =
+ if Array.length !input_seq > 0
+ then (prerr_endline "Error: Multiple input sequences specified"; exit 2)
+ else try input_seq := Util.read_sequence_file filename
+ with Util.Not_an_integer s ->
+ prerr_endline ("Error: The item '" ^ s ^ "' from file '" ^ filename ^ "' is not an integer");
+ exit 2
+
+let use_random len =
+ if Array.length !input_seq > 0
+ then (prerr_endline "Error: Multiple input sequences specified"; exit 2)
+ else input_seq := Sequence.random ~len
+
+
+
+let speclist =
+ [("-f", Arg.String input_from_file, "Read input integer sequence from file");
+ ("-i", Arg.String input_immediate, "Integer input sequence with values separated by commas");
+ ("-l", Arg.Set_int depth_limit, "Sets maximum levels to generate, default 20");
+ ("-r", Arg.Int use_random, "Use random input sequence of specified length")]
+
+
+
+module My_Algebra =
+ struct
+ type t = Z.t
+ let zero = Z.zero
+ let one = Z.one
+ let add = Z.add
+ let sub = Z.sub
+ let mul = Z.mul
+ let div = Z.div
+ let to_string = Z.to_string
+ end
+
+module My_Wall = Wall.Make (My_Algebra)
+
+
+
+let _ =
+ if Array.length Sys.argv <= 1 then (Arg.usage speclist usage_msg; exit 1);
+ Arg.parse speclist anon_fun usage_msg;
+ if Array.length !input_seq = 0 then (prerr_endline "Error: No input sequence provided"; exit 2);
+
+ let result = My_Wall.create
+ ~dimx:(Array.length !input_seq)
+ ~dimy:!depth_limit
+ ~init:!input_seq in
+ My_Wall.print result
+
+
diff --git a/src/wallsolve.ml b/src/wallsolve.ml
new file mode 100644
index 0000000..1e0e8e2
--- /dev/null
+++ b/src/wallsolve.ml
@@ -0,0 +1,162 @@
+
+
+(* Programmed by Jedidiah Barber *)
+(* Licensed under the Sunset License v1.0 *)
+
+
+let usage_msg = "Usage: wallsolve [options]"
+let input_seq = ref (Array.make 0 (Z.of_int 0))
+let depth_limit = ref 20
+
+
+
+let anon_fun item =
+ prerr_endline "Error: Don't know what to do with anonymous argument";
+ exit 2
+
+
+
+let input_immediate items =
+ if Array.length !input_seq > 0
+ then (prerr_endline "Error: Multiple input sequences specified"; exit 2)
+ else try input_seq := Util.read_integer_list items
+ with Util.Not_an_integer s ->
+ prerr_endline ("Error: The argument '" ^ s ^ "' is not an integer");
+ exit 2
+
+let input_from_file filename =
+ if Array.length !input_seq > 0
+ then (prerr_endline "Error: Multiple input sequences specified"; exit 2)
+ else try input_seq := Util.read_sequence_file filename
+ with Util.Not_an_integer s ->
+ prerr_endline ("Error: The item '" ^ s ^ "' from file '" ^ filename ^ "' is not an integer");
+ exit 2
+
+
+
+let speclist =
+ [("-f", Arg.String input_from_file, "Read input integer sequence from file");
+ ("-i", Arg.String input_immediate, "Integer input sequence with values separated by commas");
+ ("-l", Arg.Set_int depth_limit, "Sets maximum levels to generate, default 20")]
+
+
+
+module My_Algebra =
+ struct
+ type t = Poly.t
+ let zero = Poly.zero
+ let one = Poly.one
+ let add = Poly.add
+ let sub = Poly.sub
+ let mul = Poly.mul
+ let div = Poly.div
+ let to_string = Poly.to_string
+ end
+
+module My_Wall = Wall.Make (My_Algebra)
+
+
+
+let _ =
+ if Array.length Sys.argv <= 1 then (Arg.usage speclist usage_msg; exit 1);
+ Arg.parse speclist anon_fun usage_msg;
+ if Array.length !input_seq = 0 then (prerr_endline "Error: No input sequence provided"; exit 2);
+
+ if Array.length !input_seq < 2
+ then (prerr_endline "Error: Input sequence not long enough to solve"; exit 2);
+
+ (* Convert input integer sequence into the right sort of input polynomial sequence. *)
+ let poly_offset n =
+ Poly.of_list [(Z.neg !input_seq.(n), Z.one); (!input_seq.(n+1), Z.zero)] in
+ let poly_inputs = Array.init (Array.length !input_seq - 1) poly_offset in
+
+ (* Generate the number wall and look at the bottom most non-zero row. *)
+ let scribbles = My_Wall.create
+ ~dimx:(Array.length poly_inputs)
+ ~dimy:!depth_limit
+ ~init:poly_inputs in
+ let row_of_interest =
+ My_Wall.last_nonzero_row scribbles in
+
+ (* Test whether there is actually a row of zeros under the row of interest. *)
+ if row_of_interest = My_Wall.y_length scribbles - 1
+ then (prerr_endline "Error: Depth limit too low to solve sequence"; exit 2);
+
+ (* Test whether the row of interest has enough entries to be sure of them all being the same. *)
+ let exists_pred x =
+ if My_Wall.has_value scribbles ~x ~y:row_of_interest then 1 else 0 in
+ let array_sum =
+ Array.fold_left ( + ) 0 in
+ let naturals_to n =
+ Array.init n (fun n -> n) in
+ if array_sum (Array.map exists_pred (naturals_to (My_Wall.x_length scribbles))) < 3
+ then (prerr_endline "Error: Input sequence not long enough to solve"; exit 2);
+
+ (* Need to normalize the polynomials in the row as if they were set equal to zero. *)
+ let rec gcd a b =
+ Z.(if b = zero then a else gcd b (a mod b)) in
+ let normalize p =
+ let sign_correct x =
+ if Z.(Poly.get_coeff ~expnt:(Poly.degree x) x < zero) then Poly.neg x else x in
+ let coefficients =
+ List.map fst (Poly.to_list p) in
+ let divisor =
+ List.fold_left gcd Z.zero coefficients in
+ sign_correct (Poly.div p (Poly.of_zint divisor)) in
+
+ (* Now let's look at the polynomial right in the middle. *)
+ let key_poly =
+ normalize (My_Wall.get scribbles
+ ~x:((Array.length poly_inputs) / 2 + 1)
+ ~y:row_of_interest) in
+
+ (* Check for degeneracy. *)
+ if Poly.degree key_poly = Z.zero
+ then (prerr_endline "Error: Result is degree zero, not sure how that happened"; exit 2);
+
+ (* Check that all the polynomials in the row of interest are equal once normalized. *)
+ let equals_pred x =
+ let value = My_Wall.get_opt scribbles ~x ~y:row_of_interest in
+ if Option.is_some value
+ then normalize (Option.get value) = key_poly
+ else true in
+ if not (Array.for_all equals_pred (naturals_to (My_Wall.x_length scribbles)))
+ then (prerr_endline "Error: Polynomials differ, maybe try a longer sequence?"; exit 2);
+
+ (* Print out the polynomial as if it was a recurrence relation with the highest degree term *)
+ (* being on the left hand side and the remaining terms being on the right hand side. *)
+ let print_relation p =
+ let lhs = List.hd p in
+ let rhs = List.map (fun (x,y) -> (Z.neg x, y)) (List.tl p) in
+ if fst lhs <> Z.one then print_string (Z.to_string (fst lhs)); print_string "S(n) =";
+ let print_rhs_item minus plus (c,e) =
+ print_string " ";
+ if c < Z.zero then print_string minus else print_string plus;
+ if (Z.abs c) <> Z.one then print_string (Z.to_string (Z.abs c));
+ print_string ("S(n-" ^ Z.to_string Z.(snd lhs - e) ^ ")") in
+ print_rhs_item "-" "" (List.hd rhs);
+ List.iter (print_rhs_item "- " "+ ") (List.tl rhs) in
+
+ print_relation (Poly.to_list key_poly);
+ print_newline ();
+
+ (* This... should probably never trigger, but might as well check anyway. *)
+ if Poly.degree key_poly > Z.of_int (Array.length !input_seq)
+ then (prerr_endline "Error: Recurrence relation longer than input sequence... wtfhow?"; exit 2);
+
+ (* Calculate and print the predicted next number of the input sequence. *)
+ let calc_next p inp =
+ let lhs = List.hd p in
+ let rhs = List.map (fun (x,y) -> (Z.neg x, y)) (List.tl p) in
+ let rec process result terms =
+ if terms = []
+ then result
+ else let current = List.hd terms in
+ let index = Array.length inp - Z.(to_int (snd lhs - snd current)) in
+ process Z.(result + fst current * inp.(index)) (List.tl terms) in
+ Z.div (process Z.zero rhs) (fst lhs) in
+
+ print_string ("Next: " ^ Z.to_string (calc_next (Poly.to_list key_poly) !input_seq));
+ print_newline ()
+
+
generated by cgit (git 2.30.0) at 2024-12-22 19:44:32 +0000