/* BC program powerdd */
/* outputs global variables z1 and z2 */
define powerdd(a,b,d,n){
/* (a+b*sqrt(d))^n=z1+z2*sqrt(d) */
auto x1,x2,temp,temp1,temp2,temp3,temp4,temp5,y

x1=a
x2=b
y=n
z1=1
z2=0
while(y){
   while(y%2==0){
      y=y/2
      temp=x1
      temp1=x2*x2
      temp2=x1*x1
      temp3=d*temp1
      x1=temp2+temp3
      temp4=temp*x2
      x2=2*temp4
   }
   y=y-1
   temp=z1
   temp1=z2*x2
   temp2=z1*x1
   temp3=d*temp1
   z1=temp2+temp3
   temp4=temp*x2
   temp5=z2*x1
   z2=temp4+temp5
}
}