High-precision calculations

This page contains operators for multi-precision arithmetic written in PHP.
(Most of them are not optimized and perform calculations in decimal₁₀)

• mAbs()
• mSign()
• mCmp()
  

//
function mAbs($p)
{
    $p = (string)$p;
    if (!strlen($p) || $p[0] != '-') return $p;
    return substr($p, 1);
}
//
function mSign($p)
{
    return strcmp(mAbs($p), $p) ? -1:1;
}
//
function mCmp($p, $q)
{
    //
    $ps = mSign($p);
    $qs = mSign($q);
    if ($ps != $qs)
        return $ps < $qs ? -1:1;
    //
    $p = mAbs($p);
    $q = mAbs($q);
    $m = max(strlen($p), strlen($q));
    $p = str_pad($p, $m, '0', STR_PAD_LEFT);
    $q = str_pad($q, $m, '0', STR_PAD_LEFT);
    //
    for ($i = 0; $i < strlen($p); $i++)
    {
        $pn = $p[$i] - '0';
	$qn = $q[$i] - '0';
	$diff = $pn - $qn;
	if (! $diff) continue;
        return $diff < 0 ? -1:1;
    }
    return 0;
}
//

• mAdd()
• mSub()

//
function mAdd($p, $q, $base = 10)
{
    //
    $ps = mSign($p);
    $qs = mSign($q);
    $pd = $ps < 0;
    //
    if ($ps*$qs < 0)  return mSub($pd ? $q:$p, mAbs($pd ? $p:$q), $base);
    else if ($pd && $qs < 0)return '-' . mAdd(mAbs($p), mAbs($q), $base);
    //
    if (strlen($q) < strlen($p))
	list($p, $q) = [$q, $p];
    $p = str_pad($p, strlen($q), '0', STR_PAD_LEFT);
    //
    $rr = '';
    $cc = 0;
    for ($i = strlen($p) - 1; $i >= 0; $i--)
    {
	$pn = $p[$i] - '0';
	$qn = $q[$i] - '0';
	$nn = $pn + $qn + $cc;
	$cc = intval($nn / $base);
	$rr = (string)($nn % $base) . $rr;
    }
    if ($cc) $rr = $cc . $rr;
    return $rr;
}
//
function mSub($p, $q, $base = 10)
{
    $ps = mSign($p);
    $qs = mSign($q);
    if ($qs < 0)      return mAdd($p, mAbs($q), $base);
    else if ($ps < 0) return '-' . mAdd(mAbs($p), $q, $base);
    //
    $m = max(strlen($p), strlen($q));
    $p = str_pad($p, $m, '0', STR_PAD_LEFT);
    $q = str_pad($q, $m, '0', STR_PAD_LEFT);
    $s = '';
    //
    if (mCmp($p, $q) < 0)
	list($p, $q, $s) = array($q, $p, '-');
    //
    $rr = '';
    $bb = 0;
    for ($i = $m - 1; $i >= 0; $i--)
    {
	$pn = $p[$i] - '0';
	$qn = $q[$i] - '0';
	$nn = $pn - $qn - $bb;
	$bb = (int)($nn < 0);
	if ($bb) $nn += $base;
	$rr = $nn . $rr;
    }
    return $s . $rr;
}
//

• mBaseMul()
• mKaratsubaMul()

