const threads:=4;
const chunk:=10000;
const hist_digits:=6;

uses prng;
uses charset;
const hist_size:int:=ipower(3,hist_digits);

fn bit_len(x:int):int
[
	return (bsr x)+1;	
]

{ Returns 1 more than the steps necessary to get to 1
  1 returns  1
  2 returns  2
  3 returns  8
  4 returns  3
  5 returns  6
  6 returns  9
  7 returns 17
  8 returns  4 ...}
fn count_steps_to_1(input: int):int
[
	var x:=input;
	var steps:=1; { Put 1 instead of 0 so that powers of two have a ratio 1. }
	while true do [
		var zero_bits:=bsf x;
		x shr= zero_bits;
		steps += zero_bits;
		if x <= 1 then goto hit_one;
		x := x*3+1;
		steps += 1;
	]
hit_one:
	return steps;
]

fn calc_ratio~inline(x: int):real
[
	return cast(real,count_steps_to_1(x))/ bit_len(x);
]

fn collatz_scan_finite(start: int):list(real)
[
	var best_int:=start;
	var best_ratio:real:=0;
	var rlist:=fill(0.,chunk);
	
	for offset:= 0 to chunk do [
		rlist[offset]:=calc_ratio(start+2*offset);
	]
	return rlist;
]


{ Scans only odd numbers. If even given, increments to make odd. }
fn collatz_scan_best_from_chunk(implicit w: world, h:list(handle), base: int):world
[

	base or= 1;
	while true do [
		var best_ratio:real:=0;
		var best_number:=0;

		var acc_list:=empty(list(real));
		for thread := 0 to threads do acc_list+<=collatz_scan_finite~spark(base+2*thread*chunk);

		for thread := 0 to threads do [
			for offset:=0 to chunk do [
				var current_ratio:=acc_list[thread][offset];
				if current_ratio > best_ratio then [
					var current_number:=base+2*(thread*chunk+offset);
					best_ratio:=current_ratio;
					best_number:=current_number;
				]
			]
		]
		base+=2*threads*chunk;

		write(h[1],ntos(best_number)+" has ratio "+ntos(best_ratio)+nl);
		write(h[1],"Checked up to "+ntos(base-2)+nl);
		xeval w;
	]
	return w;
]

{ Scans only odd numbers. If even given, increments to make odd. }
fn collatz_scan_infinite(implicit w: world, h:list(handle), base: int):world
[
	var best_ratio:real:=0;

	base or= 1;
	while true do [
		var acc_list:=empty(list(real));
		for thread := 0 to threads do acc_list+<=collatz_scan_finite~spark(base+2*thread*chunk);

		for thread := 0 to threads do [
			for offset:=0 to chunk do [
				var current_ratio:=acc_list[thread][offset];
				if current_ratio > best_ratio then [
					var current_number:=base+2*(thread*chunk+offset);
					write(h[1],ntos(current_number)+" has ratio "+ntos(current_ratio)+nl);
					xeval w;
					best_ratio:=current_ratio;
				]
			]
		]
		base+=2*threads*chunk;
	]
	return w;
]

{ Scans only odd numbers. If even given, increments to make odd. }
fn collatz_scan_level(implicit w: world, h:list(handle), base: int, level: real):world
[
	base or= 1;
	while true do [
		var acc_list:=empty(list(real));
		for thread := 0 to threads do acc_list+<=collatz_scan_finite~spark(base+2*thread*chunk);

		for thread := 0 to threads do [
			for offset:=0 to chunk do [
				var current_ratio:=acc_list[thread][offset];
				if current_ratio >= level then [
					var current_number:=base+2*(thread*chunk+offset);
					write(h[1],ntos(current_number)+" has ratio "+ntos(current_ratio)+nl);
					xeval w;
				]
			]
		]
		base+=2*threads*chunk;
		write(h[1],"Checked up to "+ntos(base-2)+nl);
		xeval w;
	]
	return w;
]

