Karatsuba

(* Decomposition of the numbers into arrays *)
let decomposition10 a=
  let rec log10 a=
    if a<10 then 1 else 1+(log10 (a/10))
  in
  let lga=log10 a in
  let resultat=Array.create (lga) 0 in
  let rec decomp_aux a i=
    if a=0 then () else
      begin
        let b=a/10 in
          resultat.(i)<-(a-10*b);
          decomp_aux b (i+1)
      end
  in
    decomp_aux a 0;
    resultat;;

(* Power function, defined recursively *)
let rec ( ** ) nb pb=
  if pb=0 then 1 else nb*(nb**(pb-1));;

(* Gives both arrays the same length *)
let equilibre vecta vectb=
  let pluslong=max (Array.length vecta) (Array.length vectb) in
  let m=Array.make_matrix 2 pluslong 0 in
    for i=0 to pluslong-1 do
      begin
        try m.(0).(i)<-vecta.(i) with _->m.(0).(i)<-0;
      end;
      begin 
        try m.(1).(i)<-vectb.(i) with _->m.(1).(i)<-0;
      end
    done;
    m.(0),m.(1);;

(* Karatsuba multiplication function *)
let karatsuba a b=
  let decompa0=decomposition10 a
  and decompb0=decomposition10 b in
  let decompa,decompb=equilibre decompa0 decompb0 in
  let rec karatsuba_aux xi xj yi yj=
    if xi=xj && yi=yj then decompa.(xi)*decompb.(yi) else
      begin
        let x1y1=karatsuba_aux xi ((xi+xj)/2+1) yi ((yi+yj)/2+1)
        and x2y2=karatsuba_aux ((xi+xj)/2) xj ((yi+yj)/2) yj
        and x1y2=karatsuba_aux xi ((xi+xj)/2+1) ((yi+yj)/2) yj
        and x2y1=karatsuba_aux ((xi+xj)/2) xj yi ((yi+yj)/2+1)
        and m2=(xi-xj)/2+1 in
          x1y1*(10**(2*m2))+(x1y2+x2y1)*(10**m2)+x2y2
      end
  in
    karatsuba_aux (Array.length decompa -1) 0 (Array.length decompa-1) 0;;