//
function mBaseMul($p, $q)
{
    $p  = intval($p);
    $rr = '';
    $j  = 0;
    for ($i = strlen($q) - 1; $i >= 0; $i--)
    {
        $nn = $p * ($q[$i] - '0');
		$nn.= str_repeat('0', $j++);
        $rr = mAdd($nn, $rr);
    }
    return $rr;
}
//
function mKaratsubaMul($p, $q)
{
    if (! trim($p, '0'))
	return '0';
    //
    $s  = mSign($p) * mSign($q) < 0 ? '-':'';
    $p  = mAbs($p);
    $q  = mAbs($q);
    //
    $m = max(strlen($p), strlen($q));
    if ($m % 2) $m++;
	$p = str_pad($p, $m, '0', STR_PAD_LEFT);
	$q = str_pad($q, $m, '0', STR_PAD_LEFT);
    if ($m < 4)
        return $s . mBaseMul($p, $q);
    //
    $k  = intdiv($m, 2);
    $p1 = substr($p, 0, $k);
    $q1 = substr($q, 0, $k);
    $p0 = substr($p, $k);
    $q0 = substr($q, $k);
    $o0 = mSub($p0, $p1);
    $o1 = mSub($q0, $q1);
    $r0 = mKaratsubaMul($p0, $q0);
    $r1 = mKaratsubaMul($p1, $q1);
    $r2 = mKaratsubaMul(mAbs($o0), mAbs($o1));
    $ts = mSign($o0) * mSign($o1);
    $t0 = mAdd($r0, $r1);
    $t0 = $ts < 0 ? mAdd($t0, $r2) : mSub($t0, $r2);
    $t0.= str_repeat('0', $k);
    $r1.= str_repeat('0', 2*$k);
    $r  = ltrim(mAdd(mAdd($r0, $t0), $r1), '0');
    return $s . ($r ? $r:'0');
}
//

• mBaseDiv()
• mSvobodaDiv()
  
• mRecursiveDiv()
• mUnbalancedDiv()
• mWordDiv()
• mBarrettDiv()
  

