1 
    2 
    3 
    4 
    5 
    6 
    7 
    8 
    9 
   10 
   11 
   12 
   13   module Transcendentals
   14    (module ContinuedFractions, 
   15     exp1, expQ, expR, log10, logQ, logR, tanQ, tanR, cotR, atan1, atanQ, atanR)
   16   where
   17   import ContinuedFractions
   18 
   19 
   20 
   21 
   22 
   23 
   24 
   25 
   26 
   27 
   28 
   29 
   30 
   31 
   32 
   33 
   34 
   35 
   36 
   37 
   38 
   39 
   40 
   41 
   42 
   43 
   44 
   45 
   46 
   47 
   48 
   49 
   50 
   51 
   52 
   53 
   54 
   55 
   56 
   57 
   58 
   59 
   60 
   61 
   62   exp1CF :: [QRational]
   63   exp1CF = 1: f 0 where f n = (4*n+1): (-2): (-(4*n+3)): 2: f (n+1)
   64 
   65 
   66 
   67 
   68 
   69 
   70 
   71 
   72 
   73 
   74   expQ  :: QRational -> ContinuedFraction
   75   expQ  x = aaQ hom ts2
   76             where (ts1,ts2)   = splitAt 2 cf
   77                   hom         = absorbQ startH ts1
   78                   (startH,cf) = if abs x < 1 then
   79                                 (([1,0],[0,1]), oddQ x exp1CF) else
   80                                 (([qDenominator x,0],[0,qNumerator x]),
   81                                   evenQ (1/x) exp1CF) 
   82 
   83   exp1 :: ContinuedFraction
   84   exp1 = expQ (1%%1)
   85 
   86 
   87 
   88 
   89   expR  :: ContinuedFraction -> ContinuedFraction
   90   expR  x = if mightBeZero x
   91             then processQS x (splitAt 6 (oddR exp1CF))
   92             else quadraticAlgorithm ([0,1,1,0],[0,1,0,0]) x cf2
   93             where cf2 = processQS x (splitAt 6 (oddR (tail exp1CF)))
   94 
   95 
   96 
   97 
   98 
   99 
  100 
  101 
  102 
  103 
  104 
  105 
  106 
  107 
  108 
  109 
  110   logQ  :: QRational -> ContinuedFraction
  111   logQ  x = aaQ startH ts2
  112             where (ts1,ts2) = splitAt 2 (oddQ (-t*t) atan1CF)
  113                   t         = (x-1)/(x+1)
  114                   tn        = qNumerator t
  115                   td        = qDenominator t
  116                   startH    = absorbQ ([-2*td,0],[0,tn]) ts1
  117 
  118   log10 :: ContinuedFraction
  119   log10 = logQ (10%%1)
  120 
  121 
  122 
  123 
  124 
  125   logR  :: ContinuedFraction -> ContinuedFraction
  126   logR  x = quadraticAlgorithm ([2,0,0,0],[0,0,0,1]) t cf1
  127             where t   = algebraicAlgorithm ([1,(-1)],[1,1]) x
  128                   t2  = quadraticAlgorithm ([-1,0,0,0],[0,0,0,1]) t t
  129                   cf  = qsSetupAlt (drop 1 (oddR atan1CF)) t2
  130                   cf1 = quadraticAlgorithm ([0,1,0,0],[0,1,1,0]) t2 cf
  131 
  132 
  133 
  134 
  135 
  136 
  137 
  138 
  139 
  140 
  141 
  142 
  143   tan1CF :: [QRational]
  144   tan1CF = 0: f 0 where f n = (4*n+1) : (-(4*n+3)) : f (n+1)
  145 
  146   tanQ  :: QRational -> ContinuedFraction
  147   tanQ  x = if x /= 0 then aaQ startH ts2 else error "tanQ: 0 argument"
  148             where n         = qRound (abs (x/2)) + 1
  149                   (ts1,ts2) = splitAt (fromIntegral n) (bothQ x tan1CF)
  150                   startH    = absorbQ ([1,0],[0,1]) ts1
  151 
  152   tanR  :: ContinuedFraction -> ContinuedFraction
  153   tanR  x = quadraticAlgorithm ([0,1,0,0],[0,0,1,0]) cfd cfn
  154             where cf  = qsSetupAlt (drop 1 (oddR tan1CF)) x2
  155                   x2  = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  x
  156                   cfn = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  cf
  157                   cfd = quadraticAlgorithm ([0,1,1,0],[0,0,0,1]) x2 cf
  158 
  159   cotR  :: ContinuedFraction -> ContinuedFraction
  160   cotR  x = quadraticAlgorithm ([0,1,0,0],[0,0,1,0]) cfn cfd
  161             where cf  = qsSetupAlt (drop 1 (oddR tan1CF)) x2
  162                   x2  = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  x
  163                   cfn = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  cf
  164                   cfd = quadraticAlgorithm ([0,1,1,0],[0,0,0,1]) x2 cf
  165 
  166 
  167 
  168 
  169 
  170 
  171 
  172   atan1CF :: [QRational]
  173   atan1CF = 0: 1: 3: 5/4: f (8/9) 0
  174             where f :: QRational -> QRational -> [QRational]
  175                   f p n = (7/2 + 2*n)*p: (9/2 + 2*n)*p': f p'' (n+1)
  176                           where p'  = 1/(p*(n+2)*(n+2))
  177                                 p'' = 1/(p'*(n+5/2)*(n+5/2))
  178 
  179 
  180 
  181 
  182 
  183 
  184 
  185 
  186 
  187   atanQ :: QRational -> ContinuedFraction
  188   atanQ x
  189    = if x > 1  then
  190         quadraticAlgorithm ([0,-1,2,0] ,[0,0,0,1]) atan1 (cf (1/x)) else
  191      if x < -1 then
  192         quadraticAlgorithm ([0,-1,-2,0],[0,0,0,1]) atan1 (cf (1/x)) else
  193         cf x
  194      where
  195      cf x = aaQ startH ts2
  196             where
  197             (ts1,ts2) = splitAt 4 (oddQ (x*x) atan1CF)
  198             startH    = absorbQ ([qDenominator x,0],[0,qNumerator x]) ts1
  199 
  200 
  201 
  202 
  203 
  204   atan1 :: ContinuedFraction
  205   atan1 = atanQ (1%%1)
  206 
  207   atanR :: ContinuedFraction -> ContinuedFraction
  208   atanR []         = algebraicAlgorithm ([2,0],[0,1]) atan1
  209   atanR xs@(x:xs')
  210    = if x > 1 then
  211         quadraticAlgorithm ([0,-1,2,0] ,[0,0,0,1]) atan1 ys else
  212      if x < -1 then
  213         quadraticAlgorithm ([0,-1,-2,0],[0,0,0,1]) atan1 ys else
  214         atanRdirect xs
  215      where ys = atanRdirect (algebraicAlgorithm ([0,1],[1,0]) xs)
  216 
  217   atanRdirect :: ContinuedFraction -> ContinuedFraction
  218   atanRdirect x
  219    = quadraticAlgorithm ([0,0,1,0],[0,1,0,0]) cfn cfd
  220      where x2  = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x x
  221            cf  = qsSetupAlt (drop 1 (oddR atan1CF)) x2
  222            cfn = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  cf
  223            cfd = quadraticAlgorithm ([0,1,1,0],[0,0,0,1]) x2 cf
  224 
  225 
  226 
  227 
  228 
  229   mightBeZero :: ContinuedFraction -> Bool
  230   mightBeZero []     = False
  231   mightBeZero (x:xs) = x == 0
  232 
  233 
  234 
  235 
  236 
  237 
  238 
  239 
  240 
  241   aaQ :: Homography -> [QRational] -> ContinuedFraction
  242   aaQ homography@([n0,n1],[d0,d1]) xs
  243    = if n0*d1 /= n1*d0 then aaQNoCheck homography xs
  244      else if d1 /= 0   then rat2cf (n1%%d1)
  245      else if d0 /= 0   then rat2cf (n0%%d0) 
  246      else                   error "aaQ: undefined homography"
  247 
  248 
  249 
  250 
  251   aaQNoCheck homography@([n0,n1],[d0,d1]) []
  252    = rat2cf (n0%%d0)
  253   aaQNoCheck homography@([n0,n1],[d0,d1]) xs'@(x:xs)
  254    = case algebraicOutput homography x of
  255        Just o  -> o : aaQNoCheck ([d0,d1],          [n0-d0*o,n1-d1*o]) xs'
  256        Nothing ->     aaQNoCheck ([n0*n+n1*d,n0*d], [d0*n+d1*d,d0*d])  xs
  257      where (n,d) = (qNumerator x, qDenominator x)
  258 
  259 
  260 
  261 
  262   absorbQ :: Homography -> [QRational] -> Homography
  263   absorbQ homography        []     = homography
  264   absorbQ ([n0,n1],[d0,d1]) (q:qs) = absorbQ ([n0*n+n1*d,n0*d],
  265                                               [d0*n+d1*d,d0*d]) qs
  266                                      where n = qNumerator   q
  267                                            d = qDenominator q
  268 
  269 
  270 
  271 
  272   bothQ :: QRational -> [QRational] -> [QRational]
  273   bothQ x (x1:xs) = x1/x:   bothQ x xs
  274 
  275   evenQ :: QRational -> [QRational] -> [QRational]
  276   evenQ x (x1:x2:xs) = x1/x: x2:   evenQ x xs
  277 
  278   oddQ :: QRational  -> [QRational] -> [QRational]
  279   oddQ  x (x1:x2:xs) = x1:   x2/x: oddQ  x xs
  280 
  281 
  282 
  283 
  284 
  285 
  286 
  287 
  288   type State = (QRational, Homography)
  289 
  290 
  291 
  292 
  293 
  294 
  295 
  296 
  297 
  298 
  299 
  300 
  301 
  302 
  303   oddR :: [QRational] -> [State]
  304   oddR (x1:x2:xs) = (x1,([0,xn,xd,0],[xd,0,0,0])): oddR xs
  305                     where (xn,xd) = (qNumerator x2, qDenominator x2)
  306 
  307 
  308 
  309 
  310   processQS :: ContinuedFraction -> ([State],[State]) -> ContinuedFraction
  311   processQS x (hs1,hs2) = preProcess x hs1 (qsSetupAlt hs2 x)
  312 
  313 
  314 
  315 
  316 
  317   preProcess :: ContinuedFraction -> -- the argument CF
  318                 [State]           -> -- leading parts of infinite state
  319                 ContinuedFraction -> -- accumulating CF
  320                 ContinuedFraction    -- the resulting CF
  321   preProcess _ []         c = c
  322   preProcess x ((i,h):hs) c = algebraicAlgorithm ([ni,nd],[nd,0]) 
  323                               (quadraticAlgorithm h x (preProcess x hs c))
  324                               where ni = qNumerator i
  325                                     nd = qDenominator i
  326 
  327 
  328 
  329   qsSetupAlt :: [State] -> ContinuedFraction -> ContinuedFraction
  330   qsSetupAlt hs x = aaQ ([1,0],[0,1]) (qsProcess (hs, x))
  331 
  332 
  333 
  334 
  335 
  336 
  337 
  338 
  339 
  340 
  341 
  342 
  343   qsProcess :: ([State], ContinuedFraction) -> [QRational]
  344   qsProcess    (ihs, [])    = asProcess ihs
  345   qsProcess st@(ihs, (_:_)) = qsNext 1 st
  346 
  347 
  348 
  349 
  350 
  351 
  352 
  353 
  354 
  355   qsNext :: Int -> ([State], ContinuedFraction) -> [QRational]
  356   qsNext n (ihs, []) = asProcess (qsEnd ihs)
  357   qsNext n st@(((i,h):ihs), xs'@(x:xs))
  358    = if length outs > 1
  359         then init outs ++ qsNext (n+1) (ihs',xs')
  360         else qsNext (n+1) (qsIn x ((i,h): ihs), xs)
  361      where (outs, (ihs',_)) = qsTerm n st
  362 
  363 
  364 
  365 
  366 
  367 
  368   qsTerm :: Int -> ([State], ContinuedFraction)
  369            -> ([QRational],([State], ContinuedFraction))
  370   qsTerm n (ihs,xs)
  371    = (outs, (ihs',xs))
  372      where (front,back@((i,_):_)) = splitAt n ihs
  373            (outs,ihs') = foldr (qsTerm' (head xs)) ([i],back) front
  374 
  375 
  376 
  377   qsTerm' :: Integer -> State -> ([QRational], [State]) -> 
  378                                  ([QRational], [State])
  379   qsTerm' x ih@(i,h) (is,ihs)
  380    = (is',((last is',h''):ihs))
  381      where h' = qsAbsorb (init is) h
  382            i' = last is
  383            (is',h'') = qsEmit x i' ([i],h')
  384 
  385 
  386 
  387 
  388   qsAbsorb :: [QRational] -> Homography -> Homography
  389   qsAbsorb [] h = h
  390   qsAbsorb (i:is) h@([n0,n1,n2,n3],[d0,d1,d2,d3])
  391    = qsAbsorb is ([n0*n+n2*d, n1*n+n3*d, n0*d, n1*d],
  392                   [d0*n+d2*d, d1*n+d3*d, d0*d, d1*d])
  393      where n = qNumerator i
  394            d = qDenominator i
  395 
  396 
  397 
  398 
  399   qsEmit :: Integer -> QRational -> ([QRational], Homography) ->
  400                                   ([QRational], Homography)
  401   qsEmit x i ish@(is,h@([n0,n1,n2,n3],[d0,d1,d2,d3]))
  402    = case quadraticOutput h [x%%1,i] of
  403        Just o  -> qsEmit x i (is++[o%%1],
  404                               ([d0,d1,d2,d3],
  405                                [n0-d0*o, n1-d1*o, n2-d2*o, n3-d3*o]))
  406        Nothing -> ish
  407 
  408 
  409 
  410 
  411   qsIn :: Integer -> [State] -> [State]
  412   qsIn = map . input
  413          where input x (i,([n0,n1,n2,n3],[d0,d1,d2,d3]))
  414                  = (i, ([n0*x+n1, n0, n2*x+n3, n2],
  415                         [d0*x+d1, d0, d2*x+d3, d2]))
  416 
  417 
  418 
  419 
  420 
  421 
  422 
  423 
  424 
  425 
  426   qsEnd :: [State] -> [State]
  427   qsEnd = map input
  428           where input (i,([n0,_,n2,_],[d0,_,d2,_]))
  429                   = (i,([n0,n2],[d0,d2]))
  430 
  431 
  432 
  433   asProcess :: [State] -> [QRational]
  434   asProcess ihs = fst (head ihs) : asProcess (asNext ihs)
  435 
  436   asNext :: [State] -> [State]
  437   asNext ((i,h@([n0,n1],[d0,d1])):ihs'@((i',_):ihs))
  438    = case algebraicOutput h i' of
  439        Just o  -> (o%%1, ([d0,d1], [n0-d0*o,n1-d1*o])): ihs'
  440        Nothing -> asNext ((i,hIn i' h): asNext ihs')
  441      where hIn x ([n0,n1],[d0,d1]) = ([n0*xn+n1*xd,n0*xd],[d0*xn+d1*xd,d0*xd])
  442                                      where xn = qNumerator x
  443                                            xd = qDenominator x
  444 
  445 
  446 
  447 
  448