# restart; read `G:/My Drive/Aek/Arith/Sum.txt`;

with(numtheory):
with(plots):

# First version: June 7, 2022 

# The program follows the paper
# GAME VALUES OF ARITHMETIC FUNCTIONS
# of Iannucci and Larsson

# I also added some functions to section
# Saliquant and Non-Totient that are in
# our paper.


#########################
# Section 1: 
# Basic and some functions
#
# Functions:
# Mex(S), NimSum(S),
# PosOf(A,k),
#
##########################

# Input: The set of non-negative
# integers S
# Output: Mex function of S, i.e.
# the first non-negative integer
# missing
# Try: Mex({0,1,3,6,2,4,12});
Mex := proc(S) option remember;
local s;

s := 0;

while member(s,S) do 
	s := s+1;
od:

return(s);
end:

# Input: The list of 
# non-negative integers S
# Output: NimSum of numbers in S
# Try: NimSum([7,5,3]);
NimSum := proc(S) option remember;
local i,s,T,t,m;

T := [seq(convert(s,base,2), s in S)];
m := max(seq( nops(t), t in T));

T := [seq( [op(t),0$(m-nops(t))] , t in T)];
T := [seq(add( t[i], t in T) mod 2,i=1..m)];

add( T[i]*2^(i-1),i=1..nops(T));
end:


# Input: List A and integer k 
# Output: Set of the positions 
# of k in list A
# Try:
# A := [seq(n-Saliquant(2*n),n=1..100)]:
# PosOf(A,1);
PosOf := proc(A,k) option remember;
local i,S;

S := {};

for i from 1 to nops(A) do
	if A[i]=k then
		S := S union {i};
	fi:	
od:
return(S);
end:

# Calculate density of A
# Input: List A
# Output: List of the densities of value k in A
# Try: 
# A := [seq(Saliquant(2*n),n=1..100)];
# Den(A);
Den := proc(A) option remember;
local i,S;

S := [0$max(A)];

for i from 2 to nops(A) do
	S[A[i]] := S[A[i]]+1;	
od:

evalf(S/nops(A));
end:


#####################
# Section 2.1
# Maliquot, Saliquot
#####################

# Divisors game 1a)
# Input: integer n
# Output: SG(n) of Maliquot
# Try:
# Maliquot(15);
# [seq(Maliquot(n),n=1..15)];
Maliquot := proc(n) option remember;
local s,S,T;

S := divisors(n) minus {n};

T := {seq( Maliquot(s),s in S)};
return(Mex(T));
end:


# Subtract divisors game 1b)
# Input: integer n
# Output: SG(n) of Saliquot
# Try:
# Saliquot(15);
# [seq(Saliquot(n),n=1..15)];
Saliquot := proc(n) option remember;
local s,S,T;

S := divisors(n);
T := {seq( Saliquot(n-s),s in S)};
return(Mex(T));
end:


# 2-valuation
# Input: integer n
# Output: Value of 2-valuation
# i.e. number of 2 as a factor +1
# Try:
# [seq(TwoValS(n),n=0..18)];
# [seq(Saliquot(n),n=0..18)];
TwoValS := proc(n) option remember;
local t,m;

if n = 0 then
	return(0);
fi:	

t := 1:
m := n;
while m mod 2 = 0 do
	m := m/2;
	t := t+1:
od:	
	
t;
end:


#######################
# Section 2.2
# Maliquant, Saliquant
#######################

# Move to a non-divisor
# Input: non-neg integer n
# Output: SG(n) of Maliquant
# Try:
# Maliquant(12);
# [seq(Maliquant(n),n=0..20)];
Maliquant := proc(n) option remember;
local s,S,T;

if n=0 then
	return(0);
fi:	

S := {seq(i,i=0..n)} minus  divisors(n);
T := {seq( Maliquant(s) ,s in S)};
return(Mex(T));
end:


# Half odd
# Input: non-neg integer n
# Output: the index function
# for largest odd factor of n
# Try: 
# [seq(HalfOdd(n),n=0..20)];
# [seq(Maliquant(n),n=0..20)];
HalfOdd := proc(n) option remember;
local r;

if n = 0 then
	return(0);
fi:

r := n/2^(TwoValS(n)-1); 
(r+1)/2;
end:


# Subtract a non-divisor
# Input: non-neg integer n
# Output: SG(n) of Saliquant
# Try:
# Saliquant(12);
# [seq(Saliquant(n),n=0..30)];
# B :=[seq(Saliquant(4*k+2)/((17*k+4)/9),k=1..2000)]:
# sort(evalf(B));
Saliquant := proc(n) option remember;
local s,S,T;