//
function mBaseDiv($p, $q)
{
    //
    if (mCmp($p, $q) < 0) return ['0', $p];
    if (! trim($p, '0'))  return ['0','0'];
    if (! mCmp($q, '1'))) return [$p, '0'];
    //
    $k = 1;
    while (($l = $q[0] - '0') < 5)
    {
	$k *= 2;
	$p = mKaratsubaMul($p, '2');
	$q = mKaratsubaMul($q, '2');
    }
    //
    $n = strlen($q);
    $m = strlen($p) - $n;
    $l = $q[0] - '0';
    $a = array();
    //
    $qq = $q . str_repeat('0', $m);
    $a[$m] = 0;
    if (mCmp($p, $qq) >= 0)
        list($a[$m], $p) = [1, mSub($p, $qq)];
    //
    for ($i = $m - 1; $i >= 0; $i--)
    {
        $ii = str_repeat('0', $i);
        $pn = $p[$m - $i - 1] - '0';
        $pp = $p[$m - $i]     - '0';
        $nn = $pn * 10 + $pp;
        //
        $a[$i] = min(intval($nn / $l), 9);
        $t = mKaratsubaMul($q, $a[$i]) . $ii;
        $p = mSub($p, $t);
        //
        while (mSign($p) < 0)
            list($a[$i], $p) = [$a[$i] - 1, mAdd($p, $q . $ii)];
    }
    //
    $r = mBaseDiv($p, (string) $k)[0];
    return [ltrim(implode('', $a), '0'), $r];
}
//
function mSvobodaDiv($p, $q)
{
    //
    $n = strlen($q); 
    $m = strlen($p) - $n;
    //
    $kk= mBaseDiv('1' . str_repeat('0', $n + 1), $q);
    $k = $kk[0];
    if (! trim($kk[1], '0'))
    	$k = mAdd($k, '1');
    $b = mKaratsubaMul($k, $q);
    $a = array();
    //
    for ($i = $m - 1; $i > 0; $i--)
    {
        $ii = str_repeat('0', $i - 1);
        $qq = $b . $ii;
        $a[$i] = $p[$m - $i - 1];
        $p = mSub($p, mKaratsubaMul($a[$i], $qq));
        if (mSign($p) < 0)
            list($a[$i], $p) = [$a[$i] - 1, mAdd($p, $qq)];
    }
    //
    $qq = implode('', $a);
    $pp = mBaseDiv($p, $q);
    return [mAdd(mKaratsubaMul($k, $qq), $pp[0]), $pp[1]];
}
//
function mRecursiveDiv($p, $q)
{
    //
    $n = strlen($q);
    $m = strlen($p) - $n;
    if ($m < 2)
	return mBaseDiv($p, $q);
    //
    $k  = intval($m / 2);
    $kk = str_repeat('0', $k);
    $q1 = substr($q, 0, $k);
    $q0 = substr($q, $k);
    //
    [$qq1, $rr1] = mRecursiveDiv(ltrim(substr($p, 0, -2*$k), '0'), $q1);
    $a1 = mAdd($rr1 . str_repeat('0', 2*$k), substr($p, -2*$k));
    $a1 = mSub($a1, mKaratsubaMul($qq1, $q0) . $kk);
    while (mSign($a1) < 0)
	list($qq1, $a1) = [$qq1 - 1, mAdd($a1, $q . $kk)];
    //
    [$qq0, $rr0] = mRecursiveDiv(ltrim(substr($a1, 0, -$k), '0'), $q1);
    $a0 = mAdd($rr0 . $kk, substr($a1, -$k));
    $a0 = mSub($a0, mKaratsubaMul($qq0, $q0));
    while (mSign($a0) < 0)
	list($qq0, $a0) = [$qq0 - 1, mAdd($a0, $q)];
    //
    $a = mAdd($qq1 . $kk, $qq0);
    return [$a, $a0];
}
//
function mUnbalancedDiv($p, $q)
{
    $n  = strlen($q);
    $m  = strlen($p) - $n;
    $a  = '0';
    //
    while ($m > $n)
    {
	$m = $m - $n;
	[$qq, $rr] = mRecursiveDiv(ltrim(substr($p, 0, -$m), '0'), $q);
	$a = mAdd($a . str_repeat('0', $n), $qq);
	$p = mAdd(substr($p, -$m), $rr . str_repeat('0', $m));
    }
    //
    [$qq, $rr] = mRecursiveDiv($p, $q);
    return [mAdd($a . str_repeat('0', $m), $qq), $rr];
}
//
function mWordDiv($p, $w)
{
    $a = array();
    $n = strlen($p);
    $d = fmod(1 / $w, 10);
    $b = 0;
    //
    for ($i = $n - 1; $i >= 0; $i--)
    {
	$fi = $n - 1 - $i;
	$ii = $p[$fi] - '0';
	$x  = $ii - $b;
	$b0 = 0;
	//
	if ($x < 0)
	    $x += 10 + $b0++;
	//
	$a[] = fmod($d * $x, 10);
	$b1  = ($a[$fi] * $w - $x) / 10;
	$b   = $b0 + $b1;
    }
    //
    return [$a, $b];
}
//
function mBarettDiv($p, $q)
{
    $b  = strlen(ltrim($q, '0'));
    $m  = '1' . str_repeat('0', $b * 2);
    $m  = mBaseDiv($m, $q)[0];
    $p  = str_pad($p, ($b * 2), '0', STR_PAD_LEFT);
    $a1 = substr($p, 0, $b);
    $a0 = substr($p, $b);
    //
    $qq = substr(mKaratsubaMul($a1, $m), 0, -$b);
    $rr = mSub($p, mKaratsubaMul($qq, $q));
    while (mCmp($q, $rr) <= 0)
	[$qq, $rr] = [mAdd($qq, '1'), mSub($rr, $q)];
    //
    return [$qq, $rr];
}
//

• mPow()
• mRoot()
• mLog()
• mBitRev()
  

//
function mPow($p, $n)
{
    //
    $r = $p;
    for ($i = 1; $i < $n; $i++)
        $r = mKaratsubaMul($r, $p);
    return $r;
}
//
function mRoot($p, $root = 2)
{
    $uu = $p;
    $rr = (string)($root - 1);
    //
    do
    {
	$s = $uu;
	$t = mKaratsubaMul($rr, $s);
	$t = mAdd($t, mBaseDiv($p, mPow($s, $rr))[0]);
	$uu= mSvobodaDiv($t, (string)$root)[0];
    } while (mCmp($uu, $s) < 0);
    //
    return $s;
}
//
function mLog($p, $b = '2')
{
    $n = 0;
    $r = '';
    //
    while (mCmp('0', $p))
    {
	$n++;
        $p = mBaseDiv($p, $b);
	$r = $p[1] . $r;
	$p = $p[0];
    }
    //
    return [$n, $r];
}
//
function mBitRev($j, $n)
{
    [,$o] = mLog($j);
    $o = str_pad($o, $n, '0', STR_PAD_LEFT);
    $p = 1; 
    $r = 0;
    //
    for ($i = 0; $i < strlen($o); $i++)
    {
	$r += $p * ($o[$i] - '0');
	$p *= 2;
    }
    //
    return $r;
}
//