fn collatz_forward_calc(implicit w: world, h:list(handle), input: int):world
[
	var steps:=count_steps_to_1(input);
	var ratio:real;
	ratio:=steps;
	ratio:=ratio / bit_len(input);
	write (h[1],"Input "+ntos(input)+" "+ntos_base(input,2)+" steps: "+ntos(steps)+", bit length "+ntos(bit_len(input))+", ratio="+ntos(ratio)+nl);
]

fn inc_list_position(l: list(int), pos: int):list(int)
[
	if (pos >= len(l)) then
		l += fill(0,pos-len(l)+1);
	l[pos]+=1;
	return l;
]

fn add_to_hist(hist: array(list(int),[hist_size]), x: int):array(list(int),[hist_size])
[
	{Here x is even or odd}

	{ First reduce to an odd number }
	var zero_bits:=bsf x;
	x shr= zero_bits;
	if x <= 1 then goto hit_one;

	{Here x is odd}

	{If the result is divisible by 3 we are in a dead end street and we need to get out of the branch and reduce once more, we won't end up divisible by 3 because
         divisible by 3 branches don't have any subbranches, so we couldn't have been in one.}
	if x mod 3  = 0 then [
		x := x*3+1;
		{Here x is even}
		zero_bits:=bsf x;
		x shr= zero_bits;
		if x <= 1 then goto hit_one;
		{Here x is odd}
	]

	{Here x is odd}

	while true do [
		x := x*3+1;
		{Here x is even}
		zero_bits:=bsf x;
		x shr= zero_bits;
		{Here x is odd}
		var remainder:=x mod hist_size;
		hist[remainder]:=inc_list_position(hist[remainder],zero_bits);
		if x <= 1 then goto hit_one;
	]
hit_one:

	return hist;
]

fn itemize(l: list(int)):list(int)
[
	var out:=empty(int);
	for i:=0 to len(l) do out+=sparse(i,l[i]);
	return out;
]

fn hist_row2bytes(input: list(int)):bytes
[
	var out:="";
	for i:=0 to len(input) do [
		if input[i]>0 then [
			if len(out)>0 then out+=",";
			out+=ntos(input[i])+"×"+ntos(i);
		]
	]
	return out;
]

fn hist_itemize(hist:array(list(int),[hist_size])):array(list(int),[hist_size])
[
	var hist_itemized:=hist; { Prevent false accusation of unitialized variable }

	for i := 0 to hist_size do [
		hist_itemized[i]:=itemize(hist[i]);
	]
	return hist_itemized;
]

fn train_on_internal_numbers():array(list(int),[hist_size])
[
	var hist:=array_fill(empty(int),[hist_size]);


	{ Long-running numbers:
	  https://oeis.org/A284668/list

157348943716351847 has ratio 38.8448275862
931386509544713451 has ratio 38.06667
93571393692802302 has ratio 37.33929
942488749153153 has ratio 37.26000
7579309213675935 has ratio 36.96226
63728127 has ratio 36.53846
7887663552367 has ratio 36.37209
13371194527 has ratio 35.6176470588
80867137596217 has ratio 35.38298
127456255 has ratio 35.2222222222
95592191 has ratio 35.1111111111
14500812391 has ratio 35.0882352941
31694683323 has ratio 34.8571428571
12235060455 has ratio 34.8529411765
17828259369 has ratio 34.6857142857
26742389055 has ratio 34.6285714286
20056791791 has ratio 34.5428571429
30085187687 has ratio 34.4857142857
9780657630 has ratio 34.30303
4890328815 has ratio 34.303030
7335493223 has ratio 34.2424242424
226588897 has ratio 34.1785714286
19334416521 has ratio 34.1714285714
29001624783 has ratio 34.1142857143
169941673 has ratio 34.0714285714
21751218587 has ratio 34.0285714286
254912509 has ratio 34.0000000000
6189322407 has ratio 34.0000000000
12212032815 has ratio 33.9411764706
24470120911 has ratio 33.8857142857
18352590683 has ratio 33.8000000000
989345275647 has ratio 33.72500
20646664519 has ratio 33.6571428571
13040876841 has ratio 33.4117647059
25730702889 has ratio 33.3714285714
26130934783 has ratio 33.3714285714
9780657631 has ratio 33.3235294118
19298027167 has ratio 33.2857142857
14670986447 has ratio 33.2647058824
75128138247 has ratio 33.21622
11003239835 has ratio 33.1764705882
27528886025 has ratio 33.7428571429
8140184737 has ratio 33.1515151515
16504859753 has ratio 33.1176470588
12378644815 has ratio 33.0294117647
24424065631 has ratio 33.0000000000
	}

	var training_list:=list(int).[
157348943716351847,
931386509544713451,
93571393692802302,
942488749153153,
7579309213675935,
63728127,
7887663552367,
13371194527,
80867137596217,
127456255,
95592191,
14500812391,
31694683323,
12235060455,
17828259369,
26742389055,
20056791791,
30085187687,
9780657630,
4890328815,
7335493223,
226588897,
19334416521,
29001624783,
169941673,
21751218587,
254912509,
6189322407,
12212032815,
24470120911,
18352590683,
989345275647,
20646664519,
13040876841,
25730702889,
26130934783,
9780657631,
19298027167,
14670986447,
75128138247,
11003239835,
27528886025,
8140184737,
16504859753,
12378644815,
24424065631,
	];
	
	for training_number in training_list do hist:=add_to_hist(hist, training_number);	
	return hist;
]

