version 1.2

ChangeLog

@@ -1 +1,3 @@
-* Rev 1.0 Initial Release
+* Rev 1.0 Initial Release
+* Rev 1.1 Improved input validation, added convergents to output, added Pell's eqn. Minor bugfixes
+* Rev 1.2 Fixed bug; improved capabilities of 'cfrac2num' function

DESCRIPTION

@@ -1,14 +1,14 @@
 Package: contFracR
 Type: Package
 Title: Continued Fraction Generators and Evaluators
-Version: 1.0
-Date: 2022-03-24
+Version: 1.2
+Date: 2022-11-19
 Author: Carl Witthoft
-Maintainer: Carl Witthoft <carl@witthoft.com>
-Description: This library provides functions to convert numbers to continued fractions and back again.  
+Maintainer: Carl Witthoft <cellocgw@gmail.com>
+Description: Converts numbers to continued fractions and back again. A solver for Pell's Equation is provided.  The method for calculating roots in continued fraction form is provided without published attribution in such places as Professor Emeritus Jonathan Lubin, <http://www.math.brown.edu/jlubin/> and his post to StackOverflow, <https://math.stackexchange.com/questions/2215918> , or Professor Ron Knott, e.g., <https://r-knott.surrey.ac.uk/Fibonacci/cfINTRO.html> .
 License: LGPL-3
 Imports: Rmpfr, go2bigq, gmp, methods
 NeedsCompilation: no
-Packaged: 2022-04-18 20:04:24 UTC; cgw
+Packaged: 2022-12-03 02:53:11 UTC; cgw
 Repository: CRAN
-Date/Publication: 2022-04-20 08:32:29 UTC
+Date/Publication: 2022-12-03 08:00:02 UTC

MD5

@@ -1,19 +1,21 @@
-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
+2539525118cfbb1684063ac454c3df22 *ChangeLog
+65a4bf981f90c6a9dae2930f3222b66f *DESCRIPTION
+dbef16da3dfd743c2445cd2403e250da *NAMESPACE
+4d66295c10a4f33a8a3ef1932466f40d *R/cfrac.R
+59dcaee4d048da1279313a8ba9f9ddbb *R/fracHelpers.r
+dcaf23bc27f8b1fb2c25da92e191bdd0 *R/pell.r
+797e04a9c8753dd2a21ae9efabee417d *R/root2cfrac.R
+cee3cdb0dde176e55679e4aeaf0475b4 *R/sqrt2periodicCfrac.R
+6a955e3ded07a477db5a392e34f2ea78 *R/transcNumFracs.r
 87d945ce1990b788fce1cddc01cad3c0 *man/cfrac2latex.Rd
-92bbe514e6a99798fb77908b345913a9 *man/cfrac2num.Rd
+38104a7f2a9224c9e0fa22c7a7981051 *man/cfrac2num.Rd
 527d800f2060fb05f28065bfa7e525ed *man/cfrac2simple.Rd
 8e0d93bc392fa417507d9b5fd1db9715 *man/contFracR-package.Rd
 8bbd5e8f2ed5553392e5e7c681485345 *man/e2cfrac.Rd
-340dfad89cdbe601bb622111c4bb73d2 *man/num2cfrac.Rd
+f2d1775ec692815f18733fdaa512ff86 *man/num2cfrac.Rd
+b638a525538b6e45ff483ed8dedbb333 *man/pell.Rd
 2f2918ec98dbf29f9e16e32c9a7f7ecd *man/phi2cfrac.Rd
 fa38fde739e78ae14ddc8e8860c3ea4c *man/pi2cfrac.Rd
 69d16dc2fb911b0eaafa0cb1b09eebd2 *man/recipCfrac.Rd
 5fa33bd50c47a0250a0de7163dad1b8d *man/root2cfrac.Rd
-5633e3bb4f104ff6d7a4ffcc4ce0a0d0 *man/sqrt2periodicCfrac.Rd
+f3df8f21ce4c5d7d3d9670af365ca004 *man/sqrt2periodicCfrac.Rd

NAMESPACE

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

R/cfrac.R

@@ -5,7 +5,19 @@
 # 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