• mBaseGcd()
• mBezout()
• mDoubleGcd()
• mExtendedGcd()
  

//
function mBaseGcd($p, $q)
{
    //
    if (! ltrim($q, '0'))   return $p;
    if (mCmp($p, '1') <= 0) return '1';
    //
    while (true)
    {
    	list($qq, $rr) = mBaseDiv($p, $q);
	if (! ltrim($rr, '0'))
	    break;
	$p = $q;
	$q = $rr;
    }
    //
    return $q;
}
//
function mBezout($p, $q, $shortcut = true)
{
    // a = t*a+u*b (t=1, u=0)
    // b = v*a+w*b (v=0, w=1)
    //
    [$t, $u, $v, $w] = ['1', '0', '0', '1'];
    $b = false;
    $h = intdiv(max(strlen($p), strlen($q)), 2);
    //
    while (true)
    {
        [$qq, $rr] = mBaseDiv($p, $q);
	if (! ($rr = ltrim($rr, '0'))) break;
	//
	[$tt, $tv] = [$t, $v];
	$t = $u;
	$v = $w;
	$u = mSub($tt, mKaratsubaMul($u, $qq));
	$w = mSub($tv, mKaratsubaMul($w, $qq));
	//
	if ($shortcut && $b) break;
	$b = strlen($rr) <= $h;
	//
	$p = $q;
	$q = $rr;
    }
    //
    return array($q, $t, $u, $v, $w);
}
//
function mDoubleGcd($p, $q)
{
    $p = ltrim($p, '0');
    $q = ltrim($q, '0');
    $n = strlen($p);
    $m = strlen($q);
    //
    if (! ltrim($q, '0')) return $p;
    if ($m < 2)		  return mBaseGcd($p, $q);
    if (mCmp($p, $q) < 0 || $m < $n)
	return mDoubleGcd($q, mBaseDiv($p, $q)[1]);
    //
    [$_, $t, $u, $v, $w] = mBezout(
	substr($p, 0, 2),
	substr($q, 0, 2)
    );
    //
    $pp = mAdd(mKaratsubaMul($t, $p), mKaratsubaMul($v, $q));
    $qq = mAdd(mKaratsubaMul($u, $p), mKaratsubaMul($w, $q));
    return mDoubleGcd(mAbs($pp), mAbs($qq));
}
//
function mExtendedGcd($p, $q)
{
    [$g,,
     $u,,
     $v] = mBezout($p, $q, false);
    return array($g, $u, $v);
}
//

• mMod()
• mModVector()
• mModPointwise()

//
function mMod($p, $m)
{
    $q = mBaseDiv(mAbs($p), $m);
    if (mSign($p) < 0)
    {
	$q = mAdd($q[0], $q[1] != 0 ? '1':'0');
	$q = mKaratsubaMul($q, $m);
	if (! trim($q, '0')) $q = $m;
	$p = mAdd($p, $q);
    }
    //
    else if (mCmp($m, $p) <= 0)
    {
    	$q = mKaratsubaMul($q[0], $m);
	$p = mSub($p, $q);
    }
    return $p;
}
//
function mModVector($v, $m)
{
    $n = count($v);
    //
    for ($i = 0; $i < $n; $i++)
    {
    	$v[$i] = mMod($v[$i], $m);
    	$v[$i] = ltrim($v[$i], '0');
	if (! $v[$i]) $v[$i] = '0';
    }
    //
    return $v;
}
//
function mModPointwise($p, $q, $m)
{
    $pn = count($p);
    $qn = count($q);
    $rr = array();
    //
    if ($pn != $qn)
	return $rr;
    //
    for ($i = 0; $i < $pn; $i++)
    {
	$rr[] = mKaratsubaMul($p[$i], $q[$i]);
	while (mCmp($rr[$i], $m) >= 0)
	    $rr[$i] = mSub($rr[$i], $m);
    }
    //
    return $rr;
}
//

• mForwardFFT()
• mBackwardFFT()
• mSSMul()
  

//
function mForwardFFT($w, ...$vector)
{
    $k = count($vector);
    //
    if ($k <= 2) return [
	mAdd($vector[0], $vector[1]),
	mSub($vector[0], $vector[1])
    ];
    //
    $ww = mKaratsubaMul($w, $w);
    $pp = ['1', $w, $ww];
    $even = array();
    $odd  = array();
    $curr = true;
    //
    for ($i = 0, $j = $k - 4; $i < $k; $i++, $j--)
    {
	if ($i < $j)
	    $pp[] = mKaratsubaMul($pp[$i + 2], $w);
	//
	if ($curr) $even[] = $vector[$i];
	else	   $odd[] = $vector[$i];
	$curr ^= true;
    }
    //
    $rr   = [];
    $kk   = count($even);
    $even = mForwardFFT($ww, ...$even);
    $odd  = mForwardFFT($ww, ...$odd);
    [$n,] = mLog($kk - 1);
    //
    for ($j = 0; $j < $kk; $j++)
    {
        $m = mBitRev($j, $n);
	$w = $pp[$m];
	$m = mKaratsubaMul($w, $odd[$j]);
	//
	$rr[] = mAdd($even[$j], $m);
	$rr[] = mSub($even[$j], $m);
    }
    return $rr;
}
//
function mBackwardFFT($w, ...$vector)
{
    $k = count($vector);
    //
    if ($k <= 2) return [
	mAdd($vector[0], $vector[1]),
	mSub($vector[0], $vector[1])
    ];
    //
    $kk = ($k / 2);
    $pp = ['1', $w];
    for ($i = 0; $i < $k - 1; $i++)
	$pp[] = mKaratsubaMul($pp[$i + 1], $w);
    $vv    = array_chunk($vector, $kk);
    $vv[0] = mBackwardFFT($pp[2], ...$vv[0]);
    $vv[1] = mBackwardFFT($pp[2], ...$vv[1]);
    //
    for ($i = 0; $i < $kk; $i++)
    {
	$ww = $pp[$k - $i];
	$mm = mKaratsubaMul($ww, $vv[1][$i]);
	$vv[1][$i] = mSub($vv[0][$i], $mm);
	$vv[0][$i] = mAdd($vv[0][$i], $mm);
    }
    return array_merge(...$vv);
}
//
function mSSMul($p, $q)
{
    $m = mCmp($p, $q) < 0 ? $q:$p;
    $n = mLog($m)[0];
    $n = pow(2, round(log($n, 2)));
    $d = mAdd(mPow('2', $n), '1');
    //
    $k = 2*strlen($m) - 1;
    $k = pow(2, ceil(log($k, 2)));
    //
    $p = str_pad($p, $k, '0', STR_PAD_LEFT);
    $q = str_pad($q, $k, '0', STR_PAD_LEFT);
    $p = str_split(strrev($p));
    $q = str_split(strrev($q));
    // Get kth Root of Unity
    $w = '2';
    while (mCmp($w, $d) < 0)
    {
	$r = mMod(mPow($w, $k), $d);
	if (! mCmp($r, '1'))  break;
	$w =  mAdd($w, '1');
    }
    //
    $fp = mModVector(mForwardFFT($w, ...$p), $d);
    $fq = mModVector(mForwardFFT($w, ...$q), $d);
    $fr = mModPointwise($fp, $fq, $d);
    $rr = mModVector(mBackwardFFT($w, ...$fr), $d);
    $r  = '0';
    //
    for ($i = 0; $i < count($rr); $i++)
    {
	$rr[$i]	= mBaseDiv($rr[$i], $k)[0];
	$rr[$i].= str_repeat('0', $i);
	$r = mAdd($r, $rr[$i]);
    }
    //
    return $r;
}
//

• mReduce()
• mFastReduce()
  
• mConvert()
• mBinPow()
• mFastPow()
  

//
function mReduce($p, $u, $m)
{
    $h = strlen($m);
    $n = strlen($p);
    //
    $l = 1;
    for ($i = 0; $i < ($h / $l); $i++)
    {
        $q = substr($p, $n - $l*($i + 1), $l);
	$q = mKaratsubaMul($u, $q);
	$q = substr($q, -$l);
	$q = mKaratsubaMul($q, $m);
	$q.= str_repeat('0', $i * $l);
	$p = ltrim(mAdd($p, $q), '0');
    }
    //
    $p = substr($p, 0, -$h);
    if (strlen($p) > $h)
	$p = mSub($p, $m);
    return $p;
}
//
function mFastReduce($p, $u, $m)
{
    $l = strlen($m);
    $q = substr(mKaratsubaMul($u, $p), -$l);
    $r = ltrim(mAdd($p, mKaratsubaMul($q, $m)), '0');
    $r = substr($r, 0, -$l);
    //
    if (strlen($r) > $l)
	$r = mSub($r, $m);
    return $r;
}
//
function mConvert($p, $m)
{
    $p = mMod($p . str_repeat('0', strlen($m)), $m);
    return ltrim($p, '0');
}
//
function mBinPow($p, $q)
{
    $l = mLog($q)[1];
    $c = '1';
    //
    for ($j = 0; $j < strlen($l); ++$j)
    {
	$c = mKaratsubaMul($c, $c);
	if (($l[$j] - '0'))
	$c = mKaratsubaMul($c, $p);
    }
    //
    return ltrim($c, '0');
}
//
function mFastPow($p, $q, $m)
{
    // Calculate R^-1
    $r = '1' . str_repeat('0', strlen($m));
    $i = mModInverse($m, $r);
    $i = mSub($r, $i);
    //
    $l = mLog($q)[1];
    $p = mFastReduce($p, $i, $m);
    $c = mConvert('1', $m);
    //
    for ($j = 0; $j < strlen($l); ++$j)
    {
	$c = mFastReduce(mKaratsubaMul($c, $c), $i, $m);
	if (($l[$j] - '0'))
	$c = mFastReduce(mKaratsubaMul($c, $p), $i, $m);
    }
    //
    return mFastReduce($c, $i, $m);
}
//

• mToResidue()
• mFromResidue()

//
function mToResidue($p, ...$m)
{
    $l = count($m);
    //
    if ($l <= 2)
    {
	$a = array();
	for ($i = 0; $i < $l; $i++)
	    $a[] = mMod($p, $m[$i]);
	return $a;
    }
    //
    $n  = array_splice($m, intdiv($l, 2));
    $M1 = array_reduce($n, 'mKaratsubaMul', '1');
    $M0 = array_reduce($m, 'mKaratsubaMul', '1');
    //
    return array_merge(
	mToResidue(mMod($p, $M0), ...$m),
	mToResidue(mMod($p, $M1), ...$n)
    );
}
//
function mFromResidue($r, $m)
{
    $l = count($m);
    //
    if ($l <= 1)
        return $r[0];
    //
    $l = intdiv($l, 2);
    $s = array_splice($r, $l);
    $n = array_splice($m, $l);
    $a = mFromResidue($r, $m); 
    $b = mFromResidue($s, $n);
    //
    $M0 = array_reduce($m, 'mKaratsubaMul', '1');
    $M1 = array_reduce($n, 'mKaratsubaMul', '1');
    $M2 = mKaratsubaMul($M0, $M1);
    [,$u,$v] = mExtendedGcd($M0, $M1);
    //
    $L0 = mMod(mKaratsubaMul($u, $b), $M1);
    $L1 = mMod(mKaratsubaMul($v, $a), $M0);
    $r  = mAdd(mKaratsubaMul($L0, $M0), mKaratsubaMul($L1, $M1));
    if (mCmp($M2, $r) <= 0)
	$r = mSub($r, $M2);
    return $r;
}
//

• Decimal
• roundingMap

//
class Decimal
{
    public $mantissa;
    public $exponent;
    //
    function __construct($mantissa, $exponent)
    {
	$this->mantissa  = $mantissa;
	$this->exponent  = $exponent;
    }
}
//
define('ToZero',   0);
define('Nearest',  1);
define('FromZero', 2);
//
$roundingMap	 = array();
$roundingMap[0]  = array(0, 0, 0);
$roundingMap[1]  = array(0, 0, 1);
$roundingMap[10] = array(0, 0, 1);
$roundingMap[11] = array(0, 1, 1);
//

• mMantissa()
• mNormalize()
• mToDecimal()

//
function mMantissa($m, $precision)
{
    $j = strlen($m);
    list($a, $b) = array(true, '');
    //
    for ($i = 0; $i < $precision; $i++)
    {
	if (! $m) break;
	$m = ltrim(mKaratsubaMul($m, '2'), '0');
	//
	if ($j < strlen($m))
	{
	$m  = substr($m, 1);
	$b .= '1';
	if ($a)
	list($a, $i) = array(false, 0);
	continue;
	}
        $b .= '0';
    }
    //
    return $b;
}
//
function mNormalize($m)
{
    //
    $p = strpos($m, '1') + 1;
    return array(substr($m, $p), $p);
}
//
function mToDecimal($p, $precision, $bits = 11)
{
    list($a, $b) = explode('.', $p);
    $a = mLog($a)[1];
    $b = mMantissa($b, $precision);
    //
    list($m, $p) = mNormalize($a . $b);
    $m = str_pad($m, $precision, '0');
    $p = strlen($a) - $p;
    //
    return new Decimal($m, $p);
}

• mAddDecimal()
• mSubDecimal()
  
• mMulDecimal()
• mRound()
• mFinite()

//
function mAddDecimal($p, $q, $leading = true, $negate = false)
{
    if ($q->exponent < $p->exponent)
	list($p, $q) = [$q, $p];
    //
    $a = ($leading ? '1':'') . $p->mantissa;
    $b = ($leading ? '1':'') . $q->mantissa;
    //
    if (($delta = $q->exponent - $p->exponent))
    {
	$a = str_repeat('0', $delta) . $a;
	$b = str_pad($b, strlen($a), '0');
    }
    //
    $b = ($negate ? '-':'') . $b;
    list($m, $p) = mNormalize(mAdd($a, $b, 2));
    return new Decimal($m, $q->exponent);
}
//
function mSubDecimal($p, $q, $leading = true)
{
    return mAddDecimal($p, $q, $leading, true);
}
//
function mMulDecimal($p, $q, $leading = true)
{
    //
    $a = ($leading ? '1':'') . $p->mantissa;
    $b = ($leading ? '1':'') . $q->mantissa;
    $r = '0';
    //
    $j = strlen($b) - 1;
    for ($i = $j; $i >= 0; $i--)
    {
	if ($b[$i] != '1')
	    continue;
	//
	$a .= str_repeat('0', $j - $i);
	$r  = mAdd($r, $a, 2);
	$j  = $i;
    }
    //
    $r = @array_shift(mNormalize($r));
    return new Decimal($r, $p->exponent + $q->exponent);
}
//
function mRound(&$p, $n, $mode = ToZero)
{
    global $roundingMap;
    //
    if ($n >= strlen($p->mantissa))
        return;
    //
    $a = '1' . substr($p->mantissa, 0, $n);
    $b = substr($p->mantissa, $n);
    $b = str_pad($b, 3, '0', STR_PAD_LEFT);
    //
    $guard  = $b[0];
    $round  = $b[1];
    $sticky = substr_count(substr($b, 2), '1') ? '1':'0';
    $key    = (int) ($round . $sticky);
    //
    $roundingMap[1][1] = $guard - '0';
    if ($roundingMap[$key][$mode])
    $a = mAdd($a, '1', 2);
    //
    $p->mantissa = mNormalize($a)[0];
    return $p;
}
//
function mFinite($p, $bits = 11)
{
    //
    $bias = mSub(mBinPow(2, $bits - 1), '1');
    $e    = mLog(mAdd($bias, $p->exponent))[1];
    $e    = str_pad($e, $bits, '0', STR_PAD_LEFT);
    //
    return '0' .
	    $e .
	    $p->mantissa;
}

Reference:

Modern Computer Arithmetic - Richard Brent, Paul Zimmermann
On the Schonhage-Strassen Multiplication Algorithm - Lloyd Allison
Montgomery Reduction