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

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.

CodePudding user response：

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).

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...

Can't check with WolframAlpha, as that fails at it, computing "zero". But for exponent 30 it still shows a good 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.

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

CodePudding user response：

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;
}