// NumbThy.js
// Simple number theory functions written in JavaScript
// Author: Robert Campbell, <rcampbel@umbc.edu>
// Date: 7 Mar, 2010
// Freely available for use, abuse and modification (see license at bottom)

function addMod(a,b,n)
{
   var sum;
   sum = (a + b) % n;
   if (sum < 0) {sum += n;};
   return sum;
}

function multMod(a,b,n)
{
   var prod;
   prod = (a * b) % n;
   if (prod < 0) {prod += n;};
   return prod;
}

// Compute the greatest common divisor of integers a and b.
function gcd(a,b)
{
   var temp;
   if(a < 0) {a = -a;};
   if(b < 0) {b = -b;};
   if(b > a) {temp = a; a = b; b = temp;};
   while (true) {
      a %= b;
      if(a == 0) {return b;};
      b %= a;
      if(b == 0) {return a;};
   };
   return b;
}
      
// Compute the extended gcd - Given a and b compute gcd(a,b) and
// find x, y such that gcd(a,b) = x*a + y*b
function xgcd(a,b)
{
   var quot, a1=1, b1=0, a2=0, b2=1, retval=new Array();
   if(a < 0) {a = -a;};
   if(b < 0) {b = -b;};
   if(b > a) {temp = a; a = b; b = temp;};
   while (true) {
      quot = -Math.floor(a / b);
      a %= b;
      a1 += quot*a2; b1 += quot*b2;
      if(a == 0) {retval[0]=b; retval[1]=a2; retval[2]=b2; return retval;};
      quot = -Math.floor(b / a);
      b %= a;
      a2 += quot*a1; b2 += quot*b1;
     if(b == 0) {retval[0]=a; retval[1]=a1; retval[2]=b1; return retval;};
   };
   retval[0]=b; retval[1]=a2; retval[2]=b2; return retval;
}

// Compute the value of a^e(mod n) using the Russian Peasant algorithm.
function powmod(base,exp,modulus)
{
   var accum=1, i=0, basepow2=base;
   while ((exp>>i)>0) {
      if(((exp>>i) & 1) == 1){accum = (accum*basepow2) % modulus;};
      basepow2 = (basepow2*basepow2) % modulus;
      i++;
   };
   return accum;
}

// Test for primality by checking for gcd with small factors and then Miller-Rabin
// 2*3*5*7*11*13*17*19 = 9699690
function isprime(numb)
{
	return (([2,3,5,7,11,13,17,19].indexOf(numb)!=-1) || 
		((gcd(9699690,numb)==1) && isprime_mr(numb)));
}

// Test for primality using the converse of Fermat's Little Thm
// Fails (falsely declares composite to be prime) in some cases,
// including Carmichael numbers
function isprime_fermat(numb)
{
   var base;
   // Simple Fermat test - fails on Carmichael Numbers
   for(base=2; base<=5; base++){
      if(powmod(base, numb-1, numb) != 1) {
         return false;};
   };
   return true;
}

// Pseudoprime test using the Miller-Rabin algorithm
// Note: Fails for numb = 2
// Test cases: 1237, 19793 (prime), 2047, 15841 (spp base 2), 121, 703 (spp base 3), 561, 1105 (Carmichael)
function isprime_mr(numb)
{
   	var i, r, d, pow2, base, x;
   	// Compute largest pow2 so that (numb-1) = d*2^pow2
   	d = numb-1;
   	for(pow2=0; (((1<<pow2) & d)==0); pow2++){};
   	d = (d >> pow2);
   	for(base = 2; base <= 7; base++)
   	{
   		if((base>2) && ((base & 1)==0)) {continue;};  // Composite bases give redundant tests
      	x = powmod(base,d,numb);
      	if((x!=1) && (x!=(numb-1)))
      	{
         	for(r=1; r<=(pow2-1); r++)
	 	 	{
	    		x = multMod(x,x,numb);
	    		if(x==1){return false;};   // Previous value of x was sqrt(1) other than -1
	    		if(x==(numb-1)){break;};
	 		};
	 		if(x!=(numb-1)){return false;};  // base^((numb-1)/2) != 1 or -1, equiv to Fermat test
      	};
   	};
   	return true;
}

// Factor by first checking for factors of 2*3*5*7*11*13*17*19 = 9699690
// then use Pollard Rho
function factor(numb)
{
	var fact, smallfacts=[2,3,5,7,11,13,17,19];
	fact = gcd(9699690,numb);
	if(fact != 1){
		if(fact != numb)
			{return fact;}
		else 
			{
			for(var i=0; i<smallfacts.length; i++)
				{
				if((numb % smallfacts[i])==0)
					{return smallfacts[i];};
				};
			};
		}
	else
		{return factor_rho(numb);};
}

// Extract a factor from numb using Pollard's Rho method
function factor_rho(numb)
{
   var base=2, bound=3*Math.floor(Math.pow(numb,1/4)), i, fact;
   var slow=base;
   var fast = (slow*slow+1) % numb;
   for(i=1; i<bound; i++){
      slow = (slow*slow+1) % numb;
      fast = (fast*fast+1) % numb;
      fast = (fast*fast+1) % numb;
      if(fast - slow == 0){
         i=1; base++;
      } else {
         if((fact = gcd(fast-slow,numb)) != 1) return fact;
      };
   };
   return 1;
}

// Extract a factor from numb using Pollard's P-1 (aka Lambda) method
function factor_pm1(numb)
{
	var bound=Math.floor(Math.exp(Math.log(numb)/3))+6, i=2, thegcd, base=2;
	for (i=2; i<bound; i++) {
		base = powmod(base,i,numb);
		if((thegcd=gcd(base-1,numb)) != 1) {return thegcd;};
	};
	return 1;
}

// Return a list of all prime factors, in the form [[p1,e1],[p2,e2],...],
// where the prime factorization is numb = (p1^e1) * (p2^e2) * ...
function factorlist(numb)
{
	function factlist_core(numb)
	{
		var thefact;
		if(isprime(numb)) {return [numb];};
		thefact = factor(numb);
		return factlist_core(thefact).concat(factlist_core(numb/thefact));
	};
	function collect_facts(factlist)
	{
		var j=-1, outlist=[];
		for (var i=0; i<factlist.length; i++){
			if((i==factlist.length-1) || (factlist[i] != factlist[i+1])){
				outlist.push([factlist[i],i-j]);
				j=i;
			};
		};
		return outlist;
	};
	return collect_facts(factlist_core(numb).sort(function(a,b){return a - b}));
}

// (this is the Simplified BSD License, aka FreeBSD license)
// Copyright 1998-2010 Robert Campbell. All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification, 
// are permitted provided that the following conditions are met:
//
//    1. Redistributions of source code must retain the above copyright notice,
//    this list of conditions and the following disclaimer.
//
//    2. Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the 
//    documentation and/or other materials provided with the distribution.