-
+
+#TODO revision Nov 2022: think about any other possible input validation to add
+#	 can I calc convergents along the way?  
+# to get convergents WHILE calculating the cfrac, use recursion formulas. bj are denoms with b0 the integer part,  aj are nums.
+
+# "starter values A[-1] = 1, B[-1] = 0
+
+#  A0 = b0 
+#  B0 = 1
+# Aj = bj*A(j-1) + aj*A[j-2]
+# Bj =  bj*B[n-1] + aj*B[n-2] 
+#
+#
 num2cfrac <-function(num, denom = 1, ...) {
 if(length(num) >1 ) {
 	warning('only using first element of num')
@@ -37,19 +49,32 @@ 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
+# intvec = as.bigz(NULL)  #NOT a R-vector, but NULL will be "replaced" with first value assigned.  
+jlev = 2   # start vector index, with offset for A[-1] stuff
+
+# initialize convergents. Notice that first elements are the "minus 1" indices
+# since all nums are 1 at this point, can substitute that in for a[j]
+#  browser()
+ intvec <- bzero <- floor(numdenom[1]/numdenom[2] )
+
+	num = numdenom[1] - (intvec * numdenom[2] )
+	numdenom = c(numdenom[2],num) # reciprocalling
+A <- c( as.bigz(1) , bzero )  #otherwise coercion screws up bigz value 
+B <- as.bigz(c(0 ,1)) 
 # 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){
+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
+	numdenom = c(numdenom[2],num) # reciprocalling
+# numdenom is "internal" calc. The next denom is latest intvec
+	A[jlev + 1] <- intvec[jlev] * A[jlev ] +  A[jlev -1]
+	B[jlev + 1] <- intvec[jlev] * B[jlev ] +  B[jlev -1]
 	 jlev = jlev + 1
 	}
 intvec <- intvec * thesign
-return(invisible(list(denom = intvec , numdenom = ND  )) )
+return(invisible(list(denom = intvec , numdenom = ND ,convA = A, convB = B )) )
 }
 
 
@@ -103,15 +128,29 @@ return(invisible(list(eqn=eqn, tex2copy = tex2copy, texeqn = texeqn)) )
 }
 
 
-#calculating numeric value of a continued fraction representation ,
+# 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. 
+# if the particular source is that value repeated 'numterms' times. 
+# Input 'nreps' is how many times to repeat a single-value denom input
+#
 # 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)
+# TODO: check for list inputs and either ban them or unlist(carefully)
+# Revision TODO Nov 2022.  Per Hans Borcher,  
+# cfrac2num(c(0, 1,-1,1)); x
+# Error in `/.bigq`(num[[kk]], spart) : division by zero
+
+#  One can certainly argue whether negative entries shall be allowed for continued
+#  fractions, but I think it will be better if the function will check the input instead of 
+#  the system.
+
+#  store every convergent , not just the final one (cgmp)  for people to play with  
+
+#BUGFIX: nreps fouls up when length of denom is bigger, as it will be for a periodic 
+#  contfrac resulting from a square root, for example. 
 cfrac2num  <- function(denom, num = 1, init = 0, nreps= 1, msgrate = NULL ) {
 if (length(denom) == 1){
 #  "exit strategy" when denoms == 0.  
@@ -125,7 +164,11 @@ if (length(denom) == 1){
         if (length(denom) == 1) {
                 return(invisible(list(cgmp = as.bigq(denom), mdec=mpfr(denom,10),denom = denom, num = num)))
                 }
-        }
+        } else {
+      # Now apply nreps to a putatively periodic denominator
+   			denom <- c(denom[1], rep(denom[-1],times=nreps))
+        }
+
 #  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')){
@@ -138,25 +181,30 @@ if (!is(denom[[1]],'bigq')) {
 }else{
 	Kq <- denom
 }
-# just this?
-#spart <- init
-spart <- ( tail(Kq, 1) )[[1]] + init 
+# initalize
+spart <- ( tail(Kq, 1) )[[1]] + init 
+# new: save all convergents
+#  NOPE:  convergents have to be calc'd from start, not end 
+#  conv <- spart
 # 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
+        spart <- Kq[kk] + num[kk] /spart  #1 / spart
+ #       conv <- c(conv, spart)
         }
 
 	} else {
 		for (kk in loopvals) {
-		spart <- Kq[[kk]] + num[[kk]] /spart  #1 / spart
+			spart <- Kq[[kk]] + num[[kk]] /spart  #1 / spart
+#		conv <- c(conv, spart)
+        }
 		}
-	}
-mdec <- .bigq2mpfr(spart) #no reason to apply nonoptimal nbr bits
+
+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

@@ -3,6 +3,9 @@
 #   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 
