How to calculate `sin(a^b)` where `a^b` can be a very large double number?

Question

I wonder if there is an algorithm for doing that.

Consider sin(60.9^100) for example,

How to make 60.9^100 - 2pi *N in 2 pi range.

Answer 1

Here are two semi-automatic ways. They still need manual configuration, depending on the input. Could perhaps be made fully automatic, but at least they provide some way at all (and results to check further methods with). Edit: posted a fully automatic on by now.

Using Python's decimal module and its pi recipe, apparently 60.9^100 % 2pi is about 0.4826 (which can then be given to sin). Results for computing with 180 to 290 digits precision (code is at the end):

180 0.52113386128181643243087541425218797675113893601959695702254815952665...
190 0.48262316221535366629016856570286348468626721388098119253199769818223...
200 0.48262316221828443267196371207773451732899712100881145938907295835606...
210 0.48262316221828443267246775208563277802202330286500415343966588161647...
220 0.48262316221828443267246775208566344687793590859019274697600998645752...
230 0.48262316221828443267246775208566344687793479998648411772237563268709...
240 0.48262316221828443267246775208566344687793479998648362494580989872984...
250 0.48262316221828443267246775208566344687793479998648362494580991864728...
260 0.48262316221828443267246775208566344687793479998648362494580991864728...
270 0.48262316221828443267246775208566344687793479998648362494580991864728...
280 0.48262316221828443267246775208566344687793479998648362494580991864728...
290 0.48262316221828443267246775208566344687793479998648362494580991864728...

Wolfram Alpha fails at it, computing "zero". But for exponent 30 it still shows a valid result, and we match:

WolframAlpha: 6.0148312092022347033088447399833343520115646793565705028401966310...
Mine:         6.01483120920223470330884473998333435201156467935657050284019663107410...

Another way, using the first 1001 digits of pi copied from some site, and using integers until the very end, gives 0.48262316221828444 (Try it online!):

a, b = 609, 100
pi = 31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989
print(a**b * 10**(1000-b) % (2*pi) / 10**1000)

This does the operations with large integers, "scaled up" by 10^1000, until the final downscaling division gives a float.

Third way, using Python's fraction, also resulting in 0.48262316221828444 (Try it online!):

from fractions import Fraction

a, b = Fraction('60.9'), 100
pi = Fraction('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989')
print(float(a**b % (2*pi)))

Code using decimal (Try it online!):

from decimal import *

a, b = '60.9', 100

def pi():
    """Compute Pi to the current precision.

    >>> print(pi())
    3.141592653589793238462643383

    """
    getcontext().prec += 2  # extra digits for intermediate steps
    three = Decimal(3)      # substitute "three=3.0" for regular floats
    lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
    while s != lasts:
        lasts = s
        n, na = n+na, na+8
        d, da = d+da, da+32
        t = (t * n) / d
        s += t
    getcontext().prec -= 2
    return +s               # unary plus applies the new precision

for prec in range(100, 300, 10):
    setcontext(Context(prec=prec, Emax=MAX_EMAX, Emin=MIN_EMIN))

    x = Decimal(a) ** b

    try:
        print(prec, str(x % (2*pi()))[:70] + '...')
    except:
        pass

Answer 2

Using arbitrary precision floats and 600 fractional bits for division precision leads to:

 a           = 60.89999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999996385120202345673823
 b           = 100
 pi2         = 6.28318530717958647692528676655900576839433879875021164194988918461563281257241799725606965068423413596429617302656461329418768921910116446345071881625696223490056820540387704221111928924589790986076393
 a^b         = 28955379630372944287405172428351524098232652637160696453610571103951984391040712391370032514984003111097712417114477560530145388673931838194276269928567117043004481189909961875857.75649514478978553782864231122006591641737915748789056749200619756952612910126893035405495135890665595449781429534911495188619709183484163296347386211939610629876033802394440562845818548327528930373773
 a^b mod pi2 = 0.48090443944555460103493387524923864059291388194583637380891303988368325373623569292130225575175455300500465326465846619858936543877561381898189967500531981888897878580530179993696579483550649680773777

Which will serve as reference for double computation:

pow(60.9,100) = 2.895537963037287827000000000000000000000e+178

As you can see the rounding error is at magnitude ~10^160 which is hugely bigger than the wanted modulo leading to completely off results...

Using power by squaring along with modulo (modpow) to keep the subresults small like in Blackgaurd's answer:

//---------------------------------------------------------------------------
double mod(double a,double p)
    {
    a/=p;
    a-=floor(a);
    a*=p;
    return a;
    }
//---------------------------------------------------------------------------
double modpow(double a,unsigned int b,double p)
    {
    int i;
    double d;
    a=mod(a,p);
    d=mod(1.0,p);
    for (i=0;i<32;i++)
        {
        d=mod(d*d,p);
        if (DWORD(b&0x80000000))
            {
            d=mod(d*a,p);
            }
        b<<=1;
        }
    return mod(d,p);
    }
//---------------------------------------------------------------------------

Looks slightly better but still way off as during each iteration we lost precision bits of mantissa:

modpow(60.9,100,2.0*M_PI) = 0.000000000004526729

However as Kelly Bundy points out modulo on floats can't be used inside or before the pow iterations as it breaks the math (unlike on integers) making this approach unusable.

So the only way (unless there is some nice math identity to be exploited for this) to compute this is to use much bigger bitwidth numbers (either fixed point or float). The result of 60.9^100 has 178 decadic digits so you need to use at least:

178/log10(2) = 591.3 bits

I chose 600 bits for mantissa in the computations above (using power by squaring without mod inside iteration). Beware you need to have PI constant in similar bitwidth in order to make this work.

So the answer to your question is: No I see no way to compute this on double and you have to use big numbers instead.

Which depends on language you want to implement this in...

[Edit1] here the stuff computed on different accuracy bitwidth:

  600 bits: a^b mod p = 0.48090443944543787059955763581639220817922063926527443008616966832526276298682799066248021870151726073022025792461359895510004764246892419475030025711547915571989225146813845497084456179528624932006095
 1000 bits: a^b mod p = 0.48262316221828443267246775208566344687793479998648362494580991864728550767313153551910251537428629715050683896930935124670316833678465659261329197846231478482857003648065425433322518882397051713086293
10000 bits: a^b mod p = 0.48262316221828443267246775208566344687793479998648362494580991864728550767313153551910251537428629715050683896930935124670316834694077218601842581288592371259021844778135739430850760935223120750400005

Answer 3

Full Python implementation combining my first answer's semi-autimatic decimal method with Spektre's observation for how to calculate the required precision.

Our input to sin is a^b % 2π, which lies in [0, 2π), so its integer part has one digit. Beyond that one digit, I want 15 fractional digits, as Python's float have about 16 digits precision. So I use 30 digits more than the integer part of a^b in order to get 30 fractional digits (around 15 would suffice, but a few more don't hurt, and I don't want to try to compute exactly how many are needed).

from math import log10, sin
from decimal import setcontext, Context, getcontext, Decimal

def bigsin(a, b):
    prec = int(log10(a)*b) + 30
    setcontext(Context(prec=prec, Emax=prec + 10))
    a = Decimal(str(a))
    return sin(float(a**b % (2*pi())))

def pi():
    """Compute Pi to the current precision.

    >>> print(pi())
    3.141592653589793238462643383

    """
    getcontext().prec += 2  # extra digits for intermediate steps
    three = Decimal(3)      # substitute "three=3.0" for regular floats
    lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
    while s != lasts:
        lasts = s
        n, na = n+na, na+8
        d, da = d+da, da+32
        t = (t * n) / d
        s += t
    getcontext().prec -= 2
    return +s               # unary plus applies the new precision

print(bigsin(60.9, 100))

Output (Try it online!):

0.4641043156966329

Note this differs from my first answer because here I computed the sine value, not just the moduloed power.

Also note that while I start with 60.9 as a float (which doesn't represent 60.9 exactly), the str turns that into the string 60.9, so that the Decimal then does represent 60.9 exactly.

Answer 4

Here is an example of the comment by @Marat in c++:

double dmod(double x, double y) {
    return x - (int)(x/y) * y;
}

double bin_pow(double a, double exp, double mod) {
    if (!exp) return 1.0;
    double ret = 1;
    while (exp > 0) {
        if (exp & 1) ret = dmod(ret * a, mod);
        a = dmod(a * a, mod);
        exp >>= 1;
    }
    return ret;
}