fn sig2bytes(x: int):bytes
[
	return list_left_pad(ntos_base(x mod hist_size,3),hist_digits,'0');
]

fn rnd_from_list(l: list(int), state: prng_state):(prng_state,int)
[
	var rn:int;
	state,rn := prng_get_uint32(state);
	return state,l[rn mod len(l)];
]

fn number_summary(x: int, hist:array(list(int),[hist_size]), steps backtracking_step:int):bytes
[
	var ratio:real:=steps;
	var bits:=bit_len(x);
	ratio /= bits;
	return "Class "+sig2bytes(x)+" number "+ntos_base(x,2)+" "+ntos(x)+nl
              +" Ratio "+ntos_base_precision(ratio, 10, 5)+", "+ntos(steps)+" steps, "+ntos(bits)+" bits, backtracking step "+ntos(backtracking_step)+nl
              +" Selection: "+hist_row2bytes(hist[x mod hist_size]);
]

fn number_summary_nohist(x: int, steps:int):bytes
[
	var ratio:real:=steps;
	var bits:=bit_len(x);
	ratio /= bits;
	return "Class "+sig2bytes(x)+" number "+ntos_base(x,2)+" "+ntos(x)+nl
              +" Ratio "+ntos_base_precision(ratio, 10, 5)+", "+ntos(steps)+" steps, "+ntos(bits)+" bits"+nl;
]

fn print_histogram(implicit w: world, h:list (handle)):world
[
	var hist:=train_on_internal_numbers();
	write (h[1],ntos(hist_digits)+" histogram base 3 digits, histogram size 3^"+ntos(hist_digits)+"="+ntos(hist_size)+nl);
	for i := 0 to hist_size do write(h[1],"["+sig2bytes(i)+"]:"+hist_row2bytes(hist[i])+nl);
]

fn collatz_stats(implicit w: world, h:list(handle), x:int, n_steps:int, seed:int):world
[
	var hist:=train_on_internal_numbers();
	var hist_itemized:=hist_itemize(hist);
	var state := prng_init(seed);
	var bit_shift:int; { Random Number }

	x shr= bsf x; { Normalize to odd }
	var steps_to_1:=count_steps_to_1(x);
	
	for step:=0 to n_steps do [
		write(h[1],number_summary(x, hist, steps_to_1,step+1));
		var itemized_list:=hist_itemized[x mod hist_size];
		if len(itemized_list) <= 0 then [
			write(h[1],"no selection, don't know what to do in that branch, getting out of the branch and going into the next branch (x=4x+1)."+nl);
			x:=4*x+1;
			steps_to_1+=2;
		]else[
			state,bit_shift:=rnd_from_list(itemized_list, state);
			write(h[1],", selected bit shift "+ntos(bit_shift)+nl);

			x shl= bit_shift;
			steps_to_1+=bit_shift;

			eval assert(x mod 3 = 1,"Before reverse of 3x+1 must have remainder 1 mod 3");

			x:=(x-1) div 3;
			if x=1 then steps_to_1:=1;
                               else steps_to_1+=1;
		]
	]
	write(h[1],nl);
]

