subroutine longmul(x,y,z)
    parameter ichop=1000000000
    integer(8) x,y,xy,csx,csy,csxy,csz
    integer(4) x1(3),y1(3),z1(7)   
    integer(4) z(*),zz,xy1,xy2
!               csx,csy
    call chop(x,x1)     
    print *,"x1=",x1    
    call chop(y,y1)
    print *,"y1=",y1
    csx=cksum(x1,3)
    csy=cksum(y1,3)
    csxy=int(csx,8)*int(csy,8)
    print *,"csxy=",csxy
    csxy=mod(csxy,ichop-1)
    z1(1:7)=0
    do ix=1,3
       do iy=1,3
       xy=int(x1(ix),8)*y1(iy)
       print *,"ix,iy=",ix,iy," xy=",xy
       xy1=mod(xy,ichop)
       xy2=xy/ichop
       print *,"xy1,2=",xy1,xy2
       z1(ix+iy-1)  =z1(ix+iy-1)+xy1 
       if(z1(ix+iy-1)>=ichop)then
           z1(ix+iy)=z1(ix+iy)+z1(ix+iy-1)/ichop
           z1(ix+iy-1)=mod(z1(ix+iy-1),ichop)
       endif
       z1(ix+iy)    =z1(ix+iy)+xy2
       if(z1(ix+iy)>=ichop)then
           z1(ix+iy+1)=z1(ix+iy+1)+z1(ix+iy+1)/ichop
           z1(ix+iy)=mod(z1(ix+iy),ichop) 
       endif
       print 101,"z1=",z1
       enddo
    enddo
!    print 101,"z1=",z1
!  remove overflows
    do iz=1,6
        zz=z1(iz)
        z1(iz)  =mod(zz,ichop)
        z1(iz+1)=zz/ichop+z1(iz+1)
    enddo
    print 101,"z1=",z1
101 format(A,1x,8I11)
    csz=cksum(z1,6)
    print *,"csx =",csx
    print *,"csy =",csy
    print *,"csxy=",csxy
    print *,"csz =",csz
    contains
    function cksum(q,nq)
    integer(8) cksum
    integer(4) q(:nq)
    cksum=0
    do iq=1,nq
        cksum=cksum+q(iq)
    enddo
    cksum=mod(cksum,ichop-1)
    end function ! cksum
!
    subroutine chop(a,a1)
    integer(8) a,aa
    integer(4) a1(3)
    aa=a ! save
    a1(1)=mod(aa,ichop)
    aa=aa/ichop
    a1(2)=mod(aa,ichop)
    a1(3)=aa/ichop
    end subroutine
end