File:Point-doubling.R.svg: Difference between revisions

#! /usr/bin/R -f
# Point doubling illustrated

E <- y^2 ~ x^3 - x + 1
xrange <- seq(-1.5,1.5,length=4001)
point1 <- 1.0

stopifnot(E[[1]] == "~", E[[2]] == quote(y^2))

f <- call("sqrt", E[[3]])
g <- D(f,"x")
slope <- eval(g,list(x=point1))
intercept <- eval(f,list(x=point1)) - point1 * slope

r <- function(x){((intercept + slope * x)^2
        - eval(E[[3]],list(x=x)))/(x-point1)^2 }

point2 <- uniroot(r,lower=min(xrange),upper=max(xrange))$root

yrange <- eval(f,list(x=xrange))

xvals <- c(rev(xrange),xrange)
yvals <- c(rev(yrange),-yrange)

svg("point-doubling.R.svg")
plot(xvals,yvals,xlab="x",ylab="y",
        type="l",asp=1,lwd=3,lty="22",col="red",
        main=paste("point-doubling on the elliptic curve ",
                as.character(as.expression(E[[2]])),
                " = ", as.character(as.expression(E[[3]]))))
abline(0,0,col="blue")
abline(intercept,slope)
lines(c(point2,point2),
      c(-intercept-slope*point2,intercept+slope*point2),
      col=c("green","black"),lwd=2)
points(c(point1,point2,point2),
       c(eval(f,list(x=point1)),
         -intercept-slope*point2,
         intercept+slope*point2))
text(point1,eval(f,list(x=point1)),"P_1",pos=4)
text(point2,-intercept-slope*point2,"P_2 = 2*P_1",pos=4)