fn collatz_forward_print(implicit w: world, h:list(handle), input: int):world
[
	var x:=input;
	var steps_initial:=count_steps_to_1(input); { Put 1 instead of 0 so that powers of two have a ratio 1. }
	var steps_going_down:=steps_initial;
	while true do [
		var zero_bits:=bsf x;
		x shr= zero_bits;
		steps_going_down -= zero_bits;
		write (h[1],number_summary_nohist(x,steps_going_down));
		xeval w;
		if x <= 1 then goto hit_one;
		x := x*3+1;
		steps_going_down -= 1;
	]
hit_one:

	var ratio:=cast(real,steps_going_down)/bit_len(input);
	write (h[1],"Input "+ntos(input)+" "+ntos_base(input,2)+" steps: "+ntos(steps_initial)+", bit length "+ntos(bit_len(input))+", ratio="+ntos(ratio)+nl);
]

fn little_endian_4B(x: int):bytes
[
	var out:="";

	for i:=0 to 4 do [
		out+<=x and #FF;
		x shr=8;
	]
	return out;
]

fn collatz_forward_wav(implicit w: world, h:list(handle), input: int, finish_mode:bytes, sample_rate_Hz: int):world
[
	var sample_rate_LE:=little_endian_4B(sample_rate_Hz);
	var wav_header:="RIFF"+list(byte).[#FF,#FF,#FF,#FF]+"WAVEfmt "+list(byte).[16,0,0,0,1,0,1,0]
                        +sample_rate_LE+sample_rate_LE+list(byte).[1,0,8,0]+"data"+list(byte).[#FF,#FF,#FF,#FF];

	write(h[1],wav_header);
restart:
	var x:=input;
	while true do [
		var zero_bits:=bsf x;
		x shr= zero_bits;
		write (h[1],map(ntos_base(x,2),lambda(b:byte):byte[return 96+64*(b-'0');]));
		xeval w;
		if x <= 1 then [
			if (finish_mode="s") then return w;
			if (finish_mode="i") then input +=1;
			goto restart;
		]
		x := x*3+1;
	]
]


fn main
[
        if len (args)<1 then [
print_usage:
                write(h[2],""
			+nl+"collatz.ajla - a program to explore Collatz conjecture"
			+nl+"======================================================"
			+nl+"Collatz conjecture Wikipedia: https://en.wikipedia.org/wiki/Collatz_conjecture"
			+nl+"                        "
			+nl+"Functions of collatz.ajla"
			+nl+"-------------------------"
			+nl+"* Discover new Collatz numbers that have high ratio of steps to bit length, using brute force search (functions b/B,l/L,s/S)"
			+nl+"* Discover new Collatz numbers that have high ratio of steps to bit length, using heuristic backtracking in the Collatz graph (function r/R)"
			+nl+"* Examine Collatz conjecture numbers, especially in regard to the ratio of steps to bit length (functions c/C,p/P)"
			+nl+"* Sonification of Collatz sequences (function w/W)"
			+nl+""
			+nl+"Definition of the ratio"
			+nl+"-----------------------"
			+nl+"I define the ratio as ratio of steps a number takes to its bit length. Bit length is the number of bits necessary to write"
			   +" the number or the number of the topmost bit plus 1."
			+nl+"I calculate number of steps so that number 1 has 1 step, because this way, the lowest ratio is in powers of 2 and is always exactly 1.000."
			+nl+""
			+nl+"Backtracking in reverse graph using heuristic (functions h,r/R)"
			+nl+"---------------------------------------------------------------"
			+nl+"collatz.ajla first learns how long-playing numbers do it that they are"
			   +" long-playing. It has a bunch of long-play numbers programmed in on which it trains a histogram heuristic."
                           +" It stores the learned information into a histogram (can be printed out using the h function) based on several last"
			   +" digits of base 3 representation of the number (can be changed by hist_digits at the beginning of collatz.ajla). It then randomly generates"
			   +" number of bit shifts (multiplications by 2) before entering a branch, based on"
			   +" local situation (last digits in base 3). This way it tries to increase number"
			   +" of bits slowly while increasing number of iterations fast to get an"
			   +" advantageous ratio.  Its success is mediocre, it's better to use the brute"
			   +" force method (l,L,s,S) to search for long-playing numbers, but it succeeds to find a better number in some cases:"
			+nl+""
			+nl+"Example backtracking success (function r,R)"
			+nl+"-------------------------------------------"
			+nl+"Example backracking (with hist_digits=6)"
			+nl+"ajla collatz.ajla s 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 8 6"
			+nl+"Produced a number with ratio 14.58772 from a number with ratio 14.57522, increases bit count from 113 to 114 and step count from 1647 to 1663"
			+nl+"ajla collatz.ajla s 101111 2 0"
			+nl+"Produced a number with ratio 21.4 from a number with ratio 17.5, decreases bit count from 6 to 5 and increased step count from 105 to 107."
			+nl+""
			+nl+"r,R function - reverse trace the Collatz graph (backtracking)"
			+nl+"-------------------------------------------------------------"
			+nl+"collatz.ajla has a bunch of long-play numbers programmed in on which it trains the histogram heuristic. "
			+nl+""
			+nl+"Then it starts with a given starting number and gets out of the branch by doing"
			   +" 3x+1 and counts how many times to divide by 2 before arriving at an odd number."
			   +" At that odd number it fingerprints the local graph situation backwards by"
			   +" taking last several digits in the base 3 representation of that number. It"
			   +" keeps a histogram for each such local graph situation how many bit shifts were"
			   +" performed before arriving to that local graph situation odd number."
			+nl+""
			+nl+"For example:"
			+nl+"[2222]:1098×1,2×5"
			+nl+"The collatz.ajla program figured out that to arrive to an odd number with fingerprint 2222,"
                           +" it arrived 1098 times after 1 division by 2 and 2 times after 5 divisions by 2."
			+nl+""
			+nl+"If histogram is empty for a given base-3 fingerprint (typically happens when last digit of base 3 is 0,"
			   +" because branches divisible by 3 have no subbranches, but also happens when"
			   +" there is no training example that would go through such local branching situation), it"
			   +" doesn't know where to go. So it gets out of the branch, goes 2 steps up"
			   +" (multiplies by 2 two times) and gets into the next branch. This corresponds to"
			   +" the operation 4x+1, or by appending ""01"" to the end of the binary number."
			+nl+""
			+nl+"When we backtrack attempting to discover new long playing number, we can use"
			   +" this knowledge to generate probable branches into which we should go to get to"
			   +" advantageous numbers."
			+nl+""
			+nl+"I tried greedy algorithm by going into lowest branch but it produced a horrible result."
			+nl+""
			+nl+"Sometimes a better, longer-playing number can be found by taking one known long-playing number (example 931386509544713451 published on"
                           +" https://oeis.org/A284668/list which has ratio of 38.06667) and using the p/P function of collatz.ajla to see how it goes down to 1"
                           +" and it goes through an even longer playing number, in this case 157348943716351847 with ratio of 38.8448275862."
			+nl+""
                        +nl+"Usage: ajla collatz.ajla b input_number_binary"
                        +nl+"------ ajla collatz.ajla B input_number_decadic"
                        +nl+"       ajla collatz.ajla c input_number_binary"
                        +nl+"       ajla collatz.ajla C input_number_decadic"
                        +nl+"       ajla collatz.ajla h"
                        +nl+"       ajla collatz.ajla l input_number_binary min_ratio_real"
                        +nl+"       ajla collatz.ajla L input_number_decadic min_ratio_real"
                        +nl+"       ajla collatz.ajla p input_number_binary"
                        +nl+"       ajla collatz.ajla P input_number_decadic"
                        +nl+"       ajla collatz.ajla r input_number_binary n_steps_decadic seed_decadic"
                        +nl+"       ajla collatz.ajla R input_number_decadic n_steps_decadic seed_decadic"
                        +nl+"       ajla collatz.ajla s input_number_binary"
                        +nl+"       ajla collatz.ajla S input_number_decadic"
                        +nl+"       ajla collatz.ajla w input_number_binary finish_mode sample_rate_Hz"
                        +nl+"       ajla collatz.ajla W input_number_decadic finish_mode sample_rate_Hz"
                        +nl+""
                        +nl+"b,B - start brute force scan at given number and print the ones reaching best ratio from each computation chunk. "
                           +"Scans odd numbers only. If even given, increments to get odd."
                        +nl+"c,C - calculate forward, don't print intermediate results, calculate statistics"
                        +nl+"h   - print histogram result from the builtin training numbers."
                        +nl+"l,L - start scan at given number and print the ones reaching at least given ratio. Scans odd numbers only. If even given, increments to get odd."
                        +nl+"p,P - forward calculation, print intermediate odd results, calculate statistics."
                        +nl+"r,R - reverse trace the Collatz graph using histogram-based random heuristic."
                        +nl+"s,S - start brute force scan at given number and print the ones reaching best ratio so far. Scans odd numbers only. If even given, increments to get odd."
                        +nl+"w,W - WAV. Produce WAV to the standard output to listen to the bits of the forward calculation."
                        +nl+"      Example usage: ajla collatz.ajla w 111111111111111111111111111111111110111 i 11025 | mpv -"
                        +nl+"      ""finish_mode"" determines what to do when Collatz "
                           +"calculation finishes:"
		        +nl+"        i - increment - increments the input number and runs again, indefinitely."
		        +nl+"        r - repeat    - runs the same number again,                 indefinitely."
		        +nl+"        s - stop      - stops, returns to commandline. WAV has however still set length to infinite (0xFFFFFFFF)."
		        +nl+"      Funny noises: changes at      4.8 s:   111111111111111111111111111111111110111 i 11025"
		        +nl+"                    changes around 10   s:   11111111111111111111111111100000 i 11025"
		        +nl+"                    alarm-like:              1010101010101010101010101010101010101010101001 r 11025"
		        +nl+"                    machine gun-like:        11111111111111111111 r 8000"
		        +nl+"                    another sound:           10000000000000000000000000000000 i 11025"
		        +nl+"                    noise with rep. ringing: 11011011011011011011011011011011011011011011011011011011011011011011011011011011 r 8000"
		        +nl+"                    ""alien"" signal:          110110110110110110110110110110110110110110110110110110110110110110110110110110110110110110110110110111111111"
                           +" r 4800"
		        +nl+"                    flying saucer takeoff:   11011011011011011011011011011011011011011011011011011011011011011011011011011011011011011011011011 s 480"
                        +nl+""
                        +nl);
                return w;
        ]
        var x:=0;

	if args[0]<>"h" then [
		if len(args) < 2 then goto print_usage;
		if bytes_upcase(args[0]) = args[0] then x:=ston(args[1]);
				 else x:=ston_base(args[1],2);
		if is_exception x then [
			eval forced_exit(33,"Error: Input number must be binary only digits 0 and 1: "+args[1]+nl);
			while true do [ ]
		]
	]

	       if bytes_locase(args[0])="b" then collatz_scan_best_from_chunk(h,x);
	  else if bytes_locase(args[0]) = "c" then collatz_forward_calc(h,x);
	  else if args[0]="h" then print_histogram(h);
	  else if bytes_locase(args[0]) = "l" then [
		if len (args) < 3 then goto print_usage;
		var level:=instance_real_number_real64.from_bytes(args[2]);
		collatz_scan_level(h,x,level);
	] else if bytes_locase(args[0]) = "p" then collatz_forward_print(h,x);
	  else if bytes_locase(args[0]) = "r" then [
		if len (args) < 4 then goto print_usage;
		var n_steps:=ston(args[2]);
		var seed:=ston(args[3]);
		collatz_stats(h,x,n_steps,seed);
	] else if bytes_locase(args[0])="s" then collatz_scan_infinite(h,x);
	  else if bytes_locase(args[0])="w" then [
		if len (args)<4 then goto print_usage;
		var finish_mode:=args[2];
		var sample_rate_Hz:=ston(args[3]);
		collatz_forward_wav(h,x,finish_mode, sample_rate_Hz);
	]
	  else goto print_usage;
]