S := {seq(i,i=1..n)} minus divisors(n);
T := {seq( Saliquant(n-s) ,s in S)};
Mex(T);
end:

# Input: integer pow and N
# Output: set of values 2k+1 
# 1<= k <=n where
# Saliquant(2^(pow+1)*k+2^pow)
# = (2^pow-1)*k+2^(pow-1)-1.
# Try:
# TestSali(1,500);
# TestSali(4,500);
TestSali := proc(pow,N) option remember;
local k,a,S,A;

S := {};

for k from 1 to N do
	if Saliquant(2^(pow+1)*k+2^pow)
	=(2^pow-1)*k+2^(pow-1)-1 then
	S := S union {2*k+1};
	fi:
od:
S;
end:

# Search for modulo that 
# TestSali contains
#  Input: integer pow and N
# Try:
# SearchMod(2,500);
SearchMod := proc(pow,N) option remember;
local a,b,A,B,S;

S := [];
A:=TestSali(pow,N);

for b from 4 to 20 do
for a from 1 to b-1 do
if gcd(a,b) = 1 then
	B :=Andrew(167,a,b);
	if max(B) > max(A) then
		return(ERROR);
	fi:	
	if B subset A then
		S := [op(S),[a,b]];
	fi:
fi: od: od:

return(S);	
end:


# Input: integer N, a, b
# Output: set of prime (up to Nth prime)
# where p = a mod b
# Try:
# Andrew(100,5,6);
Andrew := proc(N,a,b) option remember;
local i,S;

S :={};
for i from 1 to N do
	if ithprime(i) mod b = a then
		S := S  union {ithprime(i)};
	fi:
od:
S;	
end:


# Conjecture 2 in our paper
# Input: Integers a and b
# Output: m where SG((2a+1)*2^b) = f(a,b,m),
# Try:
# [seq(Conj2(3,b),b=1..10)];
Conj2 := proc(a,b) option remember;
local m;
m:=0;

while Saliquant((2*a+1)*2^b) 
<> (2*a+1)*2^(b-1)-1/2*((2*a+1)/(2*m+1)+1) do
	m := m+1;
od:
	
m;
end:


########################
# Section 2.3
# Totative, Nontotative
########################

# Moving to a relatively 
# prime residue
# Input: non-neg integer n
# Output: SG(n) of Totative funciton
# Try:
# [seq(Totative(n),n=1..14)];
Totative := proc(n) option remember;
local s,S,T;

if n =1 then
	return(0);
fi:
	
S :={};
for s from 1 to n do
	if gcd(s,n) = 1 then
		S := S union {s};
	fi:
od:

T := {seq(Totative(s), s in S)};
Mex(T);
end:

# the order of smallest 
# prime divide n
# Input: non-neg integer n
# Output: the index of the
# smallest prime divisor of n
# (Theorem 5 of Larsson)
# Try:
# [seq(i0prime(n),n=1..44)];
# [seq(Totative(n),n=1..44)];
i0prime := proc(n) option remember;
local s;

if n =1 then
	return(0);
fi:

s := 1:
while n mod ithprime(s) <> 0 do
	s := s+1;
od:	

return(s);
end:


# gcd > 1
# Input: non-neg integer n
# Output: SG(n) of Non-totative funciton
# Try:
# [seq(NonTotative(n),n=1..14)];
NonTotative := proc(n) option remember;
local s,S,T;

S :={};
for s from 0 to n-1 do
	if gcd(s,n) > 1 then
		S := S union {s};
	fi:
od:

T := {seq(NonTotative(s), s in S)};
Mex(T);
end:


###################
# Section 3.1
# MyC(n), RecC(n),
# Totient(n),
# NonTotient(n),
#
###################

# Idea of class function in my paper
# (Lemma 4.4 etc)
# Input: non-neg integer n
# Output: number of iteration
# of phi(n) until get 2
# Try:
# [seq(MyC(n),n=1..20)];
MyC := proc(n) option remember;
local i,r;

if n = 1 or n=2 then
	return(0);
fi:
	
r := phi(n);
i := 1;

while r <> 2 do
	r := phi(r); 
	i := i+1;
od:
i;
end:


# MyC with recurrence
# Input: non-neg integer n
# Output: Recurrence to get
# MyC(n)
# Try:
# RecC(20);
# [seq(MyC(n),n=1..20)];
# [seq(RecC(n),n=1..20)];
RecC := proc(n) option remember;
local s;

