(* Karatsuba *)

(*	On va prendre B= la plus grande puissance de 10 telle que (B-1)^2 soit calculé juste: pour caml c'est 4

	dans la formule de Karatsuba:
++++++++++++++++++++++++++++++++++++++++++++++++++++
Si x=U.B+V et x'=U'.B+V', 
on nomme R=V.V', Q=-(U-V)(U'-V'), P=U.U' 
alors x.x'= U.U'(B^2+B)-(U-V).(U'-V').B+V.V'.(B+1)
soit donc  xx'=P.B^2+(P+Q+R).B+R
qui sera reformatté en c3.B^3+c2.B^2+c1.B+c0	
++++++++++++++++++++++++++++++++++++++++++++++++++++

 le coefficient de 1 sera calculé juste, il fournira une retenue au plus égale à B-2, on y ajoutera U.U'-(U-V).(U'-V')+V.V' = égal à UV'+VU', d'où un somme au plus égale à (B-1)^2+B-2 (= moins de 10^8) d'où une retenue au plus de (B-1) qui ajoutée à U.U' fournira un dernier calcul exact
*)
let B=10000;; 	(* B c'est 10^4 *)
let nbch=4;;		(* nbch: nombre de chiffres, sert à l'affichage *)

(* on commence par définir 
* un utilitaire d'affichage (sans débordement de la fenêtre)
* des div et mod qui marchent, un xor
* des promotions pour les entiers devenant Multis
* un formattage de Multi à retenues
* une multi-addition
* une multi-comparaison
* un multi-produit par un petit entier
* un multi-affichage
* des opérations concernant les MultisSous: extraits de Multis
*)

(* comme 1000! a moins de 3000 chiffres en base 10, voyons quel K prendre pour que 2^K*4 dépasse 3000. On obtient 9.55...., on va donc prendre 10 *)

let N = 16;;      (* 3 ou 5 pour tester 1024 pour de vrai *)
let N1=N-1;;  (* pour les boucles de 0 à N-1, on mettra N1 *)

(*----------------------------utilitaire d'affichage-----------------------------------*)
let global = ref (pos_out stdout);;

let saute_ligne()=
	let actuel = (pos_out stdout) in
	if (actuel-(!global))>40 
	then begin global:=actuel; print_newline() end
;;

let p_s s=saute_ligne();print_string s;;
let p_i i=saute_ligne();print_int i;;
let p_f f=saute_ligne();print_float f;;

(*---------------------------------------*)
(* le Multiaffichage *)
let Maffiche p=
	let aff_ent n=
		let sn=string_of_int n in 
		let zeros=make_string (nbch-string_length sn) `0` in
		p_s (zeros^sn^" ")
	in for k=(vect_length p -1) downto 0 do aff_ent p.(k) done;;


(*--------confection de mod et div qui marchent en caml -------------------------*)
(* q est supposé positif, p est dans Z *) 
#infix "Bdiv";;
#infix "Bmod";;
#infix "xor";;

let prefix Bmod p q=if p>=0 then p mod q else (q+(p mod q)) mod q;;
let prefix Bdiv p q= (p- (p Bmod q) ) / q;;
let prefix xor a b=(a or b) & not (a & b);;

(*---------------inclusion de N dans les Multis---------------------------------------*)
let  int_a_multi k=let u=make_vect N 0 in u.(0)<-k;u;;

let zero= int_a_multi 0;; 
let un= int_a_multi 1;;
let deux= int_a_multi 2;;

(*------------formattage d'un Multi malpropre --------------------------------------*)
let formatte t=
	let retenue=ref 0 and provisoire=ref 0 in
	for k=0 to ((vect_length t)-1) do
		provisoire:= t.(k)+ !retenue;
		t.(k)<-(!provisoire Bmod B);
		retenue:=!provisoire Bdiv B 
	done;
	if (!retenue != 0) then print_string" Attention! ça déborde " 
;;


(* -------------------- toilettage rapide d'un multi propre avant k -----------------*)
let rec nettoie m k=
	if m.(k)>(B-1)
	then let retenue=m.(k) Bdiv B in 
		begin
			m.(k+1)<-m.(k+1)+retenue;
			m.(k)<-m.(k) Bmod B;
			nettoie m (k+1)
		end;;

(*------------------ les opérations standard sur les Multis --------------------------*)
(* Majoute p q r va ajouter p et q et mettre le résultat dans r *)
let Majoute p q r=for k=0 to N1 do r.(k)<-p.(k)+q.(k)done;formatte r;;

(* Minférieur p q dit si le Multi p est inférieur au Multi q *)
let Minferieur p q=	
	let k=ref N1 in 
	while (p.(!k)=q.(!k)) & (!k >0) 	do decr k done;
	p.(!k)<=q.(!k)
;; 

(* Mmoins p q r sait déjà que p est supérieur à q, r sera p-q *)
let Mmoins p q r=for k=0 to N1 do r.(k)<-p.(k)-q.(k)done;formatte r;;

(* Mproduit n p r est un "produit externe": n est un int, seul p est un Multi, le résultat est rangé dans r *) 
let Mproduit n p r=	if n>=B 
			then print_string "Attention! ça déborde "
			else for k=0 to N1 do r.(k)<-p.(k)*n done;
			formatte r
;;

(* Mfois p q est le produit bête de p et q, il pourra servir à "vérifier Karatsuba", rien n'est fait pour l'optimiser, on n'y vérifie même pas si cela va déborder *)
let Mfois p q r=
	let nbchiff m=
		let vu=ref (-1) in 
		for k=0 to N1 
		do if m.(k)<>0 then vu:=k; done;
		vu in let nbchp=nbchiff p and nbchq=nbchiff q in
		for i=0 to !nbchp 	do
			for j=0 to !nbchq 	do
				r.(i+j)<-r.(i+j)+p.(i)*q.(j);
				nettoie r (i+j)
						done
					done
;;

(* Multicalcul de grandes factorielles *)
let rec fact n=let res=int_a_multi 1 in 
	match n with
	|0 ->res |1->res
	|n->Mproduit (n) (fact (n-1)) res;
	res
;;

(* Maffiche (fact 100);; *)

(*-------------------opérations sur les sous chaînes----------------------------------*)
(* MS signifie "MultiSous" les MSplus (...) sont comme ajoute (...) mais ils prennent leurs données dans p et q à partir de deb sur la longueur long, ils les rangent dans r de 0 à (long-1) 
MSautoinférieur compare deux sous nombres d'un même Multi, de même pour MSautomoins*)

let MSplus deb long p q r=for k=0 to (long-1) do r.(k)<-p.(deb+k)+q.(deb+k)done;formatte r;;

let  MSmoins deb long p q r=for k=0 to (long-1) do r.(k)<-p.(deb+k)-q.(deb+k)done;formatte r;;

let MSinferieur deb long p q=	
	let k=ref (deb+long-1) in 
	while (p.(!k)=q.(!k)) & (!k >deb) do decr k done;
	p.(!k)<=q.(!k)
;; 

let MS_auto_inferieur deb long p =	
	let k=ref (deb+long-1) in 
	while (p.(!k)=p.(long + !k)) & (!k >deb) do decr k done;
	p.(!k)<=p.(long + !k)
;; 

let MSautomoins deb1 deb2 long p r=for k=0 to (long-1) do r.(k)<-p.(deb1+k)-p.(deb2+k)done;formatte r;;

let MScopie deb long x y=for k= 0 to (long-1) do y.(k)<-x.(deb+k) done;; 

let MSplusdecale decale long x y=for k=0 to (long-1) do y.(decale+k)<-x.(k)+y.(decale+k) done;;

(* kara k a deba b debb c debc affectue le katsuba-produit de la tranche allant de deba à deba+k dans a  par la tranche allant de debb à debb+k dans b. le résultat est stocké dans c à partir de 0*)

(* ci-dessous le stockage des intermédiaires de calcul n'est pas optimisé, ni sur le plan vitesse ni sur le plan encombrement total 
	par contre je n'ai pas voulu recopier les morceaux intermédiaires dans des variables locales; c'est ce "scrupule" qui complique tout, et je ne sais même pas si le jeu en valait la chandelle *)

(* je n'ai pas vraiment trouvé comment morceler Karatsuba pour en faire une fonction compréhensible, j'ai (un peu) essayé d'étendre mes Multis à Z en étendant les opérations (ce fut laid), j'ai essayé de définir deux fonctions Karatsuba et karakeur s'appellant l'une l'autre (karakeur ne traitant que le problème central) et ce fut re-laid *)
 
(* le fait de passer a et b à kara va amener un grand nombre de recopiages inutiles, la solution consiste à recopier ce kara comme fonction locale de karatsuba en y enlevant a et b qui deviendront les x et y de karatsuba  *)

let rec kara k a deba b debb c =
		let kk=k+k in
		let k2 = k/2 in let a2=deba+k2 in let b2=debb+k2 in
		let P=make_vect k 0 in
		let Q=make_vect k 0 in
		let R=make_vect k 0 in
		let I=make_vect k 0 in
		let J=make_vect k2 0 in
		let u_inf_v=ref true in let u'_inf_v'=ref true in

		if k=1
		then 	begin c.(0)<-a.(deba)*b.(debb); nettoie c 0 end
		else 	begin
				kara k2 a deba b debb R;
				kara k2 a a2 b b2 P;
				u_inf_v:=not MS_auto_inferieur deba k2 a;
				u'_inf_v':=not MS_auto_inferieur debb k2 b;
(* I va contenir |U-V| *)
				if !u_inf_v 
				then MSautomoins deba a2 k2 a I
				else MSautomoins a2 deba k2 a I;
(* J va contenir |U'-V'| *)
				if !u'_inf_v'
				then MSautomoins  debb b2 k2 b J
				else MSautomoins  b2 debb k2 b J;

				MSplus 0 k P R Q;
				kara k2 I 0 J 0 I;

				if !u_inf_v xor !u'_inf_v'
				then MSplus 0 k Q I Q
				else MSmoins 0 k Q I Q;

				MScopie 0 k R c;
				MSplusdecale k2 k Q c;			
				MSplusdecale k k P c;
				formatte c;
			end;;

(* karatsuba va prendre deux entiers de même longueurs inférieures à 64 blocs de longueur 4*)
let karatsuba x y=
	let produit=make_vect (2*vect_length x) 0 in				 
	begin	 kara 64 x 0 y 0 produit; 	Maffiche produit end
;;

(* include "C:/OPTION/PROGSUP/COMPLEXI/Kara.ml";; *)

let M1=fact 10;; 
let M2=int_a_multi 0;;
kara 8 M1 0 M1 0 M2;;

 
let M1=fact 11;; 
let M2=int_a_multi 0;;
kara 8 M1 0 M1 0 M2;;
M2;;
(* à partir du carré de 30! il y a souvent un chiffre faux, le nombre de chiffres utilisé n'est pas idéal *)
(* en ayant pris N=16,  et les initialisations qui suivent, on a un résultat carrément faux ! *) 
(* c'est juste qu'il faut nettoyer w d'abord! gromalin va! *)
 let u= fact 10;; let v= fact 13;; let w=int_a_multi 57;;
kara 4 u 0 v 0 w;;
Maffiche w;;
kara 5 w 0 w 0 u;; Maffiche u;;
(* d'après AHU, page 317, c'est à partir de ... 500 chiffres, que Karatsuba commence à battre l'algorithme naïf *)