diff --git a/DESCRIPTION b/DESCRIPTION index fc8db08..0130a02 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: spt -Version: 1.1 +Version: 1.1-1 Date: 2011-10-31 Title: Sierpinski Pedal Triangle Author: Bin Wang . @@ -8,6 +8,6 @@ Depends: R (>= 2.5.0), stats Suggests: MASS Description: This package collects algorithms related to Sierpinski pedal triangle (SPT). License: Unlimited -Packaged: 2011-11-02 18:08:16 UTC; bwang +Packaged: 2011-11-05 12:59:00 UTC; bwang Repository: CRAN -Date/Publication: 2011-11-02 19:03:20 +Date/Publication: 2011-11-05 18:07:33 diff --git a/MD5 b/MD5 index 44a05a8..c6c2512 100644 --- a/MD5 +++ b/MD5 @@ -1,10 +1,10 @@ -ca496720f39e42223e10686998a6b48c *DESCRIPTION +b50393939560bcc60c38ab94de897625 *DESCRIPTION d008ac6798fdfeba4fb80de7176f68c3 *LICENCE -a9e19ab5a327f27fdd1ed031501c862f *NAMESPACE -5399d2f938beb63f730286f21262c7c8 *PORTING -bf3fb080f3e8ea2900c656f018266596 *R/kernel.R +2e46d215627025af7eed4bdd944b1b6e *NAMESPACE +180388a86e4e077d32202e9b1676be21 *PORTING +87bcb8da10ceb1042466d6fa5ecb2d12 *R/kernel.R ee3784ad0fc5e1f68654fe3055a5481a *R/zzz.R -5f46927f84b4e710334ddf9942bafe34 *man/spt.Rd +4dbd0cb0e18afb5f8dc81a2cbf6b2b99 *man/spt.Rd 8290d2e9740414e315237f0d5d4024bb *src/Makevars b0683d1aff464dec94cf1417d7e6871f *src/init.c -b6b2af4c6307892d3fdb0c5cfe0f6b66 *src/kernel.c +26da8718bb7ac51499b97b06da296311 *src/kernel.c diff --git a/NAMESPACE b/NAMESPACE index f116093..16d6553 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,4 @@ -useDynLib(spt, .registration = TRUE, .fixes = "F_") +useDynLib(spt, .registration = TRUE, .fixes = ".F_") export(.sptConnect) ##export everything that does not start with a . exportPattern("^[^\\.]") diff --git a/PORTING b/PORTING index f69ee80..2a3aedf 100644 --- a/PORTING +++ b/PORTING @@ -10,4 +10,8 @@ option chaosgame=FALSE to determine how to draw the SPT. 4) iteration number will not be specified. Taking default value 10000 for "chaosgame" and 15 for regular drawing. 5) removing "main=NULL". It can be specified, but not prespecified. 6) "tol=0.00001" is removed, -and taken as default value. \ No newline at end of file +and taken as default value. + +11/4/11: In this version, (1) we fixed a bug in computing the +SPT dimention. (2) We restricted the expiration of the external C/Fortran +functions to the namespace. (3) output a class 'sped' for further study. diff --git a/R/kernel.R b/R/kernel.R index 01aa880..b3a05f4 100644 --- a/R/kernel.R +++ b/R/kernel.R @@ -1,10 +1,15 @@ -spt <- function(A,B,x0,y0,iter, +spt <- function(A,B,x0,y0,iter,plot=TRUE, spt=TRUE,chaosgame=FALSE,main=NULL,tol=0.0001,...){ if(missing(A)||missing(B)) stop("Sizes of at least two angles need to be specified.") if(A<=0 || B<=0 || A>=180 || B >=180 || A+B>=180) stop("Invaid angle size(s) specified: 090,NA, - round(.Fortran(F_SptDimFP, as.double(angles),as.double(0))[[2]],3)) - if(is.null(main)){ - main=paste(triname,"(",A1,",",B1,",",180-A1-B1,")",sep='') - if(is.na(sptdim)) - main=paste(main, ", Dim = ",sptdim,sep='') - } - if(maxA > 90 && spt){ - delta = abs(h * sin(maxA))*.25 - xmin = xmin - delta; - delta = (xmax-xmin-ymax+ymin)*.75; - ymin = ymin - delta; - ymax = ymax + delta; - } - - ## draw the initial triangle - plot(0,0, type='n', bty='n', xaxt='n',yaxt='n', - xlab='', ylab='', asp=1, main = main, - xlim=c(xmin,xmax), ylim=c(ymin,ymax)) - lines(c(ABC$xa,ABC$xb,ABC$xc,ABC$xa), - c(ABC$ya,ABC$yb,ABC$yc,ABC$ya), - lwd=2); - - if(chaosgame){ - if(missing(x0)) x0 = (xmax-xmin)*runif(1)+xmin; - if(missing(y0)) y0 = (ymax-ymin)*runif(1)+ymin; - if(missing(iter)) iter=10000; - - points(x0,y0,pch='*', col=5) - points(ABC$xa,ABC$ya,pch=20, col=2) - points(ABC$xb,ABC$yb,pch=20, col=3) - points(ABC$xc,ABC$yc,pch=20, col=4) - if(spt){ - .GameSPT(x0,y0,ABC,iter) - }else{ - .GameST(x0,y0,ABC,iter) + ymin = 0; ymax=h; + if(maxA>90) sptdim = NA + else if(maxA==90) sptdim = 2.0 + else sptdim = .Fortran(.F_SptDimFP, as.double(angles),as.double(0))[[2]] + ## sptdim = .Fortran(F_SptDimFP, as.double(angles),as.double(0))[[2]] + if(plot){ + if(is.null(main)){ + ## if(is.na(sptdim)||sptdim < 0) + cond1 = is.null(sptdim) + cond2 = is.na(sptdim) + if(cond1 || cond2) + main=paste(triname, ", Dim = NA",sep='') + else + main=paste(triname, ", Dim = ",as.character(format(sptdim,20)),sep='') } - }else{ - if(missing(iter)) iter = 20 - Tri0 = cbind(c(xa,xb,xc), c(ya,yb,yc)); - tol = tol * (xmax-xmin) * (ymax-ymin) - if(iter>0){ - if(!spt) .ptpedal(iter,Tri0,tol) - else if(angles[1]==90){ - .sptpedal2(iter,Tri0,tol) + if(maxA > 90 && spt){ + delta = abs(h * sin(maxA))*.25 + xmin = xmin - delta; + delta = (xmax-xmin-ymax+ymin)*.75; + ymin = ymin - delta; + ymax = ymax + delta; + } + + ## draw the initial triangle + plot(0,0, type='n', bty='n', xaxt='n',yaxt='n', + xlab='', ylab='', asp=1, main = main, + xlim=c(xmin,xmax), ylim=c(ymin,ymax)) + lines(c(ABC$xa,ABC$xb,ABC$xc,ABC$xa), + c(ABC$ya,ABC$yb,ABC$yc,ABC$ya), + lwd=2); + + if(chaosgame){ + if(missing(x0)) x0 = (xmax-xmin)*runif(1)+xmin; + if(missing(y0)) y0 = (ymax-ymin)*runif(1)+ymin; + if(missing(iter)) iter=10000; + + points(x0,y0,pch='*', col=5) + points(ABC$xa,ABC$ya,pch=20, col=2) + points(ABC$xb,ABC$yb,pch=20, col=3) + points(ABC$xc,ABC$yc,pch=20, col=4) + if(spt){ + .GameSPT(x0,y0,ABC,iter) }else{ - iter = ifelse(maxA>90, min(12,iter),iter) - .sptpedal3(iter,Tri0,tol) + .GameST(x0,y0,ABC,iter) + } + }else{ + if(missing(iter)) iter = 20 + Tri0 = cbind(c(xa,xb,xc), c(ya,yb,yc)); + tol = tol * (xmax-xmin) * (ymax-ymin) + if(iter>0){ + if(!spt) .ptpedal(iter,Tri0,tol) + else if(angles[1]==90){ + .sptpedal2(iter,Tri0,tol) + }else{ + iter = ifelse(maxA>90, min(12,iter),iter) + .sptpedal3(iter,Tri0,tol) + } } } } - invisible(round(sptdim,3)); + + if(spt){ + out = structure(list(angles, + data.name = triname, + Dim = sptdim, iter=iter + ), class = "spt") + cat("\n\nThe SPT dimension is", format(sptdim,10),"\n\n"); + }else out=NULL + + invisible(out) } ## Draw PT/SPT: tri1 is the initial triangle, n is the iteration number. .ptpedal <- function(n, tri1,tol){ if(n>1){ - tmp = .Fortran(F_PTChild, as.double(as.vector(tri1)), as.double(rep(0,6))); + tmp = .Fortran(.F_PTChild, as.double(as.vector(tri1)), as.double(rep(0,6))); result = matrix(tmp[[2]], ncol=2, nrow=3); polygon(result[,1],result[,2]) tri1.1 = rbind(tri1[3,], result[2,], result[1,]); @@ -89,7 +110,7 @@ spt <- function(A,B,x0,y0,iter, .sptpedal2 <- function(n, tri1,tol){ if(n>1){ - tmp = .Fortran(F_SPTChild2, as.double(as.vector(tri1)), as.double(rep(0,2))); + tmp = .Fortran(.F_SPTChild2, as.double(as.vector(tri1)), as.double(rep(0,2))); tmp = matrix(tmp[[2]],nrow=1); segments(tmp[1],tmp[2],tri1[1,1],tri1[1,2]); tri1.1 = rbind(tmp, tri1[3,], tri1[1,]); @@ -103,7 +124,7 @@ spt <- function(A,B,x0,y0,iter, .sptpedal3 <- function(n, tri1, tol){ if(n>1){ - tmp = .Fortran(F_SPTChild3, as.double(as.vector(tri1)), as.double(rep(0,6))); + tmp = .Fortran(.F_SPTChild3, as.double(as.vector(tri1)), as.double(rep(0,6))); result = matrix(tmp[[2]], ncol=2, nrow=3); polygon(result[,1],result[,2]) tri1.1 = rbind(tri1[3,], result[2,], result[1,]); @@ -119,7 +140,7 @@ spt <- function(A,B,x0,y0,iter, .Stop <- function(triangle){ axis = as.vector(triangle); - .Fortran(F_stri, as.double(axis), as.double(0))[[2]] + .Fortran(.F_stri, as.double(axis), as.double(0))[[2]] } ########################################################################## diff --git a/man/spt.Rd b/man/spt.Rd index 51e7eb7..aecdf58 100644 --- a/man/spt.Rd +++ b/man/spt.Rd @@ -3,23 +3,19 @@ \name{spt} \alias{spt} -\alias{F_PTChild} -\alias{F_SPTChild2} -\alias{F_SPTChild3} -\alias{F_SptDimFP} -\alias{F_stri} \title{Sierpinski Pedal Triangle and Sierpinski Triangle} \description{ To draw Sierpinski pedal triangles (SPTs) and Sierpinski triangles (PTs). } \usage{ - spt(A,B,x0,y0,iter, spt=TRUE,chaosgame=FALSE,main=NULL,tol=0.0001,...) + spt(A,B,x0,y0,iter,plot=TRUE, spt=TRUE,chaosgame=FALSE,main=NULL,tol=0.0001,...) } \arguments{ \item{A,B}{The degrees of two of the three angles of a triangle.} \item{x0,y0}{The initial position of the point if "chaosgame=TRUE".} \item{iter}{Iteration number in drawing the SPTs/PTs.} + \item{plot}{Draw the SPT/PT (default) or compute the SPT dimention alone.} \item{spt}{SPT ("spt=TRUE") or PT ("spt=FALSE").} \item{chaosgame}{To draw the SPTs/PTs using regular method or via a chaos game. } \item{main}{Specify the title of the plot.} @@ -53,8 +49,8 @@ spt(45,55,spt=FALSE) spt(45,55,spt=TRUE) - spt(120,30,iter=15) - spt(10,10,iter=10) + spt(120,30,iter=3) + spt(20,20,iter=7) } \keyword{stats} diff --git a/src/kernel.c b/src/kernel.c index 9b56c32..2e7dd3e 100644 --- a/src/kernel.c +++ b/src/kernel.c @@ -89,7 +89,7 @@ void PTChild(double *xy, double *result) } double fsptdim(double a, double b, double c, double d){ - return(pow(cos(a), d)+pow(cos(b),d)+pow(cos(c),d)-1.); + return(R_pow(cos(a), d)+R_pow(cos(b),d)+R_pow(cos(c),d)-1.); } void SptDimFP(double *angles, double *dim){ @@ -101,21 +101,26 @@ void SptDimFP(double *angles, double *dim){ Secant method is used. (a generalization of Secant method: Method of False position) */ - double a,b,c; - double x0=1., x1 = 3., xi, fa,fb,fc = 1.; //initial values; + double a,b,c, tol = 6.123234e-17; + double x0=log(3.0)/log(2.0), x1 = 2.0, xi, fa,fb,fc = 1.; //initial values; a = angles[0]*PI/180.; b = angles[1]*PI/180.; c = angles[2]*PI/180.; - fa = fsptdim(a,b,c,x0); fb= fsptdim(a,b,c,x1); - while(fabs(fc) >.0000000001 && fa * fb <0.){ - xi = (x1 * fa - x0 * fb)/(fa-fb); - fc = fsptdim(a,b,c,xi); - if(fa*fc<0){ - x1=xi;} - else{ - x0=xi; - } + //a suppose to be the larget angle + if(R_pow(cos(a),x0)<= tol) + dim[0] = 2.0; + else{ fa = fsptdim(a,b,c,x0); fb= fsptdim(a,b,c,x1); + while(fabs(fc) > tol){ + xi = 0.5 *(x0+x1); + fc = fsptdim(a,b,c,xi); + if(fc<0){ + x1=xi;} + else{ + x0=xi; + } + fa = fsptdim(a,b,c,x0); fb= fsptdim(a,b,c,x1); + } + dim[0] = xi; } - dim[0] = xi; }