if n=1 or n=2 then
	return(0);
fi:

s := ifactors(n)[2][1][1];

if s = n then
	RecC(n-1)+1;
elif s mod 2 = 0 and n/s mod 2 = 0 then
	RecC(s)+RecC(n/s)+1;
else
	RecC(s)+RecC(n/s);
fi:
end:


# Number of d where gcd(d,n)=1
# Input: non-neg integer n
# Output: SG(n) of Totient funciton
# Try:
# [seq(Totient(n),n=1..30)];
Totient := proc(n) option remember;
if n=1 then
	return(0);
fi:

Mex( {Totient(phi(n))});
end:


# Input: non-neg integer n
# Output: Solution to SG(n)
# of Totient function 
# (Theorem 8 of Larsson)
# Try:
# [seq(SolTotient(n),n=1..30)];
# [seq(Totient(n),n=1..30)];
SolTotient := proc(n) option remember;
if n=1 then
	return(0);
fi:

MyC(n)+1 mod 2;
end:


# Number of d where gcd(d,n) <> 1
# Input: non-neg integer n
# Output: SG(n) of Non-totient funciton
# Try:
# [seq(NonTotient(n),n=1..30)];
NonTotient := proc(n) option remember;

if n=1 then
	return(0);
fi:	

Mex( {NonTotient(n-phi(n))});
end:


# Input: non-neg integer n
# Output: Class function C for non-totient
# Try:
# [seq(MyC2(n),n=1..30)];
# seq({seq(MyC2(2*n-1),n=2^(k-2)+1..2^(k-1))},k=2..20);
# A:=[seq(MyC2(3*(2*n-1)),n=65..128)];
MyC2 := proc(n) option remember;
if n=1 then
	return(0);
fi:

1+MyC2(n-phi(n));	
end:


###############
# Section 3.2
# MTau(n), STau(n),
#
###############

# Move to number of divisors minus 1
# Input: non-neg integer n
# Output: SG(n) of MTau funciton
# Try:
# [seq( MTau(n),n=1..30)];
MTau := proc(n) option remember;
if n = 1 then
	return(0);
fi:	
Mex( {MTau(nops(divisors(n))-1)});
end:

# Input: non-neg integer n
# Output: Class function of MTau
# This is equivalent to MTau mod 2
# Try:
# [seq( CTau(n),n=1..30)];
CTau := proc(n) option remember;
if n = 1 then
	return(0);
fi:	

1+CTau(nops(divisors(n))-1);
end:

# Subtract n from number of divisors
# Input: non-neg integer n
# Output: SG(n) of STau funciton
# Try:
# [seq( STau(n),n=1..30)];
STau := proc(n) option remember;
if n = 0 then
	return(0);
fi:	
Mex( {STau(n-nops(divisors(n)))});
end:


###############
# Section 3.3
# MBigOmega(n), 
# SmallOmega(n),
#
###############

# Move to number of 
# prime divisors: 
# counting multiplicity
# Input: non-neg integer n
# Output: SG(n) of MBigOmega funciton
# Try:
# seq(MBigOmega(n),n=1..20);
# A:=[seq(MBigOmega(n),n=1..1000000)]:
# evalf(add(a, a in A)/nops(A));
MBigOmega := proc(n) option remember;
if n=1 then
	return(0);
fi:	

Mex( {MBigOmega(bigomega(n))});
end: 


# Subtract n from number of 
# prime divisors
# Input: non-neg integer n
# Output: SG(n) of SBigOmega funciton
# Try:
# A:=[seq(SBigOmega(n),n=1..20)];
# PosOf(A,1);
SBigOmega := proc(n) option remember;
if n=1 then
	return(0);
fi:	

Mex( {SBigOmega(n-bigomega(n))});
end: 


# Class function
# Input: non-neg integer n
# Output: Class function of MBigOmega funciton
# Try:
# seq(CBigOmega(n),n=1..20);
# A:=[seq(CBigOmega(n) mod 2,n=1..20)];
# B:=[seq(MBigOmega(n),n=1..20)];
# A-B;
CBigOmega := proc(n) option remember;
if n=1 then
	return(0);
fi:	

1+ CBigOmega(bigomega(n));
end:


# Empirical desnity of class k
# Input: non-neg integer v and N
# Output: density of value v of
# CBigOmega up to N
# Try: 
# CountOM(1,1000000);
CountOM := proc(v,N) option remember;
local n,S;

S := 0;

for n from 1 to N do 
if CBigOmega(n) = v then
	S := S+1;