+
+#TODO  even tho' these are primarily helper funcs, add input validation 
+
 recipCfrac <- function( denom, num = 1, ...){
 # ensure proper lengths
 num <- rep(num,length=length(denom) -1)

R/pell.r

@@ -0,0 +1,71 @@
+# pell  x^2 - d * y^2 = 1 , d a non-square integer
+#wiki info:
+# Let   h_i / k_i denote the sequence of convergents to the regular continued fraction for sqrt(n). This sequence is unique. Then the pair  solving Pell's equation and minimizing x satisfies x1 = hi and y1 = ki for some i. This pair is called the fundamental solution. Thus, the fundamental solution may be found by performing the continued fraction expansion and testing each successive convergent until a solution to Pell's equation is found.
+# h_i are nums, k_i denoms of each successive layer in contfrac  
+
+# test example  for d ==7
+# The sequence of convergents for the square root of seven are
+# h/k (convergent)	h^2   7k^2 (Pell-type approximation)
+# 2/1			 3
+# 3/1			+2
+# 5/2			 3
+# 8/3			+1
+# so 8 / 3 is the one we want
+
+
+# rev after 1.2: more input validation, done better
+#  generate more convergents if no "win" found.  
+#  allow A,B inputs to continue  a previous run
+#  make sure $repeated == TRUE before running with the results of
+#  sqrt2periodicCfrac
+
+pell <- function (d , maxrep = 10, A = NULL, B = NULL){
+	# 	wins <- NULL  # NO!  gmp really hates NULL things
+	wins <- as.bigz( matrix(0, nrow= 1 ,ncol = 2)) 
+	d <- as.bigz(d[1])
+	# get the greatest square less than x; make sure mpfr has reasonable precision
+	dm  <- .bigz2mpfr(d, max(5 * ceiling(log10(d)), 500))
+	#validate - check that sqrt(d) is not an integer
+	if (floor(sqrt(dm)) == sqrt(dm)) {
+		stop('Pell is pointless when d is a square integer')
+	}
+	foo <- sqrt2periodicCfrac(d)
+	# check for success
+	nterm = 100  # default is 50
+	while (!foo$repeated) {
+		foo <- sqrt2periodicCfrac(d, nterms = nterm)
+		nterm = nterm *2 
+	}
+	repden <- foo$denom[-1]
+	replen <- length(repden)
+	jrep = 0
+	if (is.null(A)) {
+	# build initial stuff
+		# conv{A,B}[1] is the "starter" value A[-1], so ignore
+		for (jc in 2:length(foo$convA)) {
+			if (foo$convA[jc]^2 -d * foo$convB[jc]^2 -1 == 0) {
+				wins <-	rbind(wins,c(foo$convA[jc],foo$convB[jc]) )
+				}
+		}
+		jconv = jc+1
+		A = foo$convA
+		B = foo$convB
+		# end of if isnull convA
+	}  else {
+		A <- A
+		B <- B
+		jconv = length(A) + 1
+	}
+#   (wins < 4 because of dummy first row)
+
+	while(length(wins ) < 4 && jconv <= maxrep * replen) {
+	# need more convergents. use the repeated section of foo$denom	
+		A[jconv] <- repden[jrep + 1] *A[jconv-1] + 1 * A[jconv-2]
+		B[jconv] <- repden[jrep + 1] * B[jconv-1] + 1 * B[jconv-2]
+		if (A[jconv]^2 - d * B[jconv]^2 -1 == 0) wins <- rbind(wins,c(A[jconv],B[jconv]) )
+		# increment jrep mod(length(repden))
+		jrep <- (jrep+1) %% replen
+		jconv = jconv + 1 
+	}
+	return(invisible(list(d = d, wins=wins[-1,], cfracdata = foo, A= A, B= B)))
+}

R/root2cfrac.R

@@ -1,9 +1,13 @@
 # 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,.....
+# nums are my, (n-m)b, (n+m)b, (2n-m)b, (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
 
+# Calculating  x^(m/n)  
+
+# Nov 2022: improve input validation work.  
+
 root2cfrac <- function(x, n, m = 1, nterms = 10, ...){	
 #sanitize
 n <- floor(n[1])
@@ -12,16 +16,21 @@ m <- floor(m[1])
 # 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)
+if (x < 0 ) {
+	x <- abs(x)
+	warning('negative x not allowed. Taking abs value')
+}
+if ( m * n < 0)  {
+	stop('negative exponents not allowed. Consider entering 1/x in bigq form')
+}
+if (n == 0) stop('infinite exponent not allowed (n > 0 required)')
 # want  a^n + b = x 
-# gmp doesn't allow fraction powers
+# gmp doesn't allow fraction powers (but .bigq2mpfr "guarantees" sufficient digits are used to be exact)
 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()
+denom <- as.bigz(NULL)  
+num <- as.bigz(NULL)  
 numericdenom <- vector()
 numericnum <- vector()
 denom <- c(denom,pow.bigq(a, m))