version 1.0

ChangeLog

@@ -0,0 +1 @@
+* Rev 1.0 Initial Release

DESCRIPTION

@@ -0,0 +1,14 @@
+Package: contFracR
+Type: Package
+Title: Continued Fraction Generators and Evaluators
+Version: 1.0
+Date: 2022-03-24
+Author: Carl Witthoft
+Maintainer: Carl Witthoft <carl@witthoft.com>
+Description: This library provides functions to convert numbers to continued fractions and back again.  
+License: LGPL-3
+Imports: Rmpfr, go2bigq, gmp, methods
+NeedsCompilation: no
+Packaged: 2022-04-18 20:04:24 UTC; cgw
+Repository: CRAN
+Date/Publication: 2022-04-20 08:32:29 UTC

MD5

@@ -0,0 +1,19 @@
+205f1ac10db99343d77276c80e7ea1d7 *ChangeLog
+0b7f30f752f4b6ee51abca4e7514416b *DESCRIPTION
+d76b58ff2e101fba83cafc3fa4c11df7 *NAMESPACE
+a17a81180a3185baf4ca4f358a51901a *R/cfrac.R
+d33af6dae1cb708a7fcb8f252a71ab6e *R/fracHelpers.r
+f25365c7d5ff3826445f4f4808f4505c *R/root2cfrac.R
+19bef864ab509e4ae4da8d994df3ff57 *R/sqrt2periodicCfrac.r
+5e990ff1b6de9755c37fd4676fd3fde9 *R/transcNumFracs.r
+87d945ce1990b788fce1cddc01cad3c0 *man/cfrac2latex.Rd
+92bbe514e6a99798fb77908b345913a9 *man/cfrac2num.Rd
+527d800f2060fb05f28065bfa7e525ed *man/cfrac2simple.Rd
+8e0d93bc392fa417507d9b5fd1db9715 *man/contFracR-package.Rd
+8bbd5e8f2ed5553392e5e7c681485345 *man/e2cfrac.Rd
+340dfad89cdbe601bb622111c4bb73d2 *man/num2cfrac.Rd
+2f2918ec98dbf29f9e16e32c9a7f7ecd *man/phi2cfrac.Rd
+fa38fde739e78ae14ddc8e8860c3ea4c *man/pi2cfrac.Rd
+69d16dc2fb911b0eaafa0cb1b09eebd2 *man/recipCfrac.Rd
+5fa33bd50c47a0250a0de7163dad1b8d *man/root2cfrac.Rd
+5633e3bb4f104ff6d7a4ffcc4ce0a0d0 *man/sqrt2periodicCfrac.Rd

NAMESPACE

@@ -0,0 +1,6 @@
+exportPattern("^[[:alpha:]]+")
+importFrom("gmp", "as.bigz", "as.bigq","pow.bigq", "is.bigz", "numerator" , "denominator")
+importFrom("Rmpfr", "mpfr", ".bigq2mpfr", ".mpfr2bigz")
+importFrom("go2bigq", "go2bigq")
+importFrom("methods", "is")
+importFrom("utils", "tail")

R/cfrac.R

@@ -0,0 +1,162 @@
+# requires go2bigq; I do want that to be a separate package, though. 
+# Function  num2cfrac
+
+#Inputs can be numeric, bigz, or reasonable character strings representing numbers. Or feed a bigq to "numerator"  
+# char strings can include '-' or '+' leading sign and 'eX' exponential notation
+#  you can calculate irrationals to arbitrary precision by setting the inputs to  num = (irat_num)*1eN , denom = 1eN for N some nice integer
+# mpfr values not allowed because teensy rounding errors, e.g. "3.2400000000001" make  a hash of the fraction
+
+num2cfrac <-function(num, denom = 1, ...) {
+if(length(num) >1 ) {
+	warning('only using first element of num')
+	num <- num[1]
+	}
+if(length(denom) >1) {
+	warning('only using first element of denom')
+	denom <- denom[1]
+	}
+# new sortorific bit of code to catch bigq
+if (is(num,'bigq') ) {
+	denom <- denominator(num)
+	num <- numerator(num)
+}
+numclass <- class(num)
+denomclass <- class(denom) 
+classlist <- c('numeric','character', 'bigz')
+if (!(numclass %in% classlist) || !(denomclass %in% classlist)){
+	stop('Inputs must be bigz,bigq, numeric, or character')
+	}
+# do go2bigq on both and 'flip' the division in gmp space
+if(!is.bigz(num)) num <- go2bigq(num)
+if(!is.bigz(denom)) denom <- go2bigq(denom)
+# ok, now can reduce from bigz or bigq form. In addition, this step
+# produces a bigq in reduced form. 
+thebigq <- num[[1]]/denom[[1]]
+numdenom  <- c( abs(numerator(thebigq)), denominator(thebigq))
+thesign <- sign(thebigq)
+ND <- numdenom 
+if (numdenom[2] == 0 || numdenom[1] == 0 ) stop('zeros not allowed')
+if (numdenom[1] == 1) return (invisible(list(intvec = c(0,denom), numdenom = ND) ) )
+intvec = as.bigz(NULL)  #NOT a R-vector, but NULL will be "replaced" with first value assigned.  
+jlev = 1 # start vector index
+# numerator == 0 means there was a common factor; this is easiest way to catch it.
+# rats: is.element, %in% don't like bigq
+#while (!is.element(numdenom[1], c(0,1)) && numdenom[2] != 0 ) {
+while(numdenom[1] !=1 && numdenom[1] != 0&&numdenom[2]!= 0){
+	intvec = c(intvec,  floor(numdenom[1]/numdenom[2] ) )
+	num = numdenom[1] - (intvec[jlev] * numdenom[2] )
+	numdenom = c(numdenom[2],num) # reciprocalling
+	 jlev = jlev + 1
+	}
+intvec <- intvec * thesign
+return(invisible(list(denom = intvec , numdenom = ND  )) )
+}
+
+
+# FUNCTION cf2latex -- build the latex formula and the inline equation
+# the input 'denomvals' is the 'intvec' output of num2cfrac, or any other legal Continued Fraction form vector of integers. Similar for numvals if setting up a non-simple continued fraction. 
+# Beware: R interprets backslashes in a char string as escapers in most cases. Example:
+ # foo <- "this string has \f backslashes \\ and stuff"
+# > unlist(strsplit(foo,''))
+ # [1] "t"  "h"  "i"  "s"  " "  "s"  "t"  "r"  "i"  "n"  "g"  " "  "h"  "a"  "s" 
+# [16] " "  "\f" " "  "b"  "a"  "c"  "k"  "s"  "l"  "a"  "s"  "h"  "e"  "s"  " " 
+# [31] "\\" " "  "a"  "n"  "d"  " "  "s"  "t"  "u"  "f"  "f"   
+
+# So strings must be submitted in  double-escaped format
+
+cfrac2latex <-function(denomvals, numvals = 1, denrepeats = 1, doellipsis = FALSE,  ...) {
+vrep <- denomvals[-1] # don't repeat the lead integer
+denomvals <- c(denomvals[1], rep(vrep, times= denrepeats) )
+# Do some legerdemain to get length(numvals) to match length(denomvals)
+# and some new variable PM , pm
+numvals <- rep(numvals, length = length(denomvals)-1)
+PM = c(' - ',' + ',' + ')
+pm <- PM[sign(numvals) + 2] 
+# now remove signs from numeric numvals
+numvals <- abs(numvals)
+eqn <- paste0(denomvals[1], pm[1])
+#  stop the loop one short to avoid extra "+"
+for (jj in 2:(length(denomvals) -1) ) {
+#first denomval is integer, hence teh index shift for nums
+	eqn <- paste0(eqn, numvals[jj-1],'/(',denomvals[jj], pm[jj], collapse='')
+	}
+#now the last integer, no plus sign
+if (doellipsis){
+	eqn <- paste0(eqn, numvals[jj],'/', denomvals[jj+1],' ...',collapse='')
+	}else{
+		eqn <- paste0(eqn, numvals[jj],'/', denomvals[jj+1],collapse='')
+	}
+#finish the inline version 
+clospar <-paste0(rep(')',jj-1),collapse='')
+eqn <- paste0(eqn, clospar,collapse='' )
+# Now a Q&D conversion to LaTeX, skipping markdown "$"
+# tex2copy is useful only for copy/paste, as the \ will be parsed by console 
+# replace '1' with num value.  ([0-9]{1,}) 
+tex2copy <- gsub('([0-9]{1,})/[(]', '\\\frac{\\1}{' , eqn)
+tex2copy <- gsub('[)]', '}', tex2copy)
+#cleanup
+tex2copy <- gsub('([0-9]{1,})/', '\\\frac{\\1}{', tex2copy)
+tex2copy <- paste0(tex2copy,'}',collapse='')
+# convert the single actual character "\f" to two characters "\\" and "f" 
+texeqn <- gsub('\\f', '\\\\\\f', tex2copy)
+return(invisible(list(eqn=eqn, tex2copy = tex2copy, texeqn = texeqn)) )
+}
+
+
+#calculating numeric value of a continued fraction representation ,
+# where K is either the Continued Fraction form sequence of integers, or a single value
+# if the particular source is that value repeated 'numterms' times. 
+# outputs
+#	cgmp is the gmp num&denom
+#	cmpfr is the mpfr conversion of cgmp
+#	nums, denoms returns the input vectors used
+
+#TODO: check for list inputs and either ban them or unlist(carefully)
+cfrac2num  <- function(denom, num = 1, init = 0, nreps= 1, msgrate = NULL ) {
+if (length(denom) == 1){
+#  "exit strategy" when denoms == 0.  
+# Zero is only legal in the first position of  the contfrac sequence anyway. 
+        if( denom == 0 ) {
+                return(invisible(list(cgmp = as.bigq(0), mdec=mpfr(0,10),denom = denom, num = num)))
+                }
+# bug-ish in tail(gmp), so check nterms value
+        denom <- rep(denom, max(nreps,1) )
+# One More Thing(TM) to check: if all we have is the 'a' term:
+        if (length(denom) == 1) {
+                return(invisible(list(cgmp = as.bigq(denom), mdec=mpfr(denom,10),denom = denom, num = num)))
+                }
+        }
+#  expand num to be  length(denom) -1  if it isn't already (first denom is the integer)
+#  and if length(num) == 1 , repeat it by nterms-1
+if( !is(num[[1]],'bigq')){
+	num <- as.bigq(rep(num,length=length(denom)-1 ) )
+}else{
+	num <- rep(num,length = length(denom)-1)
+}
+if (!is(denom[[1]],'bigq')) {
+	Kq <- as.bigq(denom)
+}else{
+	Kq <- denom
+}
+# just this?
+#spart <- init
+spart <- ( tail(Kq, 1) )[[1]] + init 
+# add a status msg. Also notice reverse order as befits unwinding a cont frac
+loopvals <- (length(denom) - 1):1
+if (!is.null(msgrate)){
+	msgrate <- floor(msgrate)
+	msgnums <- loopvals[!(loopvals %% msgrate)]
+	for (kk in loopvals) {
+		if ( kk %in% msgnums) message('kk is now ', kk)
+        spart <- Kq[kk] + num[kk] /spart  #1 / spart
+        }
+
+	} else {
+		for (kk in loopvals) {
+		spart <- Kq[[kk]] + num[[kk]] /spart  #1 / spart
+		}
+	}
+mdec <- .bigq2mpfr(spart) #no reason to apply nonoptimal nbr bits
+fracpart <- spart - Kq[[kk]]
+return(invisible(list(cgmp = spart, cmpfr = mdec, fracpart = fracpart, denom = denom, num = num)))
+}