fi:
od:

evalf(S/N);
end:

# This is the  formula from LanDau's paper
# Try: LanDau(1,1000);
LanDau := proc(v,M) option remember;
local n,s,S,R;

S := {};
for n from floor(0.6*M) to 1.4*M do
	if CBigOmega(n) = v then
		S := S union {n};
	fi:
od:	

evalf(add( M^(s-1)/(s-1)!,s in S)/exp(M));
end:


#####################
# Testing for the 
# density problem
#
# Functions:
# MyPi(x,k), F1(z),
# DenseC1(x),
#
#####################

# Density of Pi(X,k)
# X = log(log(x)) 
# Try:
# MyPi(300,150);
MyPi := proc(X,k) option remember;
local A;


if k > 2*X then
	ERROR(KtooBig);
fi:	

A := F1(k/X)/exp(X)*X^(k-1)/(k-1)!;
evalf(A);
end:

# Try: F1(1);
F1 := proc(z) option remember;
local i,p,P,A;

P := seq(ithprime(i),i=1..1000);
A := 1/GAMMA(z+1)*mul( (1+z/(p-1))*(1-1/p)^z  ,p in P);

evalf(A);
end:

# Try: DenseC1(300):
DenseC1 := proc(x) option remember;
MyPi(x,1);
end:


# This one uses small x
# Try: CountC1(10^10):
CountC1 := proc(x) option remember;
local i,r;

i := floor(x/ln(x));
r := 1;

while ithprime(floor(i*r)) <= x do
	r := r+0.05; print(r);
od:

evalf(i*r/x);
end:


# Try: 
# DenseC2(100):
# seq(print(i,[DenseC2(ithprime(i*10))
# ,ApproxC2(ithprime(i*10))]),i=1..60);
DenseC2 := proc(X) option remember;
local k,S;

S := 0;

for k from floor(0.85*X) to floor(1.15*X) do
	if isprime(k) then
		S := S+MyPi(X,k);
	fi:
od:
S;	
end:

# Try:
# ApproxC2(100);
# seq(print([DenseC2(ithprime(i))
# ,ApproxC2(ithprime(i))]),i=50..55);
ApproxC2 := proc(X) option remember;
#evalf(1/ln(X));
local t,A;
A :=1/X*int(1/ln(t),t=2..X);
evalf(A);
end:


########################################

# Section 3.3.3
# Move to number of 
# prime divisors
# Input: non-neg integer n
# Output: SG(n) of MSmallOmega funciton
# Try:
# A:=[seq(MSmallOmega(n),n=1..20)];
# PosOf(A,1);
MSmallOmega := proc(n) option remember;
if n=1 then
	return(0);
fi:	

Mex( {MSmallOmega(OmegaS(n))});
end: 

# Input: non-neg integer n
# Output: Class function of 
# MSmallOmega funciton
# Try:
# A:=[seq(CSmallOmega(n),n=1..20)];
# PosOf(A,1);
CSmallOmega := proc(n) option remember;
if n=1 then
	return(0);
fi:	

1 + CSmallOmega(OmegaS(n));
end: 


# Subtract n from number of 
# prime divisors, no multiplicity
# Try:
# A:=[seq(SSmallOmega(n),n=1..20)];
# PosOf(A,1);
SSmallOmega := proc(n) option remember;
if n=1 then
	return(0);
fi:	

Mex( {SSmallOmega(n-OmegaS(n))});
end: 


# Number of prime factors 
# NOT counting multiplicity
# [seq(OmegaS(n),n=1..15)];
OmegaS := proc(n) option remember;
nops(ifactors(n)[2]);
end:


###############
# Section 4.1
# Dividing(n), 
# Omega2(n),
#
###############

# Try:
# Dividing(10);
# A:=[seq(Dividing(n),n=1..130)];
# B:=[seq(Omega2(n),n=1..130)];
Dividing := proc(n) option remember;
local s,a,S,A;

if n=1 then
	return(0);
fi:
	
S := divisors(n) minus {n};
A := {seq( NimSum([Dividing(s)$(n/s)]), s in S)};

Mex(A);
end:


# Omega2 as defined in section 4.1
# Try:
# [seq(Omega2(n),n=1..10)];
Omega2 := proc(n) option remember;
local s,m;

if n mod 2 = 0 then
	s := 1;
else
	s := 0;
fi:	

m := n;
while m mod 2 = 0 do
	m := m/2;
od:	

bigomega(m)+s;
end:



