polynomials-bernstein (empty) → 1
raw patch · 6 files changed
+1780/−0 lines, 6 filesdep +basedep +vectorsetup-changed
Dependencies added: base, vector
Files
- Algebra/Polynomials/Bernstein.lhs +1009/−0
- Algebra/Polynomials/Numerical.hs +263/−0
- Algebra/Polynomials/cnumerical.c +143/−0
- LICENSE +340/−0
- Setup.hs +3/−0
- polynomials-bernstein.cabal +22/−0
+ Algebra/Polynomials/Bernstein.lhs view
@@ -0,0 +1,1009 @@+\begin{code}+{-# OPTIONS -XUnboxedTuples -XScopedTypeVariables -XFlexibleContexts -XBangPatterns -XUndecidableInstances -XMultiParamTypeClasses -XFunctionalDependencies -XFlexibleInstances -XMagicHash -XTypeFamilies #-}+-- | Various functions for manipulating polynomials, essentially when+-- represented in the Bernstein basis, in one or two variables.+module Algebra.Polynomials.Bernstein (Bernsteinp(..),solve,Bernstein(..),+ derivate,reorient) where+++import Control.Monad.ST+import Algebra.Polynomials.Numerical++import qualified Data.Vector.Unboxed as UV+import qualified Data.Vector.Unboxed.Mutable as MUV+import qualified Data.Vector as V+import Data.Vector.Generic as GV hiding ((++))+ +-- | The type for Bernstein polynomials with an arbitrary number of variables+data Bernsteinp a b=Bernsteinp { bounds::a, coefs::UV.Vector b } deriving (Eq,Show)++type family Param a b+type instance Param Int a=a+type instance Param (Int,Int) a=(a,a)+type instance Param (Int,Int,Int) a=(a,a,a)+type instance Param (Int,Int,Int,Int) a=(a,a,a,a)++type family Inter a b+type instance Inter Int a=(a,a)+type instance Inter (Int,Int) a=(a,a,a,a)+type instance Inter (Int,Int,Int) a=(a,a,a,a,a,a)+type instance Inter (Int,Int,Int,Int) a=(a,a,a,a,a,a,a,a)+++class Bernstein a where+ -- Constructeurs+ (?)::UV.Unbox b=>Bernsteinp a b->a->b+ constant::(UV.Unbox b,Num b, Fractional b)=>b->Bernsteinp a b+ scale::(Num b, Fractional b,UV.Unbox b)=>b->Bernsteinp a b->Bernsteinp a b+ scale a (Bernsteinp i b)=Bernsteinp i $ UV.map (*a) b+ promote::(UV.Unbox b,Num b, Fractional b)=>Int->Bernsteinp Int b->Bernsteinp a b+ elevate::(UV.Unbox b,Num b, Fractional b)=>a->Bernsteinp a b->Bernsteinp a b+ eval::(UV.Unbox b,Num b,Fractional b)=>Bernsteinp a b->Param a b->b+ restriction::(UV.Unbox b,Fractional b,Num b)=>Bernsteinp a b->Param a b->Param a b->Bernsteinp a b+++instance Num (Bernsteinp a Interval)=>Intervalize (Bernsteinp a) where+ intervalize (Bernsteinp i x)=Bernsteinp i $! UV.map interval x+ intersects bxf bxg=+ let (Bernsteinp _ xx)=bxf-bxg in+ UV.all (\(Interval a b)->a<=0 && b>=0) xx++binomials::(Num a, MUV.Unbox a)=>Int->UV.Vector a+binomials m=+ UV.create $ do+ v<-MUV.replicate ((m+1)*(m+1)) 0+ MUV.write v 0 1+ let fill i+ | i>=(m+1)*(m+1) = return v+ | i`mod`(m+1) == 0 = do+ MUV.write v i 1+ fill (i+1)+ | otherwise = do+ a<-MUV.read v (i-m-1)+ b<-MUV.read v (i-m-2)+ MUV.write v i (a+b)+ fill (i+1)+ fill (m+1)++ +instance Bernstein Int where+ (?) (Bernsteinp _ a) b=a!b+ constant x=Bernsteinp 1 $ GV.singleton x+ promote _=id+ elevate r (Bernsteinp n f)=+ if r<=0 then Bernsteinp n f else+ let coef j=+ let sumAll i result+ | (i>j) || (i>=n) = result+ | otherwise =+ sumAll (i+1) $ result+(f!i)*((bin i (n-1))*((bin (j-i) r))/(bin j (n+r-1)))+ in+ (sumAll 0 0)+ binomial=binomials $ n+r-1+ bin a b=binomial!(b*(n+r)+a)+ in+ Bernsteinp (n+r) $ UV.generate (n+r) coef+ eval (Bernsteinp n points) t=+ if n==0 then 0 else runST $ do+ arr<-thaw points+ let fill i s+ | s>=n = MUV.read arr 0+ | i>=n-s = fill 0 (s+1)+ | otherwise = do+ a<-MUV.read arr i+ b<-MUV.read arr (i+1)+ MUV.write arr i $ a*(1-t)+b*t+ fill (i+1) s+ fill 0 1+ + restriction (Bernsteinp n0 points) a b=+ runST $ do+ pf<-thaw points+ let casteljau bef aft nv u v i k s j+ | i>=bef = return ()+ | k>=aft= casteljau bef aft nv u v (i+1) 0 1 0+ | s>=nv = casteljau bef aft nv u v i (k+1) 1 0+ | j>=nv = casteljau bef aft nv u v i k (s+1) 0+ | otherwise=+ let idx i_ j_ k_=(i_*nv+j_)*aft + k_+ v'=(v-u)/(1-u)+ in+ if s+j<nv then do+ pfi0<-MUV.read pf (idx i j k)+ pfi1<-MUV.read pf (idx i (j+1) k)+ MUV.write pf (idx i j k) $ (1-u)*pfi0+u*pfi1+ casteljau bef aft nv u v i k s (j+1)+ else do+ pfi0<-MUV.read pf (idx i (j-1) k)+ pfi1<-MUV.read pf (idx i j k)+ MUV.write pf (idx i j k) $ (1-v')*pfi0+v'*pfi1+ casteljau bef aft nv u v i k s (j+1)+ casteljau 1 1 n0 a b 0 0 1 0+ pff<-unsafeFreeze pf+ return $ Bernsteinp n0 pff+ ++instance Bernstein (Int,Int) where+ (?) (Bernsteinp (_,b) c) (i,j)=c!(i*b+j)+ constant x=Bernsteinp (1,1) $ UV.singleton x+ + promote 1 (Bernsteinp i x)=Bernsteinp (i,1) x+ promote _ (Bernsteinp i x)=Bernsteinp (1,i) x++ elevate (ra_,rb_) (Bernsteinp (a,b) f)=+ let ra+ | ra_>0 = ra_+ | otherwise = 0+ rb+ | rb_>0 = rb_+ | otherwise = 0+ in+ if a<=0 || b<=0 then Bernsteinp (a+ra,b+rb) $ GV.replicate ((a+ra)*(b+rb)) 0 else+ let idx i j=(i*b)+j+ idx' i j=(i*(b+rb))+j+ vect=create $ do+ v<-MUV.new ((a+ra)*(b+rb))+ let coef i j+ | i>=(a+ra) = return v+ | j>=(b+rb) = coef (i+1) 0+ | otherwise=do+ let sumAll i' j' !result+ | i'>=a || i'>i = result+ | j'>=b || j'>j = sumAll (i'+1) 0 result+ | otherwise =+ let x0=(bin i' (a-1))*((bin (i-i') ra)/(bin i (a+ra-1)))+ x1=(bin j' (b-1))*((bin (j-j') rb)/(bin j (b+rb-1)))+ in+ sumAll i' (j'+1) $! result+x0*x1*(f!(idx i' j'))+ MUV.write v (idx' i j) $ sumAll 0 0 0+ coef i (j+1)+ coef 0 0++ m=max (a+ra-1) (b+rb-1)+ bin i j=binomial!(j*(m+1)+i)+ binomial=binomials m+ in+ Bernsteinp (a+ra,b+rb) vect+ eval (Bernsteinp (n0,n1) points) (a,b)=+ if n0<=0 || n1<=0 then 0 else+ runST $ do+ pf<-thaw points+ let casteljau p0 p1 u i j s+ | i>=p0 = return ()+ | s>=p1 = casteljau p0 p1 u (i+1) 0 1+ | j>=(p1-s) = casteljau p0 p1 u i 0 (s+1)+ | otherwise = do+ x0<-MUV.read pf $ i+p0*j+ x1<-MUV.read pf $ i+p0*(j+1)+ MUV.write pf (i+p0*j) $ (1-u)*x0+u*x1+ casteljau p0 p1 u i (j+1) s+ casteljau n1 n0 a 0 0 1+ casteljau 1 n1 b 0 0 1+ MUV.read pf 0++ restriction (Bernsteinp (n0,n1) points) (a,c) (b,d)=+ runST $ do+ pf<-thaw points+ let casteljau bef aft nv u v i k s j+ | i>=bef = return ()+ | k>=aft= casteljau bef aft nv u v (i+1) 0 1 0+ | s>=nv = casteljau bef aft nv u v i (k+1) 1 0+ | j>=nv = casteljau bef aft nv u v i k (s+1) 0+ | otherwise=+ let idx i_ j_ k_=(i_*nv+j_)*aft + k_+ v'=(v-u)/(1-u)+ in+ if s+j<nv then do+ pfi0<-MUV.read pf (idx i j k)+ pfi1<-MUV.read pf (idx i (j+1) k)+ MUV.write pf (idx i j k) $ (1-u)*pfi0+u*pfi1+ casteljau bef aft nv u v i k s (j+1)+ else do+ pfi0<-MUV.read pf (idx i (j-1) k)+ pfi1<-MUV.read pf (idx i j k)+ MUV.write pf (idx i j k) $ (1-v')*pfi0+v'*pfi1+ casteljau bef aft nv u v i k s (j+1)+ casteljau 1 n1 n0 a b 0 0 1 0+ casteljau n0 1 n1 c d 0 0 1 0+ pff<-unsafeFreeze pf+ return $ Bernsteinp (n0,n1) pff++instance Bernstein (Int,Int,Int) where+ + (?) (Bernsteinp (_,b,c) d) (i,j,k)=d!(((i*b+j)*c)+k)+ constant x=Bernsteinp (1,1,1) $ UV.singleton x+ + promote 1 (Bernsteinp i x)=Bernsteinp (i,1,1) x+ promote 2 (Bernsteinp i x)=Bernsteinp (1,i,1) x+ promote _ (Bernsteinp i x)=Bernsteinp (1,1,i) x++ elevate (ra_,rb_,rc_) (Bernsteinp (a,b,c) f)=+ let ra+ | ra_>0 = ra_+ | otherwise = 0+ rb+ | rb_>0 = rb_+ | otherwise = 0+ rc+ | rc_>0 = rc_+ | otherwise = 0+ in+ if a<=0 || b<=0 || c<=0 then + Bernsteinp (a+ra,b+rb,c+rc) $ GV.replicate ((a+ra)*(b+rb)*(c+rc)) 0+ else+ let idx i j k=((i*b)+j)*c+k+ idx' i j k=((i*(b+rb))+j)*(c+rc)+k+ vect=create $ do+ v<-MUV.new ((a+ra)*(b+rb)*(c+rc))+ let coef i j k+ | i>=a+ra = return v+ | j>=b+rb = coef (i+1) 0 0+ | k>=c+rc = coef i (j+1) 0+ | otherwise=do+ let sumAll i' j' k' !result+ | i'>=a || i'>i = result+ | j'>=b || j'>j = sumAll (i'+1) 0 0 result+ | k'>=c || k'>k = sumAll i' (j'+1) 0 result+ | otherwise =+ let x0=(bin i' (a-1))*((bin (i-i') ra)/(bin i (a+ra-1)))+ x1=(bin j' (b-1))*((bin (j-j') rb)/(bin j (b+rb-1)))+ x2=(bin k' (c-1))*((bin (k-k') rc)/(bin k (c+rc-1)))+ in+ sumAll i' j' (k'+1) $! result+x0*x1*x2*(f!(idx i' j' k'))+ MUV.write v (idx' i j k) $ sumAll 0 0 0 0+ coef i j (k+1)+ coef 0 0 0++ m=max (max (a+ra-1) (b+rb-1)) (c+rc-1)+ bin i j=binomial!(j*(m+1)+i)+ binomial=binomials m+ in+ Bernsteinp (a+ra,b+rb,c+rc) vect+ eval (Bernsteinp (n0,n1,n2) points) (a,b,c)=+ if n0<=0 || n1<=0 || n2<=0 then 0 else+ runST $ do+ pf<-thaw points+ let casteljau p0 p1 u i j s+ | i>=p0 = return ()+ | s>=p1 = casteljau p0 p1 u (i+1) 0 1+ | j>=(p1-s) = casteljau p0 p1 u i 0 (s+1)+ | otherwise = do+ x0<-MUV.read pf $ i+p0*j+ x1<-MUV.read pf $ i+p0*(j+1)+ MUV.write pf (i+p0*j) $ (1-u)*x0+u*x1+ casteljau p0 p1 u i (j+1) s+ casteljau (n1*n2) n0 a 0 0 1+ casteljau n2 n1 b 0 0 1+ casteljau 1 n2 c 0 0 1+ MUV.read pf 0+ restriction (Bernsteinp (n0,n1,n2) points) (a,c,e) (b,d,f)=+ runST $ do+ pf<-thaw points+ let casteljau bef aft nv u v i k s j+ | i>=bef = return ()+ | k>=aft= casteljau bef aft nv u v (i+1) 0 1 0+ | s>=nv = casteljau bef aft nv u v i (k+1) 1 0+ | j>=nv = casteljau bef aft nv u v i k (s+1) 0+ | otherwise=+ let idx i_ j_ k_=(i_*nv+j_)*aft + k_+ v'=(v-u)/(1-u)+ in+ if s+j<nv then do+ pfi0<-MUV.read pf (idx i j k)+ pfi1<-MUV.read pf (idx i (j+1) k)+ MUV.write pf (idx i j k) $ (1-u)*pfi0+u*pfi1+ casteljau bef aft nv u v i k s (j+1)+ else do+ pfi0<-MUV.read pf (idx i (j-1) k)+ pfi1<-MUV.read pf (idx i j k)+ MUV.write pf (idx i j k) $ (1-v')*pfi0+v'*pfi1+ casteljau bef aft nv u v i k s (j+1)+ casteljau 1 (n1*n2) n0 a b 0 0 1 0+ casteljau n0 n2 n1 c d 0 0 1 0+ casteljau (n0*n1) 1 n2 e f 0 0 1 0+ pff<-unsafeFreeze pf+ return $ Bernsteinp (n0,n1,n2) pff++instance Bernstein (Int,Int,Int,Int) where+ + (?) (Bernsteinp (_,b,c,d) e) (i,j,k,l)=e!((((i*b+j)*c)+k)*d+l)+ constant x=Bernsteinp (1,1,1,1) $ UV.singleton x+ + promote 1 (Bernsteinp i x)=Bernsteinp (i,1,1,1) x+ promote 2 (Bernsteinp i x)=Bernsteinp (1,i,1,1) x+ promote 3 (Bernsteinp i x)=Bernsteinp (1,1,i,1) x+ promote _ (Bernsteinp i x)=Bernsteinp (1,1,1,i) x++ elevate (ra_,rb_,rc_,rd_) (Bernsteinp (a,b,c,d) f)=+ + let ra+ | ra_>0 = ra_+ | otherwise = 0+ rb+ | rb_>0 = rb_+ | otherwise = 0+ rc+ | rc_>0 = rc_+ | otherwise = 0+ rd+ | rd_>0 = rd_+ | otherwise = 0+ in + if a<=0 || b<=0 || c<=0 || d<=0 then + Bernsteinp (a+ra,b+rb,c+rc,d+rd) $ GV.replicate ((a+ra)*(b+rb)*(c+rc)*(d+rd)) 0+ else+ let idx i j k l=(((i*b)+j)*c+k)*d+l+ idx' i j k l=(((i*(b+rb))+j)*(c+rc)+k)*(d+rd)+l+ vect=create $ do+ v<-MUV.new ((a+ra)*(b+rb)*(c+rc)*(d+rd))+ let coef i j k l+ | i>=a+ra = return v+ | j>=b+rb = coef (i+1) 0 0 0+ | k>=c+rc = coef i (j+1) 0 0+ | l>=d+rd = coef i j (k+1) 0+ | otherwise=do+ let sumAll i' j' k' l' !result+ | i'>=a || i'>i = result+ | j'>=b || j'>j = sumAll (i'+1) 0 0 0 result+ | k'>=c || k'>k = sumAll i' (j'+1) 0 0 result+ | l'>=d || l'>l = sumAll i' j' (k'+1) 0 result+ | otherwise =+ let x0=(bin i' (a-1))*((bin (i-i') ra)/(bin i (a+ra-1)))+ x1=(bin j' (b-1))*((bin (j-j') rb)/(bin j (b+rb-1)))+ x2=(bin k' (c-1))*((bin (k-k') rc)/(bin k (c+rc-1)))+ x3=(bin l' (d-1))*((bin (l-l') rd)/(bin l (d+rd-1)))+ in+ sumAll i' j' k' (l'+1) $! result+x0*x1*x2*x3*(f!(idx i' j' k' l'))+ MUV.write v (idx' i j k l) $ sumAll 0 0 0 0 0+ coef i j k (l+1)+ coef 0 0 0 0++ m=max (max (a+ra) (b+rb)) (max (c+rc) (d+rd))+ bin i j=binomial!(j*(m+1)+i)+ binomial=binomials m+ in+ Bernsteinp (a+ra,b+rb,c+rc,d+rd) vect+ + eval (Bernsteinp (n0,n1,n2,n3) points) (a,b,c,d)=+ if n0<=0 || n1<=0 || n2<=0 || n3<=0 then 0 else+ runST $ do+ pf<-thaw points+ let casteljau p0 p1 u i j s+ | i>=p0 = return ()+ | s>=p1 = casteljau p0 p1 u (i+1) 0 1+ | j>=(p1-s) = casteljau p0 p1 u i 0 (s+1)+ | otherwise = do+ x0<-MUV.read pf $ i+p0*j+ x1<-MUV.read pf $ i+p0*(j+1)+ MUV.write pf (i+p0*j) $ (1-u)*x0+u*x1+ casteljau p0 p1 u i (j+1) s+ casteljau (n1*n2*n3) n0 a 0 0 1+ casteljau (n2*n3) n1 b 0 0 1+ casteljau n3 n2 c 0 0 1+ casteljau 1 n3 d 0 0 1+ MUV.read pf 0+ + restriction (Bernsteinp (n0,n1,n2,n3) points) (a,c,e,g) (b,d,f,h)=+ runST $ do+ pf<-thaw points+ let casteljau bef aft nv u v i k s j+ | i>=bef = return ()+ | k>=aft= casteljau bef aft nv u v (i+1) 0 1 0+ | s>=nv = casteljau bef aft nv u v i (k+1) 1 0+ | j>=nv = casteljau bef aft nv u v i k (s+1) 0+ | otherwise=+ let idx i_ j_ k_=(i_*nv+j_)*aft + k_+ v'=(v-u)/(1-u)+ in+ if s+j<nv then do+ pfi0<-MUV.read pf (idx i j k)+ pfi1<-MUV.read pf (idx i (j+1) k)+ MUV.write pf (idx i j k) $ (1-u)*pfi0+u*pfi1+ casteljau bef aft nv u v i k s (j+1)+ else do+ pfi0<-MUV.read pf (idx i (j-1) k)+ pfi1<-MUV.read pf (idx i j k)+ MUV.write pf (idx i j k) $ (1-v')*pfi0+v'*pfi1+ casteljau bef aft nv u v i k s (j+1)+ casteljau 1 (n1*n2*n3) n0 a b 0 0 1 0+ casteljau n0 (n2*n3) n1 c d 0 0 1 0+ casteljau (n0*n1) n3 n2 e f 0 0 1 0+ casteljau (n0*n1*n2) 1 n3 g h 0 0 1 0+ pff<-unsafeFreeze pf+ return $ Bernsteinp (n0,n1,n2,n3) pff+ ++instance (Num a,Fractional a,MUV.Unbox a)=>Num (Bernsteinp Int a) where++ (+) bf@(Bernsteinp m _) bg@(Bernsteinp n _)=+ let (Bernsteinp m' f')=elevate (n-m) bf+ (Bernsteinp _ g')=elevate (m-n) bg+ in+ Bernsteinp m' $ UV.generate m' $ \i->f'!i+g'!i++++ (*) ff@(Bernsteinp (af) _) gg@(Bernsteinp (ag) _)=+ if af<=0 || ag<=0 then+ Bernsteinp 0 UV.empty+ else+ let mm=(af+ag)-2+ binomial=binomials mm+ bin a b=binomial!(b*(mm+1)+a)+ in+ Bernsteinp (af+ag-1) $ create $ do+ v<-MUV.new $ af+ag-1+ let fill i+ | i>=af+ag-1 = return v+ | otherwise = do+ let mCoef' i' result+ | i'>i || i'>=af = + result+ | otherwise =+ let a=((bin i' (af-1))*(bin (i-i') (ag-1)))/(bin i (af+ag-2)) in+ mCoef' (i'+1) $+ result+a*(ff?i')*(gg?(i-i'))+ + MUV.write v i $! mCoef' (max 0 (i-ag+1)) 0+ fill (i+1)+ fill 0++ (-) bf (Bernsteinp i g)= bf+(Bernsteinp i $ UV.map negate g)++ signum _=error "No signum operation for Bernstein1"+ abs _=error "No abs operation for Bernstein1"++ fromInteger x=Bernsteinp 1 $ UV.singleton $ fromIntegral x+ +instance (Fractional a, Num a,UV.Unbox a)=>Num (Bernsteinp (Int,Int) a) where++ (+) bff@(Bernsteinp (af,bf) _) bgg@(Bernsteinp (ag,bg) _)=+ let (Bernsteinp (a,b) f')=elevate (ag-af,bg-bf) bff+ (Bernsteinp _ g')=elevate (af-ag,bf-bg) bgg+ in+ Bernsteinp (a,b) $ UV.generate (a*b) $ \i->f'!i+g'!i+++ (*) ff@(Bernsteinp (af,bf) _) gg@(Bernsteinp (ag,bg) _)=+ if af<=0 || bf<=0 || ag<=0 || bg<=0 then+ Bernsteinp (0,0) UV.empty+ else+ let mm=max (af+ag) (bf+bg)-2+ binomial=binomials mm+ bin a b=binomial!(b*(mm+1)+a)+ in+ Bernsteinp (af+ag-1,bf+bg-1) $ create $ do+ v<-MUV.new $ (af+ag-1)*(bf+bg-1)+ let idx i j=i*(bf+bg-1)+j+ fill i j+ | i>=af+ag-1 = return v+ | j>=bf+bg-1 = fill (i+1) 0+ | otherwise =+ do+ let mCoef' i' j' result+ | i'>i || i'>=af = + let b=(bin i (af+ag-2))*(bin j (bf+bg-2)) in+ result/b+ | j'>j || j'>=bf = mCoef' (i'+1) (max 0 (j-bg+1)) result+ | otherwise =+ let a=(bin i' (af-1))*(bin (i-i') (ag-1))*+ (bin j' (bf-1))*(bin (j-j') (bg-1))+ in+ mCoef' i' (j'+1) $+ result+a*(ff?(i',j'))*(gg?(i-i',j-j'))+ + MUV.write v (idx i j) $!+ mCoef' (max 0 (i-ag+1)) (max 0 (j-bg+1)) 0+ fill i (j+1)+ fill 0 0+ + (-) bf bg=bf+(bg { coefs=UV.map negate $ coefs bg })++ signum _=error "No signum operation for Bernstein1"+ abs _=error "No abs operation for Bernstein1"++ fromInteger x=Bernsteinp (1,1) $ UV.singleton $ fromIntegral x+instance (Fractional a, Num a, UV.Unbox a)=>Num (Bernsteinp (Int,Int,Int) a) where++ (+) bff@(Bernsteinp (af,bf,cf) _) bgg@(Bernsteinp (ag,bg,cg) _)=+ let (Bernsteinp (a,b,c) f')=elevate (ag-af,bg-bf,cg-cf) bff+ (Bernsteinp _ g')=elevate (af-ag,bf-bg,cf-cg) bgg+ in+ Bernsteinp (a,b,c) $ UV.generate (a*b*c) $ \i->f'!i+g'!i++ (*) ff@(Bernsteinp (af,bf,cf) _) gg@(Bernsteinp (ag,bg,cg) _)=+ if af<=0 || bf<=0 || cf<=0 || ag<=0 || bg<=0 || cg<=0 then+ Bernsteinp (0,0,0) UV.empty+ else+ let mm=(max (max (af+ag) (bf+bg)) (cf+cg))-2+ binomial=binomials mm+ bin a b=binomial!(b*(mm+1)+a)+ in+ Bernsteinp (af+ag-1,bf+bg-1,cf+cg-1) $ create $ do+ v<-MUV.new $ (af+ag-1)*(bf+bg-1)*(cf+cg-1)+ let idx i j k=(i*(bf+bg-1)+j)*(cf+cg-1)+k+ fill i j k+ | i>=af+ag-1 = return v+ | j>=bf+bg-1 = fill (i+1) 0 0+ | k>=cf+cg-1 = fill i (j+1) 0+ | otherwise =+ do+ let mCoef' i' j' k' result+ | i'>i || i'>=af = + let b=(bin i (af+ag-2))*(bin j (bf+bg-2))*+ (bin k (cf+cg-2))+ in+ result/b+ | j'>j || j'>=bf = mCoef' (i'+1) (max 0 (j-bg+1)) (max 0 (k-cg+1)) result+ | k'>k || k'>=cf = mCoef' i' (j'+1) (max 0 (k-cg+1)) result+ | otherwise =+ let a=(bin i' (af-1))*(bin (i-i') (ag-1))*+ (bin j' (bf-1))*(bin (j-j') (bg-1))*+ (bin k' (cf-1))*(bin (k-k') (cg-1))+ in+ mCoef' i' j' (k'+1) $+ result+a*(ff?(i',j',k'))*(gg?(i-i',j-j',k-k'))+ + MUV.write v (idx i j k) $!+ mCoef' (max 0 (i-ag+1)) (max 0 (j-bg+1)) (max 0 (k-cg+1)) 0+ fill i j (k+1)+ fill 0 0 0++ (-) bf bg=bf+(bg { coefs=UV.map negate $ coefs bg })++ signum _=error "No signum operation for Bernstein1"+ abs _=error "No abs operation for Bernstein1"++ fromInteger x=Bernsteinp (1,1,1) $ UV.singleton $ fromIntegral x+instance (Fractional a, Num a,UV.Unbox a)=>Num (Bernsteinp (Int,Int,Int,Int) a) where++ (+) bff@(Bernsteinp (af,bf,cf,df) _) bgg@(Bernsteinp (ag,bg,cg,dg) _)=+ let (Bernsteinp (a,b,c,d) f')=elevate (ag-af,bg-bf,cg-cf,dg-df) bff+ (Bernsteinp _ g')=elevate (af-ag,bf-bg,cf-cg,df-dg) bgg+ in+ Bernsteinp (a,b,c,d) $ UV.generate (a*b*c*d) $ \i->f'!i+g'!i++ (*) ff@(Bernsteinp (af,bf,cf,df) _) gg@(Bernsteinp (ag,bg,cg,dg) _)=+ if af<=0 || bf<=0 || cf<=0 || df<=0 || ag<=0 || bg<=0 || cg<=0 || dg<=0 then+ Bernsteinp (0,0,0,0) UV.empty+ else+ let mm=(max (max (af+ag) (bf+bg)) (max (cf+cg) (df+dg)))-2+ binomial=binomials mm+ bin a b=binomial!(b*(mm+1)+a)+ in+ Bernsteinp (af+ag-1,bf+bg-1,cf+cg-1,df+dg-1) $ create $ do+ v<-MUV.new $ (af+ag-1)*(bf+bg-1)*(cf+cg-1)*(df+dg-1)+ let idx i j k l=((i*(bf+bg-1)+j)*(cf+cg-1)+k)*(df+dg-1)+l+ fill i j k l+ | i>=af+ag-1 = return v+ | j>=bf+bg-1 = fill (i+1) 0 0 0+ | k>=cf+cg-1 = fill i (j+1) 0 0+ | l>=df+dg-1 = fill i j (k+1) 0+ | otherwise =+ do+ let mCoef' i' j' k' l' result+ | i'>i || i'>=af = + let b=(bin i (af+ag-2))*(bin j (bf+bg-2))*+ (bin k (cf+cg-2))*(bin l (df+dg-2))+ in+ result/b+ | j'>j || j'>=bf = mCoef' (i'+1) (max 0 (j-bg+1)) (max 0 (k-cg+1)) (max 0 (l-dg+1)) result+ | k'>k || k'>=cf = mCoef' i' (j'+1) (max 0 (k-cg+1)) (max 0 (l-dg+1)) result+ | l'>l || l'>=df = mCoef' i' j' (k'+1) (max 0 (l-dg+1)) result+ | otherwise =+ let a=(bin i' (af-1))*(bin (i-i') (ag-1))*+ (bin j' (bf-1))*(bin (j-j') (bg-1))*+ (bin k' (cf-1))*(bin (k-k') (cg-1))*+ (bin l' (df-1))*(bin (l-l') (dg-1))+ in+ mCoef' i' j' k' (l'+1) $+ result+a*(ff?(i',j',k',l'))*(gg?(i-i',j-j',k-k',l-l'))+ + MUV.write v (idx i j k l) $!+ mCoef' (max 0 (i-ag+1)) (max 0 (j-bg+1)) (max 0 (k-cg+1)) (max 0 (l-dg+1)) 0+ fill i j k (l+1)+ fill 0 0 0 0+ (-) bf bg=bf+(bg { coefs=UV.map negate $ coefs bg })++ signum _=error "No signum operation for Bernstein1"+ abs _=error "No abs operation for Bernstein1"++ fromInteger x=Bernsteinp (1,1,1,1) $ UV.singleton $ fromIntegral x++-- | Computes the derivative of a univariate Bernstein polynomial.+derivate::(UV.Unbox a,Num a)=>Bernsteinp Int a->Bernsteinp Int a+derivate (Bernsteinp n f)+ | n<=1 = Bernsteinp 0 $ UV.empty+ | otherwise=Bernsteinp (n-1) $ UV.generate (n-1) (\i->(f!(i+1)-f!i)*(fromIntegral $ n-1))++-- | Computes @f(1-x)@ (useful when used with Bezier curves).+reorient::(UV.Unbox a)=>Bernsteinp Int a->Bernsteinp Int a+reorient (Bernsteinp n f)=Bernsteinp n (UV.reverse f)++\end{code}++\begin{code}+{-+restrict::Int->Int->Int->Bernsteinp i Interval->Double->Double->Bernsteinp i Interval+restrict bef aft nv (Bernsteinp ix poly) a b=+ --traceShow "Restrict" $ + runST $ do+ poly'<-thaw poly :: ST s (MUV.STVector s Interval)+ casteljau poly' 0 0 1 0+ unsafeFreeze poly' >>= return.(Bernsteinp ix)+ + where+ + + (# bl,bu #)=+ let (# b0,b1 #)=minus b b a a+ (# b2,b3 #)=minus 1 1 a a+ in+ over b0 b1 b2 b3++ idx i j k= --traceShow (i,j,k,d) $+ (i*nv+j)*aft + k+ + casteljau::MUV.STVector s Interval->+ Int->Int->Int->Int->ST s ()+ + casteljau pf i k s j+ | i>=bef = return ()+ | k>=aft= casteljau pf (i+1) 0 1 0+ | s>=nv = casteljau pf i (k+1) 1 0+ | j>=nv = casteljau pf i k (s+1) 0+ + -- Au boulot+ | s+j<nv = do+ (Interval l1 u1)<-MUV.read pf (idx i j k)+ (Interval l3 u3)<-MUV.read pf (idx i (j+1) k)+ let (# l0,u0 #)=minus 1 1 a a+ (# l2,u2 #)=times l0 u0 l1 u1+ (# l4,u4 #)=times a a l3 u3+ (# l5,u5 #)=plus l2 u2 l4 u4+ MUV.write pf (idx i j k) (Interval l5 u5)+ casteljau pf i k s (j+1)+ + | otherwise = do+ (Interval l1 u1)<-MUV.read pf (idx i (j-1) k)+ (Interval l3 u3)<-MUV.read pf (idx i j k)+ let (# l0,u0 #)=minus 1 1 bl bu+ (# l2,u2 #)=times l0 u0 l1 u1+ (# l4,u4 #)=times bl bu l3 u3+ (# l5,u5 #)=plus l2 u2 l4 u4+ MUV.write pf (idx i j k) (Interval l5 u5)+ casteljau pf i k s (j+1)+-}+{-+restrict::Int->Int->Int->Bernsteinp i a->a->a->Bernsteinp i a+restrict bef aft nv (Bernsteinp ix poly) a b=+ --traceShow "Restrict" $ + runST $ do+ poly'<-thaw poly+ casteljau poly' 0 0 1 0+ unsafeFreeze poly' >>= return.(Bernsteinp ix)+ + where+ + + casteljau pf bef aft nv a b i k s j+ | i>=bef = return ()+ | k>=aft= casteljau pf bef aft nv a b (i+1) 0 1 0+ | s>=nv = casteljau pf bef aft nv a b i (k+1) 1 0+ | j>=nv = casteljau pf bef aft nv a b i k (s+1) 0+ | otherwise=+ let idx i j k=(i*nv+j)*aft + k + b'=(b-a)/(1-a)+ in+ if s+j<nv then do+ pfi0<-MUV.read pf (idx i j k)+ pfi1<-MUV.read pf (idx i (j+1) k)+ MUV.write pf (idx i j k) $ (1-a)*pfi0+a*pfi1+ casteljau pf bef aft nv a b i k s (j+1)+ else do+ pfi0<-MUV.read pf (idx i (j-1) k)+ pfi1<-MUV.read pf (idx i j k)+ MUV.write pf (idx i j k) $ (1-b')*pfi0+b*pfi1+ casteljau pf bef aft nv a b i k s (j+1)+-}++-- Le booleen veut dire "tous les coefs sont nuls"+convexHull::Int->Int->Int->Bernsteinp i Interval->Double->Double->(Bool,Double,Double)+convexHull bef aft nv (Bernsteinp _ points) a b=+ let (allzero,pointsl,pointsu)=runST $ do+ let idx i j k=(i*nv+j)*aft + k+ pl<-MUV.replicate nv (1/0)+ pu<-MUV.replicate nv (-1/0)+ let fill i j k allzero_ -- a avant, b courant, c apres+ | i>=bef = return allzero_+ | j>=nv = fill (i+1) 0 0 allzero_+ | k>=aft = fill i (j+1) 0 allzero_+ | otherwise = do+ pl0<-MUV.read pl j+ pu0<-MUV.read pu j+ MUV.write pl j (min pl0 $ ilow $ points!(idx i j k))+ MUV.write pu j (max pu0 $ iup $ points!(idx i j k))+ fill i j (k+1) (allzero_ && pl0<=0 && pu0>=0)+ allzero_<-fill 0 0 0 True+ pl'<-UV.unsafeFreeze pl+ pu'<-UV.unsafeFreeze pu+ return (allzero_,pl',pu')+ inter::Int->Int->(Double,Double)+ inter i j+ | i>j = inter j i+ | otherwise =+ let yli=pointsl!i+ yui=pointsu!i+ ylj=pointsl!j+ yuj=pointsu!j+ fi=fromIntegral i+ fj=fromIntegral j+ inter0 yi yj=+ let k=fromIntegral $ i-j in+ Interval fi fi + + (Interval yi yi)*(Interval k k)/+ (Interval yj yj-Interval yi yi)+ in+ if yli<=0 then+ if yui>=0 then+ if ylj<=0 then+ if yuj>=0 then+ -- 1 les deux sont sur la ligne+ --traceShow "1" $+ (fi,fj)+ else+ -- 2 M est sur la ligne, M' est en-dessous+ --traceShow "2" $+ (fi, iup $ inter0 yui yuj)+ else+ -- 3 M est sur la ligne, M' est au-dessus+ --traceShow "3" $+ (fi,iup $ inter0 yli ylj)+ else+ -- M est en-dessous de la ligne + if ylj<=0 then+ if yuj>=0 then+ -- 4 M' est sur la ligne+ --traceShow "4" $+ (ilow $ inter0 yui yuj, fj)+ else+ -- 5 M' est en-dessous (comme M)+ --traceShow "5" $+ (1/0,-1/0)+ else+ -- 6 M' est au-dessus+ --traceShow "6" $+ (ilow $ inter0 yui yuj, iup $ inter0 yli ylj)+ else+ -- M est au-dessus de la ligne+ if ylj<=0 then+ if yuj>=0 then+ -- 7 M' est sur la ligne+ --traceShow "7" $+ (ilow $ inter0 yli ylj,fj)+ else+ -- 8 M' est en-dessous (comme M)+ --traceShow "8" $+ (ilow $ inter0 yli ylj, iup $ inter0 yui yuj)+ else+ -- 9 M' est au-dessus+ --traceShow "9" $+ (1/0,-1/0)+ + testAll i j m0 m1+ | i>=nv = + let fn=fromIntegral (nv-1)+ (# a0,b0 #)=over m0 m0 fn fn+ (# a1,b1 #)=minus b b a a+ (# a2,b2 #)=times a0 b0 a1 b1+ (# a3,_ #)=plus a a a2 b2+ + (# c0,d0 #)=over m1 m1 fn fn+ (# c2,d2 #)=times c0 d0 a1 b1+ (# _,d3 #)=plus a a c2 d2+ in+ (False,a3,d3)+ | j>=nv = testAll (i+1) (i+1) m0 m1+ | otherwise = + let (m0',m1')=inter i j in+ testAll i (j+1)+ (min m0 m0') (max m1 m1')+ in+ --traceShow allzero $+ if allzero then (True,a,b) else+ testAll 0 0 (1/0) (-1/0)+--traceShow ("convexHull",pointsl,pointsu,m) $ m++class Box a i | a->i where+ cut::Int->a->[a]+ size::Int->a->Double+ restriction#::a->Bernsteinp i Interval->a+ variables::a->Int++instance Box (Double,Double) Int where+ cut _ (a,b)=+ let m=(a+b)/2 in+ if a<m && m<b then+ [(a,m),(m,b)]+ else+ [(a,b)]+ size _ (a,b)=b-a+ restriction# (a,b) points@(Bernsteinp n0 _)=+ let restr=restriction points (Interval a a) (Interval b b)+ (allz,a',b')=convexHull 1 1 n0 restr a b+ in+ (max a a', min b b')+ variables _ = 1+instance Box (Double,Double,Double,Double) (Int,Int) where+ cut 0 x@(a,b,c,d)=+ let m=(a+b)/2 in+ if a<m && m<b then+ [(a,m,c,d),(m,b,c,d)]+ else+ [x]+ cut _ x@(a,b,c,d)=+ let m=(c+d)/2 in+ if c<m && m<d then+ [(a,b,c,m),(a,b,m,d)]+ else+ [x]+ size 0 (a,b,_,_)=b-a+ size _ (_,_,a,b)=b-a+ + restriction# (a,b,c,d) points@(Bernsteinp (n0,n1) _)=+ let restr=restriction points (Interval a a,Interval c c) (Interval b b,Interval d d)+ (allz0,a',b')+ | n0>1 = convexHull 1 n1 n0 restr a b+ | otherwise = (False,a,b)+ (allz1,c',d')+ | n1>1 = convexHull n0 1 n1 restr c d+ | otherwise = (False,c,d)+ in+ (max a a', min b b', max c c', min d d')+ variables _=2+++instance Box (Double,Double,Double,Double,Double,Double) (Int,Int,Int) where+ cut 0 x@(a,b,c,d,e,f)=+ let m=(a+b)/2 in+ if a<m && m<b then+ [(a,m,c,d,e,f),(m,b,c,d,e,f)]+ else+ [x]+ cut 1 x@(a,b,c,d,e,f)=+ let m=(c+d)/2 in+ if c<m && m<d then+ [(a,b,c,m,e,f),(a,b,m,d,e,f)]+ else+ [x]+ cut _ x@(a,b,c,d,e,f)=+ let m=(e+f)/2 in+ if e<m && m<f then+ [(a,b,c,d,e,m),(a,b,c,d,m,f)]+ else+ [x]+ size 0 (a,b,_,_,_,_)=b-a+ size 1 (_,_,a,b,_,_)=b-a+ size _ (_,_,_,_,a,b)=b-a+ + restriction# (a,b,c,d,e,f) points@(Bernsteinp (n0,n1,n2) _)=+ let restr=restriction points (Interval a a,Interval c c,Interval e e)+ (Interval b b,Interval d d,Interval f f)+ (allz0,a',b')+ | n0>1 = convexHull 1 (n1*n2) n0 restr a b+ | otherwise = (False,a,b)+ (allz1,c',d')+ | n1>1 = convexHull n0 n2 n1 restr c d+ | otherwise = (False,c,d)+ (allz2,e',f')+ | n2>1 = convexHull (n0*n1) 1 n2 restr e f+ | otherwise = (False,e,f)+ in+ (max a a', min b b', max c c', min d d', max e e', min f f')+ + variables _=3+ +instance Box (Double,Double,Double,Double,Double,Double,Double,Double) (Int,Int,Int,Int) where+ cut 0 x@(a,b,c,d,e,f,g,h)=+ let m=(a+b)/2 in+ if a<m && m<b then+ [(a,m,c,d,e,f,g,h),(m,b,c,d,e,f,g,h)]+ else+ [x]+ cut 1 x@(a,b,c,d,e,f,g,h)=+ let m=(c+d)/2 in+ if c<m && m<d then+ [(a,b,c,m,e,f,g,h),(a,b,m,d,e,f,g,h)]+ else+ [x]+ cut 2 x@(a,b,c,d,e,f,g,h)=+ let m=(e+f)/2 in+ if e<m && m<f then+ [(a,b,c,d,e,m,g,h),(a,b,c,d,m,f,g,h)]+ else+ [x]+ cut _ x@(a,b,c,d,e,f,g,h)=+ let m=(g+h)/2 in+ if g<m && m<h then+ [(a,b,c,d,e,f,g,m),(a,b,c,d,e,f,m,h)]+ else+ [x]+ + size 0 (a,b,_,_,_,_,_,_)=b-a+ size 1 (_,_,a,b,_,_,_,_)=b-a+ size 2 (_,_,_,_,a,b,_,_)=b-a+ size _ (_,_,_,_,_,_,a,b)=b-a+ + restriction# (a,b,c,d,e,f,g,h) points@(Bernsteinp (n0,n1,n2,n3) _)=+ let restr=restriction points (Interval a a,Interval c c,Interval e e,Interval g g)+ (Interval b b,Interval d d,Interval f f,Interval h h)+ (allz0,a',b')+ | n0>1 = convexHull 1 (n1*n2*n3) n0 restr a b+ | otherwise = (False,a,b)+ (allz1,c',d')+ | n1>1 = convexHull n0 (n2*n3) n1 restr c d+ | otherwise = (False,c,d)+ (allz2,e',f')+ | n2>1 = convexHull (n0*n1) n3 n2 restr e f+ | otherwise = (False,e,f)+ (allz3,g',h')+ | n3>1 = convexHull (n0*n1*n2) 1 n3 restr g h+ | otherwise = (False,g,h)+ in+ --traceShow restr $+ (max a a', min b b', max c c', min d d', + max e e', min f f', max g g', min h h')+ + variables _=4++-- | Computes the intersection of a given Bezier hypersurface, given+-- by its graph, with plane @z=0@.+solve::(Show a,Show i,Eq a,Box a i)=>Double->V.Vector (Bernsteinp i Interval)->a->[a]+solve sizeThreshold equations h= -- traceShow h $+ let splitThreshold=0.95+ restrictAll neq box+ | neq>=V.length equations = box+ | not (check 0 box) = box+ | otherwise =+ let next=restriction# box (equations!neq) in+ restrictAll (neq+1) next+ check v box=+ (v>=(variables box)) ||+ (let s=size v box in+ 0<=s && s<(1/0) && check (v+1) box)+ + h'=restrictAll 0 h+ + isSmall v box=+ (v>=variables box) ||+ ((size v box <= sizeThreshold) && (isSmall (v+1) box))+ + in+ if isSmall 0 h' then+ if check 0 (restrictAll 0 h') then+ [h']+ else+ []+ else+ if check 0 h' then+ let cutAll v boxes+ | v>=(variables h) = boxes+ | otherwise =+ cutAll (v+1) $+ Prelude.concatMap (\b->if (size v b)>=splitThreshold*(size v h) + && (size v b)>sizeThreshold+ then+ cut v b+ else [b]) boxes+ cc=cutAll 0 [h']+ in+ case cc of+ [h'']+ | h''==h -> + [h]+ | otherwise -> Prelude.concatMap (solve sizeThreshold equations) cc+ _->+ Prelude.concatMap (solve sizeThreshold equations) cc+ else+ []+\end{code}
+ Algebra/Polynomials/Numerical.hs view
@@ -0,0 +1,263 @@+{-# CFILES cnumerical.c #-}+{-# OPTIONS -XUnboxedTuples -XMagicHash -XScopedTypeVariables -XBangPatterns -cpp -XTypeFamilies -XMultiParamTypeClasses #-}+{-# LANGUAGE ForeignFunctionInterface #-}+-- | This module contains the definition of the main arithmetic tools+-- used in Metafont'.+module Algebra.Polynomials.Numerical(+ -- * Raw operations+ fromIntegral#,plus,minus,over,times,+ sqrt#,cos#,sin#,acos#,asin#,+ -- * The 'Interval' type+ Interval(..),Intervalize(..),+ interval,intersectsd, union,+ fpred,fsucc+ ) where+++import Data.Vector.Unboxed as UV+import qualified Data.Vector.Generic.Mutable as GMV+import qualified Data.Vector.Generic as GV+import Foreign.C.Types+foreign import ccall unsafe "c_succ" c_fsucc::CDouble->CDouble+foreign import ccall unsafe "c_pred" c_fpred::CDouble->CDouble++fsucc,fpred::Double->Double+fpred=realToFrac.c_fpred.realToFrac+fsucc=realToFrac.c_fsucc.realToFrac++{-# INLINE plus #-}+-- | Interval addition+plus::Double->Double->Double->Double->(# Double, Double #)+plus !a !b !c !d=+ let !x=a+c+ !y=b+d+ in+ (# fpred x,fsucc y #)+ +-- | Interval substraction+{-# INLINE minus #-}+minus::Double->Double->Double->Double->(# Double, Double #)+minus !a !b !c !d=+ let !x=a-d+ !y=b-c+ in+ (# fpred x, fsucc y #)+-- | Interval multiplication+{-# INLINE times #-}+times::Double->Double->Double->Double->(# Double, Double #)+times !a !b !c !d=+ let !w=a*c+ !x=a*d+ !y=b*c+ !z=b*d+ + (# !aa,!bb #)=if w<x then (# w,x #) else (# x,w #)+ (# !cc,!dd #)=if y<z then (# y,z #) else (# z,y #)+ !m=min aa cc+ !m'=max bb dd+ in+ (# fpred m, fsucc m' #)++-- | Interval division+{-# INLINE over #-}+over::Double->Double->Double->Double->(# Double, Double #)+over !a !b !c !d=+ if c*d<=0 then + if a>0 then (# 1/0,1/0 #) else+ if b<0 then (# (-1/0), (-1/0) #) else+ (# 0/0, 0/0 #)+ + else+ let !w=a/c+ !x=a/d+ !y=b/c+ !z=b/d+ + !(aa,bb)=if w<x then (w,x) else (x,w)+ !(cc,dd)=if y<z then (y,z) else (z,y)+ !m=min aa cc+ !m'=max bb dd+ in+ (# fpred m, fsucc m' #)++-- | Converts an 'Integral' value into an interval.+fromIntegral#::Integral x=>x->(# Double,Double #)+fromIntegral# n=+ let !n_=fromIntegral n in+ (# fpred n_,fsucc n_ #)++-- | Interval cosine+cos#::Double->Double->(# Double,Double #)+cos# !a !b=+ let (# !_m0,!_m0' #)=if cos a<=cos b then (# cos a, cos b #) else (# cos b, cos a #)+ !m0=fpred _m0+ !m0'=fsucc _m0'+ checkUp !(k::Int) !m !m'=+ let (# !ka,!kb #)=fromIntegral# k+ (# !ka0,!kb0 #)=times ka kb (fpred pi) (fsucc pi)+ in+ if ka0>b then (# m,m' #) else+ if kb0<a then+ checkUp (k+1) m m'+ else+ if k`mod`2==0 then+ checkUp (k+1) m 1+ else+ checkUp (k+1) (-1) m'+ in+ checkUp (floor $ fpred (a/pi)) m0 m0'+-- | Interval sine+sin# ::Double->Double->(# Double,Double #)+sin# !a !b=+ let (# _m0,_m0' #)+ | sin a<sin b = (# sin a, sin b #)+ | otherwise = (# sin b, sin a #)+ m0=max (-1) $ fpred _m0+ m0'=min 1 $ fsucc _m0'+ (# pa,pb #)=(# fpred pi, fsucc pi #)+ (# ka1,kb1 #)=over pa pb 2 2+ + up (k::Int) !m !m'=+ let (# ka,kb #)=fromIntegral# k+ (# ka0,kb0 #)=times ka kb pa pb+ (# ka2,kb2 #)=plus ka0 kb0 ka1 kb1 -- kpi+pi/2+ in+ if ka2>b then+ (# m,m' #)+ else+ if kb2<a then+ up (k+1) m m'+ else+ if k`mod`2 == 0 then+ up (k+1) m 1+ else+ up (k+1) (-1) m'+ in+ up (floor $ a/pi) m0 m0'+ + +sqrt#::Double->Double->(# Double,Double #)+sqrt# !a !b=+ let sa=sqrt a+ sb=sqrt b+ sa_=max 0 (fpred sa)+ sb_=fsucc sb+ in+ (# sa_, sb_ #)++acos#::Double->Double->(# Double,Double #)+acos# !a !b=+ let aca=acos $ max (-1) a+ acb=acos $ min 1 b+ in+ (# fpred (min aca acb), fsucc (max aca acb) #)++asin#::Double->Double->(# Double,Double #)+asin# !a !b=+ let aca=asin $ max (-1) a+ acb=asin $ min 1 b+ in+ (# fpred (min aca acb), fsucc (max aca acb) #)++-- | The interval type (most of its operations are calls to the raw functions)+data Interval=Interval {ilow::Double,iup::Double} deriving (Eq, Show)++instance Floating Interval where+ cos (Interval a b)=+ let (# c,d #)=cos# a b in+ Interval c d+ sin (Interval a b)=+ let (# c,d #)=sin# a b in+ Interval c d+ sqrt (Interval a b)=+ let (# a#,b# #)=sqrt# a b in+ Interval a# b#+ acos (Interval a b)=+ let (# a#,b# #)=acos# a b in+ Interval a# b#+ asin (Interval a b)=+ let (# a#,b# #)=asin# a b in+ Interval a# b#+ pi=Interval (fpred pi) (fsucc pi)+ + +-- | Intersection of two 'Interval's.+{-# INLINE intersectsd #-}+intersectsd::Interval->Interval->Bool+intersectsd (Interval a b) (Interval c d) = b>=c && a<=d++-- | Union of two intersecting intervals (undefined behaviour if they do not intersect).+{-# INLINE union #-}+union::Interval->Interval->Interval+union (Interval a b) (Interval c d) = Interval (min a c) (max b d)++-- | Two common operations on types defined with intervals.+class Intervalize a where+ intervalize::a Double->a Interval+ intersects::a Interval->a Interval->Bool+ +-- | Converts an optimal IEEE-754 representation of a number into an+-- optimal interval containing this number.+interval::Double->Interval+interval x=Interval (fpred x) (fsucc x)++instance Num Interval where+ (+) (Interval a b) (Interval c d)=+ let (# a',b' #)=plus a b c d in+ Interval a' b'+ (-) (Interval a b) (Interval c d)=+ let (# a',b' #)=minus a b c d in+ Interval a' b'+ (*) (Interval a b) (Interval c d)=+ let (# a',b' #)=times a b c d in+ Interval a' b'+ abs x@(Interval a b)=+ if b<=0 then Interval (negate b) (negate a) else+ if a<=0 then+ Interval 0 (max b $ negate a)+ else+ x+ + signum _=undefined+ + fromInteger=interval.fromInteger+ +instance Fractional Interval where++ (/) (Interval a b) (Interval c d)=+ let (# a',b' #)=over a b c d in+ Interval a' b'++ fromRational r=+ let r'=fromRational r in+ Interval (fpred r') (fsucc r')+++newtype instance UV.MVector s Interval = MV_Interval (UV.MVector s (Double,Double))+newtype instance UV.Vector Interval = V_Interval (UV.Vector (Double,Double))+instance Unbox Interval++instance GMV.MVector UV.MVector Interval where+ basicLength (MV_Interval a)=GMV.basicLength a+ basicUnsafeSlice a b (MV_Interval c)=MV_Interval $ GMV.basicUnsafeSlice a b c+ basicOverlaps (MV_Interval a) (MV_Interval b)=GMV.basicOverlaps a b+ basicUnsafeNew a=GMV.basicUnsafeNew a >>= return.MV_Interval+ basicUnsafeReplicate a (Interval b c)=GMV.basicUnsafeReplicate a (b,c)>>=return.MV_Interval+ basicUnsafeRead (MV_Interval a) b=GMV.basicUnsafeRead a b >>= (\(u,v)->return $ Interval u v)+ basicUnsafeWrite (MV_Interval a) b (Interval c d)=GMV.basicUnsafeWrite a b (c,d)+ basicClear (MV_Interval a)=GMV.basicClear a+ basicSet (MV_Interval a) (Interval b c)=GMV.basicSet a (b,c)+ basicUnsafeCopy (MV_Interval a) (MV_Interval b)=GMV.basicUnsafeCopy a b+ basicUnsafeGrow (MV_Interval a) b=GMV.basicUnsafeGrow a b >>= return.MV_Interval++instance GV.Vector UV.Vector Interval where+ basicUnsafeFreeze (MV_Interval a)=GV.basicUnsafeFreeze a >>= return.V_Interval+ basicUnsafeThaw (V_Interval a)=GV.basicUnsafeThaw a >>= return.MV_Interval+ basicLength (V_Interval a)=GV.basicLength a+ basicUnsafeSlice a b (V_Interval c)=V_Interval (GV.basicUnsafeSlice a b c)+ basicUnsafeIndexM (V_Interval a) b=GV.basicUnsafeIndexM a b >>= (\(u,v)->return $ Interval u v)++(!#)::UV.Vector Interval->Int->(# Double,Double #)+(!#) a b=+ let Interval u v=a!b in (# u,v #)+
+ Algebra/Polynomials/cnumerical.c view
@@ -0,0 +1,143 @@+#include <stdio.h>+#include <HsFFI.h>++union stg_ieee754_dbl+{+ double d;+ struct {++#if WORDS_BIGENDIAN+ unsigned int negative:1;+ unsigned int exponent:11;+ unsigned int mantissa0:20;+ unsigned int mantissa1:32;+#else+#if FLOAT_WORDS_BIGENDIAN+ unsigned int mantissa0:20;+ unsigned int exponent:11;+ unsigned int negative:1;+ unsigned int mantissa1:32;+#else+ unsigned int mantissa1:32;+ unsigned int mantissa0:20;+ unsigned int exponent:11;+ unsigned int negative:1;+#endif+#endif+ } ieee;+ /* This format makes it easier to see if a NaN is a signalling NaN. */+ struct {++#if WORDS_BIGENDIAN+ unsigned int negative:1;+ unsigned int exponent:11;+ unsigned int quiet_nan:1;+ unsigned int mantissa0:19;+ unsigned int mantissa1:32;+#else+#if FLOAT_WORDS_BIGENDIAN+ unsigned int mantissa0:19;+ unsigned int quiet_nan:1;+ unsigned int exponent:11;+ unsigned int negative:1;+ unsigned int mantissa1:32;+#else+ unsigned int mantissa1:32;+ unsigned int mantissa0:19;+ unsigned int quiet_nan:1;+ unsigned int exponent:11;+ unsigned int negative:1;+#endif+#endif+ } ieee_nan;+};++++double c_succ(double y)+{+ union stg_ieee754_dbl su;+ + su.d=y;+ if (su.ieee.negative==0) { /* y >= 0 */+ if (su.ieee.exponent!=2047 || su.ieee.mantissa0!=0 || su.ieee.mantissa1!=0)+ if (su.ieee.mantissa1==0xffffffff) { + su.ieee.mantissa1=0; + if (su.ieee.mantissa0==1048575) { + su.ieee.mantissa0=0; + su.ieee.exponent++;+ } else { + su.ieee.mantissa0++;+ }+ } else { + su.ieee.mantissa1++;+ }+ } + else { /* y < 0 */+ if (su.ieee.exponent!=2047 || su.ieee.mantissa0!=0 || su.ieee.mantissa1==0){+ if (su.ieee.negative==1 && su.ieee.exponent==0 && su.ieee.mantissa0==0 && su.ieee.mantissa1==0) {+ su.ieee.negative=0;+ su.ieee.mantissa1=1;+ } else {+ if (su.ieee.mantissa1==0) { + su.ieee.mantissa1=0xffffffff; + if (su.ieee.mantissa0==0) { + su.ieee.mantissa0=1048575; + su.ieee.exponent--;+ } else { + su.ieee.mantissa0--;+ }+ } else { + su.ieee.mantissa1--;+ }+ }+ }+ }+ return su.d;+} /* end function q_succ */+++++++double c_pred(double y)+{+ union stg_ieee754_dbl su;++ su.d=y;+ if (su.ieee.negative==1) { /* y < 0 */+ if (su.ieee.exponent!=2047 || su.ieee.mantissa0!=0 || su.ieee.mantissa1!=0 ) + if (su.ieee.mantissa1==0xffffffff) { + su.ieee.mantissa1=0; + if (su.ieee.mantissa0==1048575) { + su.ieee.mantissa0=0; + su.ieee.exponent++;+ } else { + su.ieee.mantissa0++;+ }+ } else+ su.ieee.mantissa1++;+ } else { /* y >= 0 */+ if (su.ieee.exponent!=2047 || su.ieee.mantissa0!=0 || su.ieee.mantissa1!=0) + if (su.ieee.exponent==0 && su.ieee.mantissa0==0 && su.ieee.mantissa1==0) {+ su.ieee.negative=1;+ su.ieee.mantissa1=1;+ } else {+ if (su.ieee.mantissa1==0) {+ su.ieee.mantissa1=0xffffffff; + if (su.ieee.mantissa0==0) { + su.ieee.mantissa0=1048575; + su.ieee.exponent--;+ } else { + su.ieee.mantissa0--;+ }+ } else { + su.ieee.mantissa1--;+ }+ }+ }+ + return su.d;+} /* end function q_pred */+
+ LICENSE view
@@ -0,0 +1,340 @@+ GNU GENERAL PUBLIC LICENSE+ Version 2, June 1991++ Copyright (C) 1989, 1991 Free Software Foundation, Inc.+ 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++ Preamble++ The licenses for most software are designed to take away your+freedom to share and change it. By contrast, the GNU General Public+License is intended to guarantee your freedom to share and change free+software--to make sure the software is free for all its users. This+General Public License applies to most of the Free Software+Foundation's software and to any other program whose authors commit to+using it. (Some other Free Software Foundation software is covered by+the GNU Library General Public License instead.) You can apply it to+your programs, too.++ When we speak of free software, we are referring to freedom, not+price. Our General Public Licenses are designed to make sure that you+have the freedom to distribute copies of free software (and charge for+this service if you wish), that you receive source code or can get it+if you want it, that you can change the software or use pieces of it+in new free programs; and that you know you can do these things.++ To protect your rights, we need to make restrictions that forbid+anyone to deny you these rights or to ask you to surrender the rights.+These restrictions translate to certain responsibilities for you if you+distribute copies of the software, or if you modify it.++ For example, if you distribute copies of such a program, whether+gratis or for a fee, you must give the recipients all the rights that+you have. You must make sure that they, too, receive or can get the+source code. And you must show them these terms so they know their+rights.++ We protect your rights with two steps: (1) copyright the software, and+(2) offer you this license which gives you legal permission to copy,+distribute and/or modify the software.++ Also, for each author's protection and ours, we want to make certain+that everyone understands that there is no warranty for this free+software. If the software is modified by someone else and passed on, we+want its recipients to know that what they have is not the original, so+that any problems introduced by others will not reflect on the original+authors' reputations.++ Finally, any free program is threatened constantly by software+patents. We wish to avoid the danger that redistributors of a free+program will individually obtain patent licenses, in effect making the+program proprietary. To prevent this, we have made it clear that any+patent must be licensed for everyone's free use or not licensed at all.++ The precise terms and conditions for copying, distribution and+modification follow.++ GNU GENERAL PUBLIC LICENSE+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION++ 0. This License applies to any program or other work which contains+a notice placed by the copyright holder saying it may be distributed+under the terms of this General Public License. The "Program", below,+refers to any such program or work, and a "work based on the Program"+means either the Program or any derivative work under copyright law:+that is to say, a work containing the Program or a portion of it,+either verbatim or with modifications and/or translated into another+language. (Hereinafter, translation is included without limitation in+the term "modification".) Each licensee is addressed as "you".++Activities other than copying, distribution and modification are not+covered by this License; they are outside its scope. The act of+running the Program is not restricted, and the output from the Program+is covered only if its contents constitute a work based on the+Program (independent of having been made by running the Program).+Whether that is true depends on what the Program does.++ 1. You may copy and distribute verbatim copies of the Program's+source code as you receive it, in any medium, provided that you+conspicuously and appropriately publish on each copy an appropriate+copyright notice and disclaimer of warranty; keep intact all the+notices that refer to this License and to the absence of any warranty;+and give any other recipients of the Program a copy of this License+along with the Program.++You may charge a fee for the physical act of transferring a copy, and+you may at your option offer warranty protection in exchange for a fee.++ 2. You may modify your copy or copies of the Program or any portion+of it, thus forming a work based on the Program, and copy and+distribute such modifications or work under the terms of Section 1+above, provided that you also meet all of these conditions:++ a) You must cause the modified files to carry prominent notices+ stating that you changed the files and the date of any change.++ b) You must cause any work that you distribute or publish, that in+ whole or in part contains or is derived from the Program or any+ part thereof, to be licensed as a whole at no charge to all third+ parties under the terms of this License.++ c) If the modified program normally reads commands interactively+ when run, you must cause it, when started running for such+ interactive use in the most ordinary way, to print or display an+ announcement including an appropriate copyright notice and a+ notice that there is no warranty (or else, saying that you provide+ a warranty) and that users may redistribute the program under+ these conditions, and telling the user how to view a copy of this+ License. (Exception: if the Program itself is interactive but+ does not normally print such an announcement, your work based on+ the Program is not required to print an announcement.)++These requirements apply to the modified work as a whole. If+identifiable sections of that work are not derived from the Program,+and can be reasonably considered independent and separate works in+themselves, then this License, and its terms, do not apply to those+sections when you distribute them as separate works. But when you+distribute the same sections as part of a whole which is a work based+on the Program, the distribution of the whole must be on the terms of+this License, whose permissions for other licensees extend to the+entire whole, and thus to each and every part regardless of who wrote it.++Thus, it is not the intent of this section to claim rights or contest+your rights to work written entirely by you; rather, the intent is to+exercise the right to control the distribution of derivative or+collective works based on the Program.++In addition, mere aggregation of another work not based on the Program+with the Program (or with a work based on the Program) on a volume of+a storage or distribution medium does not bring the other work under+the scope of this License.++ 3. You may copy and distribute the Program (or a work based on it,+under Section 2) in object code or executable form under the terms of+Sections 1 and 2 above provided that you also do one of the following:++ a) Accompany it with the complete corresponding machine-readable+ source code, which must be distributed under the terms of Sections+ 1 and 2 above on a medium customarily used for software interchange; or,++ b) Accompany it with a written offer, valid for at least three+ years, to give any third party, for a charge no more than your+ cost of physically performing source distribution, a complete+ machine-readable copy of the corresponding source code, to be+ distributed under the terms of Sections 1 and 2 above on a medium+ customarily used for software interchange; or,++ c) Accompany it with the information you received as to the offer+ to distribute corresponding source code. (This alternative is+ allowed only for noncommercial distribution and only if you+ received the program in object code or executable form with such+ an offer, in accord with Subsection b above.)++The source code for a work means the preferred form of the work for+making modifications to it. For an executable work, complete source+code means all the source code for all modules it contains, plus any+associated interface definition files, plus the scripts used to+control compilation and installation of the executable. However, as a+special exception, the source code distributed need not include+anything that is normally distributed (in either source or binary+form) with the major components (compiler, kernel, and so on) of the+operating system on which the executable runs, unless that component+itself accompanies the executable.++If distribution of executable or object code is made by offering+access to copy from a designated place, then offering equivalent+access to copy the source code from the same place counts as+distribution of the source code, even though third parties are not+compelled to copy the source along with the object code.++ 4. You may not copy, modify, sublicense, or distribute the Program+except as expressly provided under this License. Any attempt+otherwise to copy, modify, sublicense or distribute the Program is+void, and will automatically terminate your rights under this License.+However, parties who have received copies, or rights, from you under+this License will not have their licenses terminated so long as such+parties remain in full compliance.++ 5. You are not required to accept this License, since you have not+signed it. However, nothing else grants you permission to modify or+distribute the Program or its derivative works. These actions are+prohibited by law if you do not accept this License. Therefore, by+modifying or distributing the Program (or any work based on the+Program), you indicate your acceptance of this License to do so, and+all its terms and conditions for copying, distributing or modifying+the Program or works based on it.++ 6. Each time you redistribute the Program (or any work based on the+Program), the recipient automatically receives a license from the+original licensor to copy, distribute or modify the Program subject to+these terms and conditions. You may not impose any further+restrictions on the recipients' exercise of the rights granted herein.+You are not responsible for enforcing compliance by third parties to+this License.++ 7. If, as a consequence of a court judgment or allegation of patent+infringement or for any other reason (not limited to patent issues),+conditions are imposed on you (whether by court order, agreement or+otherwise) that contradict the conditions of this License, they do not+excuse you from the conditions of this License. If you cannot+distribute so as to satisfy simultaneously your obligations under this+License and any other pertinent obligations, then as a consequence you+may not distribute the Program at all. For example, if a patent+license would not permit royalty-free redistribution of the Program by+all those who receive copies directly or indirectly through you, then+the only way you could satisfy both it and this License would be to+refrain entirely from distribution of the Program.++If any portion of this section is held invalid or unenforceable under+any particular circumstance, the balance of the section is intended to+apply and the section as a whole is intended to apply in other+circumstances.++It is not the purpose of this section to induce you to infringe any+patents or other property right claims or to contest validity of any+such claims; this section has the sole purpose of protecting the+integrity of the free software distribution system, which is+implemented by public license practices. Many people have made+generous contributions to the wide range of software distributed+through that system in reliance on consistent application of that+system; it is up to the author/donor to decide if he or she is willing+to distribute software through any other system and a licensee cannot+impose that choice.++This section is intended to make thoroughly clear what is believed to+be a consequence of the rest of this License.++ 8. If the distribution and/or use of the Program is restricted in+certain countries either by patents or by copyrighted interfaces, the+original copyright holder who places the Program under this License+may add an explicit geographical distribution limitation excluding+those countries, so that distribution is permitted only in or among+countries not thus excluded. In such case, this License incorporates+the limitation as if written in the body of this License.++ 9. The Free Software Foundation may publish revised and/or new versions+of the General Public License from time to time. Such new versions will+be similar in spirit to the present version, but may differ in detail to+address new problems or concerns.++Each version is given a distinguishing version number. If the Program+specifies a version number of this License which applies to it and "any+later version", you have the option of following the terms and conditions+either of that version or of any later version published by the Free+Software Foundation. If the Program does not specify a version number of+this License, you may choose any version ever published by the Free Software+Foundation.++ 10. If you wish to incorporate parts of the Program into other free+programs whose distribution conditions are different, write to the author+to ask for permission. For software which is copyrighted by the Free+Software Foundation, write to the Free Software Foundation; we sometimes+make exceptions for this. Our decision will be guided by the two goals+of preserving the free status of all derivatives of our free software and+of promoting the sharing and reuse of software generally.++ NO WARRANTY++ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,+REPAIR OR CORRECTION.++ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE+POSSIBILITY OF SUCH DAMAGES.++ END OF TERMS AND CONDITIONS++ How to Apply These Terms to Your New Programs++ If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++ To do so, attach the following notices to the program. It is safest+to attach them to the start of each source file to most effectively+convey the exclusion of warranty; and each file should have at least+the "copyright" line and a pointer to where the full notice is found.++ <one line to give the program's name and a brief idea of what it does.>+ Copyright (C) <year> <name of author>++ This program is free software; you can redistribute it and/or modify+ it under the terms of the GNU General Public License as published by+ the Free Software Foundation; either version 2 of the License, or+ (at your option) any later version.++ This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License for more details.++ You should have received a copy of the GNU General Public License+ along with this program; if not, write to the Free Software+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA+++Also add information on how to contact you by electronic and paper mail.++If the program is interactive, make it output a short notice like this+when it starts in an interactive mode:++ Gnomovision version 69, Copyright (C) year name of author+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.+ This is free software, and you are welcome to redistribute it+ under certain conditions; type `show c' for details.++The hypothetical commands `show w' and `show c' should show the appropriate+parts of the General Public License. Of course, the commands you use may+be called something other than `show w' and `show c'; they could even be+mouse-clicks or menu items--whatever suits your program.++You should also get your employer (if you work as a programmer) or your+school, if any, to sign a "copyright disclaimer" for the program, if+necessary. Here is a sample; alter the names:++ Yoyodyne, Inc., hereby disclaims all copyright interest in the program+ `Gnomovision' (which makes passes at compilers) written by James Hacker.++ <signature of Ty Coon>, 1 April 1989+ Ty Coon, President of Vice++This General Public License does not permit incorporating your program into+proprietary programs. If your program is a subroutine library, you may+consider it more useful to permit linking proprietary applications with the+library. If this is what you want to do, use the GNU Library General+Public License instead of this License.
+ Setup.hs view
@@ -0,0 +1,3 @@+import Distribution.Simple++main=defaultMain
+ polynomials-bernstein.cabal view
@@ -0,0 +1,22 @@+Name: polynomials-bernstein+Version: 1+Synopsis: A solver for systems of polynomial equations in bernstein form+Description: This library defines an optimized type for representing polynomials+ in Bernstein form, as well as instances of numeric classes and other+ manipulation functions, and a solver of systems of polynomial+ equations in this form.+Category: Math+Maintainer: Pierre-Etienne Meunier <pierreetienne.meunier@gmail.com>+License: GPL+License-file: LICENSE+Build-Type: Simple+Cabal-Version: >=1.6+source-repository this+ type: darcs+ location: http://www.lama.univ-savoie.fr/~meunier/darcs/polynomials+ tag: 1.0+Library+ Build-Depends: base<5, vector+ Exposed-modules: Algebra.Polynomials.Bernstein, Algebra.Polynomials.Numerical+ ghc-options: -O2 -Wall+ c-sources: algebra/polynomials/cnumerical.c