R/fracHelpers.r

@@ -0,0 +1,21 @@
+# simple helpers
+#  see https://math.stackexchange.com/questions/86043/
+#   move original b1 inside the a1....an sequence and setting the new b1 == 0
+# and for numerator, prepend a "1" 
+# generate terms for 1/X 
+recipCfrac <- function( denom, num = 1, ...){
+# ensure proper lengths
+num <- rep(num,length=length(denom) -1)
+denom <- c(0,denom)
+num <- c(1,num)
+return(invisible( list(denom=denom, num = num)) )
+}
+
+# generate equivalent simple denom for any num,denom vectors
+cfrac2simple <- function(denom, num=1, ...){
+foo <- cfrac2num(denom = denom, num = num)$cgmp
+thenum <- numerator(foo)
+thedenom <- denominator(foo)
+denom <- num2cfrac(thenum, thedenom)$denom
+return(invisible(list(intfrac = foo, denom = denom)) )
+}

R/root2cfrac.R

@@ -0,0 +1,42 @@
+# wikipedia  suggests that making x and y big (  nth root of z^m  with arbitrary z  = x^n + y . Maybe choose largest possible x? In simplest
+# setup,let m == 1  . Then find the largest n-th cube less than  z.   
+# lead term is a^m 
+# nums are my, (n-m)b, (n+m)b, (2n-m)y, (2n+m)b, 3,3,4,4,.....
+# denoms are  1na^(n-m), 2a^m, 3na^(n-m), 2a^m, 5na^(n-m), alt 2term and 7,9,11,...term
+
+root2cfrac <- function(x, n, m = 1, nterms = 10, ...){	
+#sanitize
+n <- floor(n[1])
+m <- floor(m[1])
+#TODO: check to see if can replace x, m with  some root K(x) and m*K . Probably 
+# doesn't guarantee faster convergence. 
+if(length(x) > 1) warning('only first element of x will be used') 
+x <- as.bigq(x[1]) #(just in case x is not an integer)
+# want  a^n + b = x 
+# gmp doesn't allow fraction powers
+xm <- .bigq2mpfr(x)
+a <- go2bigq(floor( xm^(1/n)))[[1]]
+b = x - pow.bigq(a, n )
+#TODO: use denom <- as.bigz(NULL) etc. to initialize, thus returning
+# bigz vectors instead of messy lists
+# and see about simplifying the converstion to numerics
+denom <- as.bigz(NULL)  #list()
+num <- as.bigz(NULL)  #list()
+numericdenom <- vector()
+numericnum <- vector()
+denom <- c(denom,pow.bigq(a, m))
+num <- c(num,m * b)
+jnum = 1  #indexer
+for (jt in seq(1,nterms, by = 2)) {
+#that skipseq works  but take care
+	num <- c(num, (jnum *n-m)*b,(jnum *n + m)*b )
+	denom <- c(denom,jt*n* pow.bigq(a, (n-m)), 2*pow.bigq(a,m) )
+# bump num index
+	jnum = jnum + 1
+}
+#finish off denom  but watch for oddeven jt as it always ends with odd regardless of nterms
+denom <- c(denom, (jt+3) *n* pow.bigq(a, (n-m)))
+numericdenom <- as.numeric(denom)
+numericnum <- as.numeric(num)
+return(invisible(list(num = num, denom=denom, numericnum = numericnum, numericdenom = numericdenom, x=x, n= n, m= m)) ) 
+}

R/sqrt2periodicCfrac.r

@@ -0,0 +1,77 @@
+#  algebraic calculation of square root simple continued fracs
+# This function finds the periodic expansion. Use root2cfrac for nonperiodic but
+# faster convergence. 
+# input must be int or bigz, or a charstring which will be floor-ed to bigz
+# Default limit on nterms just to avoid monstrous output
+#  sqrt(n/m) =  sqrt(n*m)/m ,so can do this for any rational number. 
+# irrationals we can only approximate anyway. 
+# reference:
+#  rt(67) is 8;   5  2  1  1  7  1  1  2  5 16
+# 
+#output definition:
+#	repeated: TRUE if repeat sequence found. FALSE if exceed nterms w/o repeat
+#	input: echo the inputs
+
+sqrt2periodicCfrac <- function (num, denom = 1,  nterms = 50, ...){
+if ( length(num) >1 || length(denom) >1 ){
+	warning('Only first element of num or denom will be used')
+	num <- num[1]
+	denom <- denom[1]
+	}
+input <- c(num,denom) #save to dump into output
+#convert to a single bigz and save the denom to re-apply at end of calculation
+#first, reduce to minimum fraction.
+frac <- abs(as.bigz(num) /as.bigz(denom))
+num <- numerator(frac) * denominator(frac)
+# This is a final multiplier to get the desired sqrt value
+thedivisor <- denominator(frac)
+# this will store the cont frac denoms 
+#TODO: replact with initialization thedenom <- as.bigz(NULL)
+# and see about simplifying the converstion to numerics
+thedenom <- as.bigz(NULL)  #list()
+numericdenom <- vector()
+#follow math.stackexchange question 2215918
+# get the greatest square less than x; make sure mpfr has reasonable precision
+kbits <- max(5 * ceiling(log10(num)), 500)  
+mpnum <- mpfr(num,kbits)
+#start the engine
+m <- floor(sqrt(mpnum))  # the first term in standard contfrac notation
+thedenom <- c(thedenom, .mpfr2bigz(m) )#this is the "integer" part of contfrac form  
+pq <- c(m,1)  # this is mpfrs
+newq <- (num - pq[1]^2) / pq[2] 
+#if newq is zero, terminate, perfect square was input
+if (newq == 0 ) {
+	return(invisible(list(thefrac = thedenom, numericfrac = numericdenom,  thedivisor=thedivisor, repeated='perfect_square', input=input) ))
+}
+mpdenom <- floor( (pq[1] + m)/newq)
+thedenom <- c(thedenom, .mpfr2bigz( (pq[1] + m)/newq) )
+# numericdenom[2] <- as.integer(thedenom[[2]])
+newp <-  mpdenom * newq - pq[[1] ]
+pq <- c(newp,newq)
+#loop on that for either nterms or until  mpdenom == 2*m 
+# Since we can only feed rationals to the function, and we convert
+# to an integer value, there can't be a 'preamble' 
+idx <- 2
+#initialize
+repeated = FALSE # indicating no repeat sequence found
+while ( (idx <= nterms) && !repeated  ) {
+	newq <- (num - pq[[1]]^2) / pq[[2]] 
+	mpdenom <- floor( (pq[1] + m)/newq)
+	newp <-  mpdenom *  newq - pq[[1]]  
+	pq <- c(newp,newq)
+	thedenom <- c(thedenom, .mpfr2bigz(mpdenom) )	# watch for giant integers... 
+	if (mpdenom == 2*m) {repeated = TRUE}
+	idx <- idx + 1
+}
+# now stick that 1/m term back in
+thedenom[1] <- thedenom[1]/thedivisor
+numericdenom <- as.numeric(thedenom)
+numericnums <- vector()
+numericnums[1] <- 1/as.integer(thedivisor)
+numericnums[2:(length(numericdenom) -1)] <- 1
+thenum <- as.bigq(numericnums)
+thenum[1] <- 1/thedivisor  #removes possible precision foulup 
+return(invisible(list(denom = thedenom, num = thenum, numericdenom = numericdenom, numericnum =numericnums,  repeated=repeated, input=input) ))
+
+}
+