###############
# Section 4.2
# DivRes(n), 
# ComGrundy(n),
# DivThrowRes(n), SolDivThrowRes(n),
# ResThrowDiv(n), SolResThrowDiv(n),  
#
###############

# Try:
# DivRes(10);
# [seq(DivRes(n),n=1..20)];
DivRes := proc(n) option remember;
local d,S,T;

if n=0 or n=1 then
	return(0);
fi:

S := {};	
for d from 1 to n-1 do
	T := [DivRes(d)$floor(n/d)
	,DivRes(n-d*floor(n/d))]; 
	S := S union {NimSum(T)};
od:

Mex(S);
end:


# Compliment Grundy
# Try:
# ComGrundy(10);
# [seq(ComGrundy(n),n=2..20)];
ComGrundy := proc(n) option remember;
local d,S,T,R;

if n=1 then
	return(0);
fi:

S := {};	

for d from 1 to floor(n/2) do
	T := [ComGrundy(d)$floor(n/d)
	,ComGrundy(n-d*floor(n/d))]; 
	S := S union {NimSum(T)};
od:

Mex(S);
end:


# Try:
# DivThrowRes(10);
# [seq(DivThrowRes(n),n=1..20)];
DivThrowRes := proc(n) option remember;
local d,S,T;

if n=0 or n=1 then
	return(0);
fi:

S := {};	
for d from 1 to n-1 do
	T := [DivThrowRes(d)$floor(n/d)]; 
	S := S union {NimSum(T)};
od:

Mex(S);
end:

# Try: 
# [seq(SolDivThrowRes(n),n=2..20)];
# [seq(Maliquant(n),n=2..20)];
SolDivThrowRes := proc(n) option remember;

if n=1 then
	0;
elif n=2 then
	1;
elif n mod 2 = 0 then
	SolDivThrowRes(n/2);	
else
	(n+1)/2;	
fi:
end:


# Try:
# [seq(ResThrowDiv(n),n=1..20)];
ResThrowDiv := proc(n) option remember;
local d,S;

if n=0 or n=1 then
	return(0);
fi:

S := {seq(ResThrowDiv(n-d*floor(n/d)),d=1..n-1)};
Mex(S);
end:


# Try:
# [seq(SolResThrowDiv(n),n=1..20)];
SolResThrowDiv := proc(n) option remember;
local k;

k :=0;

while n > 3*(2^k-1)+1 do
	k := k+1;
od:

k;
end:



###############
# Section 5:
# The Factoring Games
#
# MFactor(n), Factor1(n),
# SFactor(n),
#
###############

# Try: 
# MFactor(12);
# [seq(MFactor(n),n=2..30)];
# [seq(bigomega(n)-1,n=2..30)];
MFactor := proc(n) option remember;
local f,F,S;

if n=0 or n=1 then
	return(0);
fi:

S := {seq(NimSum([seq(MFactor(f),f in F)])
,F in (Factor1(n) minus {[n]}))};
Mex(S);
end:


# All factor of n (might be a little slow)
# Try:
# Factor1(12);
Factor1 := proc(n) option remember;
local d,A,S;

if n =1 then
	return([[]]);
fi:
	
S := {};	
for d in (divisors(n) minus {1}) do
	S := S union 
	{seq( sort([d,op(A)]),A in Factor1(n/d))};
od:

return(S);
end:


# Seems no pattern so far
# This can be quite random
# Try: 
# SFactor(12);
# [seq(SFactor(n),n=2..30)];
SFactor := proc(n) option remember;
local f,F,S;

if n=0 or n=1 then
	return(0);
fi:

S := {seq(NimSum([seq(SFactor(n-f),f in F)])
,F in Factor1(n))};

Mex(S);
end:



###############
# Section 6:
# Full Set Maliquot
#
# FullSet(n),
# SqFree(n),
#
###############

# Try: 
# [seq(Fullset(n),n=1..20)];
Fullset := proc(n) option remember;
local s,S,T;

if n=1 then
	return(0);
fi:	

S := divisors(n) minus {n}; 
#print(S);
T := {NimSum([seq( Fullset(s),s in S)])};
return(Mex(T));
end:



# Try: 
# [seq(SqFree(n),n=1..20)];
# A:=[seq(Fullset(n),n=1..150)]:
# B:=[seq(SqFree(n),n=1..150)];
# {op(A-B)};
SqFree := proc(n) option remember;
local s,S,T;

S := ifactors(n)[2];
T := {seq( s[2],s in S)};

if T = {1} then
	return(1);
else
	return(0);
fi:	

end: