/* *********************************************************************

This source file is a part of the standard library of Scol
For the latest info, see http://www.scolring.org

Copyright (c) 2013 Stephane Bisaro aka Iri.

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA, or go to

http://www.gnu.org/copyleft/lesser.txt

********************************************************************* */

/*
 * Standard functions for mathematical objects
 * See http://redmine.scolring.org/projects/tutorials/wiki/Scol_usage
 * for more informations
 */

/*! \file maths.pkg
 * \author Scol team
 * \version 0.1
 * \copyright GNU Lesser General Public License 2.0 or later
 * \brief Scol Standard Library - Maths API
 *
 * \details This API provides some functions to manipulate some mathematical objects
 * 
 * Dependancies :
 * - lib/std/stdlib.pkg
 * 
 **/
 
/* When Scol will support the larger integers, these constantes will be changed */
var STD_MFIBONACCI = 43;;
var STD_M2POWER = 127.0;;
var STD_M2POWERI = 29;;

fun std_mfibonacci (n, o1, o2)=
	if n == 0 then
		o2
	else
		let _fooS strcat itoa n strcat "  " itoa o2 -> _ in
		std_mfibonacci n-1 o2 o2+o1;;
		
/*
fun std_mfibonaccif (n, o1, o2)=
	if std_fIsZero n then
		o2
	else
		let _fooS strcat ftoa n strcat "  " ftoa o2 -> _ in
		std_mfibonaccif n-.1.0 o2 o2+.o1;;
		
fun std_mFibonacciF (n)=
	let std_mfibonacci absf itof n 0.0 1.0 -> r in
	if ((n < 0) && (!mod n 2)) then
		(-.r)
	else 
		r;;
*/

/*! \brief Return the nth value of the Fibonacci series
 * (1, 1, 2, 3, 5, 8, 13, 21, ...)
 *
 * \remark Negative numbers are supported.
 * \ingroup std_math
 * <b>Prototype :</b> : fun [I] I
 *
 * \param I : an integer (-43 to 43 included)
 * \return I : the result or nil if out of range
 **/
fun std_mFibonacci (n)=
	if STD_MFIBONACCI < abs n then
		nil
	else
		let std_mfibonacci abs n 0 1 -> r in
		if ((n < 0) && (!mod n 2)) then	// negatif et paire
			(-r)
		else 
			r;;

/*! \brief Return the power of two : \f$2^{f}\f$
 * 
 * \remark Negative numbers are NOT supported.
 *
 * \ingroup std_math
 * <b>Prototype :</b> : fun [F] F
 *
 * \param F : a float number (0 to 127 included)
 * \return F : the result or nil if out of range
 **/
fun std_m2power (x)=
	if (x >. STD_M2POWER) && (x >=. 0.0) then
		nil
	else
		pow 2.0 x;;

/*! \brief Return the power of two : \f$2^{i}\f$
 * 
 * \remark Negative numbers are NOT supported.
 *
 * \ingroup std_math
 * <b>Prototype :</b> : fun [I] I
 *
 * \param I : a float number (0 to 29 included)
 * \return I : the result or nil if out of range
 **/
fun std_m2powerI (x)=
	if (x > STD_M2POWERI) && (x >= 0) then
		nil
	else
		ftoi pow 2.0 itof x;;
		

/*! \brief Convert degree to radian
 *
 * \ingroup std_math
 * <b>Prototype :</b> : fun [F] F
 *
 * \param F : a value in degree
 * \return F : the same value in radian
 **/
fun std_mDegToRad (f)=
	mulf f ((2.0 *. PIf) /. 360.0);;

/*! \brief Convert radian to degree
 *
 * \ingroup std_math
 * <b>Prototype :</b> : fun [F] F
 *
 * \param F : a value in radian
 * \return F : the same value in degree
 **/
fun std_mRadToDeg (f)=
	divf f ((2.0 *. PIf) /. 360.0);;
	
/*! \brief Return the GCD (Greatest Common Divisor, PGCD in french) of 
 * two integers
 *
 * \ingroup std_math
 * <b>Prototype :</b> : fun [I I] I
 *
 * \param I : a first integer
 * \param I : a second integer
 * \return I : the GCD
 **/	
fun std_mGCD (int1, int2)=
	if int2 == 0 then
		int1
	else
		std_mGCD int2 mod int1 int2;;
		
/*! \brief Return the LCM (Least Common Multiple, PPCM in french) of 
 * two integers
 *
 * \ingroup std_math
 * <b>Prototype :</b> : fun [I I] I
 *
 * \param I : a first integer
 * \param I : a second integer
 * \return I : the LCM
 **/
fun std_mLCM (int1, int2)=
	if (int1 == 0) || (int2 == 0) then
		0
	else
		(int1 * int2) / (std_mGCD int1 int2);;
	
/*! \brief Calculate the power of an integer by a quick recursion
 *
 * \ingroup std_math
 * <b>Prototype :</b> : fun [I I] I
 *
 * \param I : an integer
 * \param I : a power (should be >= 0 otherwise, nil will be returned)
 * \return I : the result or nil if error
 **/
fun std_mPow (i, n)=
	if n < 0 then
		nil
	else
		if n == 0 then
			1
		else if n == 1 then
			i
		else
			let std_mPow i n/2 -> t in
			if 0 == mod n 2 then
				t*t
			else
				i*t*t;;


fun STD_M_PRIME () = 1;;
fun STD_M_COMPOSITE () = 0;;

// Write (n - 1) as 2^s * d
fun std_mprimality_step1 (s, d)=
	if (mod d 2) == 0 then
		std_mprimality_step1 s+1 d/2
	else
		[s d];;
		
/*fun std_mprimality_step2 (n)=
	if (n == 2) || (n == 3) then
		STD_M_PRIME
	else if (mod n 2) || (n < 2) then
		STD_M_COMPOSITE
	else
		let std_mprimality_step1 0 n-1 -> [s d] in
		let std_random */
		
		
		
// calcul a^b%c		
fun std_mprime_mr_power (a, b, c)=
	let 1 -> out in
	(
	while (b) do
	(
		if (b&1) then
			set out = mod (out * a) c
		else
			0;
		set a = mod (a * a) c;
		set b = b>>1;
	);
	out
	);;
	
// n−1 = 2^s * d with d odd by factoring powers of 2 from n−1
fun std_mprime_mr_witness (n, s, d, a)=
	let std_mprime_mr_power a d n -> x in
	let nil -> out in
	let 0 -> y in
	(
	while s do
	(
		set y = mod (x * x) n;		
		if (y == 1) && (x != 1) && (x != n-1) then
			set out = STD_M_COMPOSITE
		else
			0;
		set x = y;
		set s = s-1;
	);
	if out != nil then
		out
	else
		if y != 1 then 
			STD_M_COMPOSITE
		 else
			STD_M_PRIME
	);;
	
fun std_mprime_mr_det (n)=
	 if (((!(n & 1)) && (n != 2)) || (n < 2) || ((0 == mod n 3) && (n != 3))) then
		STD_M_COMPOSITE
	else if (n <= 3) then
		STD_M_PRIME
	else
		let n / 2 -> d in
		let 1 -> s in
		(
		while (!(d & 1)) do
		(
			set d = d / 2;
			set s = s + 1
		);
		if n < 1373653 then
			(std_mprime_mr_witness n s d 2) && (std_mprime_mr_witness n s d 3)
		else if n < 9080191 then
			(std_mprime_mr_witness n s d 31) && (std_mprime_mr_witness n s d 73)
		else // for supported integers by Scol (type : I), in fact else if n < 4759123141 then
			(std_mprime_mr_witness n s d 2) && (std_mprime_mr_witness n s d 7) && (std_mprime_mr_witness n s d 71)
		/* here greater number than 4759123141
		else if n < 1122004669633 then
			(std_mprime_mr_witness n s d 2) && (std_mprime_mr_witness n s d 13) && (std_mprime_mr_witness n s d 23) && (std_mprime_mr_witness n s d 1662803)
		else if n < 2152302898747 then
			(std_mprime_mr_witness n s d 2) && (std_mprime_mr_witness n s d 3) && (std_mprime_mr_witness n s d 5) && (std_mprime_mr_witness n s d 7) && (std_mprime_mr_witness n s d 11)
		else if n < 3474749660383 then
			(std_mprime_mr_witness n s d 2) && (std_mprime_mr_witness n s d 3) && (std_mprime_mr_witness n s d 5) && (std_mprime_mr_witness n s d 7) && (std_mprime_mr_witness n s d 11) && (std_mprime_mr_witness n s d 13)
		else // until 341550071728321
			(std_mprime_mr_witness n s d 2) && (std_mprime_mr_witness n s d 3) && (std_mprime_mr_witness n s d 5) && (std_mprime_mr_witness n s d 7) && (std_mprime_mr_witness n s d 11) && (std_mprime_mr_witness n s d 13) && (std_mprime_mr_witness n s d 17)
		*/
		);;
		
/*! \brief Returns if an integer is a prime number or a composite number.
 *
 * This is a deterministic implementation of the Miller-Rabin primality
 * algorithm.
 * 
 * \ingroup std_math
 * <b>Prototype :</b> : fun [I] I
 *
 * \param I : an integer
 * \return I : STD_M_COMPOSITE if number is composite, STD_M_PRIME if
 * number is a strong probably prime or nil if error
 **/
fun std_mPrimeMRDet (iNumber)=
	if std_objIsNil iNumber then
		nil
	else
		std_mprime_mr_det iNumber;;