// CDC 6600 and IBM Stretch multiply. <- multiply.html
#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#define p1 0
typedef long int L;
typedef unsigned int ui;
typedef struct{L hi; L lo;} rp;
#define C const
char p = 1;
void prp(char * m, rp x){
  if(p1) printf("2*0x%016lx+0x%016lx = %s\n", x.hi, x.lo, m);}
long unsigned int dump;
void discard(ui x){assert(x<22); if(p1) printf("D%d ", x);
   dump = (dump>>3) + ((L)x<<(64-5));}
rp cook(L x, L y, L z){
  if(p1) printf("cook %016lx %016lx %016lx\n", x, y, z);
  return (rp){x&y|x&z|y&z, x^y^z};}
typedef unsigned long int ul;
typedef struct{ul q; ul r;} res;
res divq(res, ul);
res mulq(ul, ul);
res addq(res, res);
void prr(char * m, res x){
   printf("[%016lx %016lx] = %s\n", x.q, x.r, m);}
res llsrs(res v, int n){return n<64?(res){(L)v.q>>n, (n?(v.q<<(64-n)):0)|(v.r>>n)}
          :(res){(L)v.q>>63, (L)v.q>>(n-64)};}
res llls(res v, int n){return n<64?(res){(v.q<<n)|(n?(v.r>>64-n):0),v.r<<n}
     :(res){v.r<<(n-64), 0};}
int eqr(res a, res b){return a.q==b.q && a.r==b.r;}
res mul (L C a, L b){L C th = a+(a<<1); ui carry = 0; dump = 0;
  res C pr = mulq(a, b);
  if(p) printf("%016lx*%016lx=%016lx%016lx\n", a, b, pr.q, pr.r);
  L sel(ui n){assert(n<8); n += carry; carry = (n>3); if(p1) printf("n=%d, ", n);
    return carry?~(n>6?(n==7?a:0):(n==5?th:(a<<(n==4?2:1)))):
                 (n>2?th:!n?0:(n==1?a:a<<1));}
  {rp sum = {0, 0}; int k;
  res v1(int k){return addq(
         llsrs((res){(sum.hi<<1)+sum.lo, 0}, 64-k),
         llsrs((res){0, dump}, 62-k));}
  void rpt(char * m){printf("rpt %s k=%d:\n", m, k); prp("sum", sum);
         printf("dump=%016lx c=%d\n", dump, carry); prr("v1", v1(k+3));}
  void ck(int k){assert(labs(sum.hi)<(1L<<62));
      res t = addq(v1(k), llls(mulq(a, b+carry), k));
      if(p1) {
      prr("llsrs((res){(sum.hi<<1)+sum.lo, 0}, 64-k)",
           llsrs((res){(sum.hi<<1)+sum.lo, 0}, 64-k));
      prr("llsrs((res){0, dump}, 62-k)", llsrs((res){0, dump}, 62-k));
      prr("llls(mulq(a, b+carry), k)", llls(mulq(a, b+carry), k));
      prr("t", t); printf("carry=%d\n", carry);}
      assert (eqr(t, pr));}
  for(k=0; k<48; k+=3) {if(p1) printf("---- k=%d\n", k); if(0) ck(k);
   assert(labs(sum.hi)< (1l<<62));
   if(p1) rpt("a");
   {L C z = sel(b&7); if(p1) printf("z=%016lx\n", z);
      sum = cook((sum.hi<<1)|carry, sum.lo, z);}
   b>>=3; if(p1) rpt("b");
   if(0) {res tx = v1(48);
     discard(((7&sum.hi)<<1)+(7&sum.lo));
     sum = (rp){sum.hi>>3, sum.lo>>3}; 
       if(p1) {prp("sumq", sum); printf("dmp=%016lx\n", dump);}
     {res ty = v1(51); if(p1){prr("tx", tx); prr("ty", ty);} assert(eqr(tx, ty));}}
   else {discard(((7&sum.hi)<<1)+(7&sum.lo));
     sum = (rp){sum.hi>>3, sum.lo>>3};}
   if(p1) rpt("c");
   } ck(48);
   if(carry) sum = cook(sum.hi<<1, sum.lo, a);
   if(p1) printf("end carry = %d lo=%016lx\n", carry, sum.lo);
   return v1(48);}}

#include <stdlib.h>
void tst(L a, L b){res hp = mul(a, b); res qp = mulq(a, b);
  if(!eqr(qp, hp)) {printf("Bad!\n");
    prr("hp", hp); prr("qp", qp); exit(0);}}
L rL(){return ((L)random()<<17)^random();}
int main(){{int j=64; tst(0, 7237578L); while(j--) tst (1, j<<3);}
{tst(0, 234587237578L); tst(1, 4); tst(0, 63);
    tst(1L<<47, 1L<<2); tst(1L<<4, 7L<<44); tst(1L<<47, 7L<<44);}
 if(0) {int j = 1000000000; while(j--){
     p = !(((1<<18)-1)&random()); tst(rL(), rL());}}
 {int j = 48; while(j--) {int k = 48; while(k--) tst(1L<<j, 1L<<k);}}
 {int k=128; while(k--) prr("f", llls((res){64|(1L<<63), 65537|(1L<<63)}, k));}
printf("Good\n");
return 0;}