######################################################################
## Proj5.txt Save this file as  Proj5.txt  to use it,                #
## stay in the                                                       #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read Proj5.txt <Enter>                            #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Abrar Almahmeed and Pablo Blanco as project in         #
## Dr.Doron Zeilberger 's class, Rutgers University , Spring 2026    #
##                                                                   # 
######################################################################

print():
print(`Written by Abrar Almahmeed and Pablo Blanco as a project in the Experimental Math class`):
print(`taught by Dr.Doron Zeilberger, Rutgers University , Spring 2026`):
print():
print(`To see the major procedures, type "Help()"`):
with(NumberTheory):

Help:=proc(): 
 if args=NULL then:
	print(`The procedures are:`):
	print(`CycleLenTable, FindCycle, lpf, SubpFseq, SubpFCycles`): 
	print(`Try: Help(procedurename)`):
 elif nargs=1 and args[1]=lpf then:
	print(`lpf(n): Given an integer n, returns its smallest prime factor (positive). If n is 1,0, or -1, lpf(n) returns 1.`):
	print(`Try:`):
	print(`n:=-17*19^6*23;`):
	print(`lpf(n);`):
 elif nargs=1 and args[1]=SubpFseq then:
 	print(`SubpFseq(a0,a1,N,c0,c1): Given integers a0,a1,N (and optional parameters c0,c1, which are both 1 by default),`):
	print(`returns the first N terms of the generalized subprime Fibonacci sequence with initial values a0 and a1, with`):
	print(`coefficients c0 and c1.`):
	print(`Try:`):
	print(`SubpFseq(0,1,30);`):

 elif nargs=1 and args[1]=SubpFCycles then:
	print(`SubpFCycles(M,N,c0,c1): Inputs positive integers M,N and optional parameters c0,c1 (=1 by default). `):
	print(`Outputs a table T where T[i,j]=FindCycle(SubpFSeq(i,j,N,c0,c1)) for all -M <= i, j <= M`):
	print(`Try:`):
	print(`T:=SubpFCycles(10,100):`):
	print(`T[6,5];`):
	print(`which corresponds to`):
	print(`S:=SubpFseq(6,5,100):`):
	print(`FindCycle(S);`):
	print(`T[-5,6];`):
	
 elif nargs=1 and args[1]=FindCycle then:
	print(`FindCycle(S): Given a sequence S (as a list) whose n-th term depends only on its previous two terms, return a pair [l,i].`):
	print(`l is the length of the cycle which the sequence repeats forever and i is the index of the cycle's first term in the sequence.`):
	print(`If the procedure does not find a cycle, it returns [FAIL, FAIL]. Try:`):
	print(`S:=SubpFseq(0,1,80);`):
	print(`FindCycle(S);`):

 elif nargs=1 and args[1]=CycleLenTable then:
	print(`CycleLenTable(T): Inputs a table T (such as that from SubpFCycles) where T[a0,a1]=[l, i] where a0 and a1 are initial conditions`):
	print(`for a sequence, l is the cycle length of the sequence and i is the index of the first term inside a cycle (see Help(FindCycle)).`):
	print(`Outputs a table P where P[l] returns the set of initial conditions [a0,a1] for which T[a0,a1][1]=l. Try:`):
	print(`T:=SubpFCycles(10,300):`):
	print(`P:=CycleLenTable(T):`):
	print(`Then, to see the cycle lengths that were found, use:`):
	print(`indices(P);`):
	print(`To see the initial conditions that achieve cycle length 136, use:`):
	print(`P[136];`):
 fi:
end:

############################
#lpf(n): is the smallest prime factor of the integer n.

lpf := proc(n) local i:
 if abs(n) <= 1 then: return 1: fi:

 return(PrimeFactors(n)[1]): # returns smallest prime factor of n
end:


############################
#B(x): is the a function that input an integer and output either 1 if x is a prime number and if |x|<=1,
#or the smallest prime factor of x otherwise.

B:=proc(x) 
 option remember:
 if isprime(abs(x)) then 
 	return 1:
 else 
	return lpf(x):
 fi:
end:


###########################
#SubpFseq(a0,a1,N,c1,c2): The subprime Fibonacci sequence a[n], input an initial values a0,a1, and parameters c1,c2 
#and an integer N (length of the sequence), output a sequence of the values of a[i], i=1..N. 

SubpFseq:=proc(a0,a1,N,c1:=1,c2:=1) local a,i:
 option remember:
 a:=array(0..N):
 a[0]:=a0:
 a[1]:=a1:
 for i from 2 to N do
 	a[i]:=(c1*a[i-1]+c2*a[i-2])/B(c1*a[i-1]+c2*a[i-2]):
 od:
 [seq(a[i],i=0..N)]:
end:

#############################
#SubpFCycles(M,N,c1,c2): generate SubpFseq for initial pairs(i,j) with -M<=i, j<=M
#output the table T .

SubpFCycles:=proc(M,N,c1:=1,c2:=1) local i,j,T,S:
 option remember:	
 T:=table([]):
 for i from -M to M do
    for j from -M to M do
     S:= SubpFseq(i,j,N,c1,c2):
     T[i,j]:= FindCycle(S): 
   od:
 od:
op(T):
end:

#############################
#FindCycle(S): A procedure to find cycles in the Subprime Fibonacci sequence 

FindCycle:=proc(S) local i,S1,j,idx,n :
 n:=nops(S):
 S1:=Array([0$n]):
 for i from 1 to n-1 do
 	j:=[S[i],S[i+1]]:
 	if not( member(j,S1,'idx')) then		
 		S1[i]:=j:
 	else
 		return([i-idx,idx]):
 	fi:

 od:

 return [FAIL,FAIL]:
end:

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

# inputs a table 
CycleLenTable:=proc(T) local i,j,E,S,l,e: 
 E:={indices(T)}:
 for e in E do
   l:=T[op(e)][1]: # the cycle length of the seq with starting cond. e=[i,j]
   # store e (index) into the table entry with index l
   if assigned(S[l]) then:
    S[l]:= S[l] union {e}:
   else:
    S[l]:={e}:
   fi:
od:

 op(S): # S is a table
end:


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

#CycleLengthSeq(M,N,c1:=1,c2:=1): input integer M and N output the set of all possible cycles lengths 

CycleLengthSeq:=proc(M,N,c1:=1,c2:=1) local T,i,j,L:
 T:=SubpFCycles(M,N,c1,c2):
 L:=[]:
 for i from 0 to M do
    for j from 0 to M do
       L:= [op(L),T[i,j][1]]:
    od:
 od:
 convert(L,set):
end:

#########################################################
#IsIncreasing(S): to check whether the sequence is increasing sequence

IsIncreasing:=proc(S) local i,n:
 n:=nops(S):
 for i from 2 to n do 
     if S[i]<S[i-1] then 
         return false:
     fi:
od